domain: spi phys

text:
Marco Carini
email: marco.carini@onera.fr

Nathaniel T Baker

Benjamin Leclaire

Franco Auteri

Pressure from Tomographic PIV: the Schur Complement method

Introduction

Thanks to the recent progress made by tomographic PIV techniques [START_REF] Scarano | Tomographic piv: principles and practice[END_REF]], it is of great interest the experimental-based reconstruction of a complete pressure field throughout a spatial domain, given the non-intrusive measurement of the associated velocity field. Such measurement and reconstruction techniques offer indeed unprecedented opportunities of research, ranging from the fundamental study of turbulence and the role of the pressure-velocity correlations, to the the estimation of the loads acting on structures immersed in a flow and the investigation of aero-acoustic noise.

In order to compute a pressure field from tomo-PIV data, one must integrate the global balance of momentum for the fluid, in which the inertial and viscous terms act as a known source term. Two different approaches are usually encountered in the literature to tackle this problem [START_REF] Van Oudheusden | Piv-based pressure measurement[END_REF], [START_REF] Kat | Instantaneous planar pressure determination from piv in turbulent flow[END_REF]]. The first one consists in directly integrating the pressure gradient, while the second relies on solving a Poisson equation derived by invoking the incompressible nature of the flow. Regardless of the adopted approach, it remains necessary to impose relevant boundary conditions on the pressure field in order to correctly pose the problem.

While the direct integration of the pressure gradient may theoretically be initiated from a single point in space where the pressure is known (for instance from a local measurement, or by invoking the Bernoulli equation in an irrotational part of the flow [START_REF] Van Oudheusden | Principles and application of velocimetry-based planar pressure imaging in compressible flows with shocks[END_REF]]), such approach is known to rapidly propagate errors at each incremental step, and requires omnidirectionnal marching schemes to achieve a reliable solution [START_REF] Liu | Instantaneous pressure and material acceleration measurements using a fourexposure piv system[END_REF]]. On the other hand, solving the Poisson equation requires knowing the pressure field on the boundary of the domain, which is seldom available in an experimental context. Indeed, Pitot probes only yield punctual and intrusive measurements, while Tomo-PIV volumes are usually too small to intersect an irrotational flow everywhere. Imposing Neumann boundary conditions stemming from the global balance of momentum is neither an option, as doing so would lead to an ill-posed mathematical problem [START_REF] Rempfer | On boundary conditions for incompressible navier-stokes problems[END_REF]]. Within this context, [START_REF] Logg | Error propagation dynamics of piv-based pressure field calculations: How well does the pressure poisson solver perform inherently?[END_REF] showed in particular that the overall accuracy of a reconstructed pressure field is inherently linked to the type of boundary condition used (i.e. Dirichlet vs. Neumann).

As a matter of fact, in the uncoupled pressure-velocity formulation of the incompressible Navier-Stokes equations, the strict mathematical formulation of the boundary conditions associated to the Poisson equation for pressure is of integral form [START_REF] Quartapelle | Integral conditions for the pressure in the computation of incompressible viscous flows[END_REF]]. Dirichlet conditions for the pressure field, fully equivalent to the integral ones, can be obtained through the Glowinski-Pironneau method for pressure-velocity decoupling. [START_REF] Glowinski | On a mixed finite element approximation of the stokes problem (i)[END_REF]]. This approach has been successfully employed by [START_REF] Auteri | A novel approach for reconstructing pressure from piv velocity measurements[END_REF] for pressure reconstruction in the 2D cases. The method basically involves the solution of a set of elliptic problems as many as the total number of pressure degrees of freedom on the boundary of the domain. Therefore its extension to 3D cases, although possible, will necessarily results less computationally efficient. The method proposed here relies on the same mathematical framework introduced above while exploiting a fully discrete (or algebraic) approach for the derivation of the Poisson pressure problem, which results in a straightforward mathematical formulation and a more efficient numerical implementation for large three-dimensional PIV datasets.

The paper is structured as follows. We will start by presenting the mathematical formulation and the numerical implementation of the method in sections 2 and 3, respectively. We will then validate the algorithm on an exact analytical solution of the Navier-Stokes equations in section 4, before applying the method to synthetic and experimental PIV data in sections 5 and 6. Finally conclusions are given in section 7.

Principle of the method

The so-called Pressure Schur Complement (PSC) method introduced herein builds on the approach recently introduced in planar PIV by [START_REF] Auteri | A novel approach for reconstructing pressure from piv velocity measurements[END_REF]. Its major specificity is its discrete formulation, that leads to a simplified implementation and better computational performances. We consider the unsteady Reynolds-Averaged Navier-Stokes equations, with first-order time discretization denoted by superscripts:

u n -u n-1 ∆t + (u n • ∇) u n -ν∇ 2 u n + ∇p n = -∇ • u n ⊗ u n and ∇ • u n = 0 on Ω, u n | ∂Ω = b n , (1)
that can be recast into the following Stokes problem:

-∇ 2 w + ∇ p = g and ∇ • w = 0 on Ω,

w| ∂Ω = b (2)
upon introducing the notations p = p/ν, w = u n , b = b n and g = -

1 ν u n -u n-1 ∆t + (u n • ∇) u n + ∇ • u n ⊗ u n .
Right-hand side quantities of problem (2), b and g, can be computed directly from Tomo-PIV velocity fields.

It is then either possible to reconstruct unsteady pressure fields from time-resolved data, upon choosing u n = 0, or time-averaged pressure fields, upon setting the first term in g to zero, replacing u n by the mean flow u , and u n ⊗ u n by the Reynolds stress tensor u ⊗ u . As specified in (2), only Dirichlet conditions on u n or u , directly provided by the measurements, need to be prescribed at the domain boundaries. Upon introducing a finite-element discretisation of the velocity and pressure fields, problem (2) reduces to a linear system of the form:

L u G G T 0 w p = g 0 . (3)
Now, terms w, p and g correspond to the degrees of freedom (dof) of their associated physical variables, and matrices L u , G and G T to the discretized weak forms of the Laplacian, gradient and divergence operators. Note that, depending on the specific technique adopted for the implementation of Dirichlet boundary condition on the velocity, one can have GT = G T . It is then straightforward to derive a problem pertaining only to the dofs for pressure by computing the associated Schur complement, as well known in numerical methods for incompressible flows:

G T L -1 u G p = G T L -1 u g. (4)

Numerical implementation

The numerical implementation of the PSC method has been performed by using the open-source finite element library FeNiCs [START_REF] Logg | Error propagation dynamics of piv-based pressure field calculations: How well does the pressure poisson solver perform inherently?[END_REF]] and the open-source library of linear algebraic solvers PETSc [START_REF] Balay | PETSc Web[END_REF]]. The Stokes system (2) has been discretised on an unstructured tetrahedric mesh (directly obtained from the original PIV Cartesian mesh) by using a classical Taylor-Hood P2-P1 formulation, thus fully satisfying the well-known inf-sup condition associated with the incompressibility constraint.

While the sparse matrices L u , G and G T are thus assembled, the explicit construction and storage of the PSC matrix G T L -1 u G is never performed due to its dense structure. Instead, the Schur linear system is solved by an iterative approach, through the GMRES algorithm available within PETSc, thus only requiring matrix-vector multiplications to be performed. The same algorithm, with an additional algebraic multigrid preconditionner, has also been employed for the efficient numerical inversion of the large-scale matrix L u , as required at each iteration of the external PSC GMRES loop. Therefore a two-level nested GMRES approach is used to solve (4). Both the finite element assemblage and the inversion of the involved large-scale matrices have been performed in parallel thanks to the inherent parallelism of the employed numerical libraries and the available ONERA HPC facilities.

Validation and performance assessment on an analytical solution

The PSC method and its convergence properties are here validated with respect to an analytical solution. For this purpose, we consider the following pressure field mimicking a vortical flow:

p ref = ρU 2 0 - 1 8 sin 2 (x) sin 2 (y) + sin 2 (z) + 1 4 sin 2 (y) sin 2 (z) e -6νt . (5)
The above pressure field is an exact solution to the unsteady Navier-Stokes equations subject to the external forcing f

∂ u ∂t + u • ∇ u = - ∇p ρ + ν∆ u + f ρ , (6)
with

u = U 0   -cos(x) sin(y) sin(z) 1 2 sin(x) cos(y) sin(z) 1 2 sin(x) sin(y) cos(z)   e -3νt , (7)
and

f = - 3 4 ρU 2 0 cos(x) sin(x) sin 2 (y) + sin 2 (z) 0 0 e -6νt . (8)
Figure 1 shows the magnitude of the velocity field (left) and the pressure field (right) given before. These fields are sampled across a three-dimensional cubic domain of size L = 2π, with N = 32 points in each spatial direction.

Sensitivity to spatial resolution

First of all, let us assess the convergence of the reconstructed pressure field with respect to the spatial resolution h = L/N, where N represents the number of points used to discretise the domain along each spatial direction. The relative error between the exact solution p ref (Eq. 5) and the pressure p numerically reconstructed using the Schur complement method fed with the velocity field u (Eq. 7) is quantified using both the classical 2 and ∞ norms for the nodal pressure values, respectively defined by

E 2 = ∑ Ω (p ref -p) 2 ∑ Ω p 2 ref , (9)
and

E ∞ = max |p ref -p| max |p ref | . ( 10
)
Results are reported in figure 2, clearly confirming that the pressure field reconstructed using the PSC method converges towards the exact solution p ref with second-order spatial accuracy, as expected.

Sensitivity to noise

Since experimental data is inherently corrupted by some level of noise, the robustness of the algorithm to the presence of the noise is here briefly investigated. The adopted methodology consists in adding a random noise with a Gaussian distribution to each component of the input velocity field, and at each vertex in the computational domain. The noise is assumed to be spatially uncorrelated, but time correlated (indeed assuming a white noise also in time will lead to very large errors when computing its derivative). The noise level is characterized by the standard deviation σ, with σ = 1, 5 and 10% of U 0 . Figure 3 displays the pressure mean error distribution (left) and the standard deviation distribution (right) over the entire domain. These statistics have been computed over 500 trials, which yield a convergence level better than 10 -4 . This figure shows that the spatial distribution of the mean error is not uniform throughout space. On the one hand, comparing the spatial distribution of the mean error on the reconstructed field with the velocity field (Figure 1) suggests that the mean error is maximum where the velocity gradient is maximum. On the other hand, comparing the spatial distribution of the error's standard deviation suggests that it scales as the amplitude of the velocity field, in contrast to what observed by [START_REF] Auteri | A novel approach for reconstructing pressure from piv velocity measurements[END_REF].

