Thomas Peyret

Gwenolé Corre

Mathieu Thevenin

Kevin Martin

Philippe Coussy

An Automated Design Approach to Map Applications on CGRAs

Keywords: B.5.2 [Register-Transfer-Level Implementation]: Design Aids -Automatic synthesis CGRA, Mapping, Scheduling, Binding

Coarse-Grained Reconfigurable Architectures (CGRAs) are promising high-performance and power-efficient platforms. However, their uses are still limited by the capability of mapping tools. This abstract paper outlines a new automated design flow to map applications on CGRAs. The interest of our method is shown through comparison with state of the art approaches.

INTRODUCTION

For the last two decades, Coarse-Grained Reconfigurable Architectures (CGRAs) have been mainly proposed for accelerating multimedia applications. CGRA are indeed an interesting trade-off between FPGAs and many-core architectures thanks to their power efficiency and programmability [START_REF] Taylor | Is dark silicon useful?: harnessing the four horsemen of the coming dark silicon apocalypse[END_REF]. The literature is very rich in CGRAs architectures, which distinguish by different features such as the granularity of the Processing Elements (PE) named tile, homogeneity or heterogeneity of PE, type of operators, absence/presence of Register Files (RF) or interconnection network topologies. Figure 1 presents an example of CGRA.

The result of the "compilation" of an application on a CGRA (named mapping) is the scheduling and the binding of its operations on operators and registers. This NP-complete process [START_REF] Lee | Mapping multidomain applications onto coarse-grained reconfigurable architectures[END_REF] must be automated to allow efficient mapping of complex applications. Several methods have been proposed to tackle this problem. They are split in two categories i.e. [START_REF] Brenner | Optimal Simultaneous Scheduling, Binding and Routing for Processor-like Reconfigurable Architectures[END_REF] approaches that solve scheduling and binding separately with heuristics or meta-heuristics [START_REF] Friedman | SPR: an architecture-adaptive CGRA mapping tool[END_REF][START_REF] Lee | Mapping multidomain applications onto coarse-grained reconfigurable architectures[END_REF][START_REF] Park | Edge-centric modulo scheduling for coarsegrained reconfigurable architectures[END_REF] or by combining an heuristic and an exact method [START_REF] Hamzeh | REGIMap: register-aware application mapping on coarsegrained reconfigurable architectures (CGRAs)[END_REF] and (2) approaches that solve the whole problem entirely with exact method [START_REF] Brenner | Optimal Simultaneous Scheduling, Binding and Routing for Processor-like Reconfigurable Architectures[END_REF] or meta-heuristics [START_REF] Mei | DRESC: A retargetable compiler for coarse-grained reconfigurable architectures[END_REF][START_REF] De Sutter | Placement-and-routing-based register allocation for coarsegrained reconfigurable arrays[END_REF]. This paper presents a unified approach that maps application on CGRAs. The proposed mapping flow relies on simultaneous scheduling and binding steps respectively based on a heuristic and an exact method followed by a pruning step. The graph of the application is backward traversed and dynamically transformed allowing a better exploration of the design space. This extended abstract paper is organized as follows. Section 2 depicts proposed method. Section 3 presents the experiments and discusses obtained results. Conclusion is given in Section 4.

PROPOSED METHOD

Our design flow is presented in Figure 2. Inputs are a C/C++ application code compiled to obtain a formal Control Data Flow Graph (CDFG) and the targeted CGRA's model. Objective of the method is to minimize latency under resource constraint. The proposed mapping approach allows exploring the design space while keeping computation time low.

The key idea is to combine the advantages of exact and heuristic methods while minimizing as much as possible their respective drawbacks. CDFG is mapped by processing each Data Flow Graph (DFG) of basic bloc sequentially. A list-scheduling based algorithm schedules nodes of each DFG. As it is a local greedy method, the binding is made simultaneously to ensure that at least one solution exists, hence avoiding dead-ends, and is realized incrementally by using an exact method derived from Levi's algorithm [START_REF] Levi | A note on the derivation of maximal common subgraphs of two directed or undirected graphs[END_REF]. However, as exact methods do not scale up [START_REF] Lee | Mapping multidomain applications onto coarse-grained reconfigurable architectures[END_REF], a wise pruning step is executed at the end of each scheduling cycle to remove redundant partial mappings and thus keep a reasonable number of solutions during mapping process. Besides, DFGs are dynamically transformed as needed when no mapping (i.e. during scheduling or binding) solution is found. DFGs are also backward traversed to allow for using more different graph transformations.

EXPERIMENTS AND RESULTS

The proposed synthesis flow has been fully automated using Java. GCC has been used to generate CDFGs from applications. Five applications from signal processing domain have been used for our experiments: DC filter, Elliptic filter, Moving Exponential Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage, and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). Copyright is held by the author/owner(s). Average Filter (EMA), Moving Window Deconvolution (MWD) and unsharp mask. To obtain a large spectrum of results, several constraints have been considered: CGRA size, RF size and the number of tiles the final mapping is allowed to use leading to 16 different set of constraints per application and per method.

GLSVLSI'14,

The proposed approach is compared with two approaches from state of the art. The first, named "Method 1", solves the scheduling and the binding problem separately as the initial step of [START_REF] Lee | Mapping multidomain applications onto coarse-grained reconfigurable architectures[END_REF]. It uses a forward list scheduling algorithm and binding is made by using Levi's algorithm. "Method 2" forward traverses the graph, schedules nodes by applying statically graph transformations and tries to find a mapping by using Levi's algorithm as proposed in [START_REF] Hamzeh | REGIMap: register-aware application mapping on coarsegrained reconfigurable architectures (CGRAs)[END_REF] (that have been shown to provide better results than [START_REF] De Sutter | Placement-and-routing-based register allocation for coarsegrained reconfigurable arrays[END_REF]).

Two metrics were considered: (1) success rate (percentage of time the method finds a solution when at least one of the compared methods succeeds) and ( 2) percentage of time the method gives the best latency between the compared methods. Figure 3 and Figure 4 give the comparisons between the three methods for the previously defined metrics.

Figure 3 shows that Method 1, which solves scheduling and binding totally separately, leads to the lowest success rate (~56%). Method 2, which transforms the graph a priori, provides better results (~67%) but is not as good as the proposed approach (~98%). Figure 4 shows the percentage of time each method found the best latency and shows that the Proposed Method finds it most of the time (~82%) even if it relies on a heuristic-based scheduling algorithm, while the Methods 1 and 2 find it for respectively 57% and 63% of the benchmark.

CONCLUSION

In this paper, a generic method to map applications written in high level language on CGRA architectures has been presented. Experimental results show that this method finds 82% of time the best latency, has the highest success rate and achieves 2.2 times better mappings throughput compared to the other methods and thus achieves a very good exploration of the solution space.

Figure 4 GeneralFigure 3

 43 Figure 4 General Flow and Algorithm Core

4×4 CGRA with 2D mesh torus and RF in each tile

  May 21-23, 2014, Houston, Texas, USA.

					From Neighbours & Memory
	PE PE PE Figure 1 A PE	PE PE PE PE	PE PE PE PE	PE PE PE PE	RF To Neighbours & Memory FU
	ACM 978-1-4503-2816-6/14/05.
	http://dx.doi.org/10.1145/2591513.2591552


Sylvain Lavernhe
email: lavernhe@lurpa.ens-cachan.fr

Yann Quinsat

Claire Lartigue

Christopher Brown
email: brown@wpi.edu

-Y Lartigue Quinsat -C

Lurpa

Realistic simulation of surface defects in 5-axis milling using the measured geometry of the tool

Keywords: Surface topography, Defects, Realistic modeling and simulation, Surface analysis, Surface machining

Realistic simulation of surface defects in 5-axis milling using the measured geometry of the tool

Introduction

Surface topography characterization resulting from 5-axis milling is a current and essential issue. Some studies have enhanced the importance of 3D surface topography in engineering applications. 3D surface topography influences mechanical and physical properties of contacting parts and plays a major role in surface integrity, which affects fatigue life [START_REF] Souto-Lebel | Characterization and influence of defect size distribution induced by ball-end finishing milling on fatigue life[END_REF][START_REF] Novovic | The effect of machined topography and integrity on fatigue life[END_REF]. The defect size distribution, as well as the defect orientation, are of great interest for a better understanding of the mechanical behavior under loads.

