Aude Bolopion
email: aude.bolopion@femto-st.fr.

Christian Dahmen

Christian Stolle
email: christian.stolle@uni-oldenburg.de.

Sinan Haliyo
email: haliyo@isir.upmc.fr.

Stéphane Régnier
email: regnier@isir.upmc.fr.

Sergej Fatikow

* A Bolopion

S T Éphane

R Égnier

Vision-based haptic

HAL is a

VISION BASED HAPTIC FEEDBACK FOR REMOTE MICROMANIPULATION IN A SEM ENVIRONMENT

AUDE BOLOPION * ,1 , CHRISTIAN DAHMEN * ,2 , CHRISTIAN STOLLE * ,2 , SINAN HALIYO 3 , ST ÉPHANE R ÉGNIER 3 AND SERGEJ FATIKOW 2 ABSTRACT. This paper presents an intuitive environment for remote micromanipulation composed of both haptic feedback and virtual reconstruction of the scene. To enable non expert users to perform complex teleoperated micromanipulation tasks it is of utmost importance to provide them with information about the 3D relative positions of the objects and the tools. Haptic feedback is an intuitive way to transmit such information. Since position sensors are not available at this scale, visual feedback is used to derive information about the scene. In this work, three different techniques are implemented, evaluated and compared to derive the object positions from scanning electron microscope images. The modified correlation matching with generated template algorithm is accurate and provides reliable detection of objects. To track the tool, a marker based approach is chosen since fast detection is required for stable haptic feedback. Information derived from these algorithms is used to propose an intuitive remote manipulation system, that enables users situated in geographically distant sites to benefit from specific equipments such as SEMs. Stability of the haptic feedback is ensured by the minimization of the delays, the computational efficiency of vision algorithms and the proper tuning of the haptic coupling. Virtual guides are proposed to avoid any involuntary collisions between the tool and the objects. This approach is validated by a teleoperation involving melamine microspheres with a diameter of less than 2 µm between Paris, France and Oldenburg, Germany.

INTRODUCTION

Versatile micromanipulation systems would open a wide range of applications, that would benefit different areas, such as medicine, electronics or physics. To enable non expert users to perform micromanipulation tasks intuitiveness of the system is a major requirement. In particular it is necessary to provide users with information about the 3D relative positions of the tools and the objects. So that intuitiveness is ensured it is possible to enhance the visual feedback by reconstructing virtual scenes [START_REF] Kim | Task-based and stable tele-nanomanipulation in a nanoscale virtual environment[END_REF][START_REF] Ferreira | Virtual reality and haptics for nanorobotics[END_REF]. In addition to visual display haptic feedback can be provided. Operators manipulate the tools through a joystick, and receive force feedback [START_REF] Bukusoglu | Haptic manipulation of microspheres using optical tweezers under the guidance of artificial force fields[END_REF][START_REF] Hui Xie | Atomic force microscopy based nanorobotics, modelling, simulation[END_REF]. To provide such assistance sensing capabilities are necessary to get information about the scene. Specific sensors, adapted to microscale requirements, have been developped [START_REF] Fahlbusch | Force sensing in microrobotic systems-an overview[END_REF][START_REF] Beyeler | A Six-Axis MEMS Force-Torque sensor with Micro-Newton and Nano-Newtonmeter resolution[END_REF]. However, their integration into the dedicated tools is complex, and greatly increases the difficulty of their design and their conception.

Vision is thus a promising solution to compensate for this lack of sensing [START_REF] Greminger | Vision-based force measurement[END_REF][START_REF] Pawashe | Two-dimensional vision-based autonomous microparticle assembly using nanoprobes[END_REF][START_REF] Narayana Reddy | Miniature compliant grippers with Vision-Based force sensing[END_REF][START_REF] Chawda | Vision-Based force sensing for nanomanipulation[END_REF].

Visual feedback obtained from classical optical microscopes is limited to objects whose size is greater than a few micrometers. To manipulate smaller objects specific equipments, Scanning Electron Microscopes (SEM), must be used to visualize the scene. For detection and tracking of objects and tools in SEM images, common algorithms include rigid-body based matching [START_REF] Kratochvil | Real-time rigid-body visual tracking in a scanning electron microscope[END_REF], active contours [START_REF] Sievers | Real-time object tracking inside an sem, Automated Nanohandling by Microrobots[END_REF] and template matching using cross correlation [START_REF] Sievers | Real-time object tracking for the robot-based nanohandling in a scanning electron microscope[END_REF]. Though these algorithms are flexible and robust, not all of them are suitable for the simultaneous detection of a high number of simple objects, like melamine spheres. They must also be computationally efficient to provide haptic feedback and real time reconstruction of the virtual scene.

In this work three different detection algorithms are implemented, evaluated and compared in terms of computation time, accuracy and reliability. To provide haptic feedback and real time reconstruction of the scene the modified correlation matching with generated template algorithm is used to get information about object positions. To track the tool, a marker based approach is chosen since fast detection is required. As the proposed detection algorithms are computationally efficient, they can be used to provide an intuitive remote teleoperation environment with stable haptic feedback. To validate this approach, a teleoperated task involving microspheres is performed between France and Germany. The haptic feedback teleoperation system is detailed. A special care has been taken to ensure stability by minimizing the communication delays and the vision detection computation time. Virtual guides are proposed to avoid any involuntary collision between the objects and the tool, by the transmission of adapted haptic feedback. A paper describing the software architecture of this system, that ensures the minimization of time delays as well as modularity, has been published in [START_REF] Bolopion | Remote microscale teleoperation through virtual reality and haptic feedback[END_REF]. The work presented here concentrates on computationally efficient and robust detection algorithms for SEM images in order to propose vision-based haptic feedback and virtual reconstruction of the scene. Another approach is presented in [START_REF] Bolopion | Haptic teleoperation for 3D microassembly of spherical objects[END_REF], where the relative position of the tool (a cantilever) and the objects is derived from direct measurement: the cantilever is excited at its resonant frequency, and the amplitude of the oscillations is monitored. It decreases while the cantilever approaches an object. However, it requires bringing the cantilever in the close vicinity of the object, and it does not give indication about the whole scene. In particular, it is not possible to provide an updated virtual reconstruction of the scene during the manipulation task based on direct measurements. This paper is organized as follows. Section 2 describes the micromanipulation system.

Vision algorithms, adapted to the specificities of SEM images, are proposed in Section 3 to track both the tool and the objects. Based on the position information derived from these algorithms, an intuitive manipulation system composed of both a virtual scene reconstruction and haptic feedback is proposed in Section 4. This approach is validated by experiments conducted on microspheres using teleoperation between France and Germany.

MICROMANIPULATION PLATFORM

2.1. Setup. Cantilevers are widely used inside Atomic Force Microscopes (AFM) to manipulate micron-sized objects since they are commercially available and rather cheap. Several strategies have been developed, such as pushing, rolling, or pick-and-place by adhesion, or by using two tips as a microgripper [START_REF] Hui Xie | Atomic force microscopy based nanorobotics, modelling, simulation[END_REF]. To derive interaction forces involved while manipulating objects the deformations applied on the cantilever are measured , either by using optical systems composed of lasers reflected on photodiodes, or by using piezoresistive sensors. However, these measurements are mostly limited to the vertical deflection of the cantilever, which is related to a vertical force applied on the tool. Due to the design of the cantilever, the resolution of torsional forces is too low to be used as force feedback for a manipulation task. In this work a commercially available piezoresistive contact type cantilever 1 with a stiffness of 4 N.m -1 is used. Vertical forces applied on the cantilever are translated into a change in resistance and measured as voltage change.

Voltage measurements are made using an amplifier 2 and a 16-bit analog digital converter 3 .

The cantilever is used to manipulate commercially available melanin microspheres4 (diameter: 1 -1.6 µm) loosely deployed on a silicon plate (5 mm x 5 mm). It is mounted on the fine positioning part of a nanorobotic setup [START_REF] Mick | Combined nanorobotic afm/sem system as novel toolbox for automated hybrid analysis and manipulation of nanoscale objects[END_REF]. This positioning part is piezo driven with a build-in capacitive position sensors for closed-loop movement control and high accuracy positioning 5 . Its lateral stroke is up to 100 µm and the z-range is up to 50 µm and the positioning accuracy is at the lower nanometer scale. The specimen holder has been mounted on a coarse positioning unit. This unit consists of three slip stick driven, orthogonal mounted linear axes 6 . These positioners are equipped with optical positioning sensors allowing for travel ranges of several centimeters with a repositioning accuracy up to 50 nm.

The whole setup is placed inside a high-resolution scanning electron microscope with a Schottky-emitter 7 that is also equipped with a FIB column 8 and a gas injection system.

Such microscopes provides images of objects smaller than a few micrometers, sufficiently precise to determine the relative position of the tool and the objects for manipulation tasks. This information is used in this work to set up an intuitive remote micromanipulation platform. The teleoperation platform is composed of an Omega haptic device 9 used to control the position of the cantilever. A haptic feedback is sent through the device to assist operators while manipulating objects. This master arm is a 3 degrees of freedom device for both displacement and haptic feedback. To avoid damages of the haptic interface forces higher than 3 N are truncated. Users are also immersed in a virtual reality room which comprises a 3D display module. A reconstructed image of the scene is projected in stereoscopic 3D on a 1.80 × 2.4 m 2 screen.