A more quantitative description of the impact of noise on the results is shown in Figure 4, which displays the statistical distribution of the error on the pressure field (ε P ). As a reference, the noise distribution added to all three velocity components is shown (ε U , ε V , ε W ), which confirms that it is indeed Gaussian. Interestingly, the error distribution of the reconstructed pressure dramatically departs from a Gaussian distribution, and is much flatter than the latter. This means that the expected errors made on the reconstructed pressure field are larger than those on the velocity. Yet, the amplitude of the error on the reconstructed pressure still scales according to the amplitude of the noise added to the velocity field.

Demonstration on synthetic PIV data

Let us now apply the method to synthetic PIV data numerically generated by an IDDES simulation performed by [START_REF] Gomit | Comparaison de méthodes dévaluation de la pression partir de mesures piv. 16ème Congrès Francophone de Techniques Laser Figure 8: Results of PSC pressure reconstruction from Tomo PIV measurements of the SACCON test case[END_REF]. The geometry adopted here is a NACA0015 wing of 80 mm chord and 140 mm span, placed at a 30 o angle of attack in a 230 mm × 230 mm square hydraulic tunnel. The freestream value of the incoming flow is fixed at U = 1.25 m/s, which corresponds to a Reynolds number based on the chord of Re c = 100, 000. The computed flow is then extracted along three sub-volumes of the computational domain,which mimic typical Tomo-PIV volumes (cf. Figure 5). The data within these volumes has been filtered to simulate the smoothing effect of PIV, and the resulting spatial resolution is 1 mm. The pressure field was reconstructed in volumes 2 and 4 of figure 5. Doing so makes it possible to assess the behavior of the method in two different situations: in the wake of an airfoil (volume 2), and in the presence of a rigid obstacle (volume 4). The pressure reconstruction was performed on one of ONERA's clusters. The most computational intensive case (i.e. the wake in volume 4) required to solve a linear system of size 689778 × 689778 for the pressure and 15843927 × 15843927 for the velocity, whose convergence was reached in 83 GMRES iterations. As a matter of fact, the computation took about 20 minutes, thanks to the parallelization of the job on 250 processors. The maximum memory required during the job was 152 Gb.

Figure 6 provides 2D views of the reconstructed pressure field along the mid-plane of volumes 4 and 2 respectively. The top graphs show the actual pressure field, while the bottom graphs give the relative error between the reconstructed pressure field and the reference pressure field given by the simulation. Note that these plots use log colorbars. As with the analytical case studied earlier, this figure confirms that the error distribution is not uniform throughout the PIV domain, and is strongly linked to the presence of velocity gradients. Indeed, the error appears to be highest at the location of rigid obstacles where the noslip boundary condition applies, and in the presence of strong vortical structures. A global estimate for the accuracy of the reconstructed pressure using the 2 norm (Eq. 9) yields E 2 = 10% in volume 2 and E 2 = 3% in volume 4.

Demonstration on experimental data

To assess the method's performance in real conditions, the flow around a generic flying wing model (SAC-CON configuration, see e.g. [START_REF] Schütte | Flow physics analyses of a generic unmanned combat aerial vehicle configuration[END_REF], at a chord Reynolds number of 1.1 × 10 6 and angle-ofattack of 20 o , has been considered. Tomographic PIV experiments have been performed in ONERA L1 wind-tunnel of the Lille ONERA centre. A sketch of the model and a picture of the experimental setup are reported in Figure 7. The volume illumination has been placed so as to capture the wingtip vortex, and to match with a generatrix of pressure taps located on the wing surface, yielding a reference for the reconstruction. Further details and results are described in our companion paper by Baker et al.

