I Gammoudi
email: ibtisem.gammoudi@u-bordeaux1.fr

V Raimbault

H Tarbague

F Moroté

C Grauby-Heywang

A Othmane

R Kalfat

D Moynet

D Rebière

C Dejous

T Cohen-Bouhacina

Enhanced bio-inspired microsensor based on microfluidic/bacteria/love wave hybrid structure for continuous control of heavy metals toxicity in liquid medium

Keywords:

Environmental pollution by toxic heavy metals (HM) presents a real threat for aquatic medium and human health. Therefore aquatic ecosystem management requires early warning systems for on line monitoring. Microtechnologies can give rise to innovative bio-inspired hybrid microsensors, likely to meet this need and providing cost reductions by reducing reagents consumption and manufacturing cost.

This work deals with a bacteria-based Love wave sensor, with enhanced properties provided by integration of a polydimethylsiloxane (PDMS) microfluidic network for a better control of the sample flow, and devoted to in situ monitoring of Cd(II) and Hg(II). Whole Escherichia coli (E. coli) bacteria are used as bioreceptor, mimicking in vivo enzymatic activity. They were immobilized on polyelectrolyte multilayer (PEM) films realized using layer by layer technique (LbL) with alternatively adsorption of positive and negative chains. The acoustic delay line was inserted into an electronic oscillation loop for real time monitoring.

Compared to previous work, this paper deepens the results obtained with two types of microfluidic chips (measurements in static and dynamic modes), including analysis in terms of reproducibility. These results are analyzed and interpreted thoroughly leading to assumptions about the phenomena involved in the detection mechanisms. These hypotheses are validated through a pioneering study with atomic force microscopy (AFM), performed to characterize bacteria adhesion and to establish the relationship between bacteria morphological evolution and mechanical properties. AFM was chosen for its ability to maintain the bacteria alive during the experience without inducing irreversible damage.

The resulting microsystem led to efficient HM detection, characterized by a reduced response-time (less than 60 s) and a detection limit inferior to 10 -12 M. AFM measurements have demonstrated a high bacterial attachment and the stressing effect of toxic HM on bacterial morphological state. These results are consistent with those obtained from Love wave measurements.

Introduction

Heavy metal (HM) pollution is posing significant threats to the environment and serious public health problems because of its toxicity (Bertin et al., 2006, [1]), non-biodegradability and bioaccumulation [2,3]. Indeed, HM may enter the human body through food, water, air or by skin absorption, and they become toxic when they are not metabolized by the body and accumulate in the soft tissues. Therefore, it is important to be aware of risks, to be able to track HM and to understand how they act in order to take protective measures against excessive exposure.

Many techniques as X-ray fluorescence [4], atomic emission [5] or absorption spectrometry [6] or chromatography [7] are employed for metal detection but these methods are poorly suited to in situ analysis.

Sensor is seen as a benefit compared to other means which are heavy and cumbersome. In recent years, biosensor has become an integral part of all technological advances [8,[START_REF] Hepel | Theory, Experiment and Applications[END_REF]. Several sensors have been developed for toxic compounds detection including Cadmium HM [START_REF] Sartore | [END_REF]11]. The aim of this study is to develop portable and disposable sensors. The use of bacteria as a bioreceptor mimicking in vivo enzymatic activity, can allows us to move from a multi-enzyme to a cellular bioreceptor, from in vitro to in vivo biochemical recognition. It is also a good way to study effects of HM on microorganisms, as HM generally exert an inhibitory action on microorganisms by blocking essential functional groups, displacing essential metal ions, or modifying the active conformations of biological molecules [12,13]. Exposure to HM ions can also lead to cellular disfunctioning like DNA damage [14], enzyme alteration and cell membrane disruption [15].

Biochemistry analysis systems used with microelectronic techniques are widely used in biosensor's field. Therefore the study of microfluidic sensor has become an important aspect of the analytical instrumentation. One of the polymers commonly used in developing micro-biosensor operating in liquid medium is poly-dimethylsiloxane (PDMS) [16][17][18], due to its versatile and favorable properties, which include non-toxicity, biocompatibility, flexibility, low cost, and easy fabrication [19]. Microfluidic devices obtained by photolithography and molding are now widely used because of their numerous advantages such as sample holding, reagent mixing, separation and detection [20,21].

In this paper, we present a hybrid structure combining PDMS microfluidic chip for flow control and bacteria-based Love wave transducer, devoted to in situ monitoring of HM for environmental control applications. Materials and methods involved in the device are presented in Section 2. Then, results obtained during PEM and bacteria deposition, as well as during detection of Mercury and Cadmium, are exposed and analyzed. They are compared to previous results obtained with PDMS cell for static measurements [22], bringing to light enhanced response time and sensitivity allowed by dynamical measurements with microfluidic pattern. Monitoring of bacterial stress with AFM is also used for a better understanding of mechanisms involved in the interaction between bacteria and HM.

Materials and methods

Love wave sensor

The surface acoustic wave component has been chosen for its numerous advantages, specially high sensitivity, real-time working, robustness and versatility due to mechanical transducing effect. Details on the device can be found in ref. [23]. Briefly, it consists of a dual delay line deposited on an AT cut quartz substrate (Euler angles: 0 • , 121.5 • , and 90 • ) used as the piezoelectric material. For each delay line, both Interdigital Transducers (IDTs) used for piezoelectric emission and reception of acoustic wave are made by sputtering 70 nm of gold on a 30 nm titanium layer to achieve a good are made by sputtering 70 nm of gold on a 30 nm titanium layer to achieve a good surface adhesion. IDTs are composed of 44 Ti/Au splitted-finger pairs with a periodicity ( ) of 40 m. The acoustic path length is 164 (between both IDTs), and the IDTs aperture (W) is 39 [24]. A 4 m plasma enhanced chemical vapor deposited (PECVD) SiO 2 layer is used as the guiding layer. Finally, the SiO 2 was etched upon electrical contacts. These characteristics lead to a guided shear horizontal surface acoustic wave (guided SH-SAW), also called Love wave, with synchronous frequency f 0 . This synchronous frequency is real time monitored with oscillating electronics [21] and represents the output signal of the sensor. It decreases typically upon added surfacic mass, and is also sensitive to other mechanical effects such as viscoelastic modifications [25].

PDMS chips

PDMS chips have been developed within previous works [26], in order to confine and control the fluid samples on the sensor surface. PDMS is a biocompatible material with interesting properties such as optical transparency, fast prototyping and good sealing properties on glass. The PDMS cell for static and dynamic measurements are schematically represented on Fig. 1b and can be briefly described. Both have chambers localizing fluids samples on the sensitive acoustic path between IDTs and air cavities upon IDTs to avoid capacitive coupling with conductive liquids. Both types of PDMS chips were pressed on the Love-wave component, instead of glued on it, in order to preserve easy surface accessibility. In the static configuration, fabrication process can be found in ref. [26], fluid samples were introduced with a micropipette in the open chambers. In the dynamic configuration micro channels were realized using standard soft photolithography methods [27] by means of silicon wafer with spin-coating of a SU-8 resin at 2000 rpm to obtain a negative mold with 200 m membranes as polymer template, together with micromachining for external 3D shapes. The PDMS chip was then formed from a liquid PDMS oligomer and a cross-linking agent (weight ratio of 9:1). Both components were mixed and put with a syringe pump at flow equal to 1 ml/min to obtain a good homogenization of PDMS in the mold. After that, PDMS was cured for 20 min at 95 • C and the resulting chip with the microfluidic network pattern was peeled of from the mold. Using a home-made test cell maintaining the PDMS chip with proper alignment and sealing on the Love wave device, the fluid sample was pulled through the acoustic path with programmable syringe pumps.

Chemical and biological materials

The sensitive biofilm of bacteria has been immobilized onto the surface of the Love wave platform through Polyelectrolytes (PE) multilayers. Polyelectrolytes used in our studies are a polycation, poly(allylamine hydrochloride) (PAH) (MW = 56 kDa), and a polyanion, poly (styrene sulfonate) (PSS) (MW = 70 kDa), both purchased from Sigma-Aldrich. Adsorption of PE was done from TBS (Tris Buffered saline) solutions of constant ionic strength (0.15 M) with a concentration of PE at 0.5 g/l.

PE were deposited through the "layer-by-layer" (LbL) method [28] resulting in self-assembled molecular multilayers PAH-(PSS-PAH) n . Typically, PAH and PSS solutions were alternatively deposited, adsorption and rinsing steps were then repeated to obtain a desired numbers of PE layers. In this work the PAH-PSS bi-layers ended with a PAH layer (positively charged) were used to immobilize E. coli (negatively charged). This bacterium is used as bioreceptor because of its sensitivity to HM, as presented in the introduction part. E. coli, acquired from the Centre d'Etudes du Bouchet (DGA, France), were grown overnight in Petri dishes containing LB solid agar nutritive medium. They were suspended in 0.15 M TBS (buffer) to obtain a final D.O. (optical density at 600 nm) equal to 0.6 corresponding to concentration of 10 8 cells/ml, that allowed good coverage of the biofilm on the sensor surface (results not shown).

The reference solutions of cadmium and mercury were prepared from Cd(NO 3 ) 2 •(H 2 O) 4 and Hg(NO 3 ) 2 •(H 2 O) (from Sigma-Aldrich) in the same TBS buffer (0.15 M). The platform was subjected to four selected concentrations of each metal: 10 -12 M, 10 -9 M, 10 -6 M and 10 -3 M. These concentrations were chosen as representative to study the sensitivity and reproducibility of the device over a large range and to evaluate the interest of such a sensor for environmental applications.

Atomic force microscopy

AFM was used to probe bacteria morphology and behavior in a physiological environment and in real time. In this work, the effect of the Cd 2+ concentration on the bacterial cell surface behavior was examined.

AFM experiments were carried out with a Bioscope II mounted on an Olympus inverted optical microscope and operating with the NanoScope V controller (Veeco-Brucker, Santa Barbara, CA). This technique is based on the detection of small forces existing between a tip and an object (the surface sample) which is scanned along the tip with a piezo scanner (maximum XYZ scan range of 150 m × 150 m × 12 m). The tip is attached on one end of a cantilever which is "built in" at the other end. This technique can generate information such as topography and mechanical properties at the nanometer level. Indeed, images obtained with AFM provide not only information about the (2D and 3D) topography of the studied sample (bacteria morphology) but also its mechanical properties such as adhesion, elasticity, viscoelasticity or chemical composition with friction and/or phase images [29][30][31][32][33]). In the present work, experiments were achieved in liquid and with the AFM contact mode. So, first, the biosensor surface is scanned in the three directions of space (imaging mode), this allows, for example, to provide information such as bacteria density on the biosensor surface, bacteria morphology (shape and size), etc. Then, and in order to access to the bacteria mechanical properties, force curves (FC) are obtained and collected from measuring the cantilever deflection as the sample was moved toward the tip in the z direction (spectroscopic mode).

All data presented in this paper were generated with the same cantilever (the spring constant of the cantilever used was nominally of about 0.58 N/m) and an imaging scan rate at 0.5 Hz. For each experiment, four images were recorded at the same time: trace and retrace height images, trace deflection image (signal error) and friction images (mechanical properties) [35]. For more clarity, images are flattened and only the trace height ones are shown in this paper.

Results and discussion

Acoustic responses during polyelectrolytes multilayers and bacteria deposition

The real-time characterization of PE multilayer coating is based on the modification of the acoustic wave phase velocity and so the oscillator frequency, due to mass loading of the surface.

Also, the frequency variation in real time has allowed us to optimize the flow rate. Indeed, we have selected each time a flow rate involving an improved frequency variation within the shortest time.

We will discuss in the following the results obtained for different stages of development and realization of biosensor.