2.2. Software Architecture. The robotic control system is based on the distributed control system for automated micro-and nano robotics (DCAAN) [START_REF] Stolle | Distributed control architecture for automated nanohandling[END_REF]. It consists in several servers written in C++ which are connected via CORBA (Common Object Request Broker Architecture) based on a TCP/IP network. Due to the distributed nature of the tele-control task the software has been subdivided into two parts (see Fig. 1). The one located in Paris is the operator side, and the one in Oldenburg is the manipulation setup. At the remote side three different programs are running:

• the sensor provider Olvis is responsible for processing image based information such as video streams [START_REF] Wortmann | Image processing architecture for real-time micro-and nanohandling applications[END_REF]. Depending on these algorithms different kinds of information such as sphere and cantilever positions can be derived (see Section 3),

• the AFM low level control server is processing the tele-control commands by moving the cantilever tip to the target position and sending an estimation of the actual position and a measure of the force back in return,

• the SEM low level control server enables in particular the operator to turn the experimental setup towards the focus of the SEM by moving the SEM stage, and to change image parameters such as brightness and contrast. These settings are indeed of primary importance to ensure efficient vision tracking and detection. 9 Force Dimension, http://www.forcedimension.com The Omega low level control server forms the heart of the operator side software and provides a small graphical user interface which enables the operator to monitor information such as force and position of the cantilever. In addition to the numerical display a live SEM image (see Section 3) can be provided to the operator via the GUI. The Omega low level control server is also responsible for translating the force and position feedback of the AFM-cantilever and forces and positions of the haptic device into each other.

The communication between the operator and the remote sides has been realized directly via UDP packages. The communication delay has been measured by the means of the round trip time (RTT) by requesting the current position of the cantilever 100 times.

The average RTT has been determined at 37 ms, the minimum RTT at 36 ms and the maximum RTT at 39 ms. The distributed control system enables to minimize RTT, and to ensure modularity. It is thus suitable to perform a wide range of teleoperation over important distances.

SEM BASED VISION ALGORITHMS FOR TRACKING AND DETECTION

In SEM-based micro-and nanorobotic handling, visual feedback from the SEM is one of the main sources of information about the workspace. The SEM feedback is used during the registration process where the positions of parts and robot tools are determined, as well as for position feedback during the handling process itself. The SEM as sensor itself introduces some constraints on detection and tracking algorithms, which must be adapted to the specifics of SEM images. These constraints include e.g. robustness against stochastic additive noise and changes in brightness or contrast. In addition the update rate must be sufficient to meet the requirements of haptic feedback. In our setup the images are acquired with a resolution of 627 x 627 pixels. The field of view for the handling scene is in the range of 10 to 100 µm.

The cantilever is tracked with a template matching algorithm (TM) [START_REF] Goshtasby | A two-stage cross-correlation approach to template matching[END_REF], which is one of the fastest and the most common real-time marker tracking approaches. It finds the position in the image with the best correspondence with the template. To achieve good tracking results, a marker has been structured by the FIB on top of the cantilever. This marker is used as the template for the TM, and its unique features enable to precisely determine the position of the cantilever. In the initialization step the marker is tagged manually. Automated detection will be considered in future works. To calculate the position of the cantilever tip, a offset is added to the position of the tracking marker. To be used in teleoperation, the update rate of the position of the cantilever has to be as fast as possible. Therefore the process is set to the highest priority in the Olvis framework to achieve an update rate of up to 50 Hz.

The detection of the melamine spheres is more difficult in contrast to the cantilever because the appearance of the spheres changes during the handling (Fig. 2). The cantilever shadows the spheres due to position changes, which modifies the intensity and contrast of the image. Due to the view angle depended detector shadowing effects and different brightness and texture features of the spheres, it has been decided to work on edge images to generate comparable results. Algorithms applied to this scenario need to be developed to especially meet these requirements. For this purpose three different detection methods have been selected, tested and compared:

• Contour detection with circle fitting criteria's,

• Modified correlation matching with generated templates,

• Hough circle detection. The first approach uses contour detection with circle fitting criteria. After the edge detection is applied, the contours are detected using the OpenCV [START_REF] Bradski | Programmers tool chest: The opencv library[END_REF] contour detection (cvContour) and ellipse fitting (cvFitEllipse2) functions. As a result, a number of found contours with different parameters are returned. To find the contours, which are mostly fitting a sphere, different criteria have to be calculated: • Area (A) of the bounding box, to limit the diameter search range of the spheres,

• Number of segments.

These four criteria are used to calculate a score. After applying a threshold to filter the scores, a list of validly detected sphere contours is given and the center points are calculated.

The second approach uses modified cross correlation (CC) [START_REF] Goshtasby | A two-stage cross-correlation approach to template matching[END_REF] with generated templates, similar to the ones applied during cantilever tracking. After applying the edge detection filter a CC is used to find the spheres. The formula beneath shows the mathematical basics of the used normalized cross correlation function (NCCF). In the formula w(u, v) are pixel values in the correlation window and G(x, y) are the pixel values from the search region with l and k smaller then the image size m x n of the CC function (CCF): The normalization term (NT) is

(3.3) CCF(x, y) = l-1 ∑ u=0 k-1 ∑ v=0 w(u, v) * G(x + u, y + v).
(3.4) NT (x, y) = L-1 ∑ u=0 K-1 ∑ v=0 w(u, v) 2 * L-1 ∑ u=0 K-1 ∑ v=0 G(x + u, y + v) 2 .
Finally, (3.5) NCCF(x, y) = CCF(x, y) NT (x, y) .

In contrast to normal CC templates, binary ring templates are used. These ring templates are specially generated to fit to the contour of a circle, independently of the texture values of the sphere. Therefore, a list of ring templates with different radii and ring widths are automatically generated and passed to the CC. Finally, the result list is clustered and the weighted center points are calculated. The clustering is necessary, because the same sphere can be multiply detected on a slightly different position in each CC cycle.

The third approach uses a Hough circle detection [START_REF] Gonzales | Wood, Digital image processing[END_REF] on an edge image with different radii. Therefore a modified version was implemented. For the same reason as described in the CC approach, the result list has to be clustered and finally the weighted center points of the detected spheres are calculated. The Clustering of all three methods was done by using the OpenCV function (cvSeqPartition).

All three methods have been evaluated with respect to the special requirements of AFMbased sphere handling inside the SEM. The contour detection method offers fast computation time which is necessary for teleoperation (otherwise haptic feedback is unstable [START_REF] Colgate | Passivity of a class of sampled-data systems: Application to haptic interfaces[END_REF]).

The disadvantages of this method are high false detection rates, due to increasing image noise and comparatively small spheres.

The CC with ring templates and the Hough circle detection both have a higher computation time in comparison to the contour detection method. Their advantage is higher accuracy and better reliability. On the one hand, the CC with ring templates has a higher reliability in blurry images and images with many small spheres but on the other hand, the Hough circle detection has the better computation time. Currently the computation time is in the range of 0.4 s to 2 s, depending on the density of input edge images. Another advantage of these two methods is that their algorithms can be scheduled in parallel. Currently both algorithms are implemented in OpenCl [START_REF] Trevett | The open standard for heterogeneous parallel programming[END_REF] to work on graphic card units and to optimize the calculation time.

In Figure 4, results of sphere detection of the scenes of Fig. 2 using hough circle detection can be seen. Marked are true positives in green, false positives in blue, and false negatives in red. The successful and failed detection have been determined by hand. The results of the evaluation can be seen in Table 1. As can be seen, the false positives are relatively low count. This can be still improved by further optimization of the parameters. The main errors are false negatives, which clearly depend on the various factors influencing the image as stated before.

To meet the teleoperation requirements of fast and accurate detection, as well as reliability, the CC based approach has been used for the handling experiments. In particular, the 1. Results of the evaluation from the three scenes using hough circle approach (see Figure 4). algorithm integrates SEM parameters like view field to achieve a magnification invariant detection. However, in order to avoid a slowdown of the cantilever tracking, the positions of the melamine spheres are calculated with a low priority in the Olvis framework. This choice was made since the cantilever is moved more often than the spheres. The sphere detection has an update rate of 3 to 8 Hz for searching in the sphere diameters in the range of 1 to 1.6 µm. This detection algorithm, as well as the one for the cantilever detection, have a sub-pixel accuracy.

During this application the position of the involved objects was only determined in x and y directions of the SEM image plane. The relative height of the cantilever with respect to the substrate is thus unknown. For this work, this lack of information is compensated by the force measurement of the piezoresistive sensor, which estimates the vertical force applied on the cantilever and enables to determine if the tip is or not in contact with the substrate or an object. However, this method does not compensate for the unknown orientation of the substrate surface and the image plane. This may lead in turn to an ellipsoidic shaped force gradient while the cantilever is in contact with the surface. One solution would be to include a calibration step, or image processing based on the distance between the tip and its shadow on the substrate.

INTUITIVE TELEOPERATION ENVIRONMENT