Volumic mean pressure computation using the PSC method in RANS mode has been performed using mean flows and Reynolds stresses resulting from tomographic PIV. In practice, for the case presented here, this led to solve a 58378737 × 58378737 system, which was run on 840 cores, and required up to 1.05 TB of memory. Figure 8 shows results of this reconstruction in the plane coinciding with the pressure tap generatrix at the model wall. Quantity displayed are contours of pressure coefficient C p = 2(pp ∞ )/ρU ∞ in the whole plane, and extracted closest to the wall for comparison with the taps. In the case of PSC, no particular tap was used as a reference for the absolute pressure level. Instead, it was preferred to choose this reference value so that the integral of both C p curves (PSC and pressure taps) coincide, in order not to favor a particular Figure 6: Performance assessment on two volumes extracted from an IDDES numerical simulation around a finite span NACA0015 wing (Re c = 100, 000, see [START_REF] Gomit | Comparaison de méthodes dévaluation de la pression partir de mesures piv. 16ème Congrès Francophone de Techniques Laser Figure 8: Results of PSC pressure reconstruction from Tomo PIV measurements of the SACCON test case[END_REF]. Top: Reconstructed pressure fields. Bottom: Relative error between reconstruction and ground truth (note the use of a log scale). Left column corresponds to a flow zone close to the wingtip, right column to a zone within the wing wake.

sensor. As seen in Figure 8, the PSC method appears to capture truthfully the signature of the leading edge vortex, which is characterized by a strong depression at the location of the vortex core. The pressure profile extracted closest to the wing's upper surface compares well with the pressure tap measurements (cf. Figure 8, right) along almost the whole profile, for -0.16 ≤ y ≤ -0.10. For y ≥ -0.8, i.e. closest to the vortex core, reconstructed C p however deviates from the pressure tap measurements. This is probably to be ascribed to the fact that, due to remaining light reflection, the lower edge of the reconstructed volume does not exactly coincide with the wall, but is rather located a few mm above. As it is well known, curved streamlines (such as in the case of a vortex) are associated with a pressure gradient normal to the curvature, so that extrapolation of pressure directly below a vortex core can be particularly error-prone.

Conclusion

This paper has presented a new method to reconstruct a pressure field from the tomo-PIV measurements of an incompressible flow. This method essentially relies on re-casting the Navier-Stokes equations into a forced Stokes problem, whose forcing term is completely determined by experimental data. Under this form, boundary conditions have to be prescribed only for the velocity field (which is known everywhere from the PIV measurements) by contrast to other Poisson-based approaches available in the literature, requiring explicit pressure boundary data which are rarely accessible experimentally.

First, the soundness of the algorithm has been validated against an analytical solution of the unsteady incompressible Navier-Stokes equations in order to assess the convergence of the method with respect to spatial, as well as its filtering behaviour w.r.t. to spatial noise. In particular, we have shown that the amplification of noise in the velocity data (which is unvoidable experimentally) roughly scales according to the amplitude of the noise on the velocity field.

We have then assessed the quality of the reconstructed pressure field on a more realistic application, namely a time-resolved numerical simulation of the flow around a NACA0015 airfoil, and a time-averaged actual experiment of the flow over the upper surface of a flying wing. As observed for the analytical solution, the error on the reconstructed pressure in the complex flow around the airfoil is not uniform and strongly depends on the flow gradients as encountered approaching the walls and in the core region of vortical structures.

Last but not least, the promising results obtained for the SACCON flying wing configuration have proved that the method is indeed applicable even to industrial-size facilities, which paves the way for advanced applications in the aerospace industry.

Figure 1 :Figure 2 :

 12 Figure 1: (Left) Visualization of the velocity field u given by Eq. 7. (Right) Reference pressure field given by Eq. 5

Figure 3 :

 3 Figure 3: (Left) Spatial distribution of the mean error. (Right) Spatial distribution of the error standard deviation.

Figure 4 :

 4 Figure 4: Probability distribution of the error on the reconstructed pressure p (colored symbols) for a Gaussian distribution of the error on the velocity field (empty symbols, one for each spatial coordinate).

Figure 5 :

 5 Figure 5: Visualization of the PIV volumes extracted from the numerical database.

Figure 7 :

 7 Figure 7: SACCON Experimental: sketch of the flying wing model for PIV measurement (left); picture of the PIV experimental setup (right).

Acknowledgements

The authors acknowledge the financial support of the Agence National pour la Recherche and the DGA under the grant ANR-16-ASTR-0005-01, as well as the CPER-FEDER of the Hauts-de-France and Nouvelle-Aquitaine regions.


+ P Monmousseau

A Marzuoli

E Feron

D Delahaye

[EN-A-46] Passengers on social media: A real-time estimator of the state of the US air transportation system (EIWAC 2019)

Keywords: delay estimation, ATM performance measurement, big data, machine learning

This paper aims at investigating further the use of the social media Twitter as a real-time estimator of the US Air Transportation system. Two different machine learning regressors have been trained on this 2017 passenger-centric dataset and tested on the first two months of 2018 for the estimation of air traffic delays at departure and arrival at 34 different US airports. Using three different levels of content-related features created from the flow of social media posts led to the extraction of useful information about the current state of the air traffic system. The resulting methods yield higher estimation performances than traditional state-of-the-art and off-the-shelf time-series forecasting techniques performed on flight-centric data for more than 28 airports. Moreover the features extracted can also be used to start a passenger-centric analysis of the Air Transportation system. This paper is the continuation of previous works focusing on estimating air traffic delays leveraging a real-time publicly available passenger-centered data source. The results of this study suggest a method to use passenger-centric datasources as an estimator of the current state of the different actors of the air transportation system in real-time.

INTRODUCTION

The Air Transportation System is a complex interconnected system that carried more than 631 million passengers on domestic flights in the United States in 2010 according to the Bureau of Transportation Statistics (BTS) [1]. Flight delays are still a major issue both in the United States with 27.0% of departing flights and 27.8% of arriving flights experiencing delays in 2017 [1].

Most previous studies aimed at predicting or classifying flight delays were centered on flight-centric information coming from a variety of sources with different levels of public availability, and using only very little passenger-centric data. Mueller and Chatterji [2] created a probabilistic model of delays by fitting Poisson and Normal distributions to the historic delay data from 10 airports. Rebollo and Balakrishnan [START_REF] Rebollo | Characterization and prediction of air traffic delays[END_REF] implemented a network model to classify and predict future delays on specific links or specific airports using two years of flight-centric and weather-related data. Klein et al. [START_REF] Klein | Airport delay prediction using weather-impacted traffic index (WITI) model[END_REF] and [START_REF] Sridhar | Short-term national airspace system delay prediction using weather impacted traffic index[END_REF] focused on predicting short-term weather-related delays using only past and current weather information. Aljubairy et al. [START_REF] Aljubairy | Real-time investigation of flight delays based on the Internet of Things data[END_REF] used Internet of Things in order to analyze flight-related sensors in real-time and classify the delay of an upcoming flight.

Over the past few years, NextGen [START_REF]NextGen Integration and Implementation Office[END_REF] in the United States has been advocating a shift from flight-centric metrics to passenger-centric metrics to evaluate the performance of the Air Transportation System. The failures and inefficiencies of the air transportation system not only have a significant economic impact but they also stress the importance of putting the passenger at the core of the system [START_REF]Connected World : Transforming Travel, Transportation and Supply Chains[END_REF] [START_REF]Smart travel: Unlocking economic growth and development through travel facilitation[END_REF]. Several studies have highlighted the disproportionate impact of airside disruptions on passenger door-to-door journeys. Flight delays do not accurately reflect the delays imposed upon passengers' full multi-modal itinerary. Cook et al. [START_REF] Cook | Passenger-Oriented Enhanced Metrics[END_REF] designed propagation-centric and passenger-centric performance metrics, and compare them with existing flight-centric metrics. Wang [START_REF] Wang | Methods for analysis of passenger trip performance in a complex networked transportation system[END_REF] showed that high passenger trip delays are disproportionately generated by canceled flights and missed connections. Nine of the busiest thirty-five airports cause 50% of total passenger trip delays. Congestion, flight delay, load factor, flight cancellation time and airline cooperation policy are the most significant factors affecting total passenger trip delay. NextGen intends to not only improve the predictability and resilience of the US Air Transportation System, but also to reduce door-to-door travel time for passengers.

Passengers are at the core of this system and, yet, limited quantitative information about passenger movements is publicly shared. Each aviation stakeholder only has access to a partial view of the passenger-side of air transportation operations making a system-wide datadriven picture of passenger behavior difficult to implement. The BTS provides aggregated passenger data per market but no granular information. Passenger surveys conducted by airports or airlines, while very detailed, remain limited to small samples of passengers and short time periods, and may not be representative.

Precursor work was made by Marzuoli et al. in [START_REF] Marzuoli | Implementing and validating air passenger-centric metrics using mobile phone data[END_REF] and [START_REF] Marzuoli | Passengercentric metrics for Air Transportation leveraging mobile phone and Twitter data[END_REF] using mobile phone data in order to analyze the performances of airports from the passengers' perspective. These studies validated the use of passengercentric data to better assess the overall health of the Air Transportation System. However mobile phone data is proprietary data and is not often publicly available. In order to operate in real-time, it is thus necessary to also look into other sources of passenger data available on a national scale.

Another popular source of data previously used for studying large-scale behaviors with real time availability is social media, in particular Twitter. With more than 68 millions active users in the United States [START_REF] Statista | Monthly active Twitter users in the united states[END_REF], Twitter is an important pool of user-created data that is still not fully leveraged. Twitter has already been the main focus of many studies focused on its real-time availability, especially during natural disasters with multiple works by Palen et al. [START_REF] Palen | Twitterbased information distribution during the 2009 Red River Valley flood threat[END_REF], [START_REF] Vieweg | Microblogging during two natural hazards events: What twitter may contribute to situational awareness[END_REF], [START_REF] Kireyev | Applications of Topics Models to Analysis of Disaster-Related Twitter Data[END_REF]. Terpstra et al. also studied how a real time Twitter analysis could have provided valuable information for the operational response of a natural disaster crisis management with the case of the storm hitting a festival in Belgium [START_REF] Terpstra | Towards a realtime Twitter analysis during crises for operational crisis management[END_REF]. Regarding the air transportation field, most works mining Twitter data focus on how airlines are perceived by passengers by means of sentiment analysis [START_REF] Breen | Mining twitter for airline consumer senti-ment[END_REF] or sentiment classification [START_REF] Wan | An ensemble sentiment classification system of Twitter data for airline services analysis[END_REF]. Though these works give a good insight on how passengers perceive the state of some specific actors within the air transportation system, it does not give a global idea of its health. Monmousseau et al. in [START_REF] Monmousseau | Predicting and Analyzing US Air Traffic Delays using Passenger-centric Data-sources[END_REF] used publicly available social media data created by passengers to accurately estimate and predict the hourly aggregated status of the US air transportation system.

This paper proposes to build on this previous work in order to estimate the state of the air transportation system to a finer level. Rather than predicting the number of delays across all the United States, the proposed passenger-centric models are improved and tuned to accurately estimate the state of delays for each of the 35 major airports within the United States. The created models are based on three different levels of contentrelated features created from the flow of social media posts. First results indicate that these new models can estimate the number of hourly delays with a mean absolute error of less than 3 flights for 26 of the considered airports, and of less than 6 flights for the 9 remaining airports.

The rest of the paper is structured as follows: Section 2 describes the datasets and the feature extraction process. The methodology and results of the training process are shown in Section 3, before being analyzed and exploited in Section 4. Section 5 concludes this study and discusses possible future steps.

DATASET DESCRIPTION AND FEATURE SELECTION

Dataset description

Following the initial work performed in [START_REF] Monmousseau | Predicting and Analyzing US Air Traffic Delays using Passenger-centric Data-sources[END_REF], the goal here is to use passengers behavior on social media -in particular on Twitter -in order to analyze and estimate the flight-centric health of the US air-transportation system at an airport level. In this study, the flightcentric health of an airport is described by delay related information contained within BTS data. This data is publicly available usually with a two to three month delay and this study limits itself with the BTS data from January 2017 to February 2018.

The Twitter dataset available for this study is the same as in [START_REF] Monmousseau | Predicting and Analyzing US Air Traffic Delays using Passenger-centric Data-sources[END_REF] and consists of all the tweets found using a basic search for each handle of 7 major US airlines as well as 34 major US airports (one of them having two Twitter handles). The full list of handles can be found in Table 1. Each entry consists of a timestamp, a user id, the content of the tweet and the handle used to retrieve the tweet. This dataset spans the entire period from January 1st 2017 to February 28th 2018. The extraction of features from this dataset has been improved since the previous study and is described in Section 2.2.

Figure 1 plots the total number of tweets related to each handle over the year 2017 against the total number of flights flown by each airline or to and from each airport. Airlines tend to gather more tweets than airports, and the number of tweets is not necessarily correlated to the number of flights flown per airline. The handle "@Delta" gathered the most tweets over 2017 even though Southwest Airlines carried out the most flights in 2017. Most airports are regrouped around a cluster of 10k tweets and 200k flights over 2017, with

Category

Twitter handles Airlines @united, @Delta, @AmericanAir, @SouthwestAir, @SpiritAirlines, @VirginAmerica, @JetBlue Airports @JFKairport, @ATLairport, @flyLAXairport, @fly2ohare, @DFWAirport, @DENAirport, @CLTAirport, @LASairport, @PHXSkyHarbor, @MiamiAirportMIA, @iah, @EWRairport, @MCOAirport, @Official_MCO, @SeaTacAirport, @mspairport, @DTWeetin, @BostonLogan, @PHLAirport, @LGAairport, @FLLFlyer, @BWI_Airport, @Dulles_Airport, @MidwayAirport, @Reagan_Airport, @slcairport, @SanDiegoAirport, @flyTPA, @flypdx, @flystl, @flySFO,@HobbyAirport, @flynashville, @AUStinAirport, @KCIAirport Los Angeles International airport (LAX) and Hartsfield-Jackson Atlanta International airport (ATL) being exceptions due to their higher number of tweets. In order to estimate the flight-centric health of each considered airport, this information first needs to be extracted from the BTS dataset for each airport. Only two types of delayed flights are considered here from a passenger's perspective: Flights departing with any amount of delay, and flights arriving with a delay greater than 15 minutes. Once all the flights departing an airport and all the flights arriving at the same airport are selected, the following values can be aggregated per hour:

• NumDepDelay: Number of flights departing with a delay • NumArrDelay15: Number of flights arriving with a delay greater than 15 minutes The aim of this study is to accurately estimate these two values for each airport at every hour using a single passenger-centric dataset.

Feature selection on Twitter data

Volume features

Features were extracted identically for all search handles presented in Table 1, for the exception of @Miami-AirportMIA, which does not gather enough tweets. In addition to the raw number of tweets per hour per search handle, keyword related information is also extracted from the Twitter dataset. In order to keep all the relevant tweets without having to decline all the possible forms of the chosen keywords (e.g. "delay", "delayed", "delays", etc.), simple regular expression filters were created for each keyword: Any tweet containing a word starting with the related keyword is kept and the resulting tweets are then aggregated per hour. Five keywords were chosen for this study: 'delay', 'wait', 'cancel', 'hours', 'refund'.

Topic features

Another way of exploiting information from the content of these tweets is to perform a topic analysis of the tweet database using Latent Dirichlet Allocation [START_REF] Blei | Latent Dirichlet Allocation[END_REF] (LDA). A first step in topic analysis is to clean the documents analyzed, here the tweets. This cleaning process was already performed in [START_REF] Marzuoli | Passengercentric metrics for Air Transportation leveraging mobile phone and Twitter data[END_REF] and [START_REF] Monmousseau | Predicting and Analyzing US Air Traffic Delays using Passenger-centric Data-sources[END_REF] and consists of the following steps: Any reference to websites or pictures was replaced by a corresponding keyword. Every mention to another Twitter user within a tweet (@someone) as well as most emojis were similarly replaced. Since this database contains many replies from airlines to their customers, individual signatures of each agent were also replaced by a keyword. Dates and times were also generically replaced by keywords (e.g. "3rd Jan 2017" becomes "DATE" and "4pm" becomes "TIME"). The resulting text was then filtered from common stop-words and from words occurring only once in the whole year of 2017. A list of 100 topics was then created using the Gensim [START_REF] Sojka | Software Framework for Topic Modelling with Large Corpora[END_REF] library. The topic distribution of each tweet was then calculated before averaging these distributions per hour and per search handle. The hourly standard deviation of the distributions was also extracted.

This cleaning process introduces two additional keywords that enables a quick filtering of tweets, and therefore two additional features to add per search handle: tweets containing a picture and those containing a website link. Thus seven keywords are actually considered for feature extraction: 'delay', 'wait', 'cancel', 'hours', 'refund', 'PICTURE', 'WEBSITE'.

Sentiment features

Sentiment analysis is also used here to enhance the feature set considered. Two different datasets and cleaning method were used to train three different regressors each. The first dataset used was the labelled dataset used in a Kaggle competition [START_REF] Kaggle | Twitter US airline sentiment[END_REF] and was cleaned using the same process as for the previous LDA learning. The generic keywords from the cleaning process (e.g. 'WEBSITE', 'DATE') were removed before creating the associated dictionary, as well as words appearing in less than 20 tweets or in more than 75% of the full dataset. A second dataset and cleaning process was generated based on the work of Read [START_REF] Read | Using emoticons to reduce dependency in machine learning techniques for sentiment classification[END_REF]. Emoji filters were used to extract tweets from the initial dataset and automatically label them with a positive or negative sentiment according to Table 2. The text cleaning process is improved by merging negation words ("no", "not" and "never") with the word that follows it. The tokens used for the creation of the dictionary are the resulting bigrams, i.e. combinations of two words that follow each other in a tweet, with the same frequency filter as the first method described. )", "=)", ":-)", ";)", ";-)", ":-D", ":D", "=D" Negative ":(", ":-(", "=(", ":-@", ":'(", ":-|"

For both methods, three classifiers are trained (a random forest classifier, a naive Bayesian classifier and a logistic regressor) using the scikit-learn library [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. A sentiment score is then calculated for each tweet by averaging the output of these classifiers, 0 meaning a unanimous negative sentiment and 1 a unanimous positive sentiment. The hourly average of these scores are added to the Twitter feature set.

Summary

Given the temporal nature of the data analyzed, the following features were chosen to keep track of the date: month of the year, day of the month, day of the week and hour in the day. In summary the following 8,484 features are considered:

• Hourly volume of tweets for each search handle

ESTIMATING DELAYS

The aim of this section is to see how well it is possible to estimate per airport the number of flights departing with a delay and the number of flights arriving with a delay greater than 15 minutes using the features extracted from the Twitter dataset. The dataset was split into a training set consisting of the data from the year 2017, and a testing set with the data from January and February 2018.

Methodology

For each BTS value, two different machine learning regressors were trained on the training data set: a Random Forest regressor and a Gradient Boosting regressor. These regressors were implemented from scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF] with identical hyper-parameters. The maximum depth of each regressor was limited to ten, the minimum number of samples for a split was fixed to two and the maximum number of trees was fixed at ten.

As a comparison benchmark, we used Facebook's time-series forecasting tool Prophet [START_REF] Taylor | Forecasting at Scale[END_REF] on the 2017 BTS data to forecast the full two first months of 2018. The Prophet tool is based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality [START_REF] Facebook | Prophet -Forecasting at scale[END_REF]. It is described as robust to outliers and missing data with no parameter tuning necessary, therefore the default parameters of the Prophet tool was used for this forecasting benchmark.

Lastly, the standard deviation of the BTS values in the training set were calculated to illustrate the added value of the trained regressors. The performance measures used to compare the different regressors are presented in the upcoming section 3.2.

Estimation performance measures

In order to measure the performance of the different models, two different indicators were used: the R 2 score and the mean-absolute error (MAE).

The R 2 score, also known as the coefficient of determination, is defined as the unity minus the ratio of the residual sum of squares over the total sum of squares:

R 2 = 1 -i (y i -f i ) 2 i (y i -ȳ) 2 (1)
where y is the value to be predicted, ȳ its mean and f is the predicted value. It ranges from -∞ to 1, 1 being a perfect prediction and 0 meaning that the prediction does as well as constantly predicting the mean value for each occurence. In the case of a negative R 2 , then the model has a worse prediction than if it were predicting the mean value for each occurence and therefore yields no useful predictions. Regarding the mean-absolute error, the smaller its value is, the more accurate the prediction is. It is calculated using the following formula:

MAE = 1 n i |f i -y i | ( 2
)
where n is the number of values being predicted.

Estimation results

Fig. 2 shows a comparison per airport of the meanabsolute error of the two trained regressors along with the chosen benchmark for the estimation of the number of flights departing with a delay. The standard deviation of the number of delayed departing flights at each airport during the year 2017 is also included for comparison. The Random Forest models have the best results in this case: they outperform the Gradient Boosting models at allbut-one airports (LAX) and the Facebook Prophet tool on 31 airports out of 34. For 26 airports, the Random Forest models are able to estimate the hourly number of delayed departing flights with a mean-absolute error of three flights or less, and with an error of less than six flights for the remaining airports. Fig. 3 shows a comparison per airport of the meanabsolute error of the two trained regressors along with the chosen benchmark for the estimation of the number of flights arriving with a delay greater than 15 minutes. The standard deviation of the number of delayed arriving flights at each airport during the year 2017 is also included for comparison. The Random Forest models also have the best results in this case though their relative performance are not as important as for delayed departing flights : they outperform the Gradient Boosting models at 27 airports out of 34 and the Facebook Prophet tool on 28 airports out of 34. The absolute performance is however better than for estimating the number of delayed departing flights. For 28 airports, the Random Forest models are able to estimate the hourly number of delayed departing flights with a meanabsolute error of less than three flights, and with an error of less than five flights for the remaining airports.

Fig. 4 shows a comparison per airport of the R 2 score of the two trained regressors along with the chosen benchmark for the estimation of the number of flights departing with a delay. The Random Forest models still have the best results in this case, but the model associated with LAX airport also shows the only negative score. They outperform the Gradient Boosting models at 27 airports out of 34 and the Facebook Prophet tool on 28 airports out of 34.

ANALYSIS AND APPLICATIONS

The aim of this section is to analyze the differences between the chosen models as well as explore possible applications resulting from the extracted features.

Model analysis

Figure 5 shows the hourly prediction of the number of delayed departing flights at Atlanta airport (ATL) over the period January 12th-16th for the two trained regressors along with the benchmark and the actual values. This airport was chosen since it has the highest BTS standard deviation for the number of delayed departing and arriving flights, and the period was chosen to illustrate the high variability of the number of delays from a day to another. In this example, January 12 has more than twice as many delayed flights than any other day, as well as important hourly variations. Figure 5 illustrates the main differences between the different models: The Prophet tool predicts for each day a similar daily variation with three peaks during the day yet with amplitudes varying depending on the month and the day of the week. It also predicts negative values, which underlines some limitations of the model in this case. The added value from passenger-centric data-sources is better seen on January 12 and 13, where only the Random Forest regressor is able to estimate the higher number of delays on January 12 before correctly estimating the more usual levels of January 13. The Gradient Boosting regressor doesn't estimate outliers as well as the Random Forest regressor due to the difference in their loss functions. That difference is also illustrated by the non-zero minimum of the Gradient Boosting estimation during night time.

Other applications 4.2.1 Real-time sentiment analysis

The extracted features can be fed to the trained models for accurately estimating the number of delayed flights, but they can also be used directly in order to sense the overall passenger mood. Once the sentiment analyses are conducted on the tweets, it is possible to merge them into one score per airline and monitor their evolution.

Figure 6 shows the hourly average mood for three major airlines during the Northeastern bomb cyclone studied in [START_REF] Marzuoli | Passengercentric metrics for Air Transportation leveraging mobile phone and Twitter data[END_REF]. These three airlines have a similar passenger mood evolution at the beginning and the end of the period, yet United Airlines shows a drop in passenger mood on January 4th, the day when the bomb cyclone actually hit the East coast. Though all three airlines have hubs in New York, United Airlines is the only airline with a hub at Newark International Airport (EWR) and not John F. Kennedy International Airport (JFK) nor LaGuardia Airport (LGA), which were both closed during the bomb cyclone, meaning that United Airlines probably had more dissatisfied passengers to handle on site during these extreme weather conditions.

Airports passenger map

After training the Random Forest models, it is possible to search for the most important features within the 8,484 initial features for each airport. This is achieved by using the Mean Decrease Impurity measure defined by Breiman in [START_REF] Breiman | Classification and regression trees[END_REF] and normalizing the obtained feature importances so that the sum of all feature importances is equal to one. Table 3 shows the ten features with the highest feature importances for predicting the number of delayed departing flights in ATL. Besides date related features, four of the top ten features are related to the volume of tweets containing delay keywords. Once the features gathering 99% of the total importance for estimating the number of delayed flights are extracted, it is possible to group these features per origin in order to gain some insight on how airports are related from a passenger perspective. For example, once the most important features for estimating the number of delays at ATL are extracted, it is possible to count how many of these features are issued from tweets gathered using the handle of John F. Kennedy International Airport (JFK). Figure 7 shows how ATL is connected to the other airports from this perspective. The larger the link between ATL and another airport, the more features were kept among the features gathering 99% of the total importance for estimating the number of departing delayed flights at ATL. Interestingly, this airport graph is different from the graph built from the actual BTS values. Figure 8 shows how ATL is connected to the other airports using the number of delayed departing flights from ATL. For example, although there are many delayed flights departing to Florida, few features from Floridan airports are kept. The opposite observation can be made regarding Portland (PDX): there were less than a thousand delayed flights from ATL to PDX, yet features from PDX were kept. This example illustrates the possibility of creating a yearly review of airport relationship from a passenger point of view. Future studies should investigate more thoroughly the possible correlation and relation between the passenger connection map and the delay connection map.

CONCLUSION

This paper aimed at investigating further the use of the social media Twitter as an estimator of the US Air Transportation system. Exploiting both raw volume information as well as different levels of content information within the Twitter stream enables to accurately estimate for each airport the number of flights departing with a delay and the number of flights arriving with a delay greater than fifteen minutes. This passenger-based estimation yields a better estimation performance for a majority of airports compared to using a state-of-theart and off-the-shelf forecasting tool on the flight-centric data alone. Moreover, the methods used to extract relevant features from this passenger-centric data-source can be used to gain additional real-time insight on how passengers relate to the Air Transportation system.

This study confirmed that information contained in passenger-centric datasets are useful for a better understanding of the different stakeholders within the air transportation system, and have the added benefit of being more readily and publicly available than flight centric datasets. Future studies should focus on analyzing cases when the estimation is less accurate, implying differences between the handling of passengers and that of planes. Another direction of study considered is to validate this method to other countries or regions (e.g. the European Union) where sufficient flight-centric data is available.

Figure 1 :

 1 Figure 1: Number of tweets over the number of flights during the year 2017 for the considered airlines and airports

  (7 airlines and 33 airports giving 40 features): Num_tweets_handle • Hourly volume of keyword-related tweets for each search handle (40x7 features): Num_tweets_key-word_handle • Hourly average of tweets' sentiment (40x2 features): Mean_sent_method_handle • Hourly average of topic distribution for each search handle (40x100 features): Mean_topic_handle • Hourly standard deviation of topic distribution for each search handle (40x100 features): Std_topic_handle • Month of the year, Day of the month, Day of the week and Hour in the day (4 features)

Figure 2 :

 2 Figure 2: Comparison of the mean absolute errors per airport for the trained regressors for the estimation of the number of delayed departing flights. The standard deviation of the BTS value on the training set is included for comparison.

Figure 3 :

 3 Figure 3: Comparison of the mean absolute errors per airport for the trained regressors for the estimation of the number of flights arriving with a delay greater than 15 minutes. The standard deviation of the BTS value on the training set is included for comparison.

Figure 4 :

 4 Figure 4: Comparison of the R 2 scores per airport for the trained regressors for the estimation of the number of delayed departing flights

Figure 5 :

 5 Figure 5: Predicted number of delayed departing flights at ATL by the trained regressor over the period January 12th, 2018 to January 16th, 2018. The actual number of delayed departing flights is indicated for comparison.

Figure 6 :

 6 Figure 6: Average passenger sentiment with respect to three major airlines over the period January 2nd, 2018 to January 6th, 2018, corresponding to a bomb cyclone hitting in the North-East of the US.

Figure 7 :

 7 Figure 7: Map of feature links between Atlanta airport (ATL) and the other airports for estimating the number of delayed departing flights. The larger the link, the more features were kept among the features gathering 99% of the total importance for estimating the number of departing delayed flights at ATL.

Figure 8 :

 8 Figure 8: Map of delay links between Atlanta airport (ATL) and the other airports. The larger the link, the more flights departed with a delay during 2017 from ATL towards the connecting airport. Only links with more than 1000 delayed flights in 2017 were considered.

Table 1 :

 1 Twitter handles used for gathering tweets

Table 2 :

 2 Emoji sentiment association

	Category	Emojis
	Positive	":

Table 3 :

 3 Top ten features for predicting the number of delayed departing flights at ATL

	Rank	Feature	Rank	Feature
	1	Hour	6	DayOfMonth
	2	Month	7	delay_@SouthwestAir
	3	DayOfWeek	8	num_ATL
	4	delay_@Delta	9	delay_JFK
	5	delay_ATL	10	mean_63_BWI

ACKNOWLEDGMENTS

The authors would like to thank Nikunj Oza from NASA-Ames, the BDAI team from Verizon Media in Sunnyvale as well as the French government for their financial support.

COPYRIGHT

Copyright Statement

The authors confirm that they, and/or their company or institution, hold copyright of all original material included in their paper. They also confirm they have obtained permission, from the copyright holder of any third party material included in their paper, to publish it as part of their paper. The authors grant full permission for the publication and distribution of their paper as part of the EIWAC2019 proceedings or as individual off-prints from the proceedings.


Keywords: theorem

Optimal consumption in the discrete time 2 Optimal consumption and investment in discrete time models 3 Optimal consumption in continuous time 4 Optimal consumption and investment in continuous time models 4.1 Financial market model . . . . . . . . . . . . .

Optimal consumption in the discrete time

We consider the optimal consumption problem on the interval [0, N ], where N < ∞ is the horizon, i.e. we consider a consumer who chooses the amount of its consumption c k ≥ 0 at the time moment 0 ≤ k ≤ N . Describe the dynamics of his wealth. Let x 0 = x > 0 be the initial wealth. So, at the time k = 0 we obtain the constraint on the consumption:

0 ≤ c 0 ≤ x .
Then he places the rest x -c 0 in a risk-free asset with the interest rate r ≥ 0. So, the next day his wealth is equal to

x 1 = R(x -c 0 ) and R = 1 + r .
At the time k = 1 the consumer can consume 0 ≤ c 1 ≤ x 1 and its capital at the time k = 2 is given by the following equation

x 2 = R(x 1 -c 1 ) .
Similarly, when k ≥ 1 the capital x k is written as

x k = R(x k-1 -c k-1 ) , x 0 = x , (1.1)
where the consumption vector c = (c 0 , . . . , c N -1 ) holds the following constraints

0 ≤ c k ≤ x k , 0 ≤ k ≤ N -1 . (1.2)
Recall that the prime means transpose.

Let U : R + → R be a utility function and h : R + → R be a heritage function. It is assumed that the functions U and h are increasing and continuous. In this case the objective function is

J 0 (x, c 0 , . . . , c N -1 ) = N -1 k=0 U (c k ) + h(x N ) . (1.3)
The problem is to find a sequence c * = (c * 0 , . . . , c * N -1 ) that maximizes the objective function (1.3) under the constraints (1.1) -(1.2), i.e. to solve the following extremal problem J * 0 (x) = max (c 0 ,...,c N -1 ) ∈C 0 (x)

J 0 (x, c 0 , . . . , c N -1 ) ,

where C 0 (x) ⊂ R N is the set of the vectors c = (c 0 , . . . , c N -1 ) which satisfy the conditions (1.1) -(1.2).

To resolve this problem we will use the Bellman dynamic programming method [1]. Now we describe this approach. First define the objective functions on intervals [n, N ] for 0 ≤ n ≤ N -1. We have

J n (x, c n , . . . , c N -1 ) = N -1 k=n U (c k ) + h(x N ) . (1.5)
Here the consumptions c n , . . . , c N -1 satisfy the following conditions

0 ≤ c k ≤ x k , n ≤ k ≤ N -1 , (1.6)
where x n = x and for k > n we get

x k = R(x k-1 -c k-1 ) . (1.7)
Denote by C n (x) the set in R N -n of the vectors (c n , . . . , c N -1 ) which hold the constraints (1.6) -(1.7). For 0 ≤ n ≤ N -1 we set

J * n (x) = max (c n ,...,c N -1 ) ∈C n (x)
J n (x, c n , . . . , c N -1 ) .

(1.8)

For n = N put J * N (x) = h(x) .

(1.9)

To resolve the problems (1.8) we use the following deterministic dynamic programming principle.

Theorem 1.1. (Deterministic dynamic programming principle)

The functions (J * n (x)) 0≤n≤N -1 satisfy the following Bellman equa-tions:

J * n (x) = max 0≤c≤x U (c) + J * n+1 (φ(x, c)) , (1.10)
where φ(x, c) = R(x -c).

Proof. First, note that all the sets

C n (x) are bounded in R N -n .
Therefore all the functions (J * n (x)) 0≤n≤N -1 are fined, i.e. the problem (1.8) has a solution for any 0 ≤ n ≤ N -1. Moreover, by the definition, we obtain that

J * n (x) = max (c n ,...,c N -1 ) ∈C n (x) N -1 k=n U (c k ) + h(x N ) = max 0≤c n ≤x max (c n+1 ,...,c N -1 ) ∈Q(x,c n ) U (c n ) + N -1 k=n+1 U (c k ) + h(x N ) = max 0≤c≤x U (c) + max (c n+1 ,...,c N -1 ) ∈Q(x,c) J n+1 (φ(x, c), c n+1 , . . . , c N -1 ) , where Q(x, c) = C n+1 (φ(x, c)) ⊂ R N -n-1
. This directly implies the equation (1.10).

To find the optimal objective functions (J * (n, x)) 0≤n≤N -1 we use the back recursion. It means that we start from n = N -1. We have

J * N -1 (x) = max 0≤c≤x (U (c) + h(φ(x, c))) . (1.11)
We denote by s N -1 (x) the solution of this problem. Similarly, by the recurrence we obtain the solution s n (x) of the problem (1.10) for 0 ≤ n ≤ N -1. Therefore, finally we obtain s 0 (x), where

x > 0 is initial wealth. This means that c * 0 = s 0 (x). Then, to find the consumption for every n ≥ 1 we have to calculate the optimal wealths (x * n ) 1≤n≤N . We start with n = 1. We obtain

x * 1 = R(x -c * 0 ) .
This implies that c * 1 = s 1 (x * 1 ). This consumption gives us the wealth n = 2, i.e.

x * 2 = R(x * 1 -c * 1 )
and c * 2 = s 2 (x * 2 ) etc. So, by the recurrence

x * n = R(x * n-1 -c * n-1 ) and c * n = s n (x * n ) .
Finally, the solution of (1.4) is given by the vector

c * = (c * 0 , c * 1 , . . . , c * N -1 ) .
Let us consider now the following example.

Example. Let U (x) = h(x) = x γ and 0 < γ < 1. We start to resolve the Bellman equations (1.10) with n = N -1, i.e. the equation (1.11). In this case

J * N -1 (x) = max 0≤c≤x (c γ + (x -c) γ R γ ) .
For 0 < γ < 1 the solution is

s N -1 (x) = α N -1 x , α N -1 = 1 1 + R q and q = γ 1 -γ .
This means that

J * N -1 (x) = β N -1 x γ and β N -1 = α γ N -1 + (1 -α N -1 ) γ R γ .
(1.12)

We will resolve the Bellman equations (1.10) by induction. We assume that for some 1 ≤ n ≤ N -1 the solution has the following form:

s n (x) = α n x and J * n (x) = β n x γ (1.13)
for some positive constants α n and β n . To show that these formulas are true for all 0 ≤ n ≤ N -1 it suffices to verify this for n -1.

Using (1.13) in (1.10) one obtains that

J * n-1 (x) = max 0≤c≤x (c γ + J * (n, R(x -c))) = max 0≤c≤x (c γ + β n (x -c) γ R γ ) .
We find that s n-1 (x) = α n-1 x and

J * n-1 (x) = s n-1 (x) γ + β n R γ x -s n-1 (x) γ = β n-1 x γ ,
where

α n-1 = 1 1 + R q β 1+q n and β n-1 = α γ n-1 + β n R γ (1 -α n-1 ) γ . (1.14)
This implies that (1.13) is true for all 0 ≤ n ≤ N -1. So, to find the optimal consumptions it remains to calculate the wealth (x * n ) 0≤n≤N . We start by n = 0. Our initial wealth is x > 0. Therefore,

c * 0 = s 0 (x) = α 0 x .
Then, for n = 1 we get that

x * 1 = R(x -c * 0 ) = R(1 -α 0 )x and c * 1 = s 1 (x * 1 ) = α 1 x * 1 .
By the same way for any n ≥ 1

x * n = R(x n-1 -c * n-1 ) = R(1-α n-1 )x * n-1 and c * n = s n (x * n ) = α n x * n .
This implies that for

1 ≤ n ≤ N x * n = xR n n j=1 (1 -α j-1 ) and c * n = xα n R n n j=1
(1 -α j-1 ) .

Using here the recurrence equations (1.14) with the values β N -1

given in (1.12) we can calculate all the coefficients (α j ) 0≤j≤N -2 .

Exercises

1. Let the utility function U (x) = ln x and the heritage function

h(x) = ln x.
(a) Find optimal consumption for N = 1, N = 2 and any

N ≥ 3.
(b) For which the interest rates r ≥ 0 optimal terminal wealth will be greater than initial wealth x for the cases N = 1, N = 2 and any N ≥ 3?

2. Let U (x) = h(x) = 1 -e -x .
Find optimal consumptions for N = 1. Find optimal consumptions for N = 2 and any N ≥ 3 in the case r = 0.

3. Find the optimal consumptions for the following cases:

(a) U (x) = x and h(x) = 0; N = 1, N = 2 and N = 3.

(b) U (x) = 0 and h(x) = 1 -e -x ; N = 1, N = 2 and N = 3.

2 Optimal consumption and investment in discrete time models

Let (Ω, P) be a finite probability space with Ω = {ω 1 , . . . , ω k }.

In this space we consider a (B, S) -financial market with two assets on the time interval [0, N ]:

   (B) B n = (1 + r)B n-1 , (S) S n = (1 + ρ n )S n-1 , 1 ≤ n ≤ N , (2.1)
where B 0 > 0 and S 0 > 0 are initial values. Assume that interest rates of the risk asset (ρ n ) 1≤n≤N is sequence of the independent identically distributed ( i.i.d.) random variables with values {a 1 , . . . , a m }, where a 1 < . . . < a m and P(ρ n = a i ) = p i > 0. To avoid trivial strategies we will assume that

-1 < a 1 < r < a m . (2.2)
We consider an investor with an initial capital x 0 = x > 0. On the date n = 0 we can consume a quantity 0 ≤ c 0 ≤ x. Moreover, it distributes the remaining x -c 0 between the assets (B) and (S)

as follows. The part θ 0 (x -c) is placed in the asset (S) and (1 -

θ 0 )(x -c) in (B)
. This means that on the date n = 1 (the next day)

the wealth x 1 is equal to

x 1 = R(1 -θ 0 )(x 0 -c 0 ) + (1 + ρ 1 ) θ 0 (x 0 -c 0 ) = τ 1 (θ 0 ) (x 0 -c 0 ) , where τ 1 (θ) = R + θ(ρ 1 -r) and R = 1 + r.
Let now x n the wealth on the date n. We denote by 0 ≤ c n ≤ x n the consumption on the date n and by θ n the part of the wealth invested in risk asset (S) on the same date. Similarly, we obtain

x n+1 = τ n+1 (θ n ) (x n -c n ) , (2.3)
where τ n (θ) = R + θ(ρ n -r) and R = 1 + r.

We need to introduce now an information structure, i.e. we need to describe all possible random variables which can be used

for financial strategies (v n ) 0≤n≤N , where v n = (c n , θ n ). To this end we set F 0 = {∅, Ω} and F n = σ{ρ 1 , . . . , ρ n } for 1 ≤ n ≤ N . The sequence (F n ) 0≤n≤N is called filtration. Assume that the random variables v n are measured with respect to F n for any 0 ≤ n ≤ N -1.
In this case we say that the sequence of random variables

(v n ) 0≤n≤N -1 is adapted. Moreover, the variables (x n ) 1≤n≤N should be positive. It means that for any 0 ≤ k ≤ N -1 θ k ∈ Θ and 0 ≤ c k ≤ x k , (2.4)
where

Θ = {θ ∈ R : inf k≥1 τ k (θ) ≥ 0 a.s.}
The variables (x k ) 1≤k≤N are defined by (2.3) with x 0 = x. We

denote by V 0 (x) the set of all adapted sequences v = (v n ) 0≤n≤N -1
which satisfy the constraints (2.4). Generally, we denote by

V n (x)
the set of all adapted sequences (v k ) n≤k≤N -1 which satisfy the constraints (2.4) for n ≤ k ≤ N -1 with x n = x. We define now the cost function on the time interval [n, N ] as

J n (v n , . . . , v N -1 ) = N -1 k=n U (c k ) + h(x N ) . (2.5)
Our goal is to find a financial strategy v * = (v * k ) 0≤k≤N -1 in V 0 (x) which maximizes the conditional expectation of the function

J 0 (x, v 0 , . . . , v N -1 ), i.e. J * 0 (x) = max (v 0 ,...,v N -1 ) ∈V 0 (x) E J 0 (v 0 , . . . , v N -1 )|x 0 = x . (2.6)
For this problem we will use the Bellman stochastic dynamical programming principle. According to this principle we need to find optimal strategies on the intervals

[n, N ] for 0 ≤ n ≤ N -1, i.e. J * n (x) = max (v n ,...,v N -1 ) ∈V n (x) E J n (v n , . . . , v N -1 )|x n = x . (2.7)
Theorem 2.1. (The stochastic dynamic programming principle)

The optimal cost functions (J * n (x)) 0≤n≤N -1 satisfy the following Bellman equations :

     J * n (x) = max 0≤c≤x , θ∈Θ (U (c) + J n+1 (x, c, θ)) , J * N (x) = h(x) , (2.8)
where

J n+1 (x, c, θ) = E J * n+1 (x n+1 )|x n = x, v n = (c, θ) .
Example. Let us resolve the Bellman equations (2.8) for U (x) = h(x) = x γ with 0 < γ < 1. First, note that

J n+1 (c, θ) = EJ * n+1 (τ 1 (θ)(x -c)) ,
where τ 1 (θ) = R + θ(ρ 1 -r). Thus, in this case the equation (2.8)

has the following form

J * n (x) = max 0≤c≤x max θ∈Θ n (U (c) + E J * n+1 (τ 1 (θ)(x -c)) .
Moreover, it is easy to see that

Θ = [θ min , θ max ] , (2.9)
where

θ min = - R a m -r and θ max = R r -a 1 .
Therefore, we obtain that for 0 ≤ n ≤ N -1

J * n (x) = max 0≤c≤x max θ min ≤θ≤θ max U (c) + E J * n+1 (τ 1 (θ)(x -c)) (2.10) and J * N (x) = h(x) = x γ .
We start to study these equations from n = N -1. We have

J * N -1 (x) = max 0≤c≤x max θ min ≤θ≤θ max (c γ + (x -c) γ L(θ)) = max 0≤c≤x (c γ + (x -c) γ ν) .
(2.11) where ν = max

θ min ≤θ≤θ max L(θ) and L(θ) = E (τ 1 (θ)) γ .
To study this problem one needs to solve the following equation

f (θ) = E ρ 1 -r (R + θ(ρ 1 -r)) 1-γ = m j=1 a j -r R + θ(a j -r) 1-γ p j = 0 .
(2.12)

One can check directly, that the solution of the problem (2.11) is given by

(s N -1 (x) , κ N -1 (x)) with s N -1 (x) = 1 1 + µ x and κ N -1 (x) = κ * ,
where µ = ν q , q = 1/(1 -γ) and κ * is the unique root of the equation (2.12) on the interval (θ min , θ max ). We set

s N -1 (x) = α N -1 x , where α N -1 = 1 1 + µ . (2.13)
Then, we obtain that

J * N -1 (x) = β N -1 x γ , β N -1 = α γ N -1 + ν(1 -α N -1 ) γ . (2.14)
Now, by induction we assume that for 0 ≤ n ≤ N -1

s n (x) = α n x , κ n (x) = κ * and J * n (x) = β n x γ . (2.15)
Using these conditions in (2.10), we find

J * n-1 (x) = max 0≤c≤x max θ min ≤θ≤θ max (c γ + β n (x -c) γ L(θ)) .
In the same way as for the case n = N -1 we get

s n-1 (x) = α n-1 x , κ n-1 (x) = κ * and J * n-1 (x) = β n-1 x γ , (2.16)
where

α n-1 = 1 + β γ 1 n µ -1 and β n-1 = α γ n + β n ν(1 -α n ) γ . (2.17)
Through these equations we obtain the solutions of the Bellman equations (2.10) for all 0 ≤ n ≤ N -1.

Now we have to calculate the optimal strategies (v * 0 , . . . , v * N -1 ), where v j = (c * j , θ * j ). Note that according to (2.16) for all 0 ≤ n ≤ N -1 we have θ * n = κ * . Then, taking into account that x 0 = x, we obtain c * 0 = s(0, x) = α 0 x. So, for n = 1 we get

x * 1 = τ 1 (κ * )(x -c * 0 ) = τ 1 (κ * )(1 -α 0 )x and c * 1 = s(1, x * 1 ) = α 1 x * 1 .
Similarly, for n > 1:

x * n = τ n (κ * )(1-α n-1 )x * n-1 and c * n = s(n, x * n ) = α n x *
n . This implies that

x * n = x n k=1 τ k (κ * )(1 -α k-1 ) and c * n = α n x n k=1 τ k (κ * )(1 -α k-1 )
.

Using now the recurrence equalities (2.17) we can calculate the optimal consumption for any n ≥ 1.

Exercises

1. Find the optimal strategies for U (x) = h(x) = ln x, N = 1, N = 2 and for any N ≥ 3.

2. Find the optimal strategy for N = 1, U (x) = x γ and h(x) = x with 0 < γ < 1.

3. Find the optimal strategy for N = 1, U (x) = x and h(x) = x γ with 0 < γ < 1.

Optimal consumption in continuous time

In this chapter we will consider a financial market in continuous time consisting of a riskless asset Let X t be a wealth at the moment t. To obtain an equation for the wealth we divide the interval [0, t] by the sequence (t k ) 0≤k≤n defined as

Ḃ(t) = r B(t) B 0 = 1 , 0 ≤ t ≤ T . ( 3
t k = kt n .
So, for k ≥ 1 by analogy with the discrete time case is obtained

X t k = 1 + r t k-1 ∆ + d k X t k-1 - t k t k-1 c t dt , (3.3)
where

∆ = t k -t k-1 = t n and d k = r t k t k-1 B t -B t k-1 dt .
It is clear that

lim n→∞ 1 ∆ max 1≤k≤n |d k | = 0 .
One can rewrite the equation (3.3) as

X t k -X t k-1 = r t k-1 X t k-1 ∆ - t k t k-1 c t dt + h k ∆ with h k = d k ∆ X t k-1 - t k t k-1 c t dt -r t k-1 t k t k-1 c t dt .
It is easy to deduce that

lim n→∞ 1 ∆ max 1≤k≤n |h k | = 0 .
So, for any n ≥ 1,

X t = x + n k=1 r t k-1 X t k-1 ∆ - t 0 c u du + n k=1 h k .
Then, passing to the limit when n → ∞, we obtain that for any t > 0

X t = x + t 0 r u X u du - t 0 c u du ,
where x > 0 is an initial endowment. This means that

Ẋt = r X t -c t , X 0 = x . (3.4)
Now we define the cost function as

J(x, c) := T 0 U (c t )dt + h(X T ) . (3.5)
Our goal is to resolve the following optimization problem

J * (x) = max c J(x, c) , (3.6)
where c = (c t ) 0≤t≤T . For this problem we will use the Bellman dynamic programming method. To this end we set

J(t, x, c) := T t U (c v )dv + h(X T )
and

J * (t, x) = max c J(t, x, c) = J(t, x, c * ) .
Then we introduce the Hamiltonien function as

H(t, x, z) := sup ϑ≥0 H 0 (t, x, z, ϑ) (3.7)
with

H 0 (t, x, z, ϑ) := z (r x -ϑ) + U (ϑ) .
To find the solution for the problem (3.6), one needs to study the following Hamilton-Jacobi-Bellman equation

       z t (t, x) + H(t, x, z x (t, x)) = 0 , t ∈ [0, T ] , z(T, x) = h(x) , x > 0 . (3.8)
Here z t is the partial derivative of z with respect to t. We use the same notation for all partial derivatives.

We assume that

C 1 ) There exists a solution z(t, x) for the equation (3.8).

C 2 ) There exists a measurable function θ :

[0, T ] × (0, ∞) → R + such that H(t, x, z x (t, x)) = H 0 (t, x, z x (t, x), θ(t, x))
for all 0 ≤ t ≤ T and x > 0.

C 3 ) There exists a solution of the following equation

dX * t = (r x * t -θ(t, X * t )) dt , X * 0 = x . (3.9)
The solution of the problem (3.6) is given in the following theorem.

Theorem 3.1. Assume that the conditions C 1 )-C 3 ) hold. Then for all 0 ≤ t ≤ T and for all x > 0