The steps of the Fig. 1a (top curve) show a fast decrease in real-time signal during the injection in micro-flow of both PE alternatively, followed by signal stabilization that needs only a few minutes after each PE indicating the end of the adsorption step. Optimal throughput was set at 20 l/min. It is observed that for each one of the three bilayers, the steady-state frequency change increases with each new injection (Fig. 1a,bottom). This curve also highlights the good reproducibility of these results. It has been shown that a cleaning step was needed at the end of each layer, unlike when using a PDMS chip for static measurements, for which rinsing steps seemed to be unnecessary, as similar frequency shifts were obtained, compared to same protocol with rinsing steps [22].

Bacteria immobilization was then monitored in real-time. An optimal flow rate was found to 15 l/min. It was observed a decrease of electronic signal, followed by signal stabilization after 20 min. A frequency variation 3-5 times greater than the variation found for the immobilization of bacteria in static mode [22] was obtained.

In parallel, a study was performed by AFM to optimize the protocol for the polyelectrolytes deposition in order to ensure full substrate coverage with homogeneous film and reduced number of aggregates.

Acoustic responses during mercury and cadmium detection

The real-time responses to different concentrations of Cd 2+ and Hg 2+ were studied (Fig. 2). Optimal throughput was set at 15 l/min for heavy metals. Note an increasing frequency in the presence of heavy metals by increasing the levels and up to 10 -3 mol l -1 , with a detection threshold better than to 10 -12 mol l -1 . Beyond 10 -3 mol l -1 , a significant decrease in the frequency of several kilohertz can be observed. It was attributed to physico-chemical interactions, as similar effect is observed with a control line with PE and without bacteria (results not shown). This effect can partly explain the degradation of reproducibility appearing at high concentrations (10 -3 M).

More generally, these phenomena can be attributed to changes in the viscoelastic properties, related to changes in the bacterial metabolism. Indeed, HM interact with the constituents of the membrane and inhibit some enzymatic activities of bacterium inducing a decrease in its viability. Qualitatively similar effects were observed previously, when compared with the static configuration where the solutions were injected with micropipette [22]. It should be noted that characteristics have been improved with microfluidic system: the response time has been significantly reduced with effective detection within 60 s, and the sensitivity has been enhanced (Fig. 2c), with detection limit still inferior to 10 -12 M. These results demonstrate the importance of conducting tests under dynamic conditions. Also, the rate of injection of liquid and route to the sensitive surface of the sensor, as well as the flows form, significantly affect the interaction effectiveness between bacteria and the sensitive layer, as well as the effect mechanically induced on the transducing platform and therefore the response time and overall sensitivity. Thereby, with smaller analysis chambers (100-300 m in height), and a better distri- bution of liquid, the response time and in some cases the sensor sensitivity, can be improved.

The reproducibility was also verified and some results are reported on Fig. 1a (bottom) and Fig. 2a (bottom). As for PE and bacteria, responses to low concentration of HM show a good reproducibility. The degradation appeared at 10 -6 M and 10 -3 M.

Monitoring of bacterial stress by AFM

Thorough investigation of bacteria morphology and mechanical properties evolution as function of heavy metal concentration was made with AFM measurements. Among results, it appeared that HM (Cd 2+ ) had a definite effect on the behavior of the bacterial population immobilized on the biosensor. Fig. 3a shows typically rod-shaped cells with smooth surfaces well attached on substrate. The addition of Cd 2+ induced considerable damages resulting from changes on bacterial morphology (Fig. 3b) and on its adhesive properties. Indeed, a detachment of bacteria of the substrate can be observed (bacteria 3) following scanning tip while other bacteria remained attached (Number 1 and 2) but one can observe that their contours became less well defined which suggests a modification of the bacteria integrity. Besides, the AFM tip motion destabilized bacteria because it was sensitive to the external strain. This is result of hydrodynamic forces large enough to remove the cell (probably become less stable and less adherent), involving its displacement and then its detachment from substrate after a few scans (results not shown).

In order to complete this AFM study, we carried out nanoindentation comparative experiments on bacteria before and after injection of heavy metals (Fig. 3c). The obtained force curves are consistent with topographic results. Firstly, unlike the results obtained before injection of heavy metals, cells responses in presence of Cd 2+ are not reproducible. In fact, approach-retract curves differ from one cell to another and from one area of the cell to another, indicating important variability in cell stiffness linked to variable states of degraded cells. Second, and from a qualitative point of view, superposition of approach curves obtained with bacteria after addition of Cd 2+ to those of the substrate (PE layers) indicates the loss of cell elasticity behavior which can be correlated with degradation of bacteria. Also by analyzing the retract curve in the presence of cadmium, one notices the presence of several negative cantilever deflections while withdrawing the probe from the bacteria surface indicating tip contamination with the bacterial degraded cell content. In order to make a quantitative study, it would be interesting to follow the cell elastic behavior [29,33] when varying gradually the cadmium concentration. We plan to carry out this study in the very near future. Fig. 4.

Conclusion

A bacterial based Love wave microsensor was presented. Associated to new hydrodynamical PDMS chips integrating microfluidic network, this platform has lead to fast heavy metals detection, with efficient detection in less than 60 s with enhanced sensitivity and a detection limit better than 10 -12 M compared to the work of Stobiecka et al. [34] with a limit of detection LOD of 19 nM. AFM has emerged as a good complementary technique in this study. We observed real-time cell degradation following heavy metals injection through AFM images and force measurements on E. coli bacteria confirming the assumptions made on the basis of the results obtained with Love waves. This overall work put to evidence the interesting of this hybrid approach mixing electronics and micro-organisms for efficient detection of heavy metals, specially at very low concentration (<10 -6 M).

realizing acoustic wave delay-lines. AFM experiments were carried on the NSI platform.

Fig. 1 .

 1 Fig. 1. (a) Schematic view of a dual surface acoustic Love wave (b) hydrostatic PDMS cell (c) microfluidic PDMS cell.

Fig. 2 .

 2 Fig. 2. (a) Typical responses for PEM deposition with a microfluidic cell, top: real-time frequency; bottom: steady state frequency shift due to one PE layer (PAHi or PSSi) and reproducibility (circles: PAH, squares: PSS), mean values and error bars have been calculated from 5 experiments with different delay-lines, (b) Comparison of real time frequency responses to bacteria immobilization for static (top) and hydrodynamic (bottom) protocols.

Fig. 3 .

 3 Fig. 3. Response of the hybrid sensor to Cd(II)and Hg(II): (a), (b) (top curves).Typical real time frequency shifts (base frequency: 116,383 MHz) during E. coli (bioreceptor) immobilization and increasing concentrations of Cd(II) (a) and Hg(II) (b). (a), (b) (bottom curves). Steady-state frequency shifts vs. Cd(II) (a) and Hg(II) (b) concentrations, mean values and error bars have been calculated from 4 experiments with different delay-lines. (c) Comparison of cumulative steady-state frequency shifts |F-F0| (kHz) to heavy metals detection with static and hydrodynamic protocols.

Fig. 4 .

 4 Fig. 4. Bacteria population before (a) and after (b) addition of cadmium metal (contact mode AFM in liquid medium), (c) longitudinal profile of the bacterium (d) Force curve due to the elastic behavior of the bacteria before and after addition and (e) Force curve Zoom, longitudinal profile of the bacterium.

Acknowledgments

This work on the bacterial-based hybrid biosensor for heavy metals detection was financially supported by the CMCU, project No. 10G1103. PDMS chips were designed in the frame of ANR project (Agence Nationale de la Recherche, PRECODD-BIOALERT program, and project ANR-06-ECOT-004). The authors thank M. Benoît (Réseau Technologique de Base, LAAS-CNRS, Toulouse), for


Dante Kalise
email: dante.kalise@oeaw.ac.at

Axel Kr Öner
email: axel.kroener@inria.fr

AND Karl Kunisch
email: karl.kunisch@uni-graz.at.

LOCAL MINIMIZATION ALGORITHMS FOR DYNAMIC PROGRAMMING EQUATIONS

Keywords: dynamic programming, Hamilton-Jacobi-Bellman equations, semi-Lagrangian schemes, first order primal-dual methods, semismooth Newton methods AMS subject classifications

The numerical realization of the dynamic programming principle for continuous-time optimal control leads to nonlinear Hamilton-Jacobi-Bellman equations which require the minimization of a nonlinear mapping over the set of admissible controls. This minimization is often performed by comparison over a finite number of elements of the control set. In this paper we demonstrate the importance of an accurate realization of these minimization problems and propose algorithms by which this can be achieved effectively. The considered class of equations includes nonsmooth control problems with ℓ 1 -penalization which lead to sparse controls.

1. Introduction. Since its introduction by Bellman in the 50's, dynamic programming has become a fundamental tool in the design of optimal control strategies for dynamical systems. It characterizes the value function of the corresponding optimal control problem in terms of functional relations, the so-called Bellman and Hamilton-Jacobi-Bellman (henceforth HJB) equations. We begin by briefly recalling this setting in the context of infinite horizon optimal control. We make the following assumptions. We equip R n for n ∈ N with the Euclidean norm • 2 . Furthermore, let (1.1)

f : R d × R m → R, |f (x, u) -f (y, u)| ≤ ω R x -y 2 , l : R d × R m → R, |l(x, u) -l(y, u)| ≤ ω R x -y 2 ,
for x, y ∈ R d with xy 2 ≤ R and modulus ω R : [0, ∞) → [0, ∞) of polynomial growth satisfying lim r→0 + ω R (r) = 0 (cf. Ishii [START_REF] Ishii | Uniqueness of unbounded viscosity solution of Hamilton-Jacobi equations[END_REF]), d, m ∈ N. Let the dynamics be given by ẏ(t) = f (y(t), u(t)), y(0) = x (1.2) for t > 0, where x ∈ R d and u ∈ U ≡ {u : R + → U measurable}, U ⊂ R m compact. We introduce the following cost functional J : U → R J(u) = ∞ 0 l(y(s), u(s))e -λs ds , λ > 0 , where y is the solution of (1.2) depending on x and u. By the application of the dynamic programming principle, the value function

v(x) ≡ inf u∈U J(u)
is characterized as the viscosity solution [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF]Chapter 3] of the HJB equation λv(x) + sup u∈U {-f (x, u) • ∇v(x)l(x, u)} = 0, x ∈ R d .

There exists an extensive literature concerning the construction of numerical schemes for static HJB equations. The spectrum of numerical techniques includes ordered upwind methods [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: Theory and algorithms[END_REF][START_REF] Alton | An ordered upwind method with precomputed stencil and monotone node acceptance for solving static convex Hamilton-Jacobi equations[END_REF], high-order schemes [START_REF] Xiong | Fast sweeping fifth order WENO scheme for static Hamilton-Jacobi equations with accurate boundary treatment[END_REF], domain decomposition techniques [START_REF] Cacace | A patchy dynamic programming scheme for a class of Hamilton-Jacobi-Bellman equations[END_REF] and geometric approaches [START_REF] Botkin | Stable numerical schemes for solving Hamilton-Jacobi-Bellman-Isaacs equations[END_REF], among many others (we refer to [START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF]Chapter 5,p.145] for a review of classical approximation methods). In this paper we follow a semi-Lagrangian approach [START_REF] Falcone | A numerical approach to the infinite horizon problem of deterministic control theory[END_REF], which is broadly used to approximate HJB equations arising in optimal control problems, see,e.g., [START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF]. We illustrate the basic steps in the formulation of a semi-Lagrangian scheme for our model problem.

The construction of a first-order semi-Lagrangian scheme begins by considering an Euler discretization of the system dynamics with time step h > 0 y n+1 = y n + hf (y n , u n ),

y 0 = x,
for n ∈ N 0 , x ∈ R d , and controls u n ∈ U . Then, the application of the dynamic programming principle on the discrete-time dynamics leads to the Bellman equation

v(x) = min u∈U {(1 -λh)v(x + hf (x, u)) + hl(x, u)}, x ∈ R d .
To discretize this equation in space we introduce a bounded domain Ω = [a, b] d ⊂ R d , a, b ∈ R, where we define a regular quadrangular mesh with N nodes and mesh parameter k. We denote the set of nodes by Ω k ⊂ Ω. Let the discrete value function be defined in all grid points, V := {v(x)} x∈Ω k . However, note that x + f (x, u) for x ∈ Ω k is not necessarily a grid point, and therefore the value function has to be evaluated by interpolation which is chosen as a linear one here. The interpolant I[•](x) is defined on the basis of the dataset V . The resulting fully discrete scheme then reads