Based on the position information of the tool with respect to the object that is derived by vision algorithms, an intuitive remote manipulation system is proposed. It is composed of both visual and haptic feedback. Intuitive remote manipulation systems will enable non expert users to get access to rare and expensive equipments, that are often situated in geographically distant sites. 4.1. The Intuitive Remote Manipulation System. To assist users while performing a manipulation, a stereoscopic virtual reality scene is reconstructed using Blender, an open source 3D content creation suite [START_REF]Blender[END_REF]. The scene includes the AFM tip used as the manipulation tool and a sphere as the manipulated object. In that work, it was decided to represent only the sphere closest to the cantilever tip. This facilitates the manipulation for untrained users by displaying only important information on the reconstructed scene. The scene is constructed setting 1 Blender unit as 1 µm. All along the operation, the user has complete freedom to modify the camera position and zoom. The drawback of this method is the requirement of prior knowledge of geometries and dimensions and an initial calibration process between real image and virtual reality frames. As in this case the cantilever's and manipulated objects' geometries are perfectly known and their virtual models are created prior to the operation, it is sufficient to map their real positions to the virtual scene. Details such as the precise shape of the objects, or their exact dimensions are not represented in the virtual scene since this information is not necessary for a manipulation task.

In order to calibrate the real and the virtual cantilever positions, the real cantilever is moved to the contact with substrate, then its vertical position is adjusted such as its flexion is null, without snapping the tip off the surface. This position is offset as the (0, 0, 0) on the virtual scene. The object position is obtained from image processing along with the tip position. They are then used to calibrate those initial positions with the reference positions of the virtual scene. During the manipulation and after the calibration, the Virtual Reality engine receives periodically positions of the cantilever and the manipulated object and updates them respectively.

Reconstructing a virtual scene using the position information derived with vision algorithm instead of directly transmitting a video of the SEM images enables to decrease the amount of data to transfer, and thus the communication delays.

In addition to 3D display of the scene, haptic feedback is proposed to provide intuitive manipulation environment. Users control the cantilever through a haptic interface, and feel forces which give them information about the scene (Fig. 5). The desired position of the tip P t is set using the position of the haptic handle P h , scaled down by a factor α p to ensure the precision of the positioning. Operators feel a haptic force F h derived from both force measurements and the vision algorithms.

Force measurements, obtained by the flexion of the AFM probe, are used to compute the vertical haptic force f z h transmitted to users. It is computed from the vertical force f t applied on the tip, scaled by a factor α f to increase it:

f z h = α f f t .
To compute the haptic force in the (x, y) plane (substrate plane), vision algorithms are used. They enable to determine the respective position of the tool (P t ) and the sphere (P s ).

This force is based only on the distance to the closest sphere. Indeed, simplifying the scene enables the user to concentrate on the manipulation task. FIGURE 5. Haptic coupling that connects the master device to the AFM manipulator. Users control the tool by using the joystick, and feel in return force feedback. The haptic feedback is computed using both information coming from vision algorithms and direct force measurement of the piezoresistive cantilever.

the objects in a close neighborhood of the tool could be proposed. A repulsive haptic force is computed. It pushes away the tool from the spheres in the horizontal plane and avoid involuntary contact. This force is based on the distance between the tool and the object: are the coordinates along the x (resp. y) axis. Based on that distance, a repulsive force field is computed. While the distance D is greater than a given threshold R int , the force is null since it is considered that the tip is far enough from the sphere. When this distance becomes less than the threshold, a non null repulsive force f r is computed:

(4.2) f r =        1.0 -D R int f max if D < R int 0 else
where f max is the maximum force that would be sent for D = 0. This value will never be reached since the minimum distance between the tip and the center of the sphere cannot be less then the radius of the object. This force enables users to know when they are entering the interaction radius, and to estimate the relative position of the tip and the sphere for spheres of a given dimension. Moreover, by setting the threshold value equal or slightly superior to the radius of the object, it would also be possible to let the user feel solely the contact in the horizontal plane between the tip and the object.

The force sent to the haptic device is computed by projecting f r along the x and y axes, proportionally to d x and d y . f z h is given by measurements of the tip deflection as explained above:

(4.3) F h =       f x h f y h f z h       =       d x D f r d y D f r α f f t       4.2.
Intuitive Localization of Microspheres using Remote Teleoperation. Experiments involving microspheres are conducted through a teleoperated task between France and Germany (Fig. 6). To assist users both visual and haptic indications are provided.

The first indication that users must get is if the tip is or not in contact with the substrate.

A repulsive haptic feedback in the vertical direction is then sent to the user in case of contact. Figure 7(a) shows the haptic force sent to the user. Before 10 s the tip is higher than the substrate, and the force felt by the user is null. Then the tip reaches the substrate.

The haptic force increases as the user applies a force on the substrate and decreases after 20 s as the user moves the tip upward. The haptic force falls to zero when the tip is moved away from the substrate. The same experiment is performed again after 25 s. The haptic feedback enables users to know when the cantilever is in contact with the substrate, and to estimate the effort applied on the substrate since the haptic force is proportional to the flexion (hence the vertical displacement) of the cantilever.

In addition, users must get information about the relative position of the objects with respect to the tool on the substrate plane. A repulsive haptic feedback as described in Eq .

4.3.