z(t, x) = J * (t, x) = J * (t, x, c * ) ,
where the optimal consumption c * = (c * t = ϑ(t, X * t )) 0≤t≤T is defined in C 3 ).

Proof. Let c = (c t ) 0≤t≤T be a nonnegative integrable function. We have

z(T, X T ) -z(t, X t ) = T t z t (v, X v ) + z x (v, X v ) Ẋv dv .
Using the equation (3.4) and the conditions X t = x and z(T, x) = h(x), we obtain that

z(t, x) = J(t, x, c) - T t (z t (s, X s ) + H 1 (s, X s , c s )) ds , (3.10)
where H 1 (s, x, ϑ) = H 0 (t, x, z x (t, x), ϑ). This means that z(t, x) ≥ J(t, x, c) for all c = (c t ) 0≤t≤T . Therefore, we get z(t, x) ≥ J * (t, x).

Moreover, replacing c in (3.9) by c * we obtain z(t, x) = J(t, x, c * ).

Hence Theorem 3.1.

Exercises

1. Write the Hamiltoniens for the following cases:

(a) U (x) = x γ with 0 < γ < 1.

(b) U (x) = ln x. (c) U (x) = 1 -e -x .
2. Find the optimal consumptions for the following problems : Here (w t ) 0≤t≤T is the Wiener process and the coefficients r, µ and σ > 0 are fixed parameters. We recall, that r is the interest rate, µ is the increasing coefficient and σ is volatility. Moreover, (F w t ) 0≤t≤T is the filtration generated by the Wiener process (w t ) 0≤t≤T , i.e.