The characterization of surface topography as a function of cutting conditions is thus a challenge, in particular for the prediction of the shape, size and orientation of the defects induced by machining. Several authors handle this issue through simulations. Most simulations concern 3-axis ball-end milling and, are based on the 3D modeling of both the tool geometry and the workpiece. Kim and Chu [START_REF] Kim | Texture prediction of milled surfaces using texture superposition method[END_REF] described the texture obtained in milling by using a texture superposition model considering the effect of the fillet radius. In this study, the authors investigate the influence of the feedrate, the cutter type, and the run-out. Bouzakis et al. focused on the motion of the cutting edge, which is segmented into elementary linear sections of constant width. These authors highlight the influence of the tool orientation, the transverse step and the feedrate on the machined surface quality [START_REF] Bouzakis | Determination of the chip geometry, cutting force and roughness in free form surfaces finishing milling, with ball end tools[END_REF]. Liu et al. detailed a method based on the solid modeling of each cutter flute [START_REF] Liu | Surface finish visualisation in high speed, ball nose milling applications[END_REF]. The workpiece surface is sampled in the (XY) plane by a set of parallel lines (called spikes) to the Z axis (tool axis). The cutter tool generates a volume which truncates the spikes. The remaining parts of the reference spikes provide a good approximation of the surface finish which can be linked with the cutting parameters. In a previous work, we proposed to link the machining strategy in 3-axis ball-end milling with a 3D surface roughness parameter and to optimize the machining direction ac-cording to this parameter [START_REF] Quinsat | Surface topography in ball end milling process: Description of a 3D surface roughness parameter[END_REF]. In their work, Arizmendi et al. [START_REF] Arizmendi | Model development for the prediction of surface topography generated by ball-end mills taking into account the tool parallel axis offset. experimental validation[END_REF] proposed a model accounting for the tool parallel axis offset. The cutting edge geometry is represented as a cylindrical helix projected onto a spherical surface. The topography prediction is thus obtained considering the equation of the cutting edge trajectory and the equation of the material swept by the tool. More recently, Buj-Corral et al. [START_REF] Buj-Corral | Surface topography in ball-end milling processes as a function of feed per tooth and radial depth of cut[END_REF] focused on the prediction of roughness generated in ball-end milling of flat surfaces. The surface is modeled as a set of horizontal planes (the vertical Z defining the tool axis), each of them sampled according to a square grid (X, Y). Each cutting tool tooth is modeled as an elliptical curve of increasing radius. The topography is built based on a set of intersections between the tool motion and the part. The effect of cutting conditions is taken into account as well as the effect of some characteristic of the tool geometry. In their approach, Denkena et al. [START_REF] Denkena | Kinematic and stochastic surface topography of machined TiAl6V4-parts by means of ball nose end milling[END_REF] considered the surface topography as the sum of the "kinematic topography", resulting from the cutting edge movement, and the "stochastic topography", including surface irregularities and cracks. The first one is obtained using Boolean operations, simulating material removal, between the CAD-model of the tool and a discrete model of the workpiece. The "stochastic topography" is generally obtained by subtracting the kinematic topography to surface measurements. Nevertheless, the authors proposed a method to evaluate this stochastic topography with an empirical model based on the process force simulations.

Few works handle the case of 5-axis milling. Kim et al. [START_REF] Kim | Texture prediction of milled surfaces using texture superposition method[END_REF] proposed to simulate the 3D topography obtained in 5axis milling by using a filleted-ball end tool. The envelope of the tool movement is modeled by successive tool locations according to the feed per tooth. Antoniadis et al. [START_REF] Antoniadis | Prediction of surface topomorphy and roughness in ball-end milling[END_REF] consider the general case of multi-axis with ball-end tools. The cutting edge of the tool is decomposed into an elementary cutting edge of trapezoidal shape. The workpiece is modeled as a set of linear segments (or needles) defining a brush. The geometrical model of the cutting tool follows the selected milling operation and intersects the set of needles. The part material remaining defines the milled surface and gives its 3D topography. Starting from the expression of the trajectory equation of the cutting edge relative to the workpiece and the tool paths, Zhang et al. introduced a new and general iterative method to simulate roughness in multi-axis ballend milling [START_REF] Gao | A new algorithm for the numerical simulation of machined surface topography in multiaxis ball-end milling[END_REF]23]. The authors investigated in particular the effects of the tool inclination and of the cutting mode on surface roughness. Indeed, their study showed that surface roughness remains mostly unchanged beyond a certain limiting inclination angle, and that the cutting mode influences mainly sculptured surface roughness.

In the case of 5-axis milling, we proposed in a first study a theoretical approach to simulate the 3D surface topography obtained in 5-axis milling with a filleted-ball end cutter tool integrating actual feedrate evolution. The material removal simulation relies on the well-known N-buffer method and requires the modeling of the tool geometry, the surface and the actual tool trajectory [START_REF] Lavernhe | Model for the prediction of 3D surface topography in 5-axis milling[END_REF][START_REF] Quinsat | Characterization of 3D surface topography in 5-axis milling[END_REF]. The surface is discretized according to a rectangular grid with lines directed along the local normal to the surface. The set of lines is truncated by the envelope of the tool trajectory leading to the 3D surface topography (Fig. 1).

The local actual feedrate evolution is obtained by a kinematical predictive model. The tool is assumed to be rigid. It is approximated by a local meshing built from the CAD theoretical model of the tool geometry. The comparison between simulations and experimental trials shows good agreement between the predicted 3D topography and the experimental one (Fig. 2).

Nevertheless, and as for most methods available in literature, the modeling of the cutting process only attempts geometrical fidelity on scales of the curvature of a tool with a smooth cutting edge. This leads to a good representation of the general pattern, but the model is not able to predict small defects appearing at the bottom of the valleys (Fig. 3).

However, the knowledge of the size and distribution of the defects is essential in fracture mechanics to establish the relationship between defect distribution and the fracture probability. Studies have enhanced the effect of microgeometry on fatigue life [START_REF] Novovic | The effect of machined topography and integrity on fatigue life[END_REF] [START_REF] Guillemot | Prediction of the endurance limit taking account of the microgeometry after finishing milling[END_REF]. Fatigue failure of a mechanical part can originate from micro-scratches on the surface. By considering the statistical defect size distribution that exists in an effective volume, it is possible to calculate the failure probability of the structure. For small defects (less than 5 µm) Guillemot and Souto-Lebel et al. [START_REF] Souto-Lebel | Characterization and influence of defect size distribution induced by ball-end finishing milling on fatigue life[END_REF] proposed to model the defect size population as a beta function and used this model to identify the defect distribution leading to experimental fatigue limits and compared it favorably to the distribution obtained by direct measurement.

In addition, the general direction of the texture lay or anisotropy is of great importance with regard to the load direction. To predict actual surface topography, in contrast with Dekena et al., one can consider that the difference between the measured topography and the topography that was simulated using classical tools (also called "kinematic topography") is not stochastic. The difference might rather be due to finer scale phenomena not taken into account in simulation, e.g., actual tool geometry, and tool deflection. In fact, as displayed in Fig. 3, local defects in the bottom of the valleys seem to repeat in the machining (or feed) direction. It can be observed a large population of defects smaller than 10 µm. This is probably related to defects of the edge of the cutting tool. Such defects are not predictable if the geometry of the cutting edge is supposed to be smooth. A change of scale is required.

Therefore, the prediction of local defects requires the modeling of the tool surface characteristics at a finer scale, more similar to that used for grinding [START_REF] Inasaki | Grinding process simulation based on the wheel topography measurement[END_REF][START_REF] Aurich | Development of a superabrasive grinding wheel with defined grain structure using kinematic simulation[END_REF]. Generally speaking, most simulations of the grinding process are based on models describing the macro-and micro-geometry of the grinding wheel, the grit size, the grit volume, and the grit distribution. These parameters are obtained from an experimental characterization of the grinding tool. Actually, the cutting edge distribution used to simulate the relative motion between the grinding wheel and the workpiece is obtained by measuring the wheel topography. As a result, surface roughness predictions are in good agreement with experimental ones.

In a similar way, we propose to improve 3D surface topography prediction in milling by characterizing the actual cutting edge of the tool. For this purpose, the tool cutting edge is measured using an Alicona infinite focus system which measures the actual topography of the cutting edge. The latter is integrated in the simulation method previously implemented for 5-axis milling surface topography prediction to account for the measured topography of the tool.