(1.3) V (x) = min u∈U {(1 -λh)I[V ](x + hf (x, u)) + hl(x, u)} =: G(V ), x ∈ Ω k ,
which can be solved by a fixed point iteration starting from an initial guess V 0 by

V i+1 = G(V i ) .
An alternative to the fixed point approach is the use of Howard's algorithm or iteration in the policy (control) space [START_REF] Bokanowski | Some convergence results for Howard's algorithm[END_REF][START_REF] Alla | An efficient policy iteration algorithm for dynamic programming equations[END_REF], which is faster but also utilizes the parametric minimization of the Hamiltonian. Finally, once the value function is computed this allows to derive a feedback control for a given state x ∈ R d by u(x) ∈ argmin u∈U {(1λh)I[V ](x + hf (x, u)) + hl(x, u)}.

A characteristic feature of the class of HJB equations arising in optimal control problem is its nonlinear Hamiltonian, sup u∈U {-f (x, u) • ∇v(x)l(x, u)} , which requires a parametric maximization (minimization) over the control set U . For its discrete analogue G(V ), a common practice in the literature is to compute the minimization by comparison, i.e., by evaluating the expression in a finite set of elements of U (see for instance [START_REF] Alla | An efficient policy iteration algorithm for dynamic programming equations[END_REF][START_REF]Numerical methods for differential games based on partial differential equations[END_REF][START_REF] Kröner | Optimal feedback control of the undamped wave equation by solving a HJB equation[END_REF] and references therein). In contrast to the comparison approach, the contribution of this paper is to demonstrate that an accurate realization of the min-operation on the right hand side of (1.3) can have an important impact on the optimal controls that are determined on the basis of the dynamic programming principle. In this respect, the reader can take a preview to Figure 4.2, where differences between optimal control fields obtained with different minimization routines can be appreciated. Previous works concerning the construction of minimization routines for this problem date back to [START_REF] Carlini | An efficient algorithm for Hamilton-Jacobi equations in high dimension[END_REF], where Brent's algorithm is proposed to solve high dimensional Hamilton-Jacobi -Bellman equations and to [START_REF] Cristiani | Fast semi-Lagrangian schemes for the Eikonal equation and applications[END_REF], where the authors consider a fast semi-Lagrangian algorithm for front propagation problems. In this latter reference, the authors determine the minimizer of a specific Hamiltonian by means of an explicit formula. Moreover, for local optimization strategies in dynamic programming we refer to [START_REF] Jarczyk | Lokale Optimierungsstrategien in der dynamischen Programmierung[END_REF] for Brent's algorithm and to [START_REF] Unrath | Nichtglatte Optimierungsstrategien in der dynamischen Programmierung[END_REF] for a Bundle Newton method.

In this article, a first-order primal-dual method (also known as Chambolle-Pock algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]) and a semismooth Newton method [START_REF] Ito | Lagrange multiplier approach to variational problems and applications[END_REF][START_REF] Ulbrich | Semismooth Newton Methods for Variational Inequalities and Constrained Optimization Problems in Function Spaces[END_REF] are proposed within the semi-Lagrangian scheme for the evaluation of the right hand side in (1.3). In contrast to the minimization by the comparison approach, the proposed algorithms leads to more accurate solutions for the same CPU time. Since we preserve the continuous nature of the control set, in some specific settings it is also possible to derive convergence results for our minimization strategies. Furthermore, it provides a solid framework to address challenging issues, such as nonsmooth optimal control problems with ℓ 1 control penalizations in the cost functional.

The paper is organized as follows. In Section 2 we begin by recasting the discretized Hamiltonian as a minimization problem explicitly depending on the control u. In Section 3 we introduce and adapt the Chambolle-Pock and semismooth Newton methods for the problem under consideration, and in Section 4 we present numerical examples assessing the accuracy and performance of the proposed schemes.

Explicitly control-dependent Hamiltonians.

In this section we study the numerical evaluation of the discretized nonlinear Hamiltonian (2.1)

min u∈U {(1 -λh)I[V ](x + hf (x, u)) + hl(x, u)} .
This is a non -standard minimization problem requiring the evaluation of the nonlinear mapping u → I[V ](x + hf (x, u)) which depends on the discrete dataset V and the system dynamics f (x, u). Therefore, a first step towards the construction of local minimization strategies is to recast (2.1) by assuming specific structures for I[V ], f (x, u), and l(x, u), leading to explicit piecewise linear or quadratic optimization problems on u. This section is split between the treatment of the interpolation I[V ](x + hf (x, u)), and the evaluation of l(x, u). For the sake of simplicity we restrict the presentation to the two dimensional case d = 2, although the presented framework can be generalized to higher dimensions.

The interpolation operator and local subdivisions of the control space.

To analyze the interpolation operator we consider for every node x = (x 1 , x 2 ) ∈ Ω k the patch of four triangles defined by the neighboring nodes, cf. Fig. 2.1. Assume that the arrival point x+hf (x, u) for x ∈ Ω k , u ∈ U , is located in a triangle defined by the points x = x 1 , x 2 , and x 3 in R 2 with associated values V 1 , V 2 , and V 3 as indicated in Fig. 2.1. The linear interpolation formula then reads for a mesh point Remark 2.1. Note that in the considered case of first order approximations the passage from the interpolation on one triangular sector to another one is continuous. For the first order interpolation schemes in a point x ∈ Ω k , one might also think of considering one interpolant for a macro-cell defined by a larger set of neighboring nodes. This could be computed by using 2 d suitable points. However, in practice this is similar to consider a grid of size 2h, which leads to less accuracy, which is, in particular, in higher dimensions a crucial issue. Alternatively one could consider one interpolant for the macro-cell which is defined piecewise over all the triangulars. However, this piecewise smooth interpolation is not convenient for numerical optimization. As a further alternative one can resort to higher order interpolation schemes.

x ∈ Ω k (2.2) I[V ](x) = cx 1 + dx 2 + e with (2.3)   x 1 1 x 1 2 1 x 2 1 x 2 2 1 x 3 1 x 3 2 1     c d e   =   V 1 V 2 V 3   . V 1 II I III IV V 2 V 3 V 4 V 5 x 1 + hf (x 1 , u)
According to (2.1), the interpolation operator I[V ](x) is evaluated at the arrival point x + hf (x, u). If this point is sufficiently close to x ∈ Ω k , it is in general possible to know a priori in which of the four triangular sectors it will be located. Our goal is to establish a consistent procedure to locally divide the control space in sets that generate arrival point related to a single triangular sector. For instance, in the case of the eikonal equation there is a simple correspondence between the computed triangulars in control and state spaces, see Section 2.3. In the following we restrict ourselves to dynamics of the form f (x, u) = g(x) + Bu, (2.4) for B ∈ R d×m , i.e. nonlinear in the state, and linear in the control. We assume that h > 0 is chosen such that

h sup x∈Ω,u∈U g(x) + Bu 2 ≤ √ 2 2 k. (2.5)
Let the rows of B be denoted by {b i } d i=1 . Further, let x D denote a grid point in Ω k , let I ⊂ {1, . . . , d} be an index set with complement I ⊂ , and define a triangle Q I ⊂ R d associated to x D for the interpolation and I by

Q I = {x = x D + x | x 1 ≤ k, x i ≥ 0 for i ∈ I, x i ≤ 0 for i ∈ I ⊂ }.
Then

U I = {u ∈ U | g(x D ) i + b i u ≥ 0 for i ∈ I, g(x D ) i + b i u ≤ 0 for i ∈ I ⊂ } has the property that x A = x D + h(g(x D ) + Bu) ∈ Q I for u ∈ U I , and 0 ≤ h ≤ h. In fact, x = x A -x D = h(g(x D ) + Bu)
, and we obtain with (2.5) that x 1 ≤ k, x i ≥ 0 for i ∈ I and x i ≤ 0 for i ∈ I ⊂ . Note also that

U = I∈W U I ,
where W is the set of all subsets of index sets in {1, . . . , d}. The associated interpolation operator on Q I is denoted by I I . Note that evaluation of I I [V ](x + f (x, u)) for x ∈ Q I and u ∈ U I leads to a linear dependence on u of the form

(2.6) c I u 1 + d I u 2 + e I ,
where the coefficients c I , d I and ǫ I are uniquely determined by the vertices of Q I .

Evaluation of the cost term.

Having approximated the interpolation term of the Hamiltonian, what is left is to provide an approximation of the running cost l(x, u). If this term is defined in a pointwise manner, this imposes no additional difficulty. For instance, in some of the examples we shall utilize

l(x, u) = x 2 2 + γ 2 u 2 2 = x 2 1 + x 2 2 + γ 2 (u 2 1 + u 2 2 ) x ∈ Ω k , u ∈ U. (2.7)
This introduces a constant and a quadratic term in u to be added to the above presented expressions for the interpolant (2.6). For l given as in (2.7) we will consider in the scheme (1.3) the explicit dependence on u. For a given node x ∈ Ω k we have

[V ] x = min I∈W {H I (x, V, u, I)} , (2.8)
where [V ] x denotes the value of the discrete value function at node x. Further, for each set I we have

H I (x, V, u, I) = min u∈U I {βI I [V ](x + hf (x, u)) + hl(x, u)} , (2.9) = min u∈U I β(c I u 1 + d I u 2 + e I ) + h x 2 + y 2 + γ 2 (u 2 1 + u 2 2 ) , = min u∈U I { ãI 2 u 2 1 + bI 2 u 2 2 + cI u 1 + dI u 2 + ẽI } , (2.10) with β = 1 -λh , ãI = hγ , bI = hγ , cI = βc I , dI = βd i , ẽI = βe I + h(x 2 + y 2 ) .
Remark 2.2. For minimum time optimal control problems, the resulting semi-Lagrangian scheme for points x ∈ Ω k reads

[V ] x = min u∈U {βI[V ](x + hf (x + u)) + 1 -β} , where β = e -h , x ∈ Ω k ,
see [START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF], and the Hamiltonian takes the simplified form

H I (x, V, u, I) = min u∈U I {c I u 1 + dI u 2 + ẽI }.
To evaluate the right hand side in (2.8) numerically, we compute (2.9) for I ∈ W by applying numerical optimization methods and then determine the minimizer on the macro-cell by comparison over the sets U I .

To determine the minimizer of the right hand side in (2.10) we define

F (u) = ãI 2 u 2 1 + bI 2 u 2 2 + cI u 1 + dI u 2 + ẽI .
Then the minimization problem is given by

min u∈U I F (u), I ∈ W. (2.11)
In particular, the following types of cost functionals are of interest for scalars a, b, c, d, e, r, l, s ∈ R, where a, c ≥ 0. Thereby we include also a bilinear term bu 1 u 2 to allow also bilinear interpolants over quadrangular cells. For the triangular interpolation used throughout this paper, we take b = 0.

1. Infinite horizon problem with quadratic control cost:

F 2 (u 1 , u 2 ) = a 2 u 2 1 + bu 1 u 2 + c 2 u 2 2 + du 1 + eu 2 + r. (2.12)
The functional is smooth and convex if acb 2 > 0. 2. Minimum time problem:

F MT (u 1 , u 2 ) = bu 1 u 2 + du 1 + eu 2 + r. (2.13)
The functional is non-convex if b = 0 but smooth. 3. Minimum time/Infinite horizon problem with ℓ 1 -control cost:

F MT/IH (u 1 , u 2 ) = a|u 1 | + bu 1 u 2 + c|u 2 | + du 1 + eu 2 + r. (2.14)
The functional is non-convex if b = 0 and non-smooth. 4. Infinite horizon problem with ℓ 2 -and ℓ 1 -control cost:

F IH (u 1 , u 2 ) = a 2 u 2 1 + l|u 1 | + bu 1 u 2 + c 2 u 2 2 + s|u 2 | + du 1 + eu 2 + r. (2.15)
The functional is convex if acb 2 > 0 and non-smooth. Remark 2.3.

Since U is compact all four functionals allow a minimum. If acb 2 > 0 it is unique for (2.12) and (2.13). For a discussion of optimal control problems of dynamical systems for functionals (2.14) and (2.15) we refer to [START_REF] Alt | Linear-quadratic control problems with l 1 -control cost[END_REF][START_REF] Maurer | On L 1 -minimization in optimal control and applications to robotics[END_REF].

2.3. Special case: eikonal dynamics. We consider the relation between the subdivision of the state and control space in the context of minimum time problems with eikonal dynamics. To set up the control problem we introduce a closed target T ⊂ R n , with int(T ) = ∅, and smooth boundary. We consider the minimum time function

t(x, ũ) = inf {t ∈ R + : y(t, ũ) ∈ T } if y(t, u(t)) ∈ T for some t, +∞ otherwise, (2.16)
where ũ ∈ U and y(•, u(•)) denotes the solution of (1.2) depending on the control. Furthermore we assume a small-time local controllability assumption as formulated in [13, p. 216]. The minimum time problem is defined as

T (x) = inf ũ∈U t(x, ũ)
and the minimum time function can be characterized as the viscosity solution of (2.17)

sup u∈U -f (x, u) T ∇v(x) = 1 in R \ T , v(x) = 0 on ∂T ,
where R are all points in the state space for which the time of arrival is finite.

We consider the dynamics

(2.18) f (x, u) = u 1 u 2 , U = {u ∈ R 2 : u 2 ≤ 1}, x ∈ R 2 ,
which leads to an equation of eikonal type. For this problem we have a direct relation between the direction of the control u and the identification of the arrival point in the state space; for instance, all the arrival points in the first quadrant will correspond to controls belonging to the set

U {1,2} = {u ∈ U | 0 ≤ u 1 , 0 ≤ u 2 }.
It is possible to express the interpolation operator (1.3) for given node x piecewise as

(2.19) I[V ](x + hf (x, u)) = I I [V ](x + hf (x, u)) with x + hf (x, u) in quadrant U I .
For every quadrant we obtain control-dependent formulas of the form (2.20)

I I [V ](x + hf (x, u)) = c I u 1 + d I u 2 + e I with coefficients c I , d I , e I ∈ R.
3. Algorithms for solving minimization problems of the form (2.11). We consider the following approaches: a) Minimization by comparison over a finite subset U finite ⊂ U I . b) First order primal-dual method: Chambolle-Pock algorithm, see [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. c) Second-order method: Two different types of semismooth Newton methods depending on the smoothness of the cost functional. d) If the functional is of type (2.13), the controls are located on the sphere or in the origin and can be found by a classical Newton method with a suitable chosen initialization over the parameterized sphere.

The minimization by comparison approach a) is broadly used in the literature. It consists in choosing a finite subset U finite ⊂ U where the cost function is evaluated and the minimum is selected among the corresponding values. Such a procedure induces a different optimization paradigm, in the sense that the continuous nature of the control space is replaced by a discrete approximation. If a parameter-dependent discretization of the control set is considered, where errors with respect to the continuous set can be estimated, the discretization has to be fine enough such that the error is negligible compared to the errors introduced in the different discretization steps of the overall scheme. Furthermore, accurate discretization of the control space are far from trivial even in simple cases, as in the three-dimensional eikonal dynamics with U = {u ∈ R 3 : u 2 ≤ 1}, where a spherical coordinate-based discretization of the control introduces a concentration of points around the poles. In contrast to this approach, we propose several numerical algorithms for the treatment of problems of the form (2.11), which preserve the continuous nature of the optimization problem at a similar computational cost.

3.1. The smooth case: Cost functionals of type (2.12) and (2.13). In this section we consider control constraints of the type U = {u ∈ R m : u 2 ≤ 1}. To solve the optimization problem for smooth cost functionals of type (2.12) and (2.13) we present a first-order primal-dual algorithm and a semismooth Newton method.

Primal-dual algorithm.

Here we assume that problem (2.11) can be reformulated as

min u∈R m F (u) + I K (u), (3.1)
where F is smooth, and convex and I K (u) is the indicator function of K = U I . This fits in the setting presented in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], where a primal-dual algorithm (also known as Chambolle-Pock algorithm) is formulated, when using the same F , and by choosing the mappings denoted by K and G in the reference by id and I K . We recall the algorithm in the following.

Algorithm 1: Chambolle-Pock algorithm

Data: Choose n = 0, τ > 0, σ > 0, η > 0, ϑ ∈ [0, 1], (u 0 , y 0 ) ∈ R m × R d , ū0 = u 0 . repeat Compute x n+1 = (y n+1 , u n+1 , ūn+1 ) by      y n+1 = (I + σ∂F * ) -1 (y n + σ ūn ), u n+1 = (I + τ ∂G) -1 (u n -τ y n+1 ), ūn+1 = u n+1 + ϑ(u n+1 -u n ). Set n = n + 1. until x n -x n-1 R d ×R m ×R m < η;
From identities from convex analysis, see [START_REF] Rockafellar | Convex analysis, Convex Analysis[END_REF], we have

(I + τ ∂G) -1 (y) = argmin u∈R m u -y 2 2 τ + I K (u) = P K (y), (3.2) u -(I + τ ∂F ) -1 (u) = τ I + 1 τ ∂F * -1 u τ , (3.3)
where P K denotes the projection on K. From (3.3) we have with σ = 1

τ (I + σ∂F * ) -1 (σu) = σu -σ I + 1 σ ∂F -1 (u),
and hence

(I + σ∂F * ) -1 (u) = u -σ I + 1 σ ∂F -1 1 σ u .
In the case I = {1, 2} the projection P K is given by

P K : R m → R m , P K (p) = max(0, p) max(1, max(0, p) 2 ) . (3.4)
This leads to Algorithm 2.

Algorithm 2: Chambolle-Pock algorithm (variant)

Data: Choose n = 0, τ > 0, σ > 0, η > 0, ϑ ∈ [0, 1], (u 0 , y 0 ) ∈ R m × R d , ū0 = u 0 . repeat Compute x n+1 = (y n+1 , u n+1 , ūn+1 ) by          y n+1 = y n + σ ūn -σ I + 1 σ ∂F -1 y n σ + ūn , u n+1 = P K (u n -τ y n+1 ), ūn+1 = u n+1 + ϑ(u n+1 -u n ). Set n = n + 1. until x n -x n-1 R d ×R m ×R m < η; If (I + 1 σ ∂F ) -1
is linear, the first step can be reformulated as

y n+1 = I + 1 σ ∂F -1 ∂F y n σ + ūn .
Remark 3.1. Since the cost functional has only quadratic and linear terms it is more convenient to use a linear interpolant rather than a bilinear one in the semi-Lagrangian scheme.This does not produce any bilinear terms. If F has a bilinear term it is very challenging to compute (I + 1 σ ∂F ) -1 in higher dimensions, since then ∂F depends nonlinearly on u.

3.2. Semismooth Newton method. The Chambolle-Pock algorithm is a first order algorithm which requires convexity of the functional. We refer to [START_REF] Valkonen | A primal-dual hybrid gradient method for non-linear operators with applications to MRI[END_REF] for an extension to a class of non-convex problems. As a second algorithm that we propose we turn to the semismooth Newton method which does not rely on global convexity. We recall some main aspects of semismooth functions, cf. [START_REF] Ulbrich | Semismooth Newton Methods for Variational Inequalities and Constrained Optimization Problems in Function Spaces[END_REF].

Definition 3.2 (Semismoothness). Let V ⊂ R µ be nonempty and open, µ ∈ N. Then function f : V → R µ is semismooth at x ∈ V if it is Lipschitz continuous near x and if the following limit exists for all s ∈ R ν , ν ∈ N: lim M ∈∂f (x+τ d) d→s,τ →0 + M d.
Furthermore, there holds the following chain rule.

Lemma 3.3 (Chain rule). Let V ⊂ R µ and W ⊂ R ν be nonempty open sets, g : V → W be semismooth at x ∈ V , and h : W → R η , η ∈ N, be semismooth at g(x) with g(V ) ⊂ W . Then the composite map

f := h • g : V → R η is semismooth at x. Moreover, f ′ (x, •) = h ′ (g(x), g ′ (x, •)).
To set up a semismooth Newton algorithm for (2.11) we proceed as follows. For simplicity, we focus on the case I = {1, 2}. The first-order optimality condition for cost functionals of type (2.12) and (2.13) can be formulated as (3.4). We introduce p = y -ϑ∇F (y) and choose ϑ such that ϑ∇F (u) is of the same scale as y and rewrite (3.5) as

u = P K (u -ϑ∇F (u)) ∀ϑ > 0 (3.5) for u ∈ R m , F = F 2 or F = F MT with projection P K as in
u -P K (p) = 0, u -ϑ∇F (u) -p = 0. Setting β = max(1, max(0, p) 2 ) we define E(u, p, β) =   βu -max(0, p) u -ϑ∇F (u) -p β -max(1, max(0, p) 2 )   .
Now the optimality condition can be formulated as E(z) = 0 with z = (u, p, β). Then we set up a semismooth Newton method as presented in Algorithm 3.

Algorithm 3: Semismooth Newton algorithm

Data: Choose n = 0, initialization z 0 ∈ U I × R m × R, η > 0. repeat Solve the Newton equation for δz ∈ U I × R m × R given by     β n I D χ pn ≥0 y n I -ϑ∇ 2 F (u n ) -I 0 0 χ mn≥1 -p T n Dχ pn ≥0 mn 1     δz = -E(z n ) (3.6)
(with matrix D χ pn ≥0 = diag(χ pn≥0 ) and m n = max(0, p n ) 2 ).

Update

z n+1 = z n + δz. (3.7) Set n = n + 1. until z n -z n-1 R m ×R m ×R < η;
We address solvability of (3.6). Lemma 3.4. Set M = I -ϑ∇ 2 F (u n ) and q + = max(0, q) for q ∈ R m . The Newton matrix is regular if one of the following two conditions is satisfied:

1. m n < 1:

β n = 0, 1 βn M D χp n ≥0 + I regular. 2. m n ≥ 1: β n = 0, I - 1 β n M regular (3.8) and -1 mn (p + ) T (I -1 βn M + ) -1 -M + u + n βn + z + 4 )
+ + 1 = 0, where the components of p = (p + , p -) are ordered with respect to active and inactive sets and M + denotes the submatrix corresponding to p + . Remark 3.5. For β n ≈ 1 condition (3.8) is closely related to the regularity of the Hessian of F . Moreover, if we set K = sup u∈U {1,2} ∇ 2 F (u) and assume that β n are bounded away from 1 (i.e. β n ≥ κ > 1) then the choice ϑ < κ-1 K implies (3.8). Proof of Lemma 3.4. The two cases are considered separately. For m n < 1, the Newton-Matrix is given by

DE(z n ) =   β n I D χ pn ≥0 u n M -I 0 0 0 1   .
The regularity of the matrix follows from the regularity of

β n I D χ pn ≥0 M -I and hence of 1 βn M D χp n ≥0 + I .
For m n ≥ 1 the Newton-matrix is given by

DE(z n ) =   β n I D χ pn ≥0 u n M -I 0 0 -(p + ) T mn 1  
and leads to the Newton equations

β n δu -δp + + δβu n = z 1 , (3.9) M δu -δp = z 2 , (3.10) - 1 m n (p + ) T δp + + δβ = z 3 (3.11) for suitable z 1 , z 2 ∈ R m , z 3 ∈ R independent of (δu, δp, δβ).
From the first and second equations we obtain

- 1 β n M (δp + -δβu n ) + δp = z 4
for z 4 ∈ R m . Without loss of generality we assume that the components of δp = (δp + , δp -) are ordered with respect to active and inactive sets. Then we have