(a) U (x) = h(x) = x γ with 0 < γ < 1. (b) U (x) = x γ with 0 < γ < 1 and h(x) = 0. (c) U (x) = h(x) = ln x.
F w t = σ{w u , 0 ≤ u ≤ t}.
Let now φ t ∈ R and ϕ t ∈ R be the amounts placed at the time moment t in the riskless and the risk assets respectively. In this case financial strategy is defined as the random process (φ t , ϕ t ) 0≤t≤T with the values in R 2 progressively measurable with respect to the filtration (F w t ) 0≤t≤T . The process

X t = φ t B(t) + ϕ t S(t) , 0 ≤ t ≤ T ,
is called wealth process (or the value of the portfolio) at the time moment t.

Moreover let for any 0 ≤ t ≤ T the function c t be a consumption rate on the time interval [0, t], i.e. it's a non-negative process adapted to the filtration (F w t ) 0≤t≤T and almost sure integrable on the interval [0, T ], i.e.

T 0 c t dt < ∞ a.s.
The integral t 0 c s ds is the consumption process during the time interval [0, t]. A financial strategy (φ t , ϕ t ) 0≤t≤T with consumption (c t ) 0≤t≤T is called self-financing strategy if the wealth process satisfies the following stochastic equation

X t = x + t 0 φ u dB(u) + t 0 ϕ u dS(u) - t 0 c u du , (4.2)
where x > 0 is the initial endowment. We will work with the relative quantities in relation to the wealth process, i.e. we set