Model for surface topography simulation with measured tool topography

The material removal simulation relies on the well-known N-buffer method and requires the modeling of the surface, the tool geometry and the actual tool trajectory [START_REF] Lavernhe | Model for the prediction of 3D surface topography in 5-axis milling[END_REF][START_REF] Quinsat | Characterization of 3D surface topography in 5-axis milling[END_REF].

The surface is discretized according to a rectangular grid. At each point of the grid, a line directed along the local normal to the surface is defined. The set of lines is truncated by the envelope of the tool trajectory, and the remaining part defines the 3D topography of the machined surface (Fig. 1).

The tool, considered as rigid, is approximated as a local mesh, expressed in the STL format. Only active cutting edges are considered to reduce the number of nodes and so the calculation time. In previous works, these cutting edges were assumed to be smooth. More details on the mesh characteristics are given in section 3. For each simulation step of the model, the mesh is transformed versus the workpiece according the kinematics of the 5-axis trajectory.

Concerning the tool trajectory, the input of the material simulation model is the actual tool trajectory expressed as a set of tool postures and corresponding feedrates. The latter are calculated based on a model of velocity prediction which gives a good estimation of the local feedrate of the tool teeth [START_REF] Lavernhe | Kinematical performance prediction in multi-axis machining for process planning optimization[END_REF].

The simulation of material removal thus consists of two main steps:

1. calculation of the spindle angular positions at each tool posture 2. trajectory sampling preceding the N-buffer computation Based on a linear interpolation, of the tool trajectory between the local postures {X p ,Y p , Z p }, the time interval ∆t separating two postures is calculated as follows:

∆ L i,i+1 = (X i+1 p -X i p ) 2 + (Y i+1 p -Y i p ) 2 + (Z i+1 p -Z i p ) 2 (1)
∆t i,i+1 = ∆ L i,i+1 V i+1 f -V i f 2 (2)
in which the local feedrate V f is linearly interpolated.

Let Ω spindle be the rotational velocity of the spindle, supposed to be constant regardless of fluctuations in the actual feedrate, angular positions of the tool axis α are given by:

α i+1 = α i + Ω spindle • ∆t i,i+1 (3)
The elementary trajectory defined between two tool postures is checked considering a fixed geometrical step dα to ensure a good discretized description of the envelop of the tool movement:

dt = dα Ω spindle (4)
Given this elementary time step dt, it is thus possible to discretize the time interval ∆t and determine the number of samples N samples for each linear interpolation by the integer part and then the current sampling number N * :

N * ∈ [[0 .. N i,i+1 samples ]] with N i,i+1 samples = ⌊ ∆t i,i+1 dt ⌋ (5)
This is used in the calculation of the sampled tool velocity V * f and locations of the components X * p , Y * p or Z * p along the elementary trajectory:

V * f = V i f + V i+1 f -V i f ∆t i,i+1 • N * • dt (6)
X * p = X i p + dX * = X i p + X i+1 p -X i p ∆ L i,i+1 • V * f +V i f 2 • N * • dt (7)
The components I * , J * and K * of the tool axis orientation u are computed so that the orientation vector of the tool runs in a plane created by the start and end vectors u i and u i+1 for each linear interpolation. The rotation angle θ * between the start vector u i and an intermediate orientation u * respects the ratio between the angle and the linear displacement covered:

θ * = θ i,i+1 • N * • dt ∆t i,i+1 (8)
The intermediate orientation vector u * is given by the Rodrigues' rotation formula, Eq. 9, with w the unit vector defining the axis of rotation (cross product of the start and end vectors).

u * = u i cos θ * + (w × u i ) sin θ * + w(w • u i )(1 -cos θ * ) (9)
The interpolated angular position of the tool α * is directly given by:

α * = α i + Ω spindle • N * • dt (10)
Finally, the simulated machined surface is obtained by computing the intersections between the normal lines representing the surface and the tool for each configuration

{X * p , Y * p , Z * p , I * , J * , K * , α * }.
The whole simulation method for material removal is summarized in Fig. 4. This model leads to simulations in good agreement with experimental results, in 3-axis as well as 5-axis milling and regardless of the tool geometry [START_REF] Lavernhe | Model for the prediction of 3D surface topography in 5-axis milling[END_REF]. However, as aforementioned, the model failed in the prediction of chaotic small defects. The main reason is likely due to the ideal modeling of the tool cutting edges. Therefore, the simulation is improved by integrating a realistic description of the cutting edge instead of the perfect geometrical description. The actual geometry of the cutting edge is only accessible through measurement. For this purpose, once the tool is selected its cutting edges have to be measured before machining. To control the process better, it is easier to consider just one cutting edge. Therefore, only the track associated to this tooth will be left on the part surface after machining. To illustrate our purpose, the machining of a plane surface with a ball ended tool inclined by an tilt angle β in the feed direction V f is proposed (Fig. 5).

The tool cutting edge is measured using an Alicona infinite focus. Data provided by the manufacturer state the best vertical resolution of the measurement system is 20 nm and the vertical repeatability is 8 nm. For the simulation, and to limit the number of points (18790 points), only the functional part around the contact point C c of the cutting edge is considered. This functional portion of the cutting edge is The analysis of noise on the mesh is not easy. Indeed, the study of the distribution of points around a theoretical sphere (Fig. 3) shows that it is difficult to differentiate the points associated to the cutting edge from those corresponding to the rake face and flank face. However, given the characteristics of the measuring system used, the resulting mesh is compatible with our need.

A specimen of high-strength steel (50CrMo4 bainitic steel see Tab. 2 for specimen mechanical properties) is machined in 5 axes with a constant tilt angle is taken egal to 20 • on a HSM milling center using a ball-ended tool (R = 5 mm). The cutting conditions (Ω spindle = 18000 rpm, feed per tooth 0.2 mm) are chosen according to the tool manufacturer recommendations. The machined specimen is then measured with a chromatic confocal Micromeasure station.

Simulation parameters are chosen in order to generate deviations with a magnitude lower than the measurement noise, and that is therefore studied. Hence, considering the tool radius of 5 mm, the angular step dα = 0.2 • generates a chordal deviation δ under 10 nm (Eq.11). As the linear interpolation of the trajectory is based on this parameter (dα), its does not add additional geometrical error on the surface topography.

δ = R • (1 -cos( α 2 )) (11)

Results and discussion

Results of the measurement and of the simulation are displayed in Fig. 7. To highlight the contributions of the proposed method, a surface compound of spherical caps is added. This surface is computed according to the models proposed by Kim et al. [START_REF] Kim | Texture prediction of milled surfaces using texture superposition method[END_REF] and Quinsat et al. [START_REF] Quinsat | Surface topography in ball end milling process: Description of a 3D surface roughness parameter[END_REF]. As it can be seen, the simulation using the measure tool profile appears to be a good estimation of the actual surface pattern (measured), although the amplitude of the peaks is slightly over-estimated. Actually, if only the track of the sphere is considered (Fig. 7right), the simulation lacks any of fine scale details evident on the measured surface. The apparent improvement in the simulation achieved by using the measured tool profile is confirmed by considering the roughness parameters: relative area area-scale complexity. These standardized parameters [START_REF]Geometrical product specifications (GPS) -Surface texture: Areal -Part 2: Terms, definitions and surface texture parameters[END_REF] relate to the inclinations on the surface and account for the complexity of the surface shapes at different scales. In fact, the idea is to exploit the chaotic, or fractal, nature of the surface, in particular that the apparent area of a surface depends on the scale of observation. In the area-scale analysis the topography of the surface is tiled using triangular patches. At any particular considered scale, all the patches have the same area (defining the scale), and the relative area at the scale is the calculated area divided by the nominal area [START_REF] Berglund | Milled die steel surface roughness correlation with steel sheet friction[END_REF][START_REF] Ungar | Quantification of dental microwear by tan[END_REF]1]. The calculated area is given by the number of patches at the considered scale multiplied by the area of the patch at the scale. The area-scale complexity is the slope of the relative area plot. It can be calculated between any two scales where the tiling is performed. If the surfaces can be considered multi-fractal with respect to scale, then this complexity is related to the fractal dimension [1]. These parameters, which were calculated using scale-sensitive fractal analysis [START_REF]SFrax: Surface Metrology and Fractal Analysis[END_REF], are displayed in figures 8 and 9. The relative area at different scales clearly highlights that the model based on the smooth sphere is not a good representation of the actual machined surface. At large scales, the simulation with the measured tool profile provides a better representation of the actual surface. On the other hand, at fine scales, the simulation leads to higher values of the relative area, which could be associated with the difference in the peak-to-valley magnitude previously observed. Actually, some effects, such as, abrasion and elastic strains of the material below the cutting edge present during the cutting process, are not taken into account in the simulation.