(I - 1 β n M + )δp + + δβM + u + n β n = z 4 (3.12)
and therefore

δp + = δβ(I - 1 β n M + ) -1 -M + u + n β n + z + 4 . (3.13) With (3.11) we obtain - 1 m n (p + ) T (I - 1 β n M + ) -1 -M + u + n β n + z + 4 ) + + 1 δβ = β. (3.14)
Theorem 3.6. Let ū be a strict local minimizer of (2.11). Under the assumptions of Lemma 3.4 the semismooth Newton method converges locally superlinearly or terminates after a finite number of steps at ū.

Proof. The operator E is semismooth and is as a composition of Lipschitz continuous functions again Lipschitz continuous. Furthermore from Lemma 3.4 we obtain the boundedness of the inverse derivative in a neighbourhood of ū.

Consequently, the assertion follows from [21, p. 29] and [15, p. 220, Thm. 8.3].

Approach for cost functionals of type (2.13).

To treat the problem (2.11) for cost functionals of type (2.13) we present an alternative approach. We show that all possible minimizers are located on the sphere or in the origin; cf. also the discussion in [START_REF] Cristiani | Fast semi-Lagrangian schemes for the Eikonal equation and applications[END_REF] where a linear functional is considered. To minimize over the sphere we parameterize the sphere by polar coordinates and consider the restriction of the functional on the sphere. There the minimizer can be found by applying a classical Newton method.

Lemma 3.7. For control problems (2.11) with cost functional (2.13) for r = 0 all minimizers are located either on the sphere or in the origin.

Proof. For simplicity we consider the case (2.11) for d = 2 and I = {1, 2}. Let

F (u) = du 1 + bu 1 u 2 + eu 2 , u ∈ R 2 , d, b, e ∈ R. Then we have ∇F (u) = d + bu 2 e + bu 1 , ∇ 2 F (u) = 0 b b 0 .
This implies, we have for b = 0 a saddle point with eigenvalues ±1 in

u 1 = - e b , u 2 = - d b .
For b = 0 the minimum is reached in u = (0, 0) T . Consequently, all minimizers are on the boundary of ∂U {1,2} . However, since the restriction of the bilinear functional to the axes is linear the assertion follows immediately.

In two dimensions the functional F restricted to the sphere is given by

F (ϕ) = a cos(ϕ) + b 2 sin(2ϕ) + c sin(ϕ), 0 ≤ ϕ ≤ π 2 .
3.4. The non-smooth case: Cost functionals of type (2.14) and (2.15). In this section we consider two types of control constraints, box constraints as well as Euclidean constraints. Finally, we also present a splitting approach.

Semismooth Newton in case of Euclidean norm constraints.

In this section we restrict the consideration to equations of eikonal type and present the approach for a two dimensional problem. Therefore the control set is given by U = {u ∈ R 2 : u 2 ≤ 1}. Let us consider the equivalent problem formulation min u∈R 2 E(u) + h(u), (3.15) where E : R 2 → R is smooth and

(3.16) h(u) = α(|u 1 | + |u 2 |), if u ∈ U I , ∞, if u ∈ U I .
The optimality condition for the nonsmooth problem (3.15) is given by

0 ∈ ∇E(u) + ∂h(u) ,
where ∂h denotes the subdifferential of the convex function h. Setting q = -∇E(u) the optimality condition can be equivalently written as (3.17) q = -∇E(u), q ∈ ∂h(u).

With the convex conjugate we obtain equivalently

(3.18) q = -∇E(u), u ∈ ∂h * (q) , with h * (q) = sup v∈U I (v • q -h(v)).

Note that on U

I v • q -h(v) = (q 1 -α 1 )v 1 + (q 2 -α 2 )v 2 with α 1 = sgn(L 1 )α, α 2 = sgn(L 2 )α (3.19)
and

L 1 = 1, if I = {1, 2} or I = {1}, -1, if I = {2} or I = ∅, , L 2 = 1, if I = {1, 2} or I = {2}, -1, if I = {1} or I = ∅. (3.20)
For fixed q and v ∈ U I we define

l(v) = (q 1 -α 1 )v 1 + (q 2 -α 2 )v 2 .
The sup of this linear function is attained on ∂U I . We distinguish four cases

sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = 0, sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = q 1 -α, sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = q 2 -α, sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = sup ϕ∈Λ (q 1 -α 1 ) cos ϕ + (q 2 -α 2 ) sin ϕ , with Λ =        0 ≤ ϕ ≤ π 2 if I = {1, 2}, π 2 < ϕ ≤ π if I = {2}, π < ϕ ≤ 3π 2 if I = ∅, 3π 2 < ϕ ≤ 2π if I = {1}. (3.21)
In the following for simplicity we only consider the problem for I = {1, 2}. From the first order condition we obtain the following relation between q and ϕ

tan ϕ = q 2 -α q 1 -α > 0. (3.22)
Summarizing we have

(3.23) h * (q) =          (q 1 -α) cos ϕ + (q 2 -α) sin ϕ, if q 1 , q 2 ≥ α, q 2 -α, if q 1 < α, q 2 > α, 0, if q 1 , q 2 < α, q 1 -α, if q 1 > α, q 2 < α ,
with ϕ given by (3.22) and for the generalized derivative

(3.24) ∂h * (q) =            (w 1 (q), w 2 (q)) T , if q 1 , q 2 > α, (0, 1) T , if q 1 < α, q 2 > α, (0, 0) T , if q 1 < α, q 2 < α,
(1, 0) T , if q 1 > α, q 2 < α where (w 1 (q), w 2 (q)) = ∇((q 1α) cos ϕ + (q 2α) sin ϕ) .

To set up a semismooth Newton method we introduce a piecewise affine, continuous regularization of the gradient in a ε-tube around the discontinuity and define

(3.25) (∂h * ) ε (q) =                                        (w 1 (q), w 2 (q)) T , if q 1 , q 2 > α + ε, (0, 1) T , if q 1 < α, q 2 > α + ε, 0, q 2 -α ε T , if q 1 < α, α < q 2 < α + ε, (0, 0) T , if q 1 < α, q 2 < α, q 1 -α ε , 0 T , if α < q 1 < α + ε, q 2 < α, (1, 0) T , if q 1 > α + ε, q 2 < α, r ε
(w 1 (q p ), w 2 (q p )) T , if α < q 1 , α < q 2 , and q -(α, α) T 2 ≤ ε, where r = q -(α, α) 2 , θ = tan -1 q 2α q 1α , and q p = ε(cos(θ), sin(θ)) + (α, α).

The semismooth Newton method is presented in Algorithm 4.

Algorithm 4: Semismooth Newton algorithm

Data: Choose n = 0, initialization z 0 = (q 0 , u 0 ) ∈ R 2 × U I , η > 0. repeat Solve the Newton equation for δz ∈ R 2 × U I given by

I ∇ 2 E(u n ) -D(∂h * ) ε (q n ) I δz = q n + ∇E(u n ) -(∂h * ) ε (q n ) + u n , (3.26)
where D(•) denotes the Newton derivative. Update

z n+1 = z n + δz. (3.27) Set n = n + 1. until z n -z n-1 R 2 ×R 2 < η;
We can conclude the following theorem. Theorem 3.8. Let u ∈ R 2 be a strict local minimum of (3.15). If

I -∇ 2 E(u n )D(∂h * ) ε (q n ) (3.28)
is regular for all k, the semismooth Newton converges locally superlinearly.

Proof. We introduce the piecewise affine function

G : R 2 × R 2 → R 2 × R 2 , G(q, u) = (q n + ∇E(u), -(∂h * ) ε (q) + u n ).
G is semismooth and together with condition (3.28) we get directly the regularity of the Hessian in the Newton equation given in (3.26). From both conditions we derive locally superlinear convergence as in Theorem 3.6.

3.4.2.

Semismooth Newton in case of box constraints. Now, let the control set be given by

U = {u ∈ R m : u a ≤ u ≤ u b } with u a ≤ 0 ≤ u b , u a , u b ∈ R 2 .
As in the previous section, we consider a cost functional of the type min u∈R 2 E(u) + h(u), (3.29) where E is smooth and

(3.30) h(u) = α(|u 1 | + |u 2 |), if u ∈ U I , ∞, if u ∈ U I .
Proceeding as in the previous section we distinguish four cases

sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = 0, sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = |q 1 -α 1 |u a , sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = |q 2 -α 2 |u b , sgn(q 1 -α 1 ) = sgn(L 1 ) ∧ sgn(q 2 -α 2 ) = sgn(L 2 ) : h * L1,L2 (q) = |q 1 -α 1 |u a + |q 2 -α 2 |u b ,
with α 1 , α 2 as in (3.19) and L 1 , L 2 as in (3.20). Further, there holds

(3.31) ∇h * (q) =            (sgn(q 1 ) u a , sgn(q 2 ) u b ) T , if q 1 > α, q 2 > α, (0, sgn(q 2 ) u b ) T , if q 1 < α, q 2 > α, (0, 0) T , if q 1 < α, q 2 < α,
(sgn(q 1 ) u a , 0) T , if q 1 > α, q 2 < α , for I = {1, 2} and accordingly for the other cases. To set up a semismooth Newton method we regularize the first component of the gradient with ramps of width ε at |q 1 | = α and the second component with ramps at |q 2 | = α. The overall algorithm has the same form as the semismooth Newton method in (3.26). By a similar consideration as in the previous section Theorem 3.8 holds also for the semismooth Newton method described above.

A splitting approach.

To treat the cost functionals (2.14) and (2.15) with ℓ 1 -terms we can reformulate the problem without non-differentiable terms by doubling the number of variables. This is illustrated for the two dimensional case with two dimensional control (u, v) ∈ U I . For a scalar z ∈ R we define z + = max(0, z), z -= min(0, z) (in particular z + z -= 0 and |z| = z +z -). Then problem (2.11) with cost functional (2.15) is given by

                           min u+,u-,v+,v- au 2 + + au 2 -+ cv 2 + + cv 2 -+ bu + v + -bu -v + -bu + v -+ bu -v - + l(u + -u -) + s(v + -v -) + d(u + + u -) + e(v + + v -) + r, s.t. g i + b i u + -u - v + -v - ≥ 0, u + u -= 0, v + v -= 0, u + ≥ 0, v + ≥ 0, u -≥ 0, v -≥ 0, u 2 2 + v 2 2 ≤ 1, u, v ∈ U I
with g chosen as in (2.4) and b i as the columns of B. This problem formulation can be simplified, by distinguishing four cases depending on I; for example, for U {1,2} we obtain the equivalent problem

(3.32)                  min u+,v+ au 2 + + cv 2 + + bu + v + + (l + d)u + + (s + e)v + + r, s.t. g i + b i u + v + ≥ 0, i = 1, 2, u + ≥ 0, v + ≥ 0, u 2 + + v 2 + ≤ 1.
Problem (3.32) can be solved by using algorithms of Section 3.1.

Numerical examples.

In this section we present a set of numerical tests aiming at studying the performance and accuracy of the algorithms presented in the previous sections. We begin by assessing the performance of the proposed numerical optimization routines as a separate building block.

U = {u ∈ R 2 | u 2 ≤ 1} or U = {u ∈ R 2 | 0 ≤ u 1 ≤ 1, 0 ≤ u 2 ≤ 1},

Algorithm

Tolerance Iterations CPU time

• 2 error Semismooth Newton 1E-4 58 1.47E-2 [s] 1.23E-3 Comparison (2E3 evaluations) - - 8.56E-4 [s] 4.18E-2 Comparison (1E5 evaluations) - - 1.20E-2 [s] 1.47E-2

4.2.