π t := ϕ t S(t) X t and v t := c t X t .
Taking this into account, we can rewrite the equation for X t as

dX t = (r + σπ t θ -v t )X t dt + σπ t X t dw t , X 0 = x , (4.3)
where θ = (µ -r)/σ.

Let us describe now the set of control processes ν = (ν t ) 0≤t≤T

with ν t = (π t , v t ). Definition 4.1. The control process ν = (ν t ) 0≤t≤T = (π t , v t ) 0≤t≤T
with the values in R × R + is called admissible if it is progressively measurable with respect to the filtration (F w t ) 0≤t≤T and such that

π 2 T = T 0 π 2 t dt < ∞ , T 0 v t dt < ∞ a.s. (4.4)
and the equation (4.3) has an unique strong positive solution on the time interval [0, T ]. We denote by V the set of all admissible control processes.

To emphasize that a wealth process corresponds to a control process ν we will write X ν . Now, for the initial wealth x > 0 and for a process of control (ν t ) t≥0 in V, we introduce the objective function as

J(x, ν) := E x T 0 U (v t X ν t )dt + h(X ν T ) , (4.5)
where E x is the conditional expectation for X ν 0 = x, U : R + → R and h : R + → R are utility functions.

The goal of investment and consumption is to maximize the objective function (4.5), i.e. one needs to find