The complexity of the simulation based on the tool profile (realistic simulation) reaches a maximum at 24.68 µm 2 and then declines sharply at finer scales This could be because of the refinement of the mesh to a liner scale of about 7µm, which corresponds to the leg of a triangle with an area of about 24.5µm 2 . The complexity of the measured surface drops of sharply at scales below 1µm 2 , corresponding to linear scales of about the square root of two. As the value of the diameter spot is equal to 2µm, the complexity falls at small scales. Once more, the topography model based on the smooth spherical tool (spherical caps) has a complexity that is as much as two orders of magnitude below the measured topography. This remark is confirmed by conventional roughness parameter evaluation (Tab. 3). To only consider conventional roughness characterization parameters, the surfaces are filtered using a Gaussian filter (λ cuto f f = 0.08mm). Results clearly highlight a better agreement of the value of the conventional characterization parameters obtained using the realistic simulation with those obtained by measurement. The simulation with only spherical caps is not sufficient to account for fine scale surface topography. Therefore, the simulation gives a good estimation of the surface roughness that could be used to link surface integrity with the fatigue-life of the part.

Conclusion

In the paper, an algorithm for a more realistic simulation of surface formation is proposed. This simulation relies on the well-known Z-buffer method for which a tool cuts a set of vertical spikes representing the part surface. To account for the real tool geometry, the tool cutting edge is measured using an optical means that gives a good representation of the actual edge. To assess the method, a steel workpiece is machined by 3-axis ball-end milling. Comparisons between the 3D topography obtained by realistic simulation to the measured one highlights that the simulation leads to a better estimation of the topography than does a simulation with a smooth tool. In particular, the multi-scale study of the relative area clearly shows that the simulation provides a better representation of the topography than a smooth, sphere-cap tool model. This gives a good estimation of the motif left on the part surface for a scale that can be useful to study the link between surface integrity and fatigue-life. However, the surface complexity obtained by simulation is lower than that obtained for the measured one likely due to calculation limits. A sensitive study of the effect of the tool meshing combined with the effect of the part meshing would give the adequate compromise calculation time/ quality of the representation. It would be interesting to complete the modeling by considering the mechanical behavior of the part under cutting forces. Accounting for the local part elasticity during tool cutting would likely lead to a better estimation of the peak magnitude.

Fig. 1 Fig. 2

 12 Fig. 1 3D surface topography simulation in 5 axis-milling

Fig. 3

 3 Fig. 3 Actual defects in the bottom of the valleys

Fig. 4 Fig. 5

 45 Fig. 4 Material removal simulation process

Fig. 6

 6 Fig. 6 Mesh point distribution

Table 2

 2 Mechanical properties of the bainitic steel used for the trial E(MPa) R m (MPa) Σ v (MPa) C(MPa) g A(%) Hardness(HV )

Fig. 7

 7 Fig. 7 3D obtained pattern (left-measured; middle-simulation; right-sphere track)

  Fig. 8 Relative area vs scale

Fig. 9

 9 Fig. 9 Complexity vs Scale

  dem scanning confocal microscopy and scale-sensitive fractal analyses. Scanning 25(4):185-193, DOI 10.1002/sca.4950250405 23. Zhang WH, Tan G, Wan M, Gao T, Bassir D (2008) A new algorithm for the numerical simulation of machined surface topography in multiaxis ball-end milling. Journal of Manufacturing Science and Engineering 130(1):011,003 1-11, DOI 10.1115/1.2815337

Table 1

 1 Mesh description

		Edge length (µm) Face area (µm 2 )
	Mean	0.79	0.25
	Standart Deviation 0.23	0.10

meshed in a STL format giving 36389 facets. The dimensional characteristics of the mesh used are given in table 1.

Table 3

 3 Amplitude roughness parameters[START_REF]Geometrical product specifications (GPS) -Surface texture: Areal -Part 2: Terms, definitions and surface texture parameters[END_REF]

		Sq(µm) Ssk	Sku	Sp(µm) Sv(µm) Sz(µm) Sa(µm)
	Measured	0.158	-0.0281 2.85 0.515	0.825	1.34	0.127
	Realistic simulation 0.326	0.488	2.90 0.855	0.611	1.47	0.258
	Spherical caps	0.0257	2.34	7.64 0.105	0.0153	0.120	0.0177


Luc Pronzato
email: pronzato@i3s.unice.fr

Éric Thierry

Sequential experimental design and response optimisation

Keywords:

We consider the situation where one wants to maximise a function f (θ, x) with respect to x, with θ unknown and estimated from observations y k . This may correspond to the case of a regression model, where one observes y k = f (θ, x k ) + ǫ k , with ǫ k some random error, or to the Bernoulli case where y k ∈ {0, 1}, with Pr[y k = 1|θ, x k ] = f (θ, x k ). Special attention is given to sequences given by x k+1 = arg max x f ( θk , x) + α k d k (x), with θk an estimated value of θ obtained from (x 1 , y 1 ), . . . , (x k , y k ) and d k (x) a penalty for poor estimation. Approximately optimal rules are suggested in the linear regression case with a finite horizon, where one wants to maximize N i=1 w i f (θ, x i ) with {w i } a weighting sequence. Various examples are presented, with a comparison with a Polya urn design and an up-and-down method for a binary response problem.

Introduction

We consider an optimisation problem, where one wants to maximise a function f (θ, x) with respect to x ∈ X ⊂ IR q , where θ is an unknown vector of parameters, θ ∈ IR p . Two cases will be considered. The first one corresponds to a regression problem, where one observes

y k = f (θ, x k ) + ǫ k , k = 1, 2, . . . (1)
with ǫ i an unobservable error such that E{ǫ k |ǫ 1 , . . . , ǫ k-1 } = 0 and E{ǫ 2 k |ǫ 1 , . . . , ǫ k-1 } < ∞, k = 1, . . . almost surely (a.s.).

In the second case, when the design variable is x k , one observes a Bernoulli response:

y k ∈ {0, 1} with Pr[y k = 1|θ, x k ] = f (θ, x k ) , k = 1, 2, . . . (2)
We shall denote F k the σ-field generated by y 1 , . . . , y k , E{•|F k } the corresponding posterior expectation and E{•|F 0 } the prior expectation, with a prior probability measure µ for θ. In the two cases above, the sequence of design points {x k } and observations {y k } is used to estimate θ. If the issue were only to determine the value x * that maximises f (θ, x), one could resort to optimum-design theory for choosing an appropriate sequence of inputs {x k }, see, e.g., [START_REF] Ford | A sequentially constructed design for estimating a nonlinear parametric function[END_REF][START_REF] Wu | Optimal design for percentile estimation of a quantal response curve[END_REF][START_REF] Chaloner | Bayesian design for estimating the turning point of a quadratic regression[END_REF][START_REF] Chaloner | Optimal Bayesian design applied to logistic regression experiments[END_REF][START_REF] Müller | Batch sequential design for a nonlinear estimation problem[END_REF][START_REF] Pronzato | Experimental design for estimating the optimum point in a response surface[END_REF][START_REF] Fedorov | Another view on optimal design for estimating the point of extremum in quadratic regression[END_REF]. However, here the objective is also to have each f (θ, x i ) as large as possible. The sequence {x k } must therefore fulfill two simultaneous (and generally opposite) objectives: (i) help to locate x * , (ii) be close to x * in order to maximise f (θ, x k ), k = 1, 2, . . . The problem is thus one of dual control, see the pioneering papers by Fel'dbaum [START_REF] Fel | Dual control theory[END_REF][START_REF] Fel | Dual control theory[END_REF][START_REF] Fel | The theory of dual control[END_REF][START_REF] Fel | The theory of dual control[END_REF]. Note that the x i 's will be chosen sequentially, contrary to the case considered in [START_REF] Verdinelli | Bayesian designs for maximizing information and outcome[END_REF]; that is, x i is F i-1 -measurable, i = 1 . . .