In Fig. 7(b) the user moves the tip around the object, while being in contact with the substrate (a haptic force similar to Fig. 7(a) is also provided but is not depicted here for The developed vision algorithms enable to provide an intuitive remote manipulation environment composed of both visual and haptic feedback. They are experimentally validated in complex micromanipulation scenes involving hundreds of objects.

3D virtual reconstruction of the scene

CONCLUSION

The manipulation of objects smaller than a few micrometers requires the use of Scanning Electron Microscopes. However the visual feedback provided by these systems is not enough to ensure a precise positioning of the tool with respect to the objects, especially when used by non expert operators. In this work an intuitive remote teleoperation environment is provided based on information derived from SEM images. The sphere detection algorithm is robust with respect to SEM image issues, such as shadowing, and detects reliably and efficiently micro spheres in complex micromanipulation scenes. To track the tool, a marker based approach is chosen since fast detection is required to provide stable haptic feedback and real time virtual reconstruction of the scene.

The proposed system enables intuitive manipulation of microscale objects from geographically distant sites. The virtual reconstruction of the scene is displayed in stereoscopic view, and haptic feedback is based on both the vision tracking, detection and force measurements. Virtual guides avoid any involuntary collision between the tool and the objects. To widen the range of applications, nanometer size objects should be considered.

Other tools, such as grippers, should also be integrated. The reduction of the false positives in vision detection can be done by parameter optimization. Further improvements and optimizations can still decrease the needed computation time. Next step will consist in moving the spheres, for example by pushing them with the cantilever. A special care will be taken to define the manipulation strategy. Several issues, such as sticking effects, must be addressed to perform this manipulation [START_REF] Onal | Visual servoing-based autonomous 2-D manipulation of microparticles using a nanoprobe[END_REF].

This work is a first step towards intuitive manipulation environment for microscale objects based on vision algorithms. These systems will offer great opportunities for researchers, giving them access to rare and geographically distant equipments.

FIGURE 1 .

 1 FIGURE 1. Software architecture and communication layout between the remote side in Oldenburg and the operator side in Paris

FIGURE 2 .

 2 FIGURE 2. Example images of melamine spheres imaged with the SEM. Three different scenes have been used for evaluation.

• 2 x

 2 Roundness deviation error (RDE) of the contour in comparison to the calculated radius = ContourPoint.x -BoundingBox.center.x y = ContourPoint.y -BoundingBox.center.y, • Width-to-height-ratio (Ratio W H ) of the bounding box (3.2)Ratio W H = BoundingBox.width BoundingBox.height ,

  (a) Contour detection based algorithm. (b) Cross correlation based algorithm. (c) HoughCircle based approach.

FIGURE 3 .

 3 FIGURE 3. The three used algorithms, contour detection based, cross correlation based, and HoughCircle based.

FIGURE 4 .

 4 FIGURE 4. Example images of the edge detected images with true positives, false positives and false negatives determined by the algorithm marked with color. Green: true positives, blue: false positives, red: false negatives.

  D = d x 2 + d y 2 is the distance between the sphere and the tip, and d x (resp. d y )

FIGURE 6 .

 6 FIGURE 6. Teleoperation of microspheres using a haptic interface, with 3D virtual reconstruction of the manipulation scene. The user manipulates the real AFM, situated in Germany, from the teleoperation facilities in Paris.

  Haptic force along the z axis (N) (a) Vertical localization of the cantilever. Scaling factors are set to: α p = 4800, α f = 1.10 6 .

6 D

 6 Position along the y axis (µm) Position along the x axis (µm) is ta n c e to th e s p h e re c e n te rx a x is (µ m ) D is ta nc e to th e sp he re ce nt er -y ax is (µ m ) In-plane localization of the spheres. Haptic coupling gains for the force computation are set to: R int = 5 µm, f max = 4 N, α f = 1.10 6 and α p = 4800.The magnitude of the force vectors is proportional to the distance between the tip and the object (note that x and y scales are different). The inset represents the simulated magnitude of the haptic force f r .

FIGURE 7 .

 7 FIGURE 7. Haptic feedback provided to assist the 3D localization of microspheres.

  work have been funded by the European project FIBLYS and by the project PROCOPE (DAAD-Egide). Special thanks to Uwe Mick for the technical and scientific support.

  Other scenario considering all

				P t
		f h x f h y		P s
	F h	g	P t	P s
	Haptic interface	f h z	α f	f t
				Tip P t	P s
	P h	α p

AUDE BOLOPION * ,1 , CHRISTIAN DAHMEN * ,2 , CHRISTIAN STOLLE * ,2 , SINAN HALIYO 3 , ST ÉPHANE R ÉGNIER[START_REF] Bolopion | Remote microscale teleoperation through virtual reality and haptic feedback[END_REF] AND SERGEJ FATIKOW[START_REF]Blender[END_REF]

AUDE BOLOPION * ,1 , CHRISTIAN DAHMEN * ,2 , CHRISTIAN STOLLE * ,2 , SINAN HALIYO 3 , ST ÉPHANE R ÉGNIER[START_REF] Bolopion | Remote microscale teleoperation through virtual reality and haptic feedback[END_REF] AND SERGEJ FATIKOW[START_REF]Blender[END_REF]

AUDE BOLOPION * ,1 , CHRISTIAN DAHMEN * ,2 , CHRISTIAN STOLLE * ,2 , SINAN HALIYO 3 , ST ÉPHANE R ÉGNIER 3 AND SERGEJ FATIKOW[START_REF]Blender[END_REF]

AUDE BOLOPION * ,1 , CHRISTIAN DAHMEN * ,2 , CHRISTIAN STOLLE * ,2 , SINAN HALIYO 3 , ST ÉPHANE R ÉGNIER[START_REF] Bolopion | Remote microscale teleoperation through virtual reality and haptic feedback[END_REF] AND SERGEJ FATIKOW[START_REF]Blender[END_REF]


Rama Cont

Arseniy Kukanov

Optimal order placement in limit order markets

Keywords: limit order markets, optimal order execution, execution risk, order routing, fragmented markets, transaction costs, financial engineering, stochastic approximation, Robbins-Monro algorithm

To execute a trade, participants in electronic equity markets may choose to submit limit orders or market orders across various exchanges where a stock is traded. This decision is influenced by the characteristics of the order flow and queue sizes in each limit order book, as well as the structure of transaction fees and rebates across exchanges. We propose a quantitative framework for studying this order placement problem by formulating it as a convex optimization problem. This formulation allows to study how the interplay between the state of order books, the fee structure, order flow properties and preferences of a trader determine the optimal placement decision. In the case of a single exchange, we derive an explicit solution for the optimal split between limit and market orders. For the general problem of order placement across multiple exchanges, we propose a stochastic algorithm for computing the optimal policy and study the sensitivity of the solution to various parameters using a numerical implementation of the algorithm.

Introduction

In todays' automated, electronic financial markets, the trading process is divided into several stages, each taking place on a different time horizon: portfolio allocation decisions are usually made on a monthly or daily basis and translate into trades that are executed over time intervals of several minutes to several days. Existing studies on optimal trade execution [START_REF] Bertsimas | Optimal control of execution costs[END_REF]Lo 1998, Almgren and[START_REF] Almgren | Optimal execution of portfolio transactions[END_REF] have investigated how the execution cost of a large trade may be reduced by splitting it into multiple orders spread in time. Once this order scheduling decision is taken, one still needs to specify how each individual order should be placed: this order placement decision involves the choice of an order type (limit order, market order), order size and destination, when multiple trading venues are available. Orders are filled over short time intervals of a few milliseconds to several minutes and the mechanism through which orders are filled in the limit order book are relevant for such order placement decisions. Market participants need to make such decisions thousands of times each day, and their outcomes have a large impact on each participant's transaction cost as well as on aggregate market dynamics.

Early work on optimal trade execution [START_REF] Bertsimas | Optimal control of execution costs[END_REF]Lo 1998, Almgren and[START_REF] Almgren | Optimal execution of portfolio transactions[END_REF] did not explicitly model the process whereby each order is filled, but more recent formulations have tried to incorporate some elements in this direction. In one stream of literature (see [START_REF] Obizhaeva | Optimal Trading Strategy and Supply/Demand Dynamics[END_REF], [START_REF] Alfonsi | Optimal execution strategies in limit order books with general shape functions[END_REF], [START_REF] Predoiu | Optimal Execution in a General One-Sided Limit-Order Book[END_REF]) a trader is restricted to using market orders whose execution costs are given by an idealized order book shape function. Another approach is to model the process through which an order is filled as a dynamic random process [START_REF] Cont | Statistical modeling of high-frequency financial data[END_REF], [START_REF] Cont | Price dynamics in a markovian limit order market[END_REF] and thus formulate the optimal execution problem as a stochastic control problem: this formulation has been studied in various setting with limit orders [START_REF] Bayraktar | Liquidation in limit order books with controlled intensity[END_REF], [START_REF] Gueant | General Intensity Shapes in Optimal Liquidation[END_REF]) or limit and market orders (Guilbaud andPham 2012, Huitema 2012) but its complexity makes it intractable unless restrictive assumptions are made on price and order book dynamics.

In the present work, we adopt a simpler, more tractable approach: assuming that the trade execution schedule has been specified, we focus on the task of filling each order. Decoupling the scheduling problems from the order placement problem leads to a more tractable approach which is closer to market practice and allows us to incorporate some realistic features which matter for order placement decisions, while conserving analytical tractability.

Individual order placement and order routing decisions play an important role in modern financial markets. Brokers are commonly obliged by law to deliver the best execution quality to their clients and empirical evidence confirms that a large percentage of market orders in the U.S. and Europe is sent to trading venues providing lower execution costs or smaller delays [START_REF] Boehmer | Public disclosure and private decisions: Equity market execution quality and order routing[END_REF]Jennings 2007, Foucault and[START_REF] Foucault | Competition for Order Flow and Smart Order Routing Systems[END_REF]. Market orders gravitate towards exchanges with larger posted quote sizes and low fees, while limit orders are submitted to exchanges with high rebates and lower execution waiting times (see [START_REF] Moallemi | Optimal Order Flow Routing, Exchange Competition, and the Effect of Make/Take Fees[END_REF]). These studies demonstrate how investors' aggregate order routing decisions have a significant influence on market dynamics, but a systematic study of the order routing problem from the investor's perspective is lacking. A reduced-form model for routing an infinitesimal limit order to a single destination is used by [START_REF] Moallemi | Optimal Order Flow Routing, Exchange Competition, and the Effect of Make/Take Fees[END_REF], while [START_REF] Ganchev | Censored exploration and the dark pool problem[END_REF] and [START_REF] Laruelle | Optimal split of orders across liquidity pools: a stochastic algorithm approach[END_REF] propose numerical algorithms to optimize order executions across multiple dark pools, where supply/demand is unobserved. To the best of our knowledge this paper is the first to provide a detailed treatment of investor's order placement decision in a multi-exchange market, unified with the market/limit order choice.

Our key contribution is a quantitative formulation of the order placement problem which takes into account multiple important factors -the size of an order to be executed, lengths of order queues across exchanges, statistical properties of order flows in these exchanges, trader's execution preferences, and the stucture of liquidity rebates across trading venues. Our problem formulation is tractable, intuitive and blends the aforementioned factors into an optimal allocation of limit orders and market orders across available trading venues. Order routing heuristics employed in practice commonly depend on past order fill rates at each exchange and are inherently backward-looking. In contrast, our approach is forward-looking -the optimal order allocation depends on current queue sizes and distributions of future trading volumes across exchanges. When only a single exchange is available for execution, this order placement problem reduces to the problem of choosing an optimal split between market orders and limit orders. We derive an explicit solution for this problem and analyze its sensitivity to the order size, the trader's urgency for filling the order and other factors.

Similar results are also established in a case of two trading venues under some assumptions on order flow distributions. Finally, we propose a stochastic approximation method for solving the order placement problem in the general case and demonstrate its efficiency through examples. Our numerical examples demonstrate that the use of our optimal order placement method allows to substantially decreases trading costs with respect to various 'naive' order placement strategies.

An important aspect of our framework is to account for execution risk, through the incorporation of a penalty for under-or over-filling an order. This penalty is high for time-sensitive executions or when it is costly to catch up on the unfilled portion of the order. Although market orders are executed at a less favorable price, it becomes optimal to use them when execution risk is a primary concern. Optimal limit order sizes are strongly influenced by total quantities of orders queueing for execution at each exchange and by distributions of order outflows from these queues. For example, if at one of the exchanges the queue size is much smaller than the expected future order outflow, it is optimal to place a larger limit order there. Finally, the total order size to be filled plays an important role -limit orders are used predominantly to execute small order sizes and market orders are used for medium and large orders. The amount that can be realistically filled with a limit order at each exchange is naturally constrained by the corresponding queue size and order outflow distribution, so the share of market orders in the optimal allocation increases as the total order size increases. We find that the optimal order allocation almost always splits the total quantity among all available exchanges, suggesting that there is a benefit in having multiple markets.

Section 2 describes our formulation of the order placement problem and shows that it has a global optimum. In Section 3 we derive an optimal split between market and limit orders for a single exchange. Section 4 analyzes the general case of order placement on multiple trading venues. Section 5 presents a numerical algorithm for solving the order placement problem in a general case and our simulation results, and Section 6 concludes. All proofs are presented in the Appendix.

The order placement problem

Consider a trader who has a mandate to buy S shares of a stock within a (short) time interval [0, T ].