max ν∈V J(x, ν) . (4.6)
To solve this problem, we will use the Bellman dynamic stochastic programming method which is based on the verification theorems tool for stochastic differential equations. A special version of such theorem is given in the next section.

Verification theorem

In this section we prove the verification theorem from [START_REF] Klüppelberg | Optimal consumption and investment with bounded Capital-at-Risk for power utility functions[END_REF] which will apply to the problem (4.6). We consider the problem of stochas-

lution (X ν t ) 0≤t≤T such that T 0 |a(t, X ν t , ν t )| + |b(t, X ν t , ν t )| 2 dt < ∞ a.s. (4.8)
and

E T 0 (U (t, X ν t , ν t )) -dt + (h(X ν T )) -< ∞ , (4.9)
where (a) -= -min(0, a).

We denote by V the set of all admissible control processes. Now let's define the objective functions, by setting

J(t, x, ν) := E t,x T t U (s, X ν s , ν s ) ds + h(X ν T ) , 0 ≤ t ≤ T , (4.10)
where E t,x is the conditional expectation for X ν t = x.

Our goal is to solve the following optimization problems

J * (t, x) := sup ν∈V J(t, x, ν) , 0 ≤ t ≤ T . (4.11)
To this end, let's introduce the Hamiltonian defined by

H(t, x, z 1 , z 2 ) := sup ϑ∈K H 0 (t, x, z 1 , z 2 , ϑ) (4.12)
and

H 0 (t, x, z 1 , z 2 , ϑ) := a(t, x, ϑ) z 1 + 1 2 |b(t, x, ϑ)| 2 z 2 + U (t, x, ϑ) .
To find a solution to the problem (4.11), we must study the Hamilton-Jacobi-Bellman (HJB) equation

     z t (t, x) + H(t, x, z x (t, x), z xx (t, x)) = 0 , t ∈ [0, T ] , z(T, x) = h(x) , x > 0 . (4.13)
Here z t denotes the partial derivative of z over t. We will use the same notations for all partial derivatives. Assume that H 1 ) There exists a classical solution of the problem (4.13), i.e. there exists [0, T ] × (0, ∞) → R function z which is continuously differentiable over the variable t and two times continuously differentiable over x and satisfies the equation (4.13).

H 2 ) There exists a [0, T ] × (0, ∞) → K mesurable function θ such that for any t ∈ [0, T ] and for any x > 0

H(t, x, z x (t, x), z xx (t, x)) = H 0 (t, x, z x (t, x), z xx (t, x), θ(t, x)) .
H 3 ) There is an unique strong positive solution for the following stochastic differential equation

dX * t = a * (t, X * t ) dt + b * (t, X * t ) dw t , X * 0 = x , (4.14)
where a * (t, x) = a(t, x, θ(t, x)) and b * (t, x) = b(t, x, θ(t, x)). Moreover, the optimal control process ν * = (ν * t ) 0≤t≤T with ν * t = θ(t, X * t ) belongs to the set V and for any 0 ≤ t ≤ T and x > 0

E t,x sup t≤s≤T |z(s, X * s )| < ∞ . (4.15)
Now we show that the process ν * solves the problems (4.11).

Theorem 4.1. Assume that V = ∅ and the conditions H 1 ) -H 3 )

hold. Moreover, assume that

z * = inf 0≤t≤T inf x>0 z(t, x) > -∞ . (4.16)
Then for any 0 ≤ t ≤ T and for any x > 0

z(t, x) = J * (t, x) = J * (t, x, ν * ) ,
where the strategy ν * is defined in the condition H 3 ).

Proof. Let now ν ∈ V and X ν be the wealth process with X ν t = x for a some fixed 0 ≤ t < T . Next, let us define the following stopping times

τ n = inf s ≥ t : s t |b(u, X ν u , ν u )| 2 z 2 x (u, X ν u ) du ≥ n ∧ T .
Note that the condition (4.8) implies that τ n → T as n → ∞ a.s. Taking into account that the function z(•, •) and the process (X ν t ) 0≤t≤T are continuous, we get

lim n→∞ z(τ n , X ν τ n ) = z(T, X ν T ) = h(X ν T ) a.s. (4.17)
Then, according to Ito's formula, we have

z(t, x) = τ n t U (s, X ν s , ν s ) ds + z(τ n , X ν τ n ) (4.18) - τ n t (z t (s, X ν s ) + H 1 (s, X ν s , ν s )) ds -M n ,
where

M n = τ n t b(u, X ν u , ν u ) z x (u, X ν u ) dw u and H 1 (s, x, ϑ) = H 0 (t, x, z x (t, x), z xx (t, x), ϑ) .
The definition of the stopping time τ n implies that E t,x M n = 0.

Therefore, we obtain the following inequality

z(t, x) ≥ E t,x τ n t U s ds + E t,x z n , (4.19)
where

U s = U (s, X ν s , ν s ) and z n = z(τ n , X ν τ n
). Consider the first term in this inequality. We have

E t,x z n ≥ E t,x lim n→∞ z(τ n , X ν τ n ) = E t,x h(X ν T ) .
Therefore, going to the limit as n → ∞ in (4.19) we get z(t, x) ≥ J * (t, x). In the same way, replacing ν in (4.18) by ν * we obtain the following equality

z(t, x) = E t,x τ n t U (s, X * s , ν * s ) ds + E t,x z(τ n , X * τ n ) .
Then in view of the condition (4.15) and the dominated convergence theorem we get

lim n→∞ E t,x z(τ n , X * τ n ) = E t,x lim n→∞ z(τ n , X * τ n ) = E t,x h(X * T ) .
This means that

z(t, x) = lim n→∞ E t,x τ n t U (s, X * s , ν * s ) ds + lim n→∞ E t,x z(τ n , X * τ n ) = E t,x T t U (s, X * s , ν * s ) ds + h(X * T ) = J(t, x, ν * ) .
So, using the inequality z(t, x) ≥ J * (t, x) we obtain z(t, x) = J * (t, x).

Hence Theorem 4.1.

In the next section we will see how to apply the verification theorem to the problem (4.6).

Optimal investment and consumption for the power utility functions

In this section we consider the problem (4.6) in the case where

U (x) = x γ 1 and h(x) = x γ 2 with 0 < γ 1 , γ 2 < 1.
Let us first define the following functions

A 1 (t) = γ q 1 1 T t e s t β 1 (u)du ds and A 2 (t) = γ q 2 2 e T t β 2 (u)du , (4.20)
where q i = (1 -γ i ) -1 and

β i (t) = (q i -1) r + q i 2 |θ| 2 .
Then for any 0 ≤ t ≤ T and x > 0, we define the function g(t, x) > 0 as

A 1 (t) g -q 1 (t, x) + A 2 (t) g -q 2 (t, x) = x . (4.21)
In addition, we set