A rather general formulation of the problem (see [START_REF] Ginebra | Response surface bandits[END_REF]) is:

maximize E{ ∞ i=1 w i y i |F 0 }/ ∞ i=1 w i (3)
with respect to x 1 , x 2 , . . ., with {w i } a weighting (discount) sequence. For instance, the choice w i = 1 for i = N + 1 and w i = 0 otherwise corresponds to a pure design problem, where emphasis is put on the estimation of x * after N observations; w i = 1, i = 1, . . . , N , w N +1 = K and w i = 0 for i > N + 1 corresponds to the case where the best guess for x * at step N is used for the next K steps; the finite horizon case with no discount corresponds to w i = 1, i = 1, . . . , N , w i = 0 for i ≥ N + 1, etc.

In the two cases [START_REF] Aström | Adaptive Control[END_REF][START_REF] Chaloner | Optimal bayesian experimental design for linear models[END_REF], E{y i |θ, x i } = f (θ, x i ), so that (3) becomes maximize E{

∞ i=1 w i f (θ, x i )|F 0 }/ ∞ i=1 w i . (4)
We shall mainly consider design sequences that are constructed as follows: at step k, x k+1 maximizes the sum of the predicted value of f , that is, f ( θk , x) with θk the current estimated value of θ, and a weighted penalty term α k d k (x), with d k (x), the penalty for poor estimation of θ, also a function of (x i , y i ), i = 1, . . . , k:

x k+1 = arg max x∈X f ( θk , x) + α k d k (x) . (5)
In Section 2 we give an asymptotic result (infinite horizon, no discount) obtained for the regression problem (1) when f (θ, x) is linear in θ and d k (x) is the variance function used in the construction of D-optimum designs. Other penalty functions related to L-optimum design are then suggested. Section 3 is devoted to the finite horizon case: approximately optimal strategies are suggested, and a particular sequence of weights {α k } and penalty function d k (x) for ( 5) are obtained through a series of approximations of the original problem (4). Extension to nonlinear problems is presented in Section 4, where the case of Bernoulli trials (2) is considered through an example. Finally, Section 5 concludes and points out some open issues and possible developments.

2 Asymptotic results for linear response optimisation

D-optimum penalty

Consider the case where the horizon is infinite (N = ∞) and there is no discount (w i = 1 for any i), and assume that the regression function is linear in θ, that is, f (θ, x) = r ⊤ (x)θ. We use the penalty given by the variance function used in the sequential construction of D-optimum designs, see [START_REF] Wynn | The sequential generation of D-optimum experimental designs[END_REF][START_REF] Fedorov | Theory of Optimal Experiments[END_REF],

d k (x) = d D k (x) = r ⊤ (x)M -1 k r(x) , (6)
with M k the design matrix

M k = k i=1 r(x i )r ⊤ (x i ) . (7)
Assume that r(x) is continuous on X compact, that the first K 0 regressors r(x 1 ), . . . , r(x K0 ) are such that M K0 is positive definite. Also assume for the moment that θ in (1) takes a deterministic (but unknown) value θ (we shall come back to Bayesian estimation later), and that r ⊤ (x) θ has a unique global maximiser x * in X ; that is: for all β > 0 , there exists ǫ > 0 such that

r ⊤ (x) θ + ǫ > r ⊤ (x * ) θ ⇒ x -x * < β . (8)
We estimate θ by least squares (LS):

θk = arg min θ∈Θ k i=1 [y i -r ⊤ (x i )θ] 2 , (9)
with Θ a compact subset of IR p such that θ ∈ Θ. The difficulty is that for a suitable weighting sequence {α k }, the sequence {x k } accumulates at the value x * (θ) that maximises f (θ, x) for some θ, but, at the same time, when p = dim θ > 1, a sequence too concentrated gives a singular design matrix, and thus does not yield consistent estimates. The following theorem is proved in [START_REF] Pronzato | Adaptive optimisation and D-optimum experimental design[END_REF], using the results in [START_REF] Lai | Least squares estimates in stochastic regression models with applications to identification and control of dynamic systems[END_REF] on almost sure convergence of LS estimates. We denote by ξ k the empirical measure of design points generated by [START_REF] Chaloner | Bayesian experimental design: a review[END_REF] and by ξ x the discrete measure that puts weight 1 at the point x.

Theorem 1 Assume that sequence of weights {α k } in ( 5) is such that (α k /k) log α k decreases monotonically and α k /(log k) 1+δ increases monotonically to ∞ for some δ > 0. Then, the sequence {x k } generated by ( 5) and ( 6) is such that θk → θ, (1/k) k i=1 r ⊤ (x i ) θ → r ⊤ (x * ) θ and ξ k w -→ ξ x * (in the sense of weak convergence of measures) a.s. as k → ∞, with x * = x * ( θ) = arg max x∈X r ⊤ (x) θ.

Assume now that θ has a prior µ. One can first note that from the Lebesgue dominated convergence theorem, Theorem 1 implies

E{(1/k) k i=1 r ⊤ (x i ) θ} → E{r ⊤ [x * ( θ)] θ}
when µ is supported on Θ compact. Also, using the Bayesian imbedding approach of [START_REF] Sternby | On consistency for the method of least squares using martingale theory[END_REF][START_REF] Kumar | Convergence of adaptive control schemes using least-squares parameter estimates[END_REF], one can obtain a.s. convergence results under weaker conditions on {α k } under the assumption that the errors ǫ i are i.i.d. and Gaussian N (0, σ 2 ). Indeed, for Θ = IR p and k > K 0 the LS estimator θk then coincides with the Bayesian estimator E{θ|F k } for the prior µ given by N ( θK0 , σ 2 M -1 K0 ) (M K0 is positive definite). Let Q denote the probability measure induced by {ǫ k }, and write (µ × Q)-a.s. for a property almost sure in the sense of the product measure µ×Q. From the martingale convergence theorem, θ∞ = lim k→∞ θk exists and is finite and r ⊤ (x) θ∞ is bounded on X , (µ×Q)-a.s., and the posterior covariance matrix C k tends to some limit C ∞ , (µ × Q)-a.s. When the smallest eigenvalue of M k satisfies λ min (M k ) → ∞, (µ × Q)-a.s., C ∞ is the null matrix and θk converges to θ, the value of θ that generates the observations, (µ × Q)-a.s. A straightforward extension of Theorem 1 is then as follows.

Corollary 1 Assume that the errors ǫ i are i.i.d. N (0, σ 2 ) and that sequence of weights {α k } in [START_REF] Chaloner | Bayesian experimental design: a review[END_REF] is such that α k → ∞ and α k /k → 0. Then, the sequence {x k } generated by [START_REF] Chaloner | Bayesian experimental design: a review[END_REF] and [START_REF] Dette | A note on Bayesian cand D-optimal designs in nonlinear regression models[END_REF]

is such that θk → θ, (1/k) k i=1 r ⊤ (x i ) θ → r ⊤ (x * ) θ and ξ k w -→ ξ x * (in the sense of weak con- vergence of measures), (µ × Q)-a.s., as k → ∞.
The condition on α k is weaker in Corollary 1 than in Theorem 1, but note that there may be a singular set (with respect to the Lebesgue measure) for θ for which θk is not consistent.

Using the results of [START_REF] Hu | Strong consistency of Bayes estimates in stochastic regression models[END_REF], the assumption of normality in Corollary 1 can be relaxed, provided (i) the errors ǫ i are i.i.d. with an almost everywhere strictly positive density h with respect to the Lebesgue measure, such that h ′′ is continuous and (log h) ′′ < 0, (ii) the prior measure µ is absolutely continuous with respect to the Lebesgue measure and (iii) the LS estimates θk are replaced by E{θ|F k }.

Both Theorem 1 and Corollary 1 indicate that a penalty of the form β det M k+1 / det M k , see, e.g., [START_REF] Aström | Adaptive Control[END_REF], which corresponds to taking α k = α constant in [START_REF] Chaloner | Bayesian experimental design: a review[END_REF], is not strong enough: the sequential design that is obtained will not guarantee convergence of the parameters to θ, and {x k } may accumulate at a value x different from x * ( θ).

L-optimum penalty

Using an idea similar to previous section, one could use a penalty function related to L-optimum design; that is,