The deadline T may be a fixed horizon (e.g. 1 minute) or a stopping time (triggered by market activity). To gain queue priority the trader may immediately submit K limit orders of sizes L k to various exchanges k = 1, . . . , K or submit one market order of size M . The trader's order placement decision is thus summarized by a vector X ∆ = (M, L 1 , . . . , L K ) ∈ R K+1 + whose components are nonnegative i.e. only buy orders are allowed. Our objective is to define a meaningful framework in which the trader may choose the various possibilities for this order placement decision.

We focus on limit order placement and execution and assume that a market order of size M can be filled immediately and with certainty 1 . Limit orders with quantities (L 1 , . . . , L K ) join queues of (Q 1 , . . . , Q K ) pre-existing limit orders at the best bids of K exchanges, where Q k ≥ 0. To simplify the notation, we make an assumption that all K available bid queues are lined up at the best bid price, but it is easily relaxed. Denote by (x) + ∆ = max(x, 0). If L k is constant within [0, T ], the amount purchased with a limit order on exchange k by time

T is equal to (ξ k -Q k ) + -(ξ k -Q k -L k ) + , where ξ k ∆ = C k + D k is an order outflow from the front of k-th bid queue, consisting of C k ∈ [0, Q k ]
cancelations of pre-existing orders from that queue and D k trades with contra-side marketable orders reaching that queue. We specifically note that limit order fill amounts are random, and we allow for partial fills. The total amount A(X, ξ) bought by the trader by time T with all of his orders is a function of the order allocation X and an overall bid queue outflow ξ = (ξ 1 , . . . , ξ K ):

1 This assumption is reasonable if S is small relative to the prevailing market depth. Under the assumption of immediate and certain market order execution it is easy to show that sending market orders to exchanges with high fees is always sub-optimal. We therefore consider a single exchange (with the smallest liquidity fee) for the purpose of sending a single market order.

A(X, ξ)

= M + K k=1 ((ξ k -Q k ) + -(ξ k -Q k -L k ) + ) (1)
The total price of this purchase is divided into a benchmark cost paid regardless of trader's decisions, which may be computed using mid-quote price level, and an execution cost given by

(s + f )M - K k=1 (s + r k )((ξ k -Q k ) + -(ξ k -Q k -L k ) + ), ( 2
)
where s is a half of the bid-ask spread at time 0, f is the lowest available liquidity fee and r k , k = 1, . . . , K are liquidity rebates for all exchanges. The trader can reduce the execution cost by sending more limit orders, but this leads to a risk of underfulfilling the target S because their fills are random. To capture this execution risk we include, in the objective function, a penalty for violations of target quantity in both directions:

λ u (S -A(X, ξ)) + + λ o (A(X, ξ) -S) + , (3)
where λ u , λ o are marginal penalties for, respectively underfulfilling or overfulfilling the execution target S. These penalties are motivated by a correlation that exists between limit order executions and price movements (so-called adverse selection). If A(X, ξ) < S, the trader has to purchase the remaining S -A(X, ξ) shares at time T with market orders. Adverse selection implies that conditionally on the event {A(X, ξ) < S} prices have likely moved up and the transaction cost of market orders at time T is higher than their cost at time 0, i.e. λ u > s + f . Alternatively, if

A(X, ξ) > S the trader experiences buyer's remorse -conditionally on this event prices have likely moved down and he could have achieved a better execution by being more patient. Besides adverse selection, parameters λ u , λ o may reflect trader's execution preferences. For example a trader with a positive forecast of short-term returns may prefer to trade early with a market order and set a larger value for λ u .

Problem 1 (Optimal order placement problem) An optimal order placement is a vector

X * ∈ R K+1 + solution of min X∈R K+1 + E[v(X, ξ)] (4) where v(X, ξ) := (s+f )M - K k=1 (s+r k )((ξ k -Q k ) + -(ξ k -Q k -L k ) + )+λ u (S -A(X, ξ))) + +λ o (A(X, ξ) -S) + (5)
is the sum of the execution cost and penalty for execution risk.

We will denote , ξ)]. We begin by assuming certain economically reasonable restrictions on parameter values.

V (X) = E[v(X

Assumptions

A1 λ u > 0, λ o > 0: the trader is penalized for under-and over-fulfilling the target size

A2 λ o > s + max k {r k } and λ o > -(s + f ):
it is suboptimal to over-fulfill the target size S regardless of fees and rebates A3 min k {r k } + s > 0: possibly negative rebates do not eliminate price improvement from limit order execution Proposition 1 below shows that it is not optimal to submit limit or market orders that are a priori too large or too small (larger than the target size S or whose sum is less than S). Proposition 2 guarantees the existence of an optimal solution.

Proposition 1 Consider C -a compact convex subset of R K+1 + defined by C ∆ = X ∈ R K+1 + 0 ≤ M ≤ S, 0 ≤ L k ≤ S -M, k = 1, . . . , K, M + K k=1 L k ≥ S Under assumptions A1-A3 for any X / ∈ C, ∃ X ∈ C with V ( X ) ≤ V ( X). Moreover, if min k {P(ξ k > Q k + S)} > 0, the inequality is strict: V ( X ) < V ( X).
The penalty function (3) implements a soft constraint for order sizes and effectively focuses the search for an optimal order allocation to the set C. Specific economic or operational considerations could also motivate hard constraints, e.g. M = 0 or K k=1 L k = S. Such constraints can be easily included in our framework but absent the aforementioned considerations we do not impose them here.

Proposition 2 Under assumptions A1-A3, V (X) is a convex function on R K+1 + , it is bounded below and has a global minimizer X ∈ C.

Choice of order type: limit orders vs market orders

To highlight the tradeoff between limit and market order executions in our optimization setup, we first consider a case when the asset is traded on a single exchange, and the trader has to choose an optimal split between limit and market orders. Since K = 1, we suppress a subscript 1 throughout this section.

Proposition 3 (Single exchange: optimal split between limit and market orders)

Assume that ξ has a continuous distribution and (A1-A3) hold. Denote

λ u ∆ = 2s + f + r F (Q + S) -(s + r), and
λ u ∆ = 2s + f + r F (Q) -(s + r). If λ u ≤ λ u , the optimal allocation is (M , L ) = (0, S). If λ u ≥ λ u , the optimal allocation is (M , L ) = (S, 0). If λ u ∈ (λ u , λ u )
, the optimal allocation is:

       M = S -F -1 2s + f + r λ u + s + r + Q, L = F -1 2s + f + r λ u + s + r -Q, (6)
where F (•) is a cumulative distribution function of the bid queue outflow ξ.

In the case of a single exchange, Proposition 1 implies that M + L = S, therefore there is no risk of overfullfilling the target size and λ o does not affect the optimal solution. The trader is only concerned with the risk of underfulfilling the target quantity, and balances this risk with the fee, rebate and other market information. The parameter λ u can be interpreted as trader's urgency to fill the orders, and higher values of λ u lead to smaller limit order sizes, as illustrated on Figure 1.

In contrast, the optimal market order size increases with λ u .

The optimal split between market and limit orders depends on the ratio 2s+f +r λu+(s+r) which balances marginal costs and savings from a market order. It also depends on the distribution F and the queue length Q -keeping all else constant, a trader would submit a larger limit order if its execution is more likely and vice versa. The optimal limit order size decreases with λ u as order underfulfillments become more expensive and increases with f as market orders become more expensive. Another interesting feature is that L is fully determined by Q, F and pricing parameters s, r, f, λ u , while M increases with S. As a consequence of this solution feature, as the order size S increases, a larger fraction M S of that order is executed with a market order. The solution (M , L ) depends on the entire distribution of ξ and not just on its mean, as illustrated on Figure 1 for a pair of exponential and Pareto distributions with equal means.

Optimal routing of limit orders across multiple exchanges

When multiple trading venues are available, dividing the target quantity among them provides better execution quality by reducing the risk of not filling the order. However, sending too many orders leads to an undesireable possibility of overfulfilling the target size. Proposition 4 gives a criterion for optimality of an order allocation X = (M , L 1 , . . . , L K ) that balances these risks.

Proposition 4 Assume (A1-A3), also assume that the distribution of ξ is continuous,

max k {F k (Q k + S)} < 1 and λ u < max k 2s + f + r k F k (Q k ) -(s + r k ) . Then: 1. It is optimal to submit a market order M > 0 if λ u ≥ 2s+f +max k {r k } P k {ξ k ≤Q k } -(s + max k {r k }).
2. It is optimal to submit a limit order

L j > 0 if P k =j {ξ k ≤ Q k } ξ j > Q j > λo-(s+r j ) λu+λo .
3. If 1 and 2 hold for all exchanges j = 1, . . . , K, a necessary and sufficient condition for optimality of an order allocation X ∈ C is that it solves the following equations:

P M + K k=1 ((ξ k -Q k ) + -(ξ k -Q k -L k ) + ) < S = s + f + λ o λ u + λ o (7) P M + K k=1 ((ξ k -Q k ) + -(ξ k -Q k -L k ) + ) < S ξ j > Q j + L j = λ o -(s + r j ) λ u + λ o , j = 1, . . . , K (8)
Equations (7,8) show that an order allocation is optimal as long as it sets the probabilities of underfullfilling the target quantity equal to specific thresholds computed with pricing parameters.

When the number of exchanges K is large, the probabilities in (7,8) are difficult to compute in closed-form. However, before turning to numerical procedures we investigate how these equations can be solved in a tractable case of two exchanges.

Corollary Consider the case of two exchanges with ξ 1 , ξ 2 independent with continuous distributions.

If

1. max k=1,2 {F k (Q k + S)} < 1, 2. λ u < max k=1,2 2s + f + r k F k (Q k ) -(s + r k ) , λ u ≥ 2s + f + max k=1,2 {r k } F 1 (Q 1 )F 2 (Q 2 ) -(s + max k=1,2
{r k }), and

3. F 1 (Q 1 ) < 1 - s + r 2 λ o , F 2 (Q 2 ) < 1 - s + r 1
λ o then there exists an optimal order allocation X = (M , L 1 , L 2 ) ∈ int{C} and it verifies

L 1 = Q 2 + S -M -F -1 2 λ o -(s + r 1 ) λ u + λ o (9a) L 2 = Q 1 + S -M -F -1 1 λ o -(s + r 2 ) λ u + λ o (9b) F1 (Q 1 + L 1 ) F2 (Q 2 + S -M -L 1 ) + Q 1 +L 1 Q 1 +S-M -L 2 F2 (Q 1 + Q 2 + S -M -x 1 )dF 1 (x 1 ) = λ u -(s + f ) λ u + λ o , (9c)
where F 1 (•), F 2 (•) are the cdf of ξ 1 , ξ 2 respectively.

In this solution optimal limit order quantities L 1 , L 2 are linear functions of an optimal market order quantity M . When (9a,9b) are substituted into (9c) we obtain a (non-linear) equation for M , which can be solved for a given distribution of (ξ 1 , ξ 2 ).

Example If ξ 1 , ξ 2 are exponentially distributed with means µ 1 , µ 2 respectively, then an optimal order allocation is given by:

             M = Q 1 + Q 2 + S -z L 1 = z -Q 1 + µ 2 log λ u + s + r 1 λ u + λ o L 2 = z -Q 2 + µ 1 log λ u + s + r 2 λ u + λ o , ( 10
)
where z is a solution of a transcedental equation:

1 + log (λ u + s + r 1 )(λ u + s + r 2 ) (λ u + λ o ) 2 + z µ = λ u -(s + f ) λ u + λ o e z µ , if µ 1 = µ 2 = µ (11)
or

µ 1 µ 1 -µ 2 e -z µ 1 λ u + s + r 1 λ u + λ o µ 1 -µ 2 µ 1 + µ 2 µ 2 -µ 1 e -z µ 2 λ u + s + r 2 λ u + λ o µ 2 -µ 1 µ 2 = λ u -(s + f ) λ u + λ o , if µ 1 = µ 2 (12)
Similarly to the case of one exchange, in this example an optimal market order size M is an increasing linear function of queue sizes Q 1 , Q 2 and the target quantity S, while optimal limit order sizes L i are decreasing functions of the corresponding queue sizes Q i . As in the case of one exchange, the optimal limit order sizes do not depend on the target quantity, but the optimal market order size increases with it. In addition we note that each L i depends on the order flow distribution on both exchanges through µ 1,2 and z.

Numerical solution of the optimization problem

Computing the objective function in the order placement problem (4) or its gradient at any point requires calculating an expectation (a multidimensional integral) which, aside from special examples, is generally not analytically tractable. Stochastic approximation methods, developed specifically for problems where the objective function is an expectation, turn out to be very useful for this problem. We propose in this section a procedure for computing the optimal order placement policy which does not require specifying an order outflow distribution.

We begin by briefly discussing the stochastic approximation approach and the specific method used here. Consider an objective function V (X) , ξ)] to be minimized and denote by g(X, ξ)