p(t, x) = q 1 A 1 (t) g -q 1 (t, x) + q 2 A 2 (t) g -q 2 (t, x) .
Now we define the functions which will be used to check the condition H 2 ), i.e. we set

π(t, x) = p(t, x) σx θ and v(t, x) = (γ 1 ) q 1 g -q 1 (t, x) x . (4.22)
The solution in this case is given in the following theorem.

Theorem 4.2. The optimal value of J(x, ν) is given by

max ν∈V J(x, ν) = J(x, ν * ) = A 1 (0) γ 1 g 1-q 1 (0, x) + A 2 (0) γ 2 g 1-q 2 (0, x) ,
where optimal control ν * t = (π * t , v * t ) has the following form π * t = π(t, X * t ) and v * t = π(t, X * t ). The optimal wealth process (X * t ) 0≤t≤T satisfies the following stochastic equation

dX * t = a * (t, X * t )dt + b * (t, X * t )dw t , (4.23)
where

X * 0 = x, b * (t, x) = p(t, x) θ and a * (t, x) = r x + p(t, x) |θ| 2 -γ q 1 1 g -q 1 (t, x) .
Proof. To prove this theorem we will use the theorem 4. In addition, the utility functions are defined as U (t, x, ϑ) = v γ 1 x γ 1 and h(x) = x γ 2 .

To check the conditions H 1 ) -H 3 ) calculate the Hamiltonian (4.12) in our case. We have

H(t, x, z 1 , z 2 ) = sup ϑ∈R d ×[0,∞)
H 0 (t, x, z 1 , z 2 , ϑ) , where H 0 (t, x, z 1 , z 2 , ϑ) = (r + σyθ)x z 1 + 1 2

x 2 σ 2 y 2 z 2 + v γ 1 x γ 1 -xv z 1 .

Therefore for z 2 ≤ 0 we obtain

H(t, x, z 1 , z 2 ) = H 0 (t, x, z 1 , z 2 , ϑ 0 ) = r x z 1 + z 2 1 θ 2 2|z 2 | + 1 q 1 z 1 γ 1 γ 1 γ 1 -1
, where ϑ 0 = ϑ 0 (t, x, z 1 , z 2 ) = (y 0 (t, x, z 1 , z 2 ), v 0 (t, x, z 1 , z 2 )) = (y 0 , v 0 ) and

y 0 = z 1 σx|z 2 |
θ and v 0 = γ 1 z 1 q 1

x -1 . (4.24)

Now we have to solve the problem (4.13) which has for our case the following form

       z t (t, x) + r x z x (t, x) + z 2 x (t, x) θ 2 2|z xx (t, x)| + 1 q 1 z x (t, x) γ 1 γ 1 γ 1 -1 = 0 , z(T, x) = x γ 2 .
We check directly that the solution for this problem (see [START_REF] Klüppelberg | Optimal consumption and investment with bounded Capital-at-Risk for power utility functions[END_REF]) for 0 ≤ t < T and x > 0 is z(t, x) = A 1 (t) γ 1 g 1-q 1 (t, x) + A 2 (t) γ 2 g 1-q 2 (t, x) , (4.25) which fulfills the condition (4.16) for z * = 0. Therefore, we have H(t, x, z x (t, x), z xx (t, x)) = H 0 (t, x, z x (t, x), z xx (t, x), θ(t, x)) ,

where the function θ(t, x) = (π(t, x), v(t, x)) is defined in (4.22). In [START_REF] Klüppelberg | Optimal consumption and investment with bounded Capital-at-Risk for power utility functions[END_REF] it is shown that the conditions H 1 ) -H 3 ) hold for this case. So, Corollary 4.1. Let U (x) = h(x) = x γ and 0 < γ < 1. Then the optimal value of J(x, ν) in (4.6) is J * (x) = max ν∈V J(x, ν) = J(x, ν * ) = x γ g γ q q,T + g q γ (T )

1/q
, where the optimal control process ν * t = (π * t , v * t ) for 0 ≤ t ≤ T has the following form

π * t = θ σ(1 -γ)
and v * t = g q γ (t) g q γ (T ) + T t g q γ (s) ds .

(4.26)

The optimal wealth process (X * t ) 0≤t≤T satisfies the following stochastic differential equation

dX * t = r -v * t + θ 2 1 -γ X * t dt + θ 1 -γ X * t dw t (4.27)
and X * 0 = x.

  ds represents the amount consumed on the interval [0, t].

(

  d) U (x) = ln x and h(x) = 0. (e) r = 0, U (x) = 1 -e -x and h(x) = x/e. 4 Optimal consumption and investment in continuous time models 4.1 Financial market model In this section we consider the financial Black-Scholes market consisting of a risk-free asset and a risk asset on the time interval [0, T ]. Their prices (B(t)) 0≤t≤T and (S(t)) 0≤t≤T for i = 1, . . . , d satisfy the following equations ) = r B(t)dt , B(0) = 1 ; dS(t) = S(t) (µdt + σ dw t ) , S(0) > 0 .

  1 for the stochastic differential equation (4.3). First we fixe ϑ = (y, v) with y ∈ R and v ∈ [0, ∞). For such ϑ the coefficients in the model (4.7) are defined as a(t, x, ϑ) = x ( r + σyθ -v) and b(t, x, ϑ) = σx y .

Theorem 4 . 1 implies

 41 Theorem 4.2. Now we consider the problem (4.6) for U (x) = h(x) = x γ with 40 0 < γ < 1. We set g γ (t) = exp γr + γθ 2 2(1 -γ) t .

Remark 4 . 1 .

 41 It should be noted that the strategy (4.26) -(4.27) is the well known Merton strategy obtained in[START_REF] Merton | Optimal consumption and portfolio rules in a continuous time model[END_REF].

tic control for the scalar Itô process on the time interval [0, T ]      dX ν t = a(t, X ν t , ν t ) dt + b(t, X ν t , ν t ) dw t ;

We assume that the control process ν = (ν t ) 0≤t≤T takes its values in the set K ⊆ R d and that the coefficients In this section we will modify the Definition 4.2 of admissibility for the equation (4.7) as follows.

Definition 4.2. The control process ν = (ν t ) 0≤t≤T is admissible if it is progressively measurable with respect to the filtration (F w t ) 0≤t≤T and the equation (4.7) has a unique strong positive so-

Conclusion

In this manual we considered the main optimization portfolio problems for the financial markets modeled by the stochastic differential and stochastic different equations. We considered the optimal consumption, the optimal investment and consumption problems in discrete time studied in [2], then we considered the same problems in continuous time using the approach developed in [START_REF] Karatzas | Methods of Mathematical Finance[END_REF]. For all these problems we used the stochastic calculus developed in [START_REF] Liptser | Statistics of Random Processes I. General Theory[END_REF][START_REF] Pchelintsev | Stochastic modelling for the financial markets. Part 1. Probabilistic tools[END_REF].

For all these problems we used the dynamical programming method proposed by Bellman for the optimization problem for the utilities functions of the integral forms [1]. For the Bellman principle in continuous time we studied the verification theorem method developed in [START_REF] Karatzas | Methods of Mathematical Finance[END_REF][START_REF] Klüppelberg | Optimal consumption and investment with bounded Capital-at-Risk for power utility functions[END_REF][START_REF] Pham | Smooth solutions to optimal investment models with stochastic volatilities and portfolio constraints[END_REF][START_REF] Touzi | Stochastic control problems, viscosity solutions and applications to finance[END_REF].


Najib Kacem
email: najib.kacem@femto-st.fr

Sébastien Hentz

Sébastien Baguet

Régis Dufour

Sensitive

Femto-St

Sensitive detection of nonlinear nanomechanical motion using capacitive signal down-mixing for resonant NEMS-based sensors

Introduction:

Nanoelectromechanical systems (NEMS) are emerging as strong candidates for a host of important applications in semiconductor-based technology and fundamental science [1]. At this level of miniaturization, the electric characterization of capacitive NEMS is a challenge. In fact, the highfrequency signals are normally cut off due to the limited bandwidth of the readout circuitry dominated by stray capacitance between the large amplifier and the high impedance tunneling junction [2]. In order to overcome such limitations, some methods have been suggested: using several amplifier stages [3], inserting an LC resonator [4], or more recently the scanning tunneling microscopy down-mixing readout technique [5]. Here, a new characterization technique for resonant NEMS based on capacitive downmixing has been developed and applied to M&NEMS accelerometers.

Manufacturing:

The M&NEMS process flow starts with the lithography and etching of the NEMS resonator and of the bulk contact. Then a 0.3 µm thick oxide deposition followed by a lithography and etching of the nanoresonator protection are proceeded. In the same step, an over etching of the oxide is led to open the silicon bulk contact. A few microns thick silicon epitaxial growth is done to realize the MEMS part. Contacts are defined by a 0.5 µm metal deposition followed by a lithography and etching of pads. A last lithography step and a DRIE of silicon thick layer is proceeded to make the MEMS structure, to isolate the bulk contact, and to open the SiO2 protective layer of the NEMS. The release of the sensor then is achieved by HF-vapor etching.

Experimental characterisation:

A new technique for the electrical characterization of capacitive NEMS resonators is performed in order to overcome the problems of parasitic impedances at high frequency and extend the usability of lowfrequency, high-input impedance measurement equipments. The measurements are performed on a M&NEMS resonant accelerometer shown in Fig. 1, which has been designed and fabricated in order to validate our characterization technique schemed in Fig. 2. It utilizes two AC signals and the intrinsic properties of the resonator to perform heterodyne down-mixing of the signal to a much lower frequency, which can then be detected by standard circuitry without significant signal loss. It also greatly reduces the unwanted background from the cross-talk between the detector and actuator circuits. As proof of principle, this down-mixing scheme has been applied to the detection of NEMS resonator with fundamental mode frequency of 13.5 MHz (see Fig. 3). As shown in Fig. 4, this technique permits tuning the signal to nonlinearity ratio and sensing the nonlinear nanomechanical motion in high-frequency resonators.

Conclusions and perspectives:

In a future work, the nonlinear dynamics of the considered accelerometer will be modeled while taking into account the down-mixing characterization scheme. Experimental and model results will be compared, which will complete the experimental interpretations and finally improve the presented set-up.

Word count: 487 words

Fig. 4

 4 Fig.4 Nonlinear resonance peaks measured using the developed capacitive down-mixing technique.

Fig. 2 :

 2 Fig.2: Schematic of the capacitive down-mixing setup using 2 AC voltages Ua and Ub.

Fig. 1 :

 1 Fig.1: (a): SEM view of an in-plane M\&NEMS accelerometer. (b): A focus on the gauge lets clearly appear the MEMS inertial mass of 2 µm thick, and the sub-µm resonator that has a section of 0.25 × 0.5 µm 2 . (c): A SEM image of a designed nanoresonator.

Fig. 3

 3 Fig.3 Linear resonance peaks measured using the developed capacitive down-mixing technique.