d k (x) = d L k (x) = r ⊤ (x)M -1 k H M -1 k r(x) , (10)
with H a non negative definite matrix. It is shown in [START_REF] Wu | The convergence of general step-length algorithms for regular optimum design criteria[END_REF] that when H is positive definite and the sequence {x k } is given by x k+1 = arg max x∈X d L k (x), the design measure ξ k converges to a L-optimum design measure ξ * L that minimises trace H I -1 (ξ), where I(ξ) = X r(x)r ⊤ (x)ξ(dx). Further work is required to check if a property similar to Theorem 1 can be obtained in this case.

A case that has retained much attention is when the design objective corresponds to the estimation of the point x * where f (θ, x) achieves its maximum, see [START_REF] Pronzato | Experimental design for estimating the optimum point in a response surface[END_REF], especially when f is a quadratic function of a scalar variable x, see, e.g., [START_REF] Chaloner | Bayesian design for estimating the turning point of a quadratic regression[END_REF][START_REF] Fedorov | Another view on optimal design for estimating the point of extremum in quadratic regression[END_REF]. Assume that

f (θ, x) = θ 0 + θ 1 x + θ 2 x 2 /2 , (11)
with x ∈ X , a compact subset of the real line. One has x * (θ) = -θ 1 /θ 2 , and, in the case where the errors ǫ i are i.i.d., the asymptotic variance of x * ( θk ), with θk given by ( 9), is proportional to

c ⊤ M -1 k c, with c = c(θ) = ∂x * (θ) ∂θ = (-1/θ 2 ) 0 1 x * ⊤ .
Choosing the x k 's to maximise the accuracy of the estimation of x * (θ) corresponds to c-optimum design, that is, to L-optimum design with H given by the rank-one matrix cc ⊤ . Note that the dependence of c in θ makes the problem nonlinear. A Bayesian approach is used in [START_REF] Chaloner | Bayesian design for estimating the turning point of a quadratic regression[END_REF][START_REF] Dette | A note on Bayesian cand D-optimal designs in nonlinear regression models[END_REF], based on a design criterion of the type E{c ⊤ M -1 k c}, where E{•} denotes the expectation with respect to θ for a given prior. Sequential approaches are considered in [START_REF] Ford | A sequentially constructed design for estimating a nonlinear parametric function[END_REF][START_REF] Müller | Batch sequential design for a nonlinear estimation problem[END_REF]. One can also refer to [START_REF] Chaloner | Optimal bayesian experimental design for linear models[END_REF] for the use of c-optimal design in the context of Bayesian estimation and to [START_REF] Chaloner | Bayesian experimental design: a review[END_REF] for a survey on Bayesian experimental design. Following [START_REF] Fel | Dual control theory[END_REF], a penalty function related to c-optimal design is thus

d k (x) = d c k (x) = [r ⊤ (x)M -1 k c( θk )] 2 . ( 12
)
Another approach, used in [START_REF] Pronzato | Experimental design for estimating the optimum point in a response surface[END_REF], is to derive the design criterion from the construction of a Bayesian risk related to the maximisation of f (θ, x). Assume that θ has a normal prior N ( θ0 , σ 2 Ω), that the errors ǫ i are i.i.d. N (0, σ 2 ) and that the discount factors w i satisfy w i = 1 for i = N + 1 and w i = 0 otherwise. When x 1 , . . . , x N are all chosen at the same time, this leads to the following (non sequential) L B -optimal design problem:

minimise trace H( θ0 ) (M k + Ω -1 ) -1 with H( θ) = ∂ 2 f [θ, x * (θ)] ∂θ∂θ ⊤ | θ .
The matrix H can easily be proved to be non negative definite when f is linear in θ, see [START_REF] Pronzato | Experimental design for estimating the optimum point in a response surface[END_REF], and

H(θ) = ∂r(x) ∂x ⊤ |x * (θ) ∂x * (θ) ∂θ ⊤ , (13)
which can be expressed analytically when f is quadratic in x ∈ IR q , q ≥ 1.

For instance, in the case where f is given by ( 11) (with θ 2 < 0 in order to have a function concave in x), one gets

H( θ0 ) = (-1/ θ0 2 )   0 0 0 0 1 x * ( θ0 ) 0 x * ( θ0 ) [x * ( θ0 )] 2   = -θ0 2 c( θ0 )c ⊤ ( θ0 ) .
This suggests substitution of H( θk ) for H in [START_REF] Fel | Dual control theory[END_REF] in the case of sequential design, with θk a Bayesian estimator of θ. Again, in the infinite horizon case with no discount, the choice of a weighting sequence {α k } in (5) ensuring convergence of x * ( θk ) to x * ( θ) and of ξ k to ξ x * ( θ) remains an open issue.

Linear response optimisation with finite horizon

Assume the the errors ǫ k in (1) are i.i.d. N (0, σ 2 ). We shall use an expansion in σ 2 to obtain an approximate solution to the problem (4). We first need to prove the following Lemma.

Lemma 1 Let g(x) and h(x) be two times continuously differentiable functions on X , a compact set of IR q . Assume that g has a unique global maximum at x * , an interior point of X , with ∂ 2 g(x)/∂x∂x ⊤ |x * negative definite, and let x denote the point where f (x) = g(x) + uh(x) reaches its maximum in X . Then, x -

x * = O(u) and |f (x * ) -f (x)| = O(u 2 ), u → 0.
Proof For u small enough, x is an interior point of X so that

∂f (x) ∂x ⊤ |x = 0 = ∂g(x) ∂x ⊤ |x * + (x -x * ) ⊤ ∂ 2 g(x) ∂x∂x ⊤ |x * + u ∂h(x) ∂x ⊤ |x + o( x -x * ) = (x -x * ) ⊤ ∂ 2 g(x) ∂x∂x ⊤ |x * + u ∂h(x) ∂x ⊤ |x + o( x -x * ) , and x -x * = O(u). Therefore, f (x * ) = f (x) + (x * -x) ⊤ ∂f (x) ∂x |x + 1 2 (x * -x) ⊤ ∂ 2 f (x) ∂x∂x ⊤ |x (x * -x) + o( x -x * 2 ) = f (x) + O(u 2 ) .
Next theorem shows how to construct an approximate solution to problem (4). Theorem 2 Assume that w i > 0, i = 1, . . . , N and w i = 0 otherwise; that f is linear in θ, that is, f (θ, x) = r ⊤ (x)θ, with r two times continuously differentiable in x; that the errors ǫ i are i.i.d. N (0, σ 2 ) and the prior distribution for θ is normal N ( θ0 , σ 2 Ω). Denote x * (θ) = arg max x∈X r ⊤ (x)θ and θk = E{θ|F k }, and assume that f ( θj , x) has a unique global maximum at x * ( θj ) which lies in the interior of X , with ∂ 2 f ( θj , x)/∂x∂x ⊤ |x * ( θj ) negative definite, j = 0, . . . , N -2. Define j k+1 (x) as the expected optimal gain to go at step k when x is applied:

j k+1 (x) = E{w k+1 f (θ, x) + max z∈X w k+2 j k+2 (z)|F k } , k = 0, . . . , N -2 , j N (x) = E{w N f (θ, x)|F N -1 } . It satisfies j k+1 (x) = J k+1 (x) + O(σ 4 ) , k = 0, . . . , N -2 , ( 14
)
where

J k+1 (x) = (w N + • • • + w k+2 ) × r ⊤ [x * ( θk )] θk + σ 2 2 trace [H( θk )(Ω -1 + M k ) -1 ] - σ 2 2 trace    H( θk ) N -k-2 j=0 w k+j+2 [Ω -1 j,k + r(x)r ⊤ (x)] -1    +w k+1 r ⊤ (x) θk , (15)
with

Ω j,k = Ω -1 + M k + jr[x * ( θk )]r ⊤ [x * ( θk )] -1
and M k , H(θ) respectively given by ( 7), [START_REF] Fel | The theory of dual control[END_REF]. Moreover, the strategy defined by x N = x * ( θN-1 ) and x k+1 = arg max x∈X J k+1 (x), k = 0, . . . , N -2, satisfies

x k+1 = arg max x∈X w k+1 r ⊤ (x) θk + σ 2 2 × N -k-2 j=0 w k+j+2 r ⊤ (x)Ω j,k H( θk )Ω j,k r(x) 1 + r ⊤ (x)Ω j,k r(x) , k = 0, . . . , N -2 , (16)
E N i=k+1 w i f (θ, x i )|F k = J k+1 (x k+1 ) + O(σ 4 ) , k = 0, . . . , N -2 , ( 17
)
and xk+1 -x k+1 = O(σ 4 ) , k = 0, . . . , N -2 , xN = x N , (18)
where xk+1 = arg max x∈X j k+1 (x) corresponds to the optimum strategy.