∆ = E[f (X
∆ = ∇f (X, ξ)
where the gradient is taken with respect to X. The [START_REF] Robbins | A stochastic approximation method[END_REF] stochastic approximation algorithm tackles the problem in the following way:

1: Choose X 0 and a sequence of 'step sizes' {γ n };

2: for n = 1, . . . , N do 3:

Simulate an independent random variable ξ n with distribution F 4:

Set X n = X n-1 -γ n g(X n-1 , ξ n ) 5: end
This algorithm produces an estimate X ∆ = X N , which converges to the optimal point X as N → ∞ under some technical assumptions (see e.g. [START_REF] Kushner | Stochastic Approximation and Recursive Algorithms and Applications[END_REF]). The sequence of constants {γ n } affects the rate of convergence. This sensitivity can be overcome by using for example the robust stochastic approximation of [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF] which follows the same steps as above with a constant step size γ n = γ and uses an average of iterates X ∆ = 1 N N n=1 X n instead of X N as an estimate of the optimal point. Under some weak assumptions this method achieves a performance bound V ( X ) -V (X ) ≤ DM √ N for a finite N , where D = max

X,X ∈C X -X 2 and M = max X∈C E [ g(X, ξ) 2 ] with • denoting L 2 norm 2 .
Multiplying an optimal step size γ by a constant 2 The method assumes that min X∈X {V (X)} is sought, where V (X) is a well-defined and finite-valued expectation for every X ∈ X and X is a non-empty bounded closed convex set. Moreover V (X) needs to be continuous and convex on X . The optimal step size is γ

= D √ N M .
θ > 0 leads to a performance bound of the same order of magnitude max{θ, θ -1 } DM √ N , i.e. the method is "robust" to step size misspecifications. For our problem we can further bound

D ≤ √ KS and M ≤ (s + f + λ u + λ o ) 2 + K k=1 (s + r k + λ u + λ o ) 2 1/2
We assume that on each iteration X n ∈ int{C} -this is enforced by rescaling X n when needed and does not affect the convergence -then the stochastic gradient in our problem is given by:

g(X n , ξ) =       s + f -λ u 1 {A(Xn,ξ)<S} + λ o 1 {A(Xn,ξ)>S} -(s + r 1 )1 {ξ1>Q1+L1,n} -λ u 1 {A(Xn,ξ)<S,ξ1>Q1+L1,n} + λ o 1 {A(Xn,ξ)>S,ξ1>Q1+L1,n} . . . -(s + r K )1 {ξK>QK+LK,n} -λ u 1 {A(Xn,ξ)<S,ξK>QK+LK,n} + λ o 1 {A(Xn,ξ)>S,ξK>QK+LK,n}      
Note that g(X n , ξ) depends on random variables ξ only through indicator functions, which have simple economic meaning. For example 1 {A(Xn,ξ)<S} = 1 if the target quantity was underfulfilled on the n-th step and 1 {ξk>Qk+Lk,n} = 1 if there was an opportunity to execute a larger limit order at exchange k on that step. This leads to a simple interpretation of numerical iterations -at each step order sizes are increased or decreased depending on whether or not the target quantity was underor overfullfilled and whether or not a larger limit order could be filled at a given exchange. If a model for ξ is available, one can use it to simulate ξ and compute a numerical solution that takes into account specific order flow assumptions (e.g. forecasts of future trading volume). Alternatively, one can use past order fill data to compute indicator functions in g(X n , ξ) and obtain a non-parametric numerical solution for the order placement problem (using the empirical distribution of past order fills instead of assuming a functional form for F ).

We analyze the numerical stability and convergence of estimates X by comparing them with an analytical solution in the case of one exchange. For this computation we use Q = 2000 shares, ξ ∼ P ois(µT ), µ = 2200 shares per minute, T = 1 minute and S = 1000 shares. The pricing parameters (in dollars per share) are s = 0.02, r = 0.002, f = 0.003 and fall in a typical value range for US equities. Finally, the penalty costs (also in dollars per share) are set to λ o = 0.024, λ u = 0.026.

According to (6) the optimal allocation (M , L ) = (728, 272) shares. Numerical estimates X were computed for five initial points X 0 and different number of iterates N in the stochastic approximation, using a step size γ

∆ = √ KS N (s + f + λ u + λ o ) 2 + N K k=1 (s + r k + λ u + λ o ) 2 -1/2
.

For each choice of X 0 and N we simulated an additional L = 1000 observations of ξ to estimate the objective values V (X) with sample averages

W (X) = 1 L L i=1 v(X, ξ i ))
. Figures 2 and3 show that estimates converge to X regardless of X 0 . When X 0 = X , estimates remain close to that point. Convergence is also quite fast -after as few as 50 samples the algorithm can be within 2% of the optimal objective value. In the worst case of initial points on the boundary it can take a few thousand samples to converge. It is also worth noting that convergence in terms of the objective value occurs significantly faster than convergence in terms of the order allocation vector.

Figure 2

Convergence of objective values to an optimal point for different inital points.

Figure 3

Convergence of order allocation vectors to an optimal point for different inital points. We also estimated savings from optimal order placement and from dividing a limit order among multiple exchanges. Denote a pure market order allocation by X M = (S, 0, . . . , 0), a single limit order allocation by X L = (0, S, 0, . . . , 0) and an equal split allocation by

K M L 1 L 2 L 3 L 4 L 5 /S W (XM ) W (XL) W (XE) W ( X ) S =
X E = ( S K + 1 , S K + 1 , . . . , S K + 1
).

Table 1 presents produced by stochastic approximation clearly outperform the naive benchmarks, especially when a target quantity S is relatively small and cost savings of limit orders can be fully captured.

Comparing W (X L ) and W (X E ) we also see that splitting a limit order across multiple exchanges can be very advantageous when limit order fills are independent. Since multiple exchanges in this example are copies of each other, the algorithm splits the total limit order amount equally among them. This is not the same as the allocation X E because the latter sets a market order size to S K + 1 , which may be too large or too small depending on S and the number of available exchanges. Another interesting feature of the numerical solution X is a tendency to oversize the total quantity of limit orders, which is clearly observed for S = 1000, 5000 and K = 4, 5. This may be a consequence of assumed independence between ξ k -by submitting large orders to multiple exchanges the algorithm reduces the probability of underfullfilling the target quanitity with a relatively low probability of overfulfilling it.

To illustrate the structure of a numerical solution we performed a sensitivity analysis with K = 2 exchanges and parameters

Q 1 = Q 2 = 2000, S = 1000, ξ 1,2 ∼ P ois(µ 1,2 T ), µ 1 = 2600, µ 2 = 2200, T =
1, s = 0.02, r 1 = r 2 = 0.002, f = 0.003, λ u = 0.26 and λ o = 0.24. Varying some of these parameters one at a time we plot the numerical solution X after N = 1000 iterations, together with an analytical solution for a single exchange. The results are presented on Figures 4 and5. Similarly to the single exchange case, limit order sizes on two exchanges L 1 , L 2 decrease and market order size M increases as the penalty λ u increases. Increasing the half-spread s, the rebate r 1 or the fee f makes a limit order on exchange number one more attractive, so L 1 increases and M decreases.

Because the penalty λ u is large in this example, execution risk is more important than fees and rebates, therefore the queue size Q 1 and the order outflow mean µ 1 have a much stronger effect on the optimal solution than r 1 . Both decreasing the Q 1 and increasing µ 1 make a limit order fill more likely at exchange number one and L 1 increases3 . Finally, as in the case of a single exchange, the target size S has a strong effect on the optimal order allocation. Only limit orders are used while S is small, but as it becomes larger it is difficult to fill that amount solely with limit orders and the optimal market order size begins to grow to limit the execution risk.

Conclusion

We formulated and solved the problem of placing a small batch of orders on multiple trading venues. In the case when there is a single exchange we derived an optimal split between a limit and a market order sizes. For more general cases, we proposed and tested a stochastic approximation algorithm that is shown to quickly converge to an optimal point. We explored the properties of an optimal order allocation policy and showed that splitting an order across multiple exchanges can lead to a substantial reduction in transaction costs.

Appendix. Proofs

Note: to avoid extra notation in the following proofs we refer to the function v(X, ξ) by v(X).

Proposition 1 Consider C -a compact convex subset of R K+1 + defined by

C ∆ = X ∈ R K+1 + 0 ≤ M ≤ S, 0 ≤ L k ≤ S -M, k = 1, . . . , K, M + K k=1 L k ≥ S Under assumptions A1-A3 for any X / ∈ C, ∃ X ∈ C with V ( X ) ≤ V ( X). Moreover, if min k {P(ξ k > Q k + S)} > 0, the inequality is strict: V ( X ) < V ( X).
Proof: First, for any allocation X that has M > S, we automatically have A( X) > S and we can show that the (random) cost and penalty of X is larger than those of X naive ∆ = (S, 0, . . . , 0) ∈ C:

v( X, ξ) -v(X naive , ξ) = (s + f )( M -S) - K k=1 (s + r k )((ξ k -Q k ) + -(ξ k -Q k -L k ) + )+ λ o M -S + K k=1 ((ξ k -Q k ) + -(ξ k -Q k -L k ) + ) = (λ o + s + f )( M -S) + K k=1 (λ o -s -r k )((ξ k -Q k ) + -(ξ k -Q k -L k ) + ) > 0,
which holds for all random ξ. Therefore, V ( X) > V (X naive ). Similarly, for any allocation X with

Lk > S -M define a different allocation X by M = M , L j = Lj , ∀j = k and L k = S -M . Then v( X, ξ) -v( X , ξ) = 0 on the event B = {ω|ξ k (ω) < Q k + S -M }. On its complementary event B c , v( X, ξ) -v( X , ξ) = -(s + r k )((ξ k -Q k -S + M ) + -(ξ k -Q k -Lk ) + ) +λ o ((ξ k -Q k -S + M ) + -(ξ k -Q k -Lk ) + ). Therefore V ( X) -V ( X ) = E v( X, ξ) -v( X , ξ)|B P(B) + E v( X, ξ) -v( X , ξ)|B c P(B c ) = 0 + E (λ o -(s + r k ))((ξ k -Q k -S + M ) + -(ξ k -Q k -Lk ) + )|B c P(B c ) ≥ 0 with a strict inequality if P(B c ) > 0. If X /
∈ C, we can continue truncating limit order sizes L j > S -M following the same argument. Each time the truncation does not increase the objective function and finally we obtain X ∈ C, such that

V ( X ) ≤ V ( X). Next, if X is such that M - K k=1 Lk < S define s = S -M - K k=1 Lk , take M = M , L k = Lk , k = 1, . . . , K -1 and L K = Lk + s. Then, on the event B = ω|ξ K (ω) < Q K + LK we have v( X, ξ) = v( X , ξ). However, on the event B c , v( X, ξ)-v( X , ξ) = (s+r K )((ξ K -Q K -LK ) + -(ξ k -Q k -LK -s) + )+λ u ((ξ K -Q K -LK ) + -(ξ k -Q k -LK -s) + ), therefore V ( X) -V ( X ) = E v( X, ξ) -v( X , ξ)|B P(B) + E v( X, ξ) -v( X , ξ)|B c P(B c ) = 0 + E (λ u + (s + r k ))((ξ K -Q K -LK ) + -(ξ k -Q k -LK -s) + )|B c P(B c ) ≥ 0 with a strict inequality if P(B c ) > 0.
Proposition 2 Under assumptions A1-A3, V (X) is a convex function on R K+1 + , bounded from below and admits a global minimizer X ∈ C.

Proof: First, note that (ξ k -Q k ) + -(ξ k -Q k -L k ) + are concave functions of L k . Therefore, A(X, ξ)
is concave as a sum of concave functions. Similarly, the cost term in v(X, ξ) is a sum of convex functions, as long as r k ≥ -s, k = 1, . . . , K and is itself a convex function. Second, since S -A(X, ξ) is a convex functon of X, and the function h(x) ∆ = λ u (x) + -λ o (-x) + is convex in x for positive λ u , λ o , so the penalty term h (S -A(X, ξ))) is also convex.