Interplay with the semi-Lagrangian scheme. Having embedded the minimization routines within a semi-Lagrangian scheme, we show in Fig. 4.1 the evolution of the average iteration count per gridpoint (at every fixed point iteration of the semi-Lagrangrian scheme), for both a minimum time and an infinite horizon optimal control problem subject to eikonal dynamics. Note the difference in the evolution of the subiteration count depending on the considered control problem. In the minimum time problem nodes that have not received information propagating from the optimality front are still minimized with irrelevant information until they are reached by the optimality front, whereas in the infinite horizon case "correct" information is available for every node from the first iteration due to the presence of a running cost. In this latter case, the impact of the available information from the previous semi-Lagrangian iteration is similar to a warm start of the optimization routine.

(u, x) = 1 2 x 2 2 + γ 2 2 u 2 2 , γ 2 > 0, λ > 0 .
As common setting, the fixed point iteration is solved until

V n+1 -V n ≤ 1 5 k 2 , n = 1, 2, 3, . . . ,
where k stands for the space discretization parameter and the stopping tolerance for the inner optimization routine is set to 10 -4 . Further common parameters are λ = 0.1 and γ 2 = 2, and specific settings for every problem can be found in Table 4.4.

Test Ω U h Test 1 [-1, 1] 2 u 2 ≤ 1 √ 2 4 k Test 2 [-1, 1] 3 u 2 ≤ 1 1 2 k Test 3 [-1, 1] 3 u ∞ ≤ 0.3 1 5 k
Table 4.4: Parameters for Tests 1-3.

Test 1: 2D eikonal dynamics. In this test we consider eikonal dynamics of the form

(4.2) f (x, u) = (u 1 , u 2 ) T , (u 1 , u 2 ) 2 ≤ 1 .
We study the accuracy and performance of the semi-Lagrangian scheme with different minimization routines: discretization of the control set and minimization by comparison, a semismooth Newton method given in Algorithm 3, and the approach by Chambolle and Pock given in Algorithm 2. In order to make a fair study of the different routines, we choose the discrete set of controls for the comparison algorithm such that the CPU time for the semismooth Newton method and the Chambolle-Pock algorithm is almost the same. Table 4.5 shows errors in both the value function and in the optimal control between the exact solutions v and u and their numerical approximations V h and U h for the different schemes. Errors are computed with respect to the exact solution of the Hamilton-Jacobi-Bellman equation, which is not readily available in the literature. Since it is useful for numerical investigations, it is provided in Appendix B. The results show that we have similar CPU time of the approaches and independently of the meshsize, the optimization-based schemes yield more accurate approximations of the value function and the associated optimal control than the approach based on comparison. Further results are shown in Fig. 4.2, where we also consider nonhomogeneous eikonal dynamics

(4.3) f (x, u) = (1 + χ x2>0.5 (x))(u 1 , u 2 ) T , (u 1 , u 2 ) 2 ≤ 1
with χ x2>0.5 (x) corresponds to the indicator function of the set {x = (x 1 , x 2 ) |x 2 > 0.5}. The figure shows that both approaches, the SL-scheme with a semismooth inner optimization block and the one with a comparison-based routine, lead to very similar value functions. This is clearly not the case for the optimal control fields depicted in rows 2 and 3 of Fig. 4.2. Even by a post-processing step it would be difficult to obtain the results in the third row from those in the second row. Test 2: 3D eikonal dynamics. In order to illustrate that our approach can be extended to higher dimensions, we consider the infinite horizon optimal control of three dimensional eikonal dynamics Errors and CPU time are shown in Table 4.6. Since the values for the semismooth Newton and the Chambolle-Pock algorithms are similar, we only include the data for this latter one. With the Chambolle-Pock algorithm we obtain more accurate solutions of the control field for a similar amount of CPU time as the comparison approach.

Algorithm CPU time v -V h 1 u -U h 1 CPU time v -V h 1 u -U h
f (x, u) = (u 1 , u 2 , u 3 ) T , (u 1 , u 2 , u 3 ) 2 ≤ 1 .
Test 3: Triple integrator with two control variables. In this test the dynamics are given by a triple integrator with two control variables subject to box constraints The purpose of this example is to stress that the minimization strategy that we have introduced can be also applied to non-eikonal dynamics, where the correspondence between octants in the state space and the control field is not trivial. For the sake of completeness the control space decomposition procedure for this example can be found in Appendix A. In this particular case, every octant will have associated a different rectangular sector in the control space for its arrival points. Results for the value function, optimal controls and trajectories are shown in Fig. 4.3. In the second row of this figure, we observe distinct differences between the optimal controls obtained from the Chambolle-Pock and comparison-based algorithms. These differences in the control lead to different approximated steady states as it can be seen in the left of the second row of Fig. 4.3. To highlight also the effect of closed-loop control we carried out an experiment with additive structural and output noise. In the third row of Fig. 4.3, the resulting states from open and closed-loop control can be compared.

f (x, u) = (x 2 , x 3 + u 1 , u 2 ) T , |u 1 | ≤ a, |u 2 | ≤ b, a, b ∈ R. k = 0.1, (20 3 DoF) k = 0.05, (40 3 DoF) Algorithm CPU time v -V h 1 u -U h 1 CPU time v -V h 1 u -U h

4.4.

Infinite horizon problems with combined ℓ 2 and ℓ 1 -costs. We present two numerical tests for infinite horizon optimal control problems with cost function and running cost given by

J(u, x) = ∞ 0 l(x(s), u(s))e -λs ds , l(u, x) = 1 2 x 2 2 + γ 2 2 u 2 2 + γ 1 u 1 , γ 2 > 0, γ 1 > 0, λ > 0 .
The stopping rule for the fixed point iteration is defined in the same way as in the previous set of examples. All the optimization-based routines have been solved with a semismooth Newton method as the inner block, with a regularization parameter ε = 1e -3. Further common parameters are λ = 0.1 and γ 2 = 2, and specific settings for every problem can be found in Table 4.7.

Test Open and closed loop solutions under state perturbation Test 4: 2D eikonal dynamics. This test considers the same two dimensional dynamics as in Test 1, the only difference being the inclusion of an ℓ 1 -term in the cost functional. For the inner optimization block we apply the semismooth Newton method presented in Algorithm 4. Results are shown in Fig. 4.4, where differences in the shape of the value function can be observed as the γ 1 weight increases. In the second row of Fig. 4.4, the effect of sparsity on the first control component can be seen from the fact that it is identically zero in a band around the origin. Moreover, the width of the band increases with γ 1 . As for the value function, introducing the ℓ 1 -term breaks its radially symmetric structure, as it is shown in the first row of Fig. 4.4.

Ω U h k Test 4 [-1, 1] 2 u 2 ≤ 1 √ 2 4 k 0.025 Test 5 [0, 2π] 3 u ∞ ≤ 0.3 0.2k 0.2
Test 5: 3D car model. In our last test, the dynamics are given by a nonlinear 3D car model with two control variables: costs. The addition of an ℓ 1 -term shrinks the region of a given isovalue. The second and third rows depict the effect of sparsity in both control variables, creating regions of zero action. A direct consequence of this can be seen at the bottom of Fig. 4.5, where the inclusion of the additional ℓ 1 -cost generates optimal trajectories which do not reach the origin due to the interplay between the control cost and the discount factor of the optimal control problem.

f (x, u) = (u 1 cos(x 3 ), u 1 sin(x 3 ), u 2 ) with u 1 ∈ [-ω 1 , ω 1 ], and u 2 ∈ [-ω 2 , ω 2 ], ω 1 > 0, ω 2 > 0.
Concluding remarks. We have presented a numerical approach for the solution of HJB equations based on a semi-Lagrangian discretization and the use of different local minimization strategies for the approximation of the corresponding numerical Hamiltonian. The numerical results show a more accurate resolution of the optimal control field at a similar computational cost as the currently used schemes. Furthermore, the proposed approach can be also adapted to treat non-differential costs such as ℓ 1 -penalizations on the control. Since the numerical approximation of the Hamiltonian constitutes one building block within the construction of approximation schemes for HJB and related equations, the idea of using local minimization techniques can be conveniently recast in similar problems, such as front propagation problems and differential games, and in related approximation techniques, like fast marching schemes, policy iteration methods and high-order approximations.

Appendix A: Decomposition of the control space for Test 3. We make an explicit presentation of the decomposition of the control space of Test 3. The where U = {u ∈ R 2 | u 2 ≤ 1}, for which we want to obtain an explicit solution. If The ansatz v 1 = Ar 2 , leads to λAr 2 + 2 γ A 2 r 2 -r 2 2 = 0 and hence A = γ 4 ( λ 2 + 4 γ -λ). The resulting expression for v 1 is the solution of (4.4), if (v 1 ) r ≤ γ which results in r ≤ 2((λ 2 + 4 γ )

1 2λ) -1 =: r. We note that v 1 (r) = γ 2 r. We turn to the case that 1 γ ∇v 2 > 1. In this case the maximum in (4.4) is achieved on the boundary of U at u * = -∇v ∇v 2 . Again we look for a radially symmetric solution in polar coordinates Note that this latter coupling yields that v(r) is a C 1 (R) function and therefore it is also a classical solution of eq. (4.4). Summarizing we have

v(r) =        γ 4 ( λ 2 + 4 γ -λ)r 2 for r ≤ r 1 2λ r 2 - 1 λ 2 r + γ 2λ + 1 λ 3 + de -λr for r > r,
where d is given in (4.6). Rows 2 and 3: effect of sparsity on the optimal control field. Bottom: optimal trajectories with and without sparsity.

Fig. 2 . 1 :

 21 Fig. 2.1: Arrival points and related control sets.

4. 1 . 2 +

 12 Preliminary tests. We consider a generic two-dimensional minimization problem of L • u + γ u 1 subject to Euclidean norm or box constraints, i.e.

Fig. 4 . 1 :

 41 Fig. 4.1: Subiteration count for 2D control problems with eikonal dynamics

k

  = 0.05, (40 2 DoF) k = 0.025, (80 2 DoF)

1

  Optimal control u1 via semismooth algorithm

Fig. 4 . 2 :

 42 Fig. 4.2: Infinite horizon control with 2D eikonal dynamics. Left: continuous dynamics. Right: discontinuous dynamics.

  Effect on the state trajectories of the minimization routine x 1 (t) Chambolle-Pock algorithm x 2 (t) Chambolle-Pock algorithm x 3 (t) Chambolle-Pock algorithm x 1 (t) Comparison algorithm x 2 (t) Comparison algorithm x 3 (t) Comparison algorithm Effect on the input signal of the minimization routine u 1 (t) Chambolle-Pock algorithm u 2 (t) Chambolle-Pock algorithm u 1 (t) Comparison algorithm u 2 (t) Comparison algorithm

Fig. 4 . 3 :

 43 Fig. 4.3:Triple integrator with two controls, numerical results with k = 0.05. Top: different isosurfaces. Middle: differences in the optimal control lead to different trajectories. Bottom: the feedback approach leads to robustness with respect to noise in the dynamics.

Fig. 4 . 4 :

 44 Fig. 4.4: Sparse control of eikonal dynamics. Top: inclusion of an ℓ1-cost breaks the radially symmetric structure of the solution. Bottom: the sparsity of the optimal control translates into a zero band around of the origin, depending on the weight γ1.

1 γ

 1 ∇v 2 < 1, then u * = -1 γ ∇v provides a maximum for the expression in brackets, and the HJB-equation becomes coordinates (r, ϕ) this equation can be expressed asλv(r, ϕ) + 1 2γ (v r (r, ϕ) 2 + 1 r 2 v ϕ (r, ϕ) 2 ) -