Proof Straightforward matrix manipulation shows that x k+1 given by ( 16) maximises [START_REF] Ginebra | Response surface bandits[END_REF]. We prove ( 14) and ( 17) by backward induction on k. For k = N -2, we have

j N -1 (x) = w N E{r ⊤ [x * ( θN-1 )]θ|F N -2 } + w N -1 r ⊤ (x) θN-2 = w N E{r ⊤ [x * ( θN-1 )] θN-1 |F N -2 } + w N -1 r ⊤ (x) θN-2 .
When x * ( θN-2 ) is an interior point of X , a second-order Taylor expansion around θN-2 similar to that used in [START_REF] Pronzato | Experimental design for estimating the optimum point in a response surface[END_REF] gives

j N -1 (x) = w N {r ⊤ [x * ( θN-2 )] θN-2 + σ 2 2 trace [H( θN-2 )(Ω -1 + M N -2 ) -1 ]} -w N σ 2 2 trace {H( θN-2 ) r(x)r ⊤ (x) + Ω -1 + M N -2 -1 } + w N -1 r ⊤ (x) θN-2 + O(σ 4 ) , which proves (14) for k = N -2. Since E{ N i=N -1 w i f (θ, x i )|F N -2 } = j N -1 (x N -1 ), it also proves (17) for k = N -2.
Assume that ( 14) is true at step k. At step k -1 we have j k (x) = E{w k f (θ, x) + max z∈X j k+1 (z)|F k-1 }. Using Lemma 1 with u = σ 2 , we get

j k (x) = w k r ⊤ (x) θk-1 + E{J k+1 [x * ( θk )]|F k-1 } + O(σ 4 ) .
Grouping the terms r ⊤ [x * ( θk )] θk and using a second-order Taylor expansion around θk-1 , we obtain

j k (x) = (w N + • • • + w k+2 + w k+1 ) r ⊤ [x * ( θk-1 )] θk-1 + σ 2 2 trace [H( θk-1 )(Ω -1 + M k-1 ) -1 ] -(w N + • • • + w k+2 + w k+1 ) σ 2 2 × trace {H( θk-1 )[Ω -1 + M k-1 + r(x)r ⊤ (x)] -1 } +(w N + • • • + w k+2 ) σ 2 2 × trace {H( θk-1 )[Ω -1 + M k-1 + r(x)r ⊤ (x)] -1 } - σ 2 2 trace    H( θk-1 ) N -k-2 j=0 w k+j+2 [Ω -1 j+1,k-1 + r(x)r ⊤ (x)] -1    +w k r ⊤ (x) θk-1 + O(σ 4 ) .
A simplification of the terms on the last three lines gives j k (x) = J k (x) + O(σ 4 ), which proves [START_REF] Ford | A sequentially constructed design for estimating a nonlinear parametric function[END_REF]. Since x k+1 maximises (15), similar arguments using ( 17) and Lemma 1 give

E{ N i=k w i f (θ, x i )|F k-1 } = J k (x k ) + O(σ 4
), and ( 17) is proved.

Finally, since x k+1 maximises J k+1 (x), ( 14) and Lemma 1 with u = σ 4 give [START_REF] Kpamegan | An optimizing up-and-down design[END_REF].

Remark 1 The assumptions on f ( θj , x) and x * ( θj ) for j = 0, . . . , N -2, in Theorem 2 are most often difficult, if not impossible, to check beforehand. Note, however, that it is always possible to apply the strategy and check the assumptions afterwards: if they are satisfied, the theorem applies and the strategy used is approximately optimal in the sense of the theorem.

The property stated in Theorem 2 suggests a simpler suboptimal strategy, as given in the following corollary.

Corollary 2 Under the same conditions, and with the same notations, as in Theorem 2, the strategy defined by x ′ N = x * ( θN-1 ) and

x ′ k+1 = arg max x∈X w k+1 r ⊤ (x) θk + σ 2 2 N -k-1 j=1 j w k+j+1 r ⊤ (x)Ω 0,k H( θk )Ω 0,k r(x) 1 + jr ⊤ (x)Ω 0,k r(x) , k = 0, . . . , N -2 , ( 19
)
is approximately optimal in the following sense:

max x∈X j k+1 (x) -E{ N i=k+1 w i f (θ, x ′ i )|F k } = O(σ 4 ) , ( 20
)
|x k+1 -x ′ k+1 | = O(σ 2 ) , (21)
for k = 0, . . . , N -1.

Proof We first prove that E{

N i=k+1 w i f (θ, x ′ i )|F k } = J k+1 (x ′ k+1 ) + O(σ 4 ), for k = N -2, . . . , 1.
For k = N -2 we have again

E{ N i=N -1 w i f (θ, x ′ i )|F N -2 } = j N -1 (x ′ N -1 ) which equals J N -1 (x ′ N -1 ) + O(σ 4
) from Theorem 2. Assume that the property is true at step k. At step k -1 we have

E{w k f (θ, x) + N i=k+1 w i f (θ, x ′ i )|F k-1 } = w k r ⊤ (x) θk-1 +E{J k+1 (x ′ k+1 )|F k-1 } + O(σ 4 ) ,
and easy matrix manipulation using ( 19) and [START_REF] Ginebra | Response surface bandits[END_REF] shows that

x ′ k+1 = arg max x∈X [J k+1 (x) + σ 2 h(x)], for some h(x). Lemma 1 with u = σ 2 gives J k+1 (x ′ k+1 ) = J k+1 [x * ( θk )] + O(σ 4
), and similarly to the proof of Theorem 2, E{w k f (θ, x)

+ N i=k+1 w i f (θ, x ′ i )|F k-1 } = J k (x) + O(σ 4
). Finally, since the optimal strategy xk+1 maximises j k+1 (x) = J k+1 (x) + O(σ 4 ), whereas x ′ k+1 maximises a function that takes the form J k+1 (x) + σ 2 h(x), Lemma 1 gives [START_REF] Lai | Asymptotic properties of nonlinear least squares estimates in stochastic regression models[END_REF][START_REF] Lai | Least squares estimates in stochastic regression models with applications to identification and control of dynamic systems[END_REF]. When w i = 1, i = 1, . . . , N and w i = 0 otherwise, with N large, the strategy ( 19) can be further approximated by

x ′′ k+1 = arg max x∈X r ⊤ (x) θk + (N -k -1) σ 2 2 r ⊤ (x)Ω 0,k H( θk )Ω 0,k r(x) r ⊤ (x)Ω 0,k r(x) , k = 0, . . . , N -2 , ( 22
)
which has the form (5).

Remark 2 It is clear from the proof of Corollary 2 that properties similar to [START_REF] Lai | Asymptotic properties of nonlinear least squares estimates in stochastic regression models[END_REF][START_REF] Lai | Least squares estimates in stochastic regression models with applications to identification and control of dynamic systems[END_REF] can be obtained for other strategies than [START_REF] Kumar | Convergence of adaptive control schemes using least-squares parameter estimates[END_REF]. What makes this strategy attractive is that it amounts to substituting r(x) for r[x( θk )] in J k+1 (x) to construct x ′ k+1 . This makes the rule ( 19) suboptimal compared to ( 16) that directly maximises J k+1 (x) (compare ( 21) to ( 18)), but, on the other hand, r(x ′ k+1 ) can be expected to be close to r[x( θk )] and the loss can be expected to be small.

Example 1

We assume that f (θ, x) is given by [START_REF] Fel | Dual control theory[END_REF], with X = [-1, 1], and that the observations y k are generated by (1) for θ = θ = (0 3.2 -8) ⊤ with the errors ǫ k i.i.d. N (0, σ 2 ), σ = 1. The prior for θ is N ( θ0 , σ 2 Ω), with θ0 = (2 -4 -1) ⊤ and Ω = 10 6 I 3 , with I 3 the 3-dimensional identity matrix. Note that these numerical values give a prior guess for x * at -4, whereas the true location is at 0.4.