If λ o > s + max k {r k } the function V (X) is also bounded from below since v(X, ξ) ≥ -(s + max k {r k })S.
Finally, since V (X) is convex, it is also continous and reaches a local minimum V min on the compact set C at some point X ∈ C. By convexity, V min is a global minimum of V (X) on C.

Moreover, since λ o > s + max k {r k }, Proposition 1 guarantees that V min < V ( X) for any X / ∈ C, so V min is also a global minimum of V (X) on R K+1 + .
Proposition 3 Consider the case of a single exchange where ξ has a continuous distribution.

Denote by λ u

∆ = 2s + f + r F (Q + S) -(s + r) and λ u ∆ = 2s + f + r F (Q) -(s + r).
Then, under Assumptions (A1-A3):

• If λ u ≤ λ u , the optimal allocation is (M , L ) = (0, S). If λ u ≥ λ u , the optimal allocation is (M , L ) = (S, 0).

• If λ u ∈ (λ u , λ u ), the optimal allocation is given by ( 6).

Proof: By Proposition 1 there exists an optimal split (M , L ) ∈ C between limit and market orders. Moreover for K = 1 the set C reduces to a line M + L = S so it is sufficient to find M .

Restricting L = S -M implies that {A(X, ξ) > S} = ∅, {A(X, ξ) < S, ξ > Q + L} = ∅, and we can rewrite the objective function as

V (M ) = E (s + f )M -(s + r)((ξ -Q) + -(ξ -Q -S + M ) + ) + λ u (S -M -((ξ -Q) + -(ξ -Q -S + M ) + )))) + . ( 13
)
For M ∈ (0, S) the expression under the expectation in ( 13) is bounded for all ξ and differentiable with respect to M for almost all ξ, so we can compute V (M ) = dV (M ) dM by interchanging the order of differentiation and integration (see e.g. [START_REF] Aliprantis | Principles of Real Analysis[END_REF], Theorem 24.5):

V (M ) = E s + f + (s + r)1 {ξ>Q+S-M } -λ u 1 {ξ<Q+S-M } = 2s + f + r -(s + r + λ u )F (Q + S -M ) (14) Note that if λ u ≤ 2s + f + r F (Q + S) -(s + r), then V (M ) ≥ 0 for M ∈ (0, S) and therefore V is non-decreasing at these points. Checking that V (S) -V (0) ≥ (s + f -λ u )S + (λ u + s + r)S(1 - F (Q + S)) ≥ 0 we conclude that M = 0. Similarly, if λ u ≥ 2s + f + r F (Q) -(s + r), then v(M ) ≤ 0 for all M ∈ (0, S) and V (M ) is non-increasing at these points. Checking that V (S) -V (0) ≤ (s + f -λ u )S + (λ u + s + r)S(1 -F (Q)) ≤ 0 we conclude that M = S. Finally, if λ u is between
these two values, ∃ > 0, such that V ( ) < 0, V (S -) > 0 and by continuity of V there is a point where V (M ) = 0. This M is optimal by convexity of V (M ) and ( 6) solves equations v(M , ξ) = 0, L = S -M .