(4. 5 ) 2 = γ 2 . 2 ( 3 + 4

 522234 λv(r) + |v r (r)| -r 2Assuming that v r ≥ 0 we make the Ansatz v 2 (r) = ar 2 + br + c + de -λr . Inserting into (4.5) and comparing coefficients we obtain v de -λr .Continuous concatenation with v 1 at r implies that (

44551Fig. 4 . 5 :

 45 Fig. 4.5: Sparse control of a 3D car model. Top: different isosurfaces with different ℓ1-cost.Rows 2 and 3: effect of sparsity on the optimal control field. Bottom: optimal trajectories with and without sparsity.

Table 4 .

 4 respectively. Results presented in Tables 4.1, 4.2 and 4.3 show the different performance scenarios found under different costs and constraints. For every setting we can observe that the minimization by the comparison algorithm (see Section 3) requires very fine discretizations of the control variable (and thus, higher CPU time) to achieve similar error levels as the optimization-based counterpart.

	Algorithm	Tolerance Iterations CPU time	• 2 error
	Chambolle-Pock	1E-4	9	7.3E-5 [s]	4.31E-5
	Semismooth Newton	1E-4	5	1.3E-2 [s]	7.74E-9
	Comparison (1E4 evaluations)	-	-	1.1E-3 [s]	1.6E-2
	Comparison (2E3 evaluations)	-	-	6.6E-4 [s]	4.1E-2
	Algorithm	Tolerance Iterations CPU time	• 2 error
	Semismooth Newton	1E-4	101	4.53E-3 [s] 1.51E-3
	Comparison (4E4 evaluations)	-	-	5.18E-3 [s] 5.77E-3
	Comparison (2E5 evaluations)	-	-	2.08E-2 [s] 2.80E-3

1: Performance tests for an ℓ2-cost (i.e. γ = 0) with Euclidean norm constraint.

Table 4 . 2

 42

: Combined ℓ2-and ℓ1-cost (γ = 0.1) with box constraint.

Table 4 . 3 :

 43 Combined ℓ2-and ℓ1-cost (γ = 0.1) with Euclidean norm constraint.

Table 4 . 5 :

 45 Infinite horizon control of 2D eikonal dynamics. CPU time and errors for a semi-Lagrangian scheme with different minimization routines. The comparison algorithm was run with a discrete set of 1280 control points in every node.

	1

Table 4 . 6 :

 46 Infinite horizon control of 3D eikonal dynamics. CPU time and errors for a semi-Lagrangian scheme with different minimization routines. The comparison algorithm was run with a discrete set of 5120 control points.

	1

Table 4 . 7 :

 47 Parameters for Tests 4-5.

dynamics is given by

For a given departure point x d = (x d 1 , x d 2 , x d 3 ) ∈ R 3 and a sufficiently small h, we want to identify a relation between subsets U I of the control space U = [-a, a] × [-b, b] and the location of the corresponding arrival points

In the three-dimensional case, the control set U can be decomposed into at most eight disjoint subsets, one per octant with

for W = {1, 2, 3}. Since every octant defines a unique linear interpolant, we solve eight different minimization problems, and then compute the global nodal minimizer by comparison. For instance, let us consider the sector

where the evaluation of the arrival point is defined by the linear interpolant I {1,2,3} depending on the points x d , x d + (k, 0, 0), x d + (0, k, 0), and x d + (0, 0, k). The subset U {1,2,3} related to this interpolant reduces to

Note that for this definition to be meaningful, we need to assume that x d 2 ≥ 0, otherwise U {1,2,3} = ∅. Also note that the condition (x 1 -x d 1 )+(x 2 -x d 2 )+(x 3 -x d 3 ) ≤ k is omitted by assuming a sufficiently small h ensuring it. The next subset U {2,3} relates to the sector

and is given by

Note that U {1,2,3} ≡ U {2,3} , however, U {2,3} is nonempty only when x d 2 ≤ 0. Therefore, for a given departure point, depending on the coordinate x d 2 , only one of these two subsets will be active with arrival points in different sectors Q {1,2,3} or Q {2,3} (with different interpolation data). In a similar way, the remaining six subsets can be obtained. By assuming a linear control term, the identification of the control subsets is simple, as it will only require the resolution of linear inequalities where the departure point enters as a fixed data.


François Guillaume
email: f.guillaume@ism.u-bordeaux1.fr

Yeny A Tobon

Sébastien Bonhommeau

Jean-François Létard

Lucie Moulet

Eric Freysz

Photoswitching of the spin crossover polymeric material [Fe(Htrz)2(trz)](BF4) under continuous laser irradiation in a Raman scattering experiment

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Transition metal complexes in which the spin state of the metal changes under external stimuli such as temperature, light irradiation or pressure have been extensively studied since a few decades due to their high relevance for potential applications in information processing, sensors and/or display devices [1][2][START_REF]Spin Crossover in Transition Metal Compounds I[END_REF], and even more recently in the photovoltaic domain [START_REF] Etrillard | [END_REF]. The spin crossover (SCO) phenomenon is a transition from a low-spin (LS) ground state electron configuration to a high-spin (HS) metastable state of the metal d atomic orbitals. Future applications of SCO materials require such bistability over a broad temperature range centered at $300 K, as in polymeric SCO compounds [START_REF]Spin Crossover in Transition Metal Compounds I[END_REF][5][6][7][START_REF] Letard | Spin Crossover Transition Metal Compound III[END_REF][START_REF] Aromí | [END_REF] where the highly cooperative spin transition process leads to a large thermal hysteresis loop. Some of the most promising polymeric SCO compounds are formed with the triazole ligands threefold bridging the Fe(II) metal ions (Figure 1). Within this family of SCO compounds, [Fe(Htrz) 2 (trz)](BF 4 ) (with Htrz = 1,2,4-4H-triazole and trz = deprotonated triazolato (À) ligand) has been investigated quite extensively over the last years [START_REF] Etrillard | [END_REF]6,[10][11][12][13][14] in particular because it displays a $40 K wide thermal hysteresis centered around 370 K.

The present work aims to investigate the effect of the laser intensity of a continuous laser on the hysteresis loops of the [Fe(Htrz) 2 (trz)](BF 4 ) compound. It is known that intense nanosecond laser light pulses applied within a thermal hysteresis loop can trigger in SCO materials the LS to HS transition and in some cases the reverse HS to LS transition, [15][16][17][18][19][20] provided an energy density threshold value is overcome [20,21]. However, the respective contribution of photo-induced and thermally-driven effects in this laser-induced process is still a matter of debate, and might depend on the material. In the case of a purely thermal origin at least, it might be possible to mimic somewhat the observed effects using a continuous laser irradiation.

Hereafter the spin state of the [Fe(Htrz) 2 (trz)](BF 4 ) complex will be probed by means of Raman microscopy, where the Raman probe laser will also provide the light excitation allowing thermal effects to be potentially generated. Since several years Raman spectroscopy is a well-recognized technique to characterize the phase transition process in SCO compounds, and especially their thermal hysteresis loops [12,14,[22][START_REF] Tuchagues | Spin Crossover Transition Metal Compound III[END_REF][START_REF] Wolny | [END_REF][25][26][27][28][29][30][31]. The spin transition phenomenon in [Fe(Htrz) 2 (trz)](BF 4 ) samples has been analyzed recently by Urakawa et al. [12] by means of an experimental approach combining Raman spectroscopy and X-ray powder diffraction measurements (XRPD). In this Letter, hysteresis loops measured with these different techniques were found to differ significantly. According to the Raman data, the spin transition was completed within 320-340 K with almost no hysteresis, while the XRPD data showed ⇑ Corresponding author. a $40 K-wide hysteresis. This result was interpreted by Urakawa et al. [12] as due to the different sensitivities of the techniques towards short-and long-range ordering of the crystal structure or to the distinct penetration depths of the external probes. However, these experimental results may also suggest, as pointed out by Wolny et al. [START_REF] Wolny | [END_REF] in a recent review, that 'heating effects due to the Raman laser occurred'. Indeed, we found in a previous study on another material belonging also to this family of SCO compounds [22] that, even at room temperature, outside the hysteresis loop, the laser could induce the LS to HS transition. The impact of the laser intensity in a Raman experiment on the spin state of a SCO material within the hysteresis loop has also been addressed by Craig et al. [32] very recently. Obviously, the laser light interacts with the sample in a conventional Raman experiment and the main objective of this Letter is to get some insight on the impact of the laser excitation on the Raman measurements.

Experimental

The powder samples of [Fe(Htrz) 2 (trz)](BF 4 ) were synthesized according to the method already published [6]. Raman spectra were collected in backscattering geometry using excitation lines centered at 633 nm and 532 nm and a LabRam HR-Evolution spectrometer (Horiba Jobin-Yvon) with 2 cm À1 spectral resolution. A few test measurements with the 785 nm excitation line were also performed in a backscattering configuration using an Xplora Raman spectrometer (Horiba Jobin-Yvon). These spectrometers were coupled to an Olympus microscope with an ultra long working distance 50Â MSPlan objective (numerical aperture NA = 0.55). A Linkam THS 600 stage positioned under the microscope objective was used to perform the heating and cooling cycles with a rate of 2 K per minute. Once the set temperature was reached, the sample was maintained at this temperature for 5 min for thermal equilibration before exposing it to the laser beam. The sample was initially warmed at 400 K in an atmosphere of dry N 2 to eliminate water. The sample was subsequently held in dry N 2 at atmospheric pressure during all experiments. The laser power at the output of the microscope objective was systematically measured prior to measurements and could be tuned using density filters that attenuate the laser intensity. It should be pointed out that the laser intensity (or power density defined as the ratio power/irradiated sample surface) will obviously depend on the specifications of the microscope objective. In our experiments, the laser spot diameter defined theoretically [33] as 1.22k/NA was approximately 1.4 lm at k = 633 nm and 1.2 lm at k = 532 nm using the 50Â objective with NA = 0.55.

Results and discussion

We have reported in Table 1 the main parameters characterizing the transition temperatures and the width of the hysteresis loop of [Fe(Htrz) 2 (trz)](BF 4 ) taken from literature data [6,10,11,18]. The LS to HS transition temperature T up obtained upon sample warming ranges from 379 to 386 K, and the HS to LS transition temperature T down obtained upon sample cooling ranges from 341 to 347 K. The width of the hysteresis loop is about 40 K. Clearly, the values determined by Urakawa et al. [12] by means of Raman scattering (T up = T down = 320 K) are not consistent with the other data listed in Table 1. In the following, we will refer to the averaged transition temperatures indicated in Table 1 as 'static' values.

The Raman spectra at 633 nm of the iron(II) SCO [Fe(Htrz) 2 (trz)](BF 4 ) polymer at 300 K in the LS state and at 393 K in the HS state are displayed in Figure 2. To minimize the interaction between the laser and the sample, very low laser intensity (0.8 lW/lm 2 at the sample) was set. As already mentioned in the paper of Urakawa et al. [12], several Raman bands can be identified as spin state markers (Figure 2). The most spectacular changes when crossing the LS to HS transition are observed below 500 cm À1 , i.e. in a spectral region where external vibrations are active, and particularly those involving the metal-ligand stretching or ligand-metal-ligand bending vibrations [22]. The LS marker bands are found at ca. 200, 213, 287 and 300 cm À1 and the HS marker bands at ca. 105, 150 and 182 cm À1 .

We have carried out preliminary Raman scattering experiments at room temperature to investigate the effect of the laser intensity and the excitation wavelength on the spin state. We found that for

Table 1

Typical features (T up , T down and hysteresis width) characterizing the thermal transition in the [Fe(Htrz) 2 (trz)](BF 4 ) SCO complex obtained using various experimental techniques. In bold are reported the average values of these characteristics to the exclusion of those obtained by Urakawa et al. [12].

Technique

T up T down Width K K K DSC [6] 379 343 39 DTA [6] 383 341 42 Reflectivity [6] 381 347 34 Reflectivity [18] 380 346 34 Magnetic susceptibility [6] 385 345 49 Magnetic susceptibility [11] 384 343 41 Magnetic susceptibility [10] 386 343 43 Mössbauer spectroscopy [6] 380 344 36 Average 382 344 40 X-ray diffraction [12] 370 330 40 Raman spectroscopy [12] 320 320 0 a 8 s sample irradiation at the 3 mW/lm 2 laser intensity, the complete LS to HS conversion was reached under 532 nm excitation, a partial conversion under 633 nm excitation and no detectable conversion under 785 nm irradiation. This may be understood as a consequence of the sample pink color at room temperature, since the maximum of light absorption due to the d-d 1 A 1 ? 1 T 1 transition [18,34] of the iron(II) complex is $500-600 nm. The maximum efficiency for the spin conversion LS ? HS, either photo-induced or thermally induced, is therefore expected for green excitation wavelengths. However, there is still weak absorption at higher wavelengths in the tail of the absorption band, i.e. in the near-infrared [35]. The absorption of the focused laser light necessarily heats the sample at the focal point and leads to a rapid increase of its local temperature, which may eventually cause strong sample damages under prolonged light exposition. This latter undesirable effect can be however limited (providing that the laser intensity is not too high) by rotating, swinging [12] the sample, or using galvanometric mirrors [14] to scan it fast. Figures 3 and4 show Raman spectra associated with a single region of the [Fe(Htrz) 2 (trz)](BF 4 ) sample excited at 633 nm along the ascending branch of its hysteresis loop, at 340 K and 365 K, respectively. The laser exposition time was systematically set to 160 s. Owing to the proportionality of the Raman scattering and incoming laser beam intensities, experimental Raman intensities were divided by the incident laser intensities to allow comparisons between spectra obtained with different laser powers. Note that Raman spectra (2) and ( 5) related to the highest laser intensities exhibit greater signal to noise ratios than the other spectra resulting from measurements at lower laser powers, because the same integration time was kept for all the spectra.

Outside the hysteresis loop at 340 K, the LS state is preserved under irradiation at the 0.2 mW/lm 2 laser intensity, as revealed by the Raman spectrum (1) in Figure 3 similar to the LS spectrum depicted in Figure 2. The consecutive excitation at the higher 1 mW/lm 2 laser intensity leads however to a mixture of LS and HS species, with a large majority of HS form. This assertion is corroborated by the presence of intense HS marker bands centered at 105, 150 and 182 cm À1 and weak LS marker bands at 200, 213, 287 and 300 cm À1 (Raman spectrum (2) in Figure 3). This LS/HS mixture relaxes back to the LS state after decreasing the laser intensity down to the 0.2 mW/lm 2 initial value, as evidenced by the sharp resemblance between Raman spectra (1) and (3) in Figure 3. In conclusion, outside the hysteresis loop, the LS to HS transition can be triggered by the continuous Raman-probe laser and it produces a stationary HS fraction as long as the irradiation is maintained. As the light-induced HS state is not stable however, the LS state is recovered as soon as the laser power density is sufficiently reduced.

Inside the hysteresis loop at 365 K, the LS state is still preserved under irradiation at the 0.06 mW/lm 2 laser intensity, as revealed by the Raman spectrum (4) in Figure 4 similar to the LS spectra displayed in Figures 2 and3. This suggests that the 0.06 mW/lm 2 laser intensity is not high enough to induce a spin transition at 365 K. Under 0.2 mW/lm 2 laser irradiation, a HS/LS mixture, as observed under 1 mW/lm 2 irradiation outside the hysteresis loop at 340 K, is seemingly produced as revealed by the striking resemblance between the Raman spectra (5) in Figure 4 and (2) in Figure 3. However, the Raman spectrum (6) obtained under excitation back to the 0.06 mW/lm 2 initial laser intensity matches almost perfectly the Raman spectrum (5). A long-lived metastable laser-induced HS state is populated within the hysteresis loop thereby. 4), increased up to 0.2 mW/lm 2 in ( 5) and finally reduced down to 0.06 mW/lm 2 in (6). The excitation wavelength was 633 nm and the integration time was 160 s. To further illustrate the effect of continuous light on the [Fe(Htrz) 2 (trz)](BF 4 ) SCO material using a conventional Raman spectrometer, the spin transition has been triggered in a planar and uniform region of the powder sample (inside the hysteresis loop at 350 K) with a 10Â microscope objective (NA = 0.25) so that the theoretical diameter of the laser spot was estimated to be $3 lm. After 2 s illumination at 0.4 mW/lm 2 power density, this value was reduced down to 0.06 mW/lm 2 in order to image the irradiated region of the sample (with a 50Â objective to have a better lateral spatial resolution) without affecting its spin state (Figure 5). For each Raman spectrum of the map the integrated intensity of the HS marker band around 100 cm À1 was measured and normalized to 1 when the sample contains 100% of HS species.

The HS fraction is then defined as the parameter c complying 0 < c < 1. In Figure 5, the yellow spot in the center (diameter <10 lm) represents the sample region where more than 85% of the species are in the HS state. It coincides approximately to the laser spot diameter. Interestingly, a 25 lm-broad region between the center in yellow and the border in blue (with more than 90% of LS species) contain both HS and LS species. The population of HS species in a region that was not directly excited by the laser may be explained in terms of heat diffusion and/or cooperative processes. The dissipation of energy following the irradiation that causes the local heating of the sample in neighboring regions of the laser spot would induce the observed partial spin conversion. This phenomenon suggests in particular that the size matching of pump and probe laser spots is crucial to evaluate correctly the light-induced HS fraction under pulsed laser excitation.

A laser-induced thermal effect is expected to perturb the thermal hysteresis loop, and shift in particular transition temperatures T up and T down . Figure 6 shows the change in HS fraction c (calculated from intensity integration of HS and LS marker bands at 100 cm À1 and 275 cm À1 , respectively) upon thermal cycling between 298 K and 420 K. In Figure 6, open symbols are related to c values obtained upon heating the [Fe(Htrz) 2 (trz)](BF 4 ) sample from 298 K to 420 K. Filled symbols correspond to c values obtained upon cooling the sample from 420 K to 298 K. These data points describe the ascending and descending branches of the hysteresis loop for 3 different laser power densities, namely $0.02 mW/lm 2 (circles), $0.06 mW/lm 2 (squares) and $0.2 mW/ lm 2 (diamonds). At low laser power density (below 0.02 mW/ lm 2 ), T up = 383 K and T down = 346 K. These temperatures are very close to their respective static values (T up = 382 K and T down = 344 K in Table 1). Thus, for such low laser intensity, the thermal hysteresis loop deduced from Raman measurements compares very well with those measured using the different techniques listed in Table 1, unlike the assertions of Urakawa et al. [12] However, upon increasing the laser intensity at the sample up to $0.2 mW/lm 2 , the apparent T up decreases significantly and the LS to HS conversion on the ascending branch becomes more gradual while neither T down nor the shape of descending branch are affected.

Surprisingly, experimental observations reported here reproduce well those published by Gallé et al. who excited the same kind of compounds using a pulsed laser irradiation [16,17]. In their case, a single laser pulse was also able to trigger the formation of stable HS species upon exciting LS species inside the SCO hysteresis loop while these HS entities relaxed back to the LS state within $1 ms upon irradiation outside the hysteresis loop [16,17], which agrees well with our own observations. Furthermore, the downshifted ascending branch upon continuous irradiation at 0.2 mW/ lm 2 (open diamonds in Figure 6) show a strong resemblance with the distribution, as function of the temperature, of the LS/HS ratio generated by the one-shot laser pulse experiment in another SCO compound [15] ([Fe(PM-BiA) 2 (NCS) 2 ] (PM-BiA = N-2 0 -pyr-idylmethylene-4-aminobiphenyl). In order to interpret their data, Freysz et al. [15] proposed a model where the partial conversion of LS iron(II) centers to the HS state is due to a laser-induced local heating, which mainly derives from the strong absorption of light by molecular entities in the LS state. According to this model we may suggest, in first approximation, that the weak absorption of HS species at 633 nm can be neglected indeed. As a consequence, the residual LS fraction gradually decreases as the number of HS particles increases owing to the reduced absorption of the 633 nm line in the HS state. Below a critical value of LS fraction there are not enough absorbing centers that would promote further the LS to HS transitions. On the contrary, the local heating is almost negligible above 360 K on the descending branch because all species are in the HS state. When the temperature of the system approaches T down , the LS fraction starts increasing (it is 20% at 345 K for example) but, it is still too small to allow the continuous laser excitation to trigger a thermal LS to HS conversion. This behavior yields the insensitivity to the laser power of the descending branch of the hysteresis loop. A more detailed model will be provided in a near future. Finally, this light-induced spin conversion mechanism should explain the reason why a photo-induced phase transition from the HS to the LS state cannot occur in the triazole-based polymeric samples, using a continuous laser at wavelength 633 nm as well as a single laser pulse.

Conclusion

We have shown, in agreement with Craig et al. [32], that a simple continuous laser with relatively small power (a class 2 laser with power 1 mW would be largely sufficient) may be exploited for inducing the photoswitching of a spin crossover material. Furthermore we found that, outside the hysteresis loop, the LS to HS transition can be triggered by the continuous Raman-probe laser and produce a stationary HS fraction as long as the irradiation is maintained. As the light-induced HS state is not stable however, the LS state is recovered as soon as the laser power density is sufficiently reduced. In contrast a long-lived HS state can be populated within the hysteresis loop even at low laser power. Thanks to the microscale spatial resolution of our Raman microscope, we have also shown that the dissipation of energy following the irradiation causes the local heating of the sample in neighboring regions of the laser spot, inducing partial LS to HS conversion. These results indicate that Raman microscopy performed on these SCO compounds offer many interesting advantages for promoting and characterizing optical data recording.

In addition, providing that very low laser intensities are used (<0.02 mW/lm 2 ), the thermal hysteresis loop deduced from Raman measurements compares very well with those measured using different techniques described in the literature. The thermal hysteresis loops are modified for laser intensities larger than 0.02 mW/lm 2 and the observed changes can be interpreted in terms of local heating of the sample. Finally, this light-induced spin conversion mechanism should explain the reason why a photoinduced phase transition from the HS to the LS state cannot occur in the triazole-based polymeric samples, using a continuous laser at wavelength 633 nm as well as a single laser pulse.

Figure 1 .

 1 Figure 1. Schematic view of the iron(II) SCO [Fe(Htrz) 2 (trz)](BF 4 ) polymer.

Figure 2 .

 2 Raman spectra of the iron(II) SCO [Fe(Htrz) 2 (trz)](BF 4 ) polymer in the LS state (top) and HS state (bottom).

Figure 3 .

 3 Figure 3. Raman spectra of the iron(II) SCO [Fe(Htrz) 2 (trz)](BF 4 ) polymer recorded outside the hysteresis loop at 340 K. The laser intensity is set to 0.2 mW/lm 2 in (1), increased up to 1 mW/lm 2 in (2) and finally reduced to 0.2 mW/lm 2 in (3). The excitation wavelength was 633 nm and the integration time was 160 s.

Figure 4 .

 4 Figure 4. Raman spectra of the iron(II) SCO [Fe(Htrz) 2 (trz)](BF 4 ) polymer recorded inside the hysteresis loop at 365 K. The laser intensity is 0.06 mW/lm 2 in (4), increased up to 0.2 mW/lm 2 in (5) and finally reduced down to 0.06 mW/lm 2 in (6). The excitation wavelength was 633 nm and the integration time was 160 s.

Figure 5 .

 5 Figure 5. Raman map (3 lm steps) of the powder sample irradiated at 633 nm with a laser intensity 0.4 mW/lm 2 and using a 10Â objective. The color scale bar in the inset indicates c values. The excitation wavelength for recording the map was 633 nm, the laser intensity was 0.06 mW/lm 2 and the integration time was 80 s for each spectrum. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Figure 6 .

 6 Figure 6. Evolution of the HS fraction c versus temperature as deduced from the Raman spectra for 3 different laser intensities at 633 nm, namely $0.02 mW/lm 2 (circles), $0.06 mW/lm 2 (squares) and $0.2 mW/lm 2 (diamonds). The temperature is either increased from 298 to 420 K (open symbols) or decreased from 420 to 298 K (filled symbols). The sigmoid functions (solid and dotted lines) were drawn to guide the eyes.

Acknowledgment

The authors thank the cluster Advanced Materials in Aquitaine (post-doctoral grant to Y. T.), the Conseil Régional d'Aquitaine and Europe (FEDER program) for funding equipment of the Vibrational Spectroscopy and Imaging (SIV) platform at ISM.