We take w i = 1, i = 1, . . . , N and w i = 0 otherwise, with N = 100, and compare four different strategies: S 1 corresponds to (5) with the penalty [START_REF] Dette | A note on Bayesian cand D-optimal designs in nonlinear regression models[END_REF] and α k = σ 2 (log k) 2 , and S 2 , S 3 and S 4 correspond respectively to [START_REF] Müller | Batch sequential design for a nonlinear estimation problem[END_REF], [START_REF] Kumar | Convergence of adaptive control schemes using least-squares parameter estimates[END_REF] and [START_REF] Hardwick | Optimizing a unimodal response function for binary variables[END_REF]. Table 1 presents the results obtained for 100 independent repetitions of the experiment (for each experiment, the same values of observations errors ǫ k are used for the four strategies). As expected, performances improve from S 1 to S 4 .

Figure 1 (resp. 2) presents a typical realization of the sequences {x k } (resp. {f ( θ, x k )}) generated by the four strategies. One can notice in Figure 1 that x k converges to x * = -θ1 / θ2 = 0.4 for S 1 (see Theorem 1), which does not seem to be the case for the other strategies. However, the performance measured in terms of ( 4) is much better for the other strategies, in particular S 4 that makes a particularly good compromise between estimation and optimisation, see Figure 2 where the optimum value f [ θ, x * ( θ)] is indicated by the dashed line.

4 Extension to nonlinear situations

Nonlinear regression

Whereas Theorem 1 relies on Corollary 3 of [START_REF] Lai | Least squares estimates in stochastic regression models with applications to identification and control of dynamic systems[END_REF], deriving similar asymptotic results in the case of a nonlinear regression model would require using Theorem 2 of [START_REF] Lai | Asymptotic properties of nonlinear least squares estimates in stochastic regression models[END_REF]. The conditions in [START_REF] Lai | Asymptotic properties of nonlinear least squares estimates in stochastic regression models[END_REF] being more stringent than those in [START_REF] Lai | Least squares estimates in stochastic regression models with applications to identification and control of dynamic systems[END_REF], one may conjecture that α k should increase faster than any power of log k when k → ∞, compare with Theorem 1. Using Bayesian estimation instead of LS would certainly require weaker conditions on α k , as it was the case for Corollary 1 in comparison with Theorem 1, but deriving precise results on this is still an open issue.

In the case of finite horizon, the approaches of Section 3 can be used through a linearisation of the model response. Notice, however, that this implies that terms of order O(σ 2 ) are neglected, so that the properties [START_REF] Hu | Strong consistency of Bayes estimates in stochastic regression models[END_REF] and [START_REF] Kpamegan | An optimizing up-and-down design[END_REF] are no longer true.

Example 2: Bernoulli trials

We consider the situation defined by (2), with

f (θ, x) = θ 1 + θ 2 θ 2 [1 -exp(-x/θ 1 )] exp(-x/θ 2 ) ,
with X = {x min , x min + δ, x min + 2δ, . . . , x max }, where x min = 2, x max = 100, δ = 2, which gives #X = 50. This may correspond to a dose-response problem, where efficiency is small at low doses x and toxic failures occur at large doses, see, e.g., [START_REF] Kpamegan | An optimizing up-and-down design[END_REF]. We assume that θ has a normal prior N ( θ0 , Ω), with θ0 = (70 70) ⊤ and Ω = 10 3 I 2 . The observations are generated with θ = θ = (3 50) ⊤ . Figure 3 presents the two probabilities f ( θ, x) (full line) and f ( θ0 , x) (dashed line) as functions of x. The prior guess x * ( θ0 ) is at 48.52, whereas the true value is at x * ( θ) = 8.615. We estimate θ by the maximum a posteriori estimator and denote θk = arg max

k i=1 log{y i f (θ, x i ) + (1 -y i )[1 -f (θ, x i )]} - 1 2 (θ -θ0 ) ⊤ Ω -1 (θ -θ0 ) .
The expression of f (θ, x) gives x * (θ) = -θ 1 log[θ 1 /(θ 1 + θ 2 )] and

c(θ) = ∂x * (θ) ∂θ = -θ 2 /(θ 1 + θ 2 ) -log[θ 1 /(θ 1 + θ 2 )] θ 1 /(θ 1 + θ 2 ) .
The design matrix M k is replaced by M k ( θk ), with that would be the Fisher information matrix evaluated at θ if the design levels x i were nonrandom constants. Note that it still has the form [START_REF] Durham | A sequential design for maximizing the probability of a favorable response[END_REF], with now

M k (θ) = k i=1 ∂f (θ, x i ) ∂θ ∂f (θ, x i ) ∂θ ⊤ 1 f (θ, x i )[1 -f (θ, x i )] , 0
r(x) = r(θ, x) = 1 {f (θ, x)[1 -f (θ, x)]} 1/2 ∂f (θ, x) ∂θ .
We shall compare four strategies. S 1 and S 2 are defined by ( 5) with the weights α k = k 1/4 , and penalty functions ( 6) and ( 12) respectively. S 3 corresponds to a randomized Polya urn design, see [START_REF] Durham | A sequential design for maximizing the probability of a favorable response[END_REF], and S 4 to an optimizing up-and-down design, see [START_REF] Kpamegan | An optimizing up-and-down design[END_REF]. S 3 is initialized by putting one ball of each label in the urn, S 4 samples by pairs, at (P kδ/2, P k + δ/2) at step k, and we take P 1 at the median of X .

Figure 4 (resp. 5) presents a typical realization of the sequences {x k } (resp. {f ( θ, x k )}) generated by the four strategies. On Figure 5, the optimum value f [ θ, x * ( θ)] is indicated by the dotted line. Figure 6 presents histograms of the values of x k generated by the four strategies.

Convergence to the maximum of f (θ, x) is faster for S 1 and S 2 than for S 3 and S 4 . The difference in performance is even stronger when δ decreases, that is, when the number of points in X increases. Whereas the value of δ has no significant effect of S 1 and S 2 , decreasing δ makes the behaviour of S 3 more erratic and increases the time S 4 requires to move away from the initial values of x k . This effect is already apparent on Figures 4 to 6. A possible compromise strategy would be to start with S 3 or S 4 and a large value of δ (that is, few point in X ), and then decrease δ and switch to S 1 or S 2 . form is known for f , that is, f (x) = f (θ, x) with θ unknown, different strategies have been suggested.

An approximately optimal strategy has been constructed in the case of a linear regression model with finite horizon. In more general situations, we considered sequences constructed according to the rule x k+1 = arg max x f ( θk , x) + α k d k (x), with θk an estimated value of θ obtained from (x 1 , y 1 ), . . . , (x k , y k ) and d k (x) a penalty for poor estimation. This is easily implemented in the case of binary responses, and numerical examples illustrate that such strategies may outperform more traditional ones, like Polya urn designs and up-and-down methods. Such a result should not be a surprise, since methods based on a parametric form of f use more information than model-free approaches. On the other hand, their asymptotic behaviour is more difficult to study, due to the intricate connection between estimation of θ and optimisation. Only the linear regression case, with a penalty related to D-optimum design, seems to have been considered so far, see [START_REF] Pronzato | Adaptive optimisation and D-optimum experimental design[END_REF]. Extension to nonlinear regression problems is currently under study. In practice, it might be reasonable to start with a model-free approach, in order to gather enough information about the form of f , and then switch to a parametric approach of the type considered here. This raises many issues (choice of model structure, choice of switching time, robustness with respect to model mis-specification, etc.) which, together with comparisons with other approaches (e.g., the bandit design of [START_REF] Hardwick | Optimizing a unimodal response function for binary variables[END_REF]), deserve further study.

Fig. 1

 1 Fig. 1 Sequences {x k } generated by strategies S1 to S4 (from left to right) in Example 1

1 Fig. 2

 12 Fig. 2 Sequences {f ( θ, x k )} generated by strategies S1 to S4 (from left to right) in Example 1. The optimum value f [ θ, x * ( θ)] is indicated by the dashed line

Fig. 3

 3 Fig. 3 Probabilities f ( θ, x) (full line) and f ( θ0 , x) (dashed line) as functions of x in Example 2

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4 Sequences {x k } generated by strategies S1 to S4 (from left to right) in Example 2

Table 1

 1 Empirical means and standard deviations (std) of (1/N )

	N i=1 f ( θ, xi)

Acknowledgments. This work was presented at a meeting in Bertinoro on December 2000. Prof. A. Giovagnoli from the University of Bologna, who organised the meeting, is gratefully acknowledged for her invitation.