Proposition 4 Assume (A1-A3), also assume that the distribution of ξ is continuous,

max k {F k (Q k + S)} < 1 and λ u < max k 2s + f + r k F k (Q k ) -(s + r k ) . Then: 1. It is optimal to submit a market order M > 0 if λ u ≥ 2s+f +max k {r k } P k {ξ k ≤Q k } -(s + max k {r k }).
2. It is optimal to submit a limit order

L j > 0 if P k =j {ξ k ≤ Q k } ξ j > Q j > λo-(s+r j ) λu+λo .
3. If 1 and 2 hold for all exchanges j = 1, . . . , K, a necessary and sufficient condition for optimality of an order allocation X ∈ C is that it solves equations (7)-( 8).

Proof: Proposition 2 implies the existence of an optimal order allocation X ∈ C. First, we define X M ∆ = (S, 0, . . . , 0) and prove that X = X M by contradiction. If X M were optimal in the problem (4) it would also be optimal in the same problem with a constraint L k = 0, k = j, for any one j. In other words, the solution (S, 0) would be optimal for any one-exchange problem, defined by using only exchange j. But by our assumption, there existsJ such that λ

u < 2s + f + r J F J (Q J ) -(s + r J ) and
Proposition 3 implies that (S, 0) is not optimal for the J-th single-exchange subproblem, leading to a contradiction.

The function v(X, ξ) is bounded for X ∈ C and for all ξ, differentiable with respect to M and L k , k = 1, . . . , K for X ∈ C\ {X M } for almost all ξ. Applying the same theorem as in the proof of Proposition 3 we conclude that V (X) is differentiable for X ∈ C\ {X M } and we can compute all of its partial derivatives by interchanging the order of differentiation and integration. The KKT conditions for problem (4) and X ∈ C\ {X M } are

s + f -λ u P(A(X , ξ) < S) + λ o P(A(X , ξ) > S) -µ 0 = 0 (15) -(s + r k )P(ξ k > Q k + L k ) -λ u P(A(X , ξ) < S, ξ k > Q k + L k )+ λ o P(A(X , ξ) > S, ξ k > Q k + L k ) -µ k = 0, k = 1, . . . , K (16)
M ≥ 0, L k ≥ 0, µ 0 ≥ 0, µ k ≥ 0, µ 0 M = 0, µ k L k = 0, k = 1, . . . , K (17)
Since the objective function V (.) is convex, conditions ( 15)-( 17) are both necessary and sufficient for optimality. The first result of this proposition follows from considering any X with M = 0:

V ( X) ≥ λ u SP k {ξ k ≤ Q k } -(s + max k {r k })SP k {ξ k ≤ Q k } ≥ (s + f )S = V (X M
) and we already argued that ∃X with V (X ) < V (X M ), so X = X and therefore M > 0. Rearranging terms in a j-th equality ( 16) we obtain P(ξ j > Q j + L j ) λ o -(s + r j ) -(λ u + λ o )P(A(X , ξ) < S|ξ j > Q j + L j ) -µ j = 0 (18)

The term in square brackets in (18) is negative for any X ∈ C\ {X M } with L j = 0, because P(A(X, ξ) < S|ξ j > Q j + L j ) > P

k =j {ξ k ≤ Q k } ξ j > Q j >
λo-(s+r j ) λu+λo

by assumption and since µ j ≥ 0 the condition ( 16) cannot be satisfied with L j = 0. We showed that M > 0, L j > 0 for all j = 1, . . . , K and therefore, µ 0 = µ 1 = • • • = µ K = 0 by complimentary slackness. Then the KKT conditions (15)-( 17) reduce to ( 7)-( 8).

Corollary Assume that two exchanges are available for execution, ξ 1 is independent of ξ 2 and the distribution of (ξ 1 , ξ 2 ) is continuous. Also assume that:

1. max k=1,2 {F k (Q k + S)} < 1 2. λ u < max k=1,2 2s + f + r k F k (Q k ) -(s + r k ) , λ u ≥ 2s + f + max k=1,2 {r k } F 1 (Q 1 )F 2 (Q 2 ) -(s + max k=1,2 {r k }) 3. F 1 (Q 1 ) < 1 - s + r 2 λ o , F 2 (Q 2 ) < 1 - s + r 1
λ o Then an optimal order allocation X = (M , L 1 , L 2 ) ∈ int{C} and it solves the equations (9a-9c). L k = S are ruled out by directly checking (16). Finally, L 1 = 0 and L 2 = 0 are also ruled out by (16). For example if L 1 = 0, then by Proposition 1 M + L 2 = S and in ( 16) µ 2 = 0 by complimentary slackness, P(A(X , ξ) < S, ξ 2 > Q 2 + L 2 ) = P(A(X , ξ) > S, ξ 2 > Q 2 + L 2 ) = 0. But then ( 16) cannot hold because P(ξ 2 > Q 2 + L 2 ) > 0.

For any X ∈ int{C}, A(X, ξ) > S if and only if all the following three inequalities are satisfied:

ξ 1 > Q 1 + S -M -L 2 (19a) ξ 2 > Q 2 + S -M -L 1 (19b) ξ 1 + ξ 2 > Q 1 + Q 2 + S -M (19c)
These inequalities give a simple characterization of the event {A(X, ξ) > S} and their equivalence is directly verified by considering subsets of (ξ 1 , ξ 2 ) forming a complete partition of R 2 + . Case 1: ξ 1 > Q 1 + L 1 , ξ 2 > Q 2 + L 2 . Since L 1 + L 2 + M > S, we have A(X, ξ) = L 1 + L 2 + M > S and at the same time all of the inequalities (19a-19c) are satisfied, so they are trivially equivalent in this case.

Case 2:

ξ 1 > Q 1 + L 1 , Q 2 ≤ ξ 2 ≤ Q 2 + L 2 . Because of the condition ξ 1 > Q 1 + L 1 , (19a) is satisfied.
We have in this case that A(X, ξ) = L 1 + ξ 2 -Q 2 + M and thus A(X, ξ) > S if and only if (19b) is satisfied. Finally, ξ 1 > Q 1 + L 1 together with (19b) imply (19c), so A(X, ξ) > S and (19a-19c) are equivalent in this case.

Case 3:

ξ 2 > Q 2 + L 2 , Q 1 ≤ ξ 1 ≤ Q 1 + L 1 .
Similarly to Case 2 we can show that inequalities (19a-19c) are satisfied if and only if A(X, ξ) > S.

Case 4:

Q 1 + S -M -L 2 < ξ 1 ≤ Q 1 + L 1 , Q 2 + S -M -L 1 < ξ 2 ≤ Q 2 + L 2 .
This set is nonempty because 0 < S -M -L 1 < L 2 and similarly for L 1 , L 2 reversed. only (19c) needs to be checked. We can write A(X, ξ) = ξ 1 -Q 1 + ξ 2 -Q 2 + M > S if and only if (19c) holds, so A(X, ξ) > S is equivalent to (19a-19c).

Figure 1

 1 Figure 1 Optimal limit order size L for one exchange. The parameters for this figure are: Q = 2000, S = 1000, s = 0.02, r = 0.002, f = 0.003. Colors correspond to different order outflow distributions -exponential with means 2200 and 2500 and Pareto with mean 2200 and a tail index 5.

  outputs from the numerical algorithm with X 0 = X E , N = 1000, L = 1000 for different order sizes S and number of exchanges K = 1, . . . , 5. The parameters s, f, r, λ u , λ o are same as in the previous simulation and exchanges are identical replicas of each other: r k = r, Q k = Q and ξ n,k ∼ P ois(µT ) are i.i.d. copies of each other, where k = 1, . . . , K, n = 1, . . . , N . Order allocations

Figure 4

 4 Figure 4Sensitivity analysis for a numerical solution X = (M, L1, L2) with two exchanges and an optimal solution (M a , L a ) with the first exchange only.

Figure 5

 5 Figure 5Sensitivity analysis for a numerical solution X = (M, L1, L2) with two exchanges and an optimal solution (M a , L a ) with the first exchange only.

Proof:

  Solutions on the boundary of C are sub-optimal: M = 0 and M = S are ruled out by assumption 2, L 1 = S -M and L 2 = S -M are ruled out by assumption 3 and (16). Solutions with M + K k=1

The observed drop in L1 for large µ1 and small Q1 is a feature of this example, we were not able to replicate it for other distributions of ξ.

Case 5: Outside of Cases 1-4, either (19a) or (19b)

completely symmetric, and it shows that neither A(X, ξ) > S nor (19a-19c) hold in this case.

Next, we use inequalities (19a-19c) to characterize the set {A(X, ξ) > S} in the first-order conditions ( 7)-( 8). We observe that in the two-exchange case

and then use the independence of ξ 1 and ξ 2 to compute

Together with ( 7) and ( 8), this leads to a pair of equations for limit orders sizes:

whose solution is given by L 1 , L 2 from (9a,9b). To obtain the equation ( 9c), we rewrite the first equation in (7,8) using the inequalities (19a-19c). Then P (A(X, ξ) > S) may be computed as the integral of the product measure F 1 ⊗ F 2 over the region defined by

This integral is given by