Pierre Robin
email: pierre.robin@ec-nantes.fr.

Alban Leroyer

Jeroen Wackers

Tackling Modern Sailing Challenges with a CFD-based Dynamic VPP

Keywords: CFD, VPP, sailing, hydrofoil, optimisation, FSI AGR Adaptive Grid Refinement AWA Apparent Wind Angle AWS Apparent Wind Speed CFD Computational Fluid Dynamics CPU Central Processing Unit DoF

A dynamic Velocity Prediction Program (VPP) integrated in a Computational Fluid Dynamics (CFD) code is described. The program computes aerodynamic forces and directly inputs them to the CFD solver, which solves both the flow and the motions of the boat. The VPP code controls the heel of the boat and optimises a sail power parameter for a given True Wind Angle and True Wind Speed. This approach allows naval architects to swiftly and precisely compare several yacht designs in real sailing configurations using only a few CFD computations. Several advanced features of this program are covered in this paper including convergence criteria, automatic grid refinement, foil fluidstructure interaction, multiple aerodynamics models and rudder control.

INTRODUCTION

Velocity Prediction Programs (VPP) have been used in the design of modern sailing ships for a few decades. For naval architects, these tools allow to assess and compare the performance of several yacht designs and eventually identify the best one, and as computer technology improves so do VPP codes. Classic static VPPs compute the equilibrium of the forces applied to the yacht and most rely on an interpolation performed in a matrix of the hydrodynamic forces. This matrix, previously obtained with towing tank trials, is now computed using CFD simulations in 80 to 200 configurations (with different combinations of pitch, trim, draft, heel, speed, etc.). Even with the relentless progresses made in this field, such calculations are very expensive in CPU-time. Moreover, with such a process, any change in the hull design means that a new matrix should be computed, since the hydrodynamic forces applied on the yacht are likely to be different. Thus, traditional VPPs are somewhat impractical for hull-shape optimisation. Additionally, the CFD is performed on an isolated hull, so the interactions with the appendages are approximate at best, if not neglected entirely. Furthermore, although ship performance in waves will be a major challenge in coming years for offshore sailing, predicting motion in a sea state with a VPP remains difficult: it can only be approximated using successive steady-state force inputs.

For these reasons, we developed a dynamic VPP in which the aerodynamic forces are directly applied to the boat during the CFD computation, while the flow solver handles the hydrodynamic forces. This approach results in dynamic changes in all degrees of freedom (DoF), such as trim, sinkage, heel, etc, leading to a coupling between aerodynamic and hydrodynamic effects. [START_REF] Roux | Strongly coupled VPP and CFD RANSE code for sailing yacht performance prediction[END_REF] already performed 5-DoF simulations using a coupling between a CFD code solving hydrodynamics and a potential-flow code computing aerodynamic forces for an AC90 yacht. Later, [START_REF] Böhm | A Velocity Prediction Procedure for Sailing Yachts with a hydrodynamic Model based on integrated fully coupled RANSE-Free-Surface Simulations[END_REF], Lindstrand [START_REF] Lindstrand Levin | Sailing yacht performance prediction based on coupled CFD and rigid body dynamics in 6 degrees of freedom[END_REF] and [START_REF] Persson | An Improved Procedure for Strongly Coupled Prediction of Sailing Yacht Performance[END_REF] presented VPPs using CFD solvers to compute hydrodynamic forces and body motions, while aerodynamics rely on an analytical model. The main innovation of our work is the use of modern techniques that have proven to be really useful when simulating a sailing yacht, such as overset grids or free-surface capture. Indeed, since our VPP program is embedded in the commercial CFD program FINE/Marine, several key advanced features already in the software can be used, such as quasi-static body motions, adaptive grid refinement or modal approach for fluid-structure interactions. In addition to this, other features more tailored to our use have been developed, including a PID-controlled rudder, a sockmesh approach for foils and VPP-tailored convergence control. By nature, the main objective of this high-fidelity tool is not to provide a complete polar for the studied design, but to compare different designs in key wind conditions all along the design process. It can also be used dynamically, to simulate a yacht in waves for instance, but this is still under development. While this paper will present some of the features developped for ship in waves simulations, it will mainly focus on the search of steady states to optimise the yacht's speed.

The development of this new tool reflects Finot-Conq's strategy to incorporate more high-fidelity simulations in their design process, while limiting the associated CPU costs. The VPP covered in this paper has been used for Finot-Conq's two new IMOCAs (Paprec-Arkéa and For People), launched in early 2023 with the goal of winning the Vendée Globe in 2024-2025.

Section 2 provides details of the VPP algorithm, while section 3 covers the aerodynamic modules. The hydrodynamic solver is discussed in section 4. Finally, section 5 presents results obtained using the VPP.

VPP ALGORITHM

Our dynamic CFD-based VPP (Cadence Design Systems and Finot-Conq, 2021), is integrated in the commercial CFD software FINE/Marine. It allows naval architects to optimise boat speed as a function of the sails trim (represented by a sail power parameter, see Section 3) in a single 5 or 6 DoF computation.

The general process of the VPP is described in Fig. 1. The VPP is a module added to FINE/Marine, thus its use requires first to set up a classical FINE/Marine computation. Placing the yacht in a configuration close to its final position (in terms of cardan angles and displacement) will help the simulation to converge and reduce CPU-time. Additionally, the VPP module requires the user to provide a True Wind Angle (TWA), a True Wind Speed (TWS) and some of the yacht's main characteristics (length over all, freeboard height, air draft, etc.).

Figure 1. VPP code general process

Once this data is read and FINE/Marine has started, the VPP code updates sail forces at each time step (covered in Section 3), adjusts the rudder angle (Section 2.1), and keeps the heel of the yacht within user-specified bounds (Section 2.3). The resulting forces and moments are then input to the flow solver ISIS-CFD, which solves the flow and the motions of the yacht. During the simulation, the VPP checks if an equilibrium state is reached using a specific set of convergence criteria (Section 2.2). If these criteria are satisfied, a new value for the sail power parameter is chosen by the optimisation algorithm and used in a new stage of the same computation. At the end of the computation, several equilibrium states have been found (and thus several (speed,sail power) doublets), and the VPP outputs an optimised yacht's speed and sail power parameter for the yacht's geometry and the doublet (TWA,TWS) provided.

Rudder Models

For course control, two rudder models have been developed. The first consists in an analytical rudder balancing the M z moment induced by the rest of the forces applied on the boat. When using this method, the user is required to provide a boat geometry without rudder, as well as rudder polars. The rudder's angle is computed from these polars to provide the desired M z, then all the rudder forces are evaluated at that angle and output to the flow solver.

However, to improve the accuracy of the boat's behavior, and in the scope of carrying out future shipin-waves simulations, a second approach using a CFD-modelled rudder with a PID controller has been developed. For this approach, the rudder's geometry needs to be included and is placed in an overset domain (see section 4.1). The rudder angle δ(t) is controlled using, as error e, the difference between the yacht's course angle and a given target bearing. The rudder angle is updated with a given period T using the following law:

δ(t) = K p * e(t) + K i * t t-T e(t) dt + K d * de(t) dt (1)
In addition to the PID-controller, a simplified model of an autopilot's actuator has been implemented with realistic angular speed limitations. This method is still a work in progress, as only linear integration and derivation schemes are implemented and the gains K p , K i , K d are, for the moment, determined by trial-and-error only and will require a lot of fine-tuning using dedicated methodologies.

Convergence Criteria

To test if the yacht has reached an equilibrium state, specific convergence criteria have been implemented. Indeed, yachts equipped with foils can induce small-scale unsteadiness in the computation, thus requiring to make sure that key physical quantites undergo limited changes around a converged value. To do so, the user has to specify:

• Which physical values (forces, moments, motions quantities) will be studied for convergence,

• After how much simulation time the convergence check is first carried out (t start ),

• The amount of time between two convergence tests (t laps ),

• The duration (in seconds) of data collected to form a convergence sample t cs ,

• For every physical value: two thresholds to compare the data collected to:

α ref and β ref ,
Indeed, for every physical value studied, two additional indicators will be computed: the average of the absolute value of the deviation from the mean value α, whose expression for a quantity x is:

α = 1 t cs tcs |x -x| (2)
where x is the mean value of the quantity x. But also the maximal absolute value of the deviations from the mean value β, whose expression for a quantity x is:

β = max tcs |x -x| (3)
If for every physical value studied, α and β are inferior to their respective thresholds αref and βref , an equilibrium state is considered to be reached. The general workflow of the method is described Fig2:

Sail Power Parameter Optimisation

To find the optimum yacht's speed at a given operating point, the sail power parameter, that represents the crew trimming the sails, is optimised. The way this parameter affects the sail forces will be covered in Section3.

In a single VPP computation, several equilibrium states are reached, each one being associated with an intermediate calculation using a different value of the sail power parameter. A minimum of 3 intermediates calculations is required to calculate an optimum for the sail power parameter. More can be performed if desired, to refine the optimum values obtained. During each intermediate calculation, the sail power parameter may be adjusted if needed to keep the yacht's heel angle between the boundaries specified by user.

The first intermediate computation is initiated with the parameter's maximum value. Once it is considered converged, the sail power parameter is set to its minimum value. Once the second intermediate computation has converged, a third one starts with the mean value of the two previously used. Then the VPP code finds an optimal value for the sail power parameter using parabolic optimisation, based on the first three iterations. If more points are required by the user, the sail power parameter value is set at the computed optimum and another computation starts. At the end of it, a new parabolic optimisation is performed and a new optimum for the sail power parameter is found.

Previous work at Finot-Conq attempted to perform sail power parameter optimisation using other methods, but parabolic optimisation was chosen as the most CPU-efficient. It has also proven to be sufficient as the relationship between speed and sail power parameter is correctly captured by a parabola. For the same reason, multi-criterion optimisation using reefing along with sail power was discarded. Indeed, it is more efficient and more realistic to perform separate VPP runs for different sail plans, due to the discontinuous nature of reefing on a real yacht. An example of results for a 40-feet fast cruiser is presented in Tab. 1 and Fig. 3. The maximal speed of the yacht is not reached for the maximum value of sail power parameter, as a reduced heeling allows the boat to sail faster. The 4th optimisation iteration eventually presented a limited interest here, since the difference in speed obtained at the end of the 4th iteration (0.8488) is really close to the optimum computed at the end of the 3rd iteration (0.8471).

AERODYNAMIC LIBRARY

An important part of the VPP code lies in the computation of aerodynamic forces. Our VPP can either use analytical models to determine sail aerodynamics or interpolate from a data matrix. The forces and moments applied by the sails to the boat are considered to be a function of AW A, AW S and of the sails' trim, characterized by the sail power parameter. The aerodynamic forces are computed at each time step to account for any change in these quantities.

Analytical Model

Aerodynamic forces and moments can be computed using analytical methods inspired by the Offshore Racing Congress (ORC) VPP methodology (ORC, 2021). To do so, formulas based on several yacht's sails dimensions are provided to compute sail areas, lift and drag coefficients and center of effort height for each sail. Compared to the ORC formulations, the main sail and the jib are divided in smaller areas in order to take into account the wind gradient. When using this method, the sail power parameter is a coefficient applied linearly to the aerodynamic lift force: if the sail power coefficient is equal to 0.95, this means the crew trimmed the sails to get 95% of the maximum lift.

Interpolation from Aerodynamics Matrices

In addition to this analytical model, the VPP can also use aerodynamics matrices obtained by external CFD computations, typically those made by the sailmakers. Currently, an interpolation module using cubic interpolations supplies the VPP code with the aerodynamic coefficients for lift and drag. Results for the surge coefficient obtained using this interpolation method for a mainsail and jib aerodynamics matrix are presented in Fig. 4.

When using aerodynamic matrices provided by sails designer, the sail power parameter is also the parameter that captures the trim of the sails (sheet and traveller settings) and is one of the variable of the interpolation alongside the apparent wind angle and the heel of the boat. The aerodynamic forces are then computed using the following formula :

F x = 1 2 ρ S C x V 2 (4)
F y = 1 2 ρ S C y V 2 (5)
with F x being the aerodynamic force along x-axis, F y the one along the y-axis, ρ the air' density, S the surface of the considered sail, C x the surge coefficient, C y the sway coefficient (both coefficients' values being interpolated from the sailmakers' matrix) and V the velocity of the flow at the sail's effort height (obtained from AWS).

However, it could prove useful to be able to interpolate using more variables, therefore an interpolation method based on radial-basis functions is being developed. For such purposes, several radial-basis functions are available [START_REF] Volpi | Development and validation of a dynamic metamodel based on stochastic radial basis functions and uncertainty quantification[END_REF], but based on our early results, the radial function ϕ(r) = r 3 appears to be an efficient choice to interpolate inside a sail plan's aerodynamic efforts matrix, and it is easy to implement.

Windage

Windage forces acting on the different parts of the sailing boat (hull, mast, rigging, crew, etc.) are also computed using analytical formulas, inspired by the ORC VPP methodology. The windage forces consider both front and side surface areas of the different elements of the boat. These forces are then added to the sail forces and the sum is input to the hydrodynamic solver as external forces applied on the boat. It is noted that, in some sailmakers' aerodynamic matrices, these forces are already taken into account, thus no additional windage forces need to be computed. This solver features a mixture-fluid formulation to model the water-air interface [START_REF] Queutey | An Interface Capturing Method for Free-Surface Hydrodynamic Flows[END_REF][START_REF] Wackers | Free-surface viscous flow solution methods for ship hydrodynamics[END_REF]. The flow equations are discretised in a finite-volume framework, using pressure-velocity coupling obtained through implicit time integration with a Rhie & Chow SIMPLE-type method. The discretisation is face-based. While all unknown state variables are cellcentered, the systems of equations used in the implicit time stepping procedure are constructed face by face. This technique poses no specific requirements on the topology of the cells. Therefore, the grids can be completely unstructured, and cells with an arbitrary number of arbitrarily-shaped faces are accepted.

HYDRODYNAMIC SOLVER

Computing Yacht Motions

For our 5 or 6-DoF VPP computations, a coupling between the flow, the external forces and yacht motions is necessary as the motion of the yacht modifies the flow, and as the flow applies forces on the yacht, modifying its kinematics.

Thus, ISIS-CFD uses non-linear iterations inside each single time step to solve both the Navier-Stokes equations and the yacht's dynamics at the same time and perform the coupling between flow and body motion [START_REF] Leroyer | Numerical Methods for RANSE Simulations of a Self-Propelled Fish-like Body[END_REF]. For every non-linear iteration the motions of the yacht are updated, the mesh is deformed accordingly and the flow is solved. Forces and moments applied on the yacht can then be computed using the newly obtained flow before a new non-linear iteration starts.

To handle large motions of bodies, ISIS-CFD relies on an overset method, with a finely meshed domain around the body overlapping a background domain and moving through it. These two domains are dynamically coupled using interpolation, either through a least squares approach or with a weighted distance scheme. The VPP optimisation algorithm changes the sail power parameter, thereby inducing major changes in the yacht's attitude, making the use of overset grids a must. When simulating an IMOCA and all its appendages in our VPP, up to 4 overset domains are used: the hull and the keel are within the yacht overset domain (blue in Fig. 5) which also contains two rudder overset domains (in red), as well as the sock-mesh overset domain for the foil (in yellow).

Adaptive Grid Refinement and Sock-mesh Approach

One of the most important features of FINE/Marine is the adaptive grid refinement (AGR) algorithm, with several available criteria that allow the mesh to be automatically refined locally, during the compu-tation [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF][START_REF] Wackers | Industrial application of adaptive grid refinement -the case of FINE/Marine[END_REF]. This powerful tool enhances the capturing of flow details, either for waves, wake flows, or in boundary layers around bodies, especially around hydrofoil edges.

When a good accuracy for the forces applied on bodies as well as a good understanding of the flow in the boundary layer is needed (which is typically the case in calculations involving hydrofoils) the Multi-Surface and Flux-Component Hessian (MS-FCH) criterion is the go-to criterion [START_REF] Wackers | Towards automatic resistance evaluation using mesh adaptation[END_REF]. When using this AGR algorithm, the refinement criterion is based not only on the free-surface position but also on the Hessian of both pressure and velocity. Above all, calculations using the MS-FCH criterion allow users to generate a coarser initial mesh on bodies, and to leave local mesh refinements to the AGR algorithm. When simulating a boat equipped with foils, box-shaped domains are typically employed in Hexpress (FINE/Marine's mesh generator). For these, the cells are not aligned with the curved foil geometry, so the mesh generator creates many very small cells of mediocre quality to capture the foil (see Fig. 6a).

To tackle this issue, a good practice consists in using an initial curved overset domain whose shape follows the foil's shape (see Fig. 6b) in addition to the yacht's box-shaped domain, as this allows for better-quality initial meshes. Furthermore, the benefits of using the MS-FCH criterion are reinforced when using a sock initial mesh for the foil, since the initial mesh is of better quality and lighter, thus leading to better precision and CPU-time savings [START_REF] Robin | Starting off the right foot with foil sock approach and AGR criterion[END_REF].

Fluid-structure Interaction

For IMOCA foils of roughly 6 meters length, the deflection at the tip can reach up to 1 meter in certain sailing conditions. These large deformations lead naval architects and engineers to take into account fluid-structure interactions for these appendages. Indeed, calculations and experience both show that fluid-structure interaction plays a major role in an IMOCA's foil's performance. Its bending (mostly rotation of sections around the length-wise axis of the boat) modifies the ratio of vertical and lateral lift. Its torsion has a major influence on the lift distribution, and as a consequence on the lift/drag ratio. The combination of both effects, depending upon the design of the foil, can have a stabilizing or a destabilizing influence on the yacht's behavior at sea. Accurate capture of the foil's deformations is therefore essential to correctly assess the performance of foiling or foils-assisted boats, both in a VPP-optimised steady equilibrium or dynamic, for example in waves.

As of now, FINE/Marine handles fluid-structure interaction through a modal approach [START_REF] Leroyer | Efficient and robust FSI RANSE simulations around elongated bodies for hydrodynamic applications, minimally intrusive for the beam solver[END_REF][START_REF] Mouton | Towards unsteady approach for future flutter calculation[END_REF]. The first eigenmodes deformations of the structure are computed in an external finite element solver and input to the flow solver. During the computation, structural deformations are obtained by computing the temporal evolution of the amplitude of these selected eigenmodes through a modal equation solved internally. This approach allows for CPU-time savings, since this FSI approach does not need to interact directly with a structural solver during the flow resolution. In addition to this, using a third-part external FEA solver allows to account for any internal design, shape or material characteristics. The behaviour of the structure is computed using an internal module, and so is the rigid motion. The sock-mesh approach simplifies the handling of the deformation of the mesh around the foil as the domain is of limited size and its curves follow the foil's shape.

However, the modal approach does not take into account non-linearities, and further works will compare this method to another approach where a coupling with a beam model structural solver is made.

For the naval architect, the beam approach is also interesting as it does not require the separate calculation step of the first deformation eigenmodes, and because the description of a foil as a beam, which is usually known, is directly input to the code.

RESULTS

40-feet Fast Cruiser

A VPP simulation has been carried out for a modern 40-feet fast cruiser, using the analytical aerodynamic model and the analytical rudder model covered previously. The parameters of the simulation are described in Tab. 2: In addition to Tab. 1 and Fig. 3, results for several physical quantities are presented in Fig. 7. The impact of the change in sail parameter power can be clearly observed. When initially set to the maximum value (1.0), the sail power parameter eventually reaches a smaller value (0.9343) due to the maximal heeling being attained. Then the sail power parameter is set to its minimal value (0.62), the boat rights itself, but at the cost of a loss of speed. The parameter is then set to an intermediate value (0.7771), at which the boat sails faster than for the maximal sail power parameter value. Using these results, a first optimised point is computed, with a sail power parameter of 0.8470. When used as the input of a fourth iteration, it brings a better yacht's speed, and the final parabolic optimisation eventually finds an optimal sail power parameter value of 0.8488.

IMOCA

The following calculations have been made during the design campaign of two IMOCAs for the Vendee Globe 2024. The simulated yacht has its canting keel and the leeward side foil. The interpolation in an aerodynamics matrix is used, as well as the analytical rudder model. Other parameters of the calculation are listed in Tab. 3. Results of the optimisation are presented in Tab. 4 and Fig. 8. For confidentiality reasons, some results have been normalized. The effect of the sail power parameter on various physical quantities can be observed in Fig. 9 and the effect on the yacht's attitude in Fig. 10. In this case, the maximum heel angle is set at 20°for external reasons. At the end of the first 3 iterations, it appears that the maximum speed is attained for the maximum power parameter, and no further iteration is necessary.

A larger maximum heel angle would likely result in a similar result. For this particular computation, it was expected to obtain the maximal sailing speed for the maximal value of the sail power parameter. Hence, fine convergence criteria were applied to the first stage of the computation, where the VPP code looks for an equilibrium state associated to a sail power parameter equal to 1. Coarser convergence criteria were applied during the following iterations, which have inferiors sail power parameters value, to reduce CPU-time as it was an industrial project.

DISCUSSION

When using our VPP, naval architects are able to obtain precise results to assess their designs' performance. For sailing yachts without foils, a single operating point, defined by a single (TWA, TWS) doublet, is typically calculated within half a day on a recent dual-AMD 7773X computer. For the purpose of comparing it to another design, a design can be correctly assessed with 5 to 10 carefully chosen operating points. In comparison, on the same machine, the 100+ computations required to build the hull's hydrodynamics matrix that feeds a traditional static VPP would take over a week, and require much more user time to setup and post-process.

For racing yachts with foils, the computation time is significantly longer, because of the larger number of cells required to capture the geometry of the foil with sufficient accuracy, and because of the occasional larger motions. The computation time of the required operating points is comparable to that of the hull's hydrodynamic matrix.

In this case (and in the previous case as well), the advantage of our VPP lies in the precision of the method: The interactions between the hull and the appendages, and between the different appendages, are entirely taken into account as are the surface-piercing effets of the foils. Calculations have shown these interactions to be very significant, leading to differences in leeway, for example, of up to 2°between our VPP and a static VPP with separate lifting line models for the appendages. This approach is therefore complementary to traditional VPPs which can give more operating points and can be calibrated using the high-fidelity results.

Being able to capture the behavior of the yacht in waves is currently a major topic of development in ocean-going racing yachts. Thanks to the developments currently being implemented, the ability to carry out high fidelity 6-DoF dynamic simulation of a foiling boat in these more realistic sailing conditions will be an added benefit of this VPP.

CONCLUSION

The aim of this work is to develop a VPP that will allow naval architects to perform high-fidelity studies of modern sailing boats using state-of-the-art CFD techniques. This paper showed through two yachts that such results can easily be obtained. Moreover, in the scope of industrial use, this VPP program permits to perform high-fidelity studies for a few operating points of a design, in a way that is complementary to traditional VPPs.

Like many VPPs, the results obtained using this program would need comparison to experimental or real values to correctly obtain validation. However, based on Finot-Conq's experience and on a few sea trials, the results provided by the VPP are very satisfying.

Future work includes beam-elements for FSI and their comparison to the modal approach. Ventilation and cavitation will definitely be addressed, as well as a convergence method better tailored for shipin-waves simulations.

Figure 2 .

 2 Figure 2. VPP convergence method diagram

Figure 3 .

 3 Figure 3. Example of sail power parameter optimisation

Figure 4 .

 4 Figure 4. Surge coefficient interpolation using cubic interpolations on AWA and heel. The power is fixed and equal to 1

Figure 5 .

 5 Figure 5. Use of overset domains for a typical IMOCA VPP simulation.

Figure 6 .

 6 Figure 6. Comparison of meshing methods for IMOCA VPP simulation.

Figure 7 .

 7 Figure 7. Results for a 40 feet fast cruiser

Figure 8 .

 8 Figure 8. Sail power parameter optimisation for an IMOCA.

  Figure 9. Results for an IMOCA 60' VPP simulation

Table 1 . Example of results of optimisation for a 40-feet fast cruiser (sail power range: 0.62 -1).

 1

	Iteration	Sail power parameter Yacht's speed (m/s) Heel angle •
	1	0.9343	3.32	-23.00
	2	0.6200	3.16	-11.86
	3	0.7772	3.33	-17.39
	4	0.8471	3.36	-20.01
	Optimum found	0.8488	3.36	not computed

  Averaged Navier Stokes (RANS) or RANS/Large Eddy Simulation hybrid turbulence models for multifluid flows which is developed by the METHRIC team of the LHEEA Laboratory. All calculations presented in this paper used RANS equations along k -ω SST turbulence model, with y + values defined using FINE/Marine best practices:

	y + = max	Å y + min , min	Å 30 +	(Re -1e 6 ) * 270 1e 9	, y + max	ãã	(6)
	with y + max = 300, y + min = 50 and Re the Reynolds Number.

FINE/Marine is a marine-oriented computing suite, distributed by Cadence Design Systems (formerly Numeca International) and used by many naval architects for various hydrodynamics applications. The suite features an unstructured hexahedral mesher called Hexpress, as well as the flow solver ISIS-CFD. The latter is an incompressible unsteady Navier-Stokes equations solver with various Reynold-

Table 2 . Computation parameters for the 40-feet fast cruiser simulation

 2

	TWS TWA Initial yacht's speed (m/s)	Sail plan	Power range
	10kts 50 •	8 kts	Mainsail + Jib	0.62-1

Table 3 . Computation parameters for the IMOCA 60' VPP simulation

 3

	TWS TWA Initial Yacht Speed (normalized)	Sail Plan	Power Range
	20kts 135 •	1	Mainsail + Code Zero	0-1

Table 4 . Results of optimisation for an IMOCA.

 4

	Iteration	Sail Power Parameter Yacht Speed (normalized) Heel Angle •
	1	0.9367	0.920	-19.99
	2	0.0000	0.762	-10.14
	3	0.4683	0.854	-14.23
	Optimum found	0.9367	0.920	-19.99


Heart Failure with Preserved ejection Fraction: a haemodynamic explanation of exertional dyspnoea.

Running title: haemodynamics in heart failure

Frederic MYON 1 , Elizabeth CURTIS 1 , Guillaume L'OFFICIAL 1 , Erwan DONAL This work has several strengths. Robust selection criteria thanks to using a combined score to define patients with HFpEF, built on echocardiographic functional and morphological parameters in addition to the natriuretic peptide value (HFA-PEFF score), as recommended by the 2021 European Heart Society guidelines [START_REF] Pieske | How to diagnose heart failure with preserved ejection fraction: the HFA-PEFF diagnostic algorithm: a consensus recommendation from the Heart Failure Association (HFA) of the European Society of Cardiology (ESC)[END_REF][START_REF] Authors | ESC Guidelines for the diagnosis and treatment of acute and chronic heart failure: Developed by the Task Force for the diagnosis and treatment of acute and chronic heart failure of the European Society of Cardiology (ESC). With the special contribution of the Heart Failure Association (HFA) of the ESC[END_REF]. The expired gas analysis confirms the degree of exercise intensity (maximal in all groups with a mean respiratory exchange ratio close to 1.10), which is avoids underestimating the maximal oxygen consumption. The authors used the mPAP/CO slope as the main parameter, rather than the mPAP alone, which allows the assessment of disproportionate rise in mPAP to changes in CO, as with rising CO there is a physiological increase in mPAP. Using lung ultrasound to confirm lung congestion during exercising was useful to confirm the findings(4, 5).

To link the pathophysiology with clinical outcomes, the authors then evaluated composite outcomes of all-cause mortality or heart failure events. At a follow-up of six hundred days, there was a significant difference in event free survival in HFpEF patients with Intra-and interobserver correlation coefficients are strong (0.82 to 0.86) but was based on a limited number of exams (n=40). Another challenge is echocardiographic CO estimation, the VTI curve is prone to several errors during exercise, as the sampling volume can be affected by respiration, and the clarity can impaired. The use of averaged automatic traces using the new echo-machines could be evaluated to validate the assessment during exercise [START_REF] Hasler | Pressure-Flow During Exercise Catheterization Predicts Survival in Pulmonary Hypertension[END_REF]. The authors do not give information about the main limiting effort factor (central, peripheral, or mixed) or patient's respiratory threshold value. This would have been interesting in this study where HFpEF patients and control only reach 55 to 60% of the predicted peak oxygen consumption [START_REF] Hasler | Pressure-Flow During Exercise Catheterization Predicts Survival in Pulmonary Hypertension[END_REF]. Additional parameters that could have been evaluated include the 6minutes' walk test distances and the quality-of-life questionnaires that reflect the impact of the disease on these patients [START_REF] Wolsk | Resting and exercise haemodynamics in relation to six-minute walk test in patients with heart failure and preserved ejection fraction[END_REF].

In the statistical analysis, the authors assessed the mPAP/CO slope mean values at peak of the exercise, and the timing of this can be brought into question as we know that diastolic dysfunction appears within the first few minutes of effort when the exercise load is quite weak, and the E and A waves are still distinct, facilitating analysis(2).

Heart failure is an epidemic, with more than 50% being due to HFpEF. HFpEF is not well understood and trials to improve prognosis have so far been disappointing (1) (figure 1). The HFpEF patients suffering from disproportionate exercise-induced pulmonary hypertension may benefit from pulmonary arterial hypertension specific therapy, and we are eager to see the results of ongoing studies [START_REF] Udelson | Effect of Praliciguat on Peak Rate of Oxygen Consumption in Patients With Heart Failure With Preserved Ejection Fraction: The CAPACITY HFpEF Randomized Clinical Trial[END_REF][START_REF] Armstrong | Effect of Vericiguat vs Placebo on Quality of Life in Patients With Heart Failure and Preserved Ejection Fraction: The VITALITY-HFpEF Randomized Clinical Trial[END_REF]. A better understanding of the syndrome will make the selection of targeted treatments easier [START_REF] Galli | Phenomapping Heart Failure with Preserved Ejection Fraction Using Machine Learning Cluster Analysis: Prognostic and Therapeutic Implications[END_REF].

1 1:

 1 University of Rennes, CHU Rennes, Inserm, LTSI -UMR 1099, F-35000 Rennes, France There is no conflict of interest related to this manuscript for any of the authors Correspondence Erwan DONAL Cardiology -CHU RENNES Hopital Pontchaillou 35000 Rennes Faxe: +33299282510 Tel: + 33299282525 erwan.donal@chu-rennes.fr The study of Y. Saito & al shares an interesting non-invasive exercise evaluation in heart failure with preserved ejection fraction (HFpEF) (1). They compared patients with HFpEF versus controls to shed light on the pulmonary hemodynamics during exercise to explain the exertional dyspnea reported by these patients. They assessed the trend of the pulmonary pressures registered during exercise estimated by echocardiography and its relationship with the associated cardiac output elevation. A retrospective analysis of 345 patients from a single centre tertiary Japanese hospital who were being evaluated for exertional dyspnea were divided into patients defined as having HFpEF (166), and controls (179) who did not reach diagnostic criteria for HFpEF nor pulmonary artery hypertension confirmed by right heart catheterization. The median value of the mean pulmonary arterial pressure (mPAP) over the cardiac output (CO) slope is higher in the HFpEF group (5.2 mmHg/L/min, IQR 3.4-8.0) than in the controls (3.7 mmHg/L/min, IQR 2.4-5.6). This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1002/ejhf.2877 This cut-off of 5.2 mmHg/L/min has been used by the authors to define the pathophysiological concept of 'disproportionate exercise-induced pulmonary hypertension (DEi-PH) in HFpEF'. With this threshold, the statistical analysis classified half of the group of HFpEF patients as affected by abnormal pulmonary hemodynamic response to exercise. This subgroup experienced impaired exercise capacity (lower peak oxygen consumption, worse ventilatory efficiency), poor left ventricular functional reserve (limited in CO augmentation) and right ventricle dysfunction at peak exercise (impaired RV-PA coupling, increased rightsided congestion)(1).

  DEi-PH compared to controls (*p=0.0003) and those without DEi-PH ( †p=0.01) Regarding limitations, there is a potential underestimation of pulmonary pressures during exercise due to noninvasive measurement. Only a subset of the patients underwent right heart catheterization during exercise to validate the echocardiographic parameters.

  Recent works have tried to subdivide this heterogenous and aging population in several phenotypes, to try identify therapeutic targets that may differ between subgroups (8-10). The work of Y. Saito & al advanced our understanding of the symptomatology of HFpEF matching symptoms to pathophysiology(1). This complements work by C. Lam et al showing a high prevalence of pulmonary hypertension which may be severe in the HFpEF population(11).They found that pulmonary venous hypertension is one of the main mechanisms, but does not alone explain the full severity pulmonary hypertension in HFpEF. These results are in favor of an additional role of pulmonary vascular remodeling leading to pulmonary arterial hypertension[START_REF] Lam | Pulmonary hypertension in heart failure with preserved ejection fraction: a community-based study[END_REF]. Trials are underway to assess the role of pulmonary hypertension medications in this group (Cadence, Endeavor…).With these results suggesting the role of pulmonary vascular remodeling, their findings should be prospectively validated[START_REF] Saito | Disproportionate Exercise-induced Pulmonary Hypertension in Relation to Cardiac Output in Heart Failure with Preserved Ejection Fraction -A Noninvasive Echocardiographic Study[END_REF]. Prospective evaluation and early enrollment may help to understand better why some HFpEF patients develop DEi while others don't. Repeated exercise echocardiography tests may help detect a slight decline of diastolic or right heart function prior to the development of pulmonary hypertension. Other clinical studies would be also useful to validate the threshold of 5.2 mmHg/L/min. This cut-off value differs from to the ones obtained by the sensitivity analysis of this work (3 mmHg/L/min) or from other teams predicting 3-year mortality in patient with pulmonary arterial hypertension or chronic thromboembolic pulmonary hypertension (4).

Figure 1 :

 1 Figure 1: summary slide about the value of the concept of 'disproportionate exercise-


Ioana Lal
email: ioanalal04@gmail.com

Irinel-Constantin Morȃrescu
email: constantin.morarescu@univ-lorraine.fr

Jamal Daafouz
email: jamal.daafouz@univ-lorraine.fr

Lucian Buşoniu
email: lucian@busoniu.net.

Optimistic planning for control of hybrid-input nonlinear systems 1

Keywords: Optimal control, planning, hybrid-input nonlinear systems, near-optimality analysis

We propose two branch-and-bound, optimistic planning algorithms for discrete-time nonlinear optimal control problems in which there is a continuous and a discrete action (input). The dynamics and rewards (negative costs) must be Lipschitz but can otherwise be general, as long as certain boundedness conditions are satisfied by the continuous action, reward, and Lipschitz constant of the dynamics. We start by investigating the structure of the space of hybrid-input sequences. Based on this structure, we propose for the first algorithm an optimistic selection rule that picks for refinement (branching) the subset with the largest upper bound on the value. At the price of a higher budget, the second method reduces the reliance on the Lipschitz constant, by refining all sets that are potentially optimistic. This effectively means that the Lipschitz constant is automatically optimized. The way to select the largest-impact action along which to refine the sets is the same for both algorithms, and still depends on the Lipschitz constant. We provide convergence rate guarantees for both methods, which link the computational budget to the near-optimality of the action sequences returned, in a way that depends on a problem complexity measure. We also give empirical results for a nonlinear problem, where the algorithms are applied in receding horizon, and depending on the budget either one or the other algorithm works better.

Introduction

We consider optimal control of nonlinear systems in which the inputs (control actions) consist of a continuous component and a discrete one; we refer to such systems as hybrid-input. These systems can be encountered in several fields, such as robotics [?], [?], industrial multiple-tanks systems [?], [?], [?], hydraulic systems [?], the automotive industry, for joint control of engine power and the transmission gear [?], [?], etc. A number of techniques were used for solving this type of problems, such as branch-and-bound approaches [?], switching control [?], or model-predictive control (MPC) [?], [?],

[?], [?]. Compared to these methods, our approach can handle problems with more general dynamics and cost functions, while focusing on infinite-horizon optimal control, rather than finite-horizon. In addition, nearoptimality bounds are provided, along with two fullyspecified, directly implementable algorithms (which is not always the case with other approaches). These bounds are then used for proving convergence rates for both algorithms. Specifically, we consider hybrid-input systems for which the dynamics can be generally nonlinear and the cost functions arbitrary, as long as both are Lipschitz with respect to the state and continuous action. This is not a greatly restrictive property, since usual dynamics and cost functions satisfy the constraint. The continuous action must be scalar, a restriction that can be relaxed at extra computational cost. For such systems, we propose two methods, called OPHIS and SOPHIS: Optimistic Planning for Hybrid-Input Systems, and Simultaneous OPHIS. Both algorithms are a nonstandard flavor of MPC, produce at each given state an open-loop sequence, and are meant to be applied in receding horizon. The two algorithms also belong to the optimistic planning (OP) class of algorithms. OP approaches explore the space of infinitely long sequences of actions, and re-fine regions in which optimal solutions may be located.

OPHIS creates a partition of the set of hybrid-input sequences iteratively, by choosing for refinement at each iteration an optimistic set, i.e. one that maximizes an upper bound on the value. For any set that is chosen for refinement, a dimension (time step) is also chosen, together with the type of split (continuous or discrete).

OPHIS depends on the Lipschitz constant in two places: set selection and dimension selection. In practice, the Lipschitz constant is difficult to estimate: a too large value makes the algorithm slow, and a too small value invalidates it. Therefore, we introduce our second algorithm, SOPHIS, which removes the dependence on the Lipschitz constant in set selection, by refining several sets per iteration: any that may be optimistic regardless of the value of the Lipschitz constant. Although dimension selection still has to depend on the constant, we expect that set selection has a larger impact, and this intuition is confirmed by experiments in Section 6: SOPHIS is less sensitive to the value of the Lipschitz constant.

Regarding experimental performance, both algorithms have their use: SOPHIS works better for larger budgets, whereas for smaller ones, the OPHIS approach of focusing this limited budget on one value of the Lipschitz constant still pays off. Analytically, the convergence rate of SOPHIS (when properly tuned) is almost as good as that of OPHIS. However, this is not the full story: the fact that SOPHIS expands sets for all possible Lipschitz constant values means in effect that the rates are those for the best possible value of the constant; in effect, SOPHIS automatically optimizes the Lipschitz constant for the set selection component.

The first method, OPHIS, has already been described in our preliminary conference paper [?], where the focus was on deriving the method and empirically validating it in two separate problems. In addition, only an a-posteriori guarantee was given for OPHIS. In this paper we introduce an extension to this algorithm, namely SOPHIS. Furthermore, we provide proof of convergence rates to the global infinite-horizon discounted optimum for both methods. We define a complexity measure of any problem to which (S)OPHIS is applied, and prove that the algorithms converge faster for smaller complexity measures. In particular, in the simplest type of problem, convergence is according to an exponential in a power of the computational budget. This is significant because it gives strong guarantees of near-optimality, which to our knowledge, are not given for other hybridinput methods in the literature [?], [?], [?].

We exemplify OPHIS and SOPHIS in simulations. Both methods succeed in controlling a robot arm where one link is continuously controlled and the other can have a brake either applied or not. OPHIS wins for small budgets, but SOPHIS is better for large budgets and less sensitive to the Lipschitz constant value.

OPHIS essentially combines optimistic planning for deterministic, discrete-input systems (OPD) [?], [?] and OP for continuous-input systems (OPC) [?]. The new algorithm SOPHIS, which does not need the Lipschitz value in the set selection, is similar to the extension of OPC to SOPC in [?]. The key novel technical challenge here is that the structure of the hybrid space of solutions is significantly more complicated than for either OPC or OPD, and consequently so are the refinement rules that we must come up with to explore this space. OPHIS in fact specializes to OPD when the continuous action is removed, and to OPC when the discrete action is removed; and SOPHIS specializes to SOPC.

A branch and bound approach related to optimistic planning is used in [?], in combination with sparse direct collocation. However, no near-optimality analysis is provided there. A key difference between usual hybrid-input control approaches and (S)OPHIS is that the former primarily concentrate on stability, such as in [?], where a switching control strategy is employed, whereas here we aim for near-optimality. In effect, by using discounting and imposing a joint condition on the discount factor and Lipschitz constant of the dynamics, the system dynamics are instead required to satisfy a certain contractiveness property. Promising guarantees of stability have been obtained both for the exactly optimal solution of general discounted optimal control [?] and for discrete-input optimistic planning [?]. However, the behavior of (S)OPHIS when refining continuous actions is much more intricate, and analyzing its stability is left for future work.

In the landscape of tree search methods [?], OP is a bestfirst method, since it aims to refine the optimistic set. Compared to Monte-Carlo tree search [?], which uses random sampling to make the search more efficient, the OP flavor that we use has the advantage of providing deterministic guarantees.

Next, Section 2 formalizes the problem and Section 3 discusses the background on OPC, SOPC, and OPD. OPHIS and SOPHIS are described in Section 4, while Section 5 provides the convergence rate analysis. Finally, Section 6 gives simulation results for a two-link robot arm, and Section 7 gives the conclusions of this paper.

Problem statement

We consider an optimal control problem for a hybridinput, nonlinear system x k+1 = f (x k , u k ), with x ∈ X ⊆ R m and u ∈ U consisting of a continuous action and a discrete one:

u k = [c k d k ] T
(1) where c k ∈ R and d k ∈ {0, 1, ..., p}, p ∈ N. We define a reward function ρ : X × U → R, that takes as input a state-action pair (x k , u k ): r k+1 = ρ(x k , u k ). Starting from an initial state x 0 , we define an infinitely-long sequence of actions u ∞ = (u 0 , u 1 , ...) and its infinitehorizon discounted value:

v(u ∞ ) = ∞ k=0 γ k ρ(x k , u k ) (2)
where γ ∈ (0, 1) is the discount factor (note that γ = 1 is excluded). We aim in principle to find the optimal value v * := sup u∞ v(u ∞ ) and a sequence that achieves it.

Assumption 1. We have (i)

r k ∈ [0, 1] and (ii) c k ∈ [0, 1].
Together with discounting, the bounded rewards ensure boundedness of the sequence values. Bounded costs are typical in e.g. reinforcement learning for control [?] and they could follow either from physical limitations in the system, or from saturating an a priori unbounded reward function, which changes the optimal solution but may be sufficient. Note that now U = ([0, 1] × {0, 1, ..., p}).

The restriction to scalar continuous actions can be relaxed, but at significant extra computational cost for the extended algorithm; see [?], supplementary material for such an extension in OPC. The growth would be exponential in the number of inputs.

Assumption 2. (i) Both the dynamics and the rewards are Lipschitz with respect to the state and the continuous action, i.e., ∃L f , L p s.t. ∀x, x ∈ X and c, c ∈ [0, 1]:

f (x, [c, d] T ) -f (x , [c , d] T ) ≤ L f ( x -x + |c -c |) |ρ(x, [c, d] T ) -ρ(x , [c , d] T )| ≤ L ρ ( x -x + |c -c |) (3) (ii)
The following inequality is satisfied:

γL f < 1 (4)
It should be noted that Lipschitz continuity (i) is only imposed w.r.t. the continuous component c of the action; whereas the variation w.r.t. d can be arbitrary. This is a useful feature because switches often induce discontinuities in the system. Note also that condition (i) allows the dynamics and rewards to be nondifferentiable as long as they are still Lipschitz. This helps to model effects like saturations, actuator dead-zones, etc. Condition (ii) means that the dynamics need not be strictly contractive on their own, but should become so when combined with a shrink rate equal to the discount factor γ.

Lemma 3. For any two sequences u ∞ , u ∞ ∈ U ∞ , we have:

|v(u ∞ )-v(u ∞ )| ≤ L ρ D-1 k=0 |c k -c k |γ k 1 -(γL f ) D-k 1 -γL f + γ D 1 -γ (5) where D is equal to the first step k at which d k = d k .
The proof of Lemma 3 can be found in [?]. The property in Lemma 3 drives our entire algorithm: actions will be prioritized for refinement according to their importance in (5). The two terms on the right-hand side of the inequality correspond to continuous and discrete actions, respectively, and are fundamentally different from each other. When there is no continuous action, the summation disappears and the formula simplifies to γ D 1-γ , the OPD metric [?]. Conversely, eliminating the discrete action is equivalent to taking D → ∞, so we get

Lρ 1-γL f ∞ k=0 |c k -c k |γ k , the OPC metric [?].
3 Previous optimistic planning algorithms OPHIS specializes to OPD when there is no continuous action, and to OPC when there is no discrete action. Under some conditions, SOPHIS also specializes to SOPC when there is no discrete action. Therefore, understanding these basic algorithms is important.

Optimistic planning for continuous actions (OPC) aims to find an infinitely-long sequence of continuous actions c k that maximize the objective function v, without discrete component d. The full explanation is given in [?]; here we will provide a short description of the algorithm. OPC refines a collection of infinitely-dimensional (hyper)boxes of the form (µ

0 ו • •×µ K-1 ×[0, 1]×[0, 1] • • • )
where µ k is the interval of actions at step k, and K is the first unrefined dimension; from there on, all the intervals are full, [0, 1]. OPC starts with the full set of sequences for K = 0, and iteratively refines it by selecting at each iteration an optimistic set with the largest upper bound on the value, and splitting it into M equal pieces along a well-chosen dimension. For each box, OPC needs to simulate the system with the sequence at the center of the set, up until step K -1, and store the resulting state and reward trajectories. Finally, OPC returns such a center sequence that maximizes the discounted value along these K steps. SOPC [?] is an extension of OPC, where all the sets that may have the largest upper bound for any Lipschitz value are expanded. Both the selection and splitting rules, through modifications, become independent from the Lipschitz constant.

Optimistic planning for discrete actions (OPD) [?] is an algorithm used for systems with discrete inputs. It works by building a tree, starting from a root node corresponding to the entire set of possible actions {0, . . . , p} ∞ . Then, at each step, an optimistic leaf node is chosen for expansion, by maximizing an upper bound. Each node is expanded by making its next discrete action definite. For instance the root node will have p + 1 children, one for each value of d 0 , and the remaining actions remain free. Expanding any level-1 node makes the action d 1 definite with p + 1 children at level 2, and so on. OPD returns a sequence on the tree with maximal sum of discounted rewards for the definite actions.

Hybrid-input algorithms

In the sequel we derive two novel algorithms that can search for sequences of hybrid inputs. The main idea is to iteratively split an optimistic set of inputs like in OPC and OPD, but unlike those algorithms, by refining either the discrete or the continuous action. For this, we look at the uncertainty on the value of the set, and choose the action with the largest impact on the uncertainty.

The key technical novelty compared to OPC and OPD is that continuous-and discrete-action refinements must be alternated, in a way that is dictated by the intricate geometry of the space of hybrid-input sequences.

A set is represented by an interval µ for each continuous action and a discrete action set σ for each refined step k:

S i = ∞ k=0 (µ i,k , σ i,k ) (6)
where by the product of sets we mean the repeated application of the cross-product ×, and notation (µ, σ) means a set in which c ∈ µ and d ∈ σ. For clarity, from now on we will refer to the set per step k, (µ i,k , σ i,k ), as a pair, and the infinite-horizon S i as a set.

A set also has two characteristics: D i and C i , representing the number of refined discrete and continuous dimensions, respectively. For all

k ≥ C i , µ i,k = [0, 1]. For all k < D i , σ i,k = d i,k
, a single, definite value, and for all k ≥ D i , σ i,k = {0, 1, ..., p}. A sequence of actions corresponding to each set is then (u i,0 , u i,1 , u i,2 , ...), where:

u i,k = c i,k , d i,k T (7) and c i,k ∈ µ i,k , d i,k ∈ σ i,k .
A continuous split can be done along any dimension k ≤ C i , by dividing the interval µ i,k into M equal pieces and thus generating M new sets. A discrete split is always done for dimension k = D i , by adding p + 1 new sets that make discrete action d k definite, one set for each discrete possibility. Note that the ways of splitting continuously and discretely, respectively, are fundamentally different. We provide example splits of each type below.

OPHIS

In this subsection, we focus on the OPHIS algorithm.

In order to decide which set to split, let us first consider reward r i,k+1 = ρ(x i,k , u i,k ), where with a slight abuse of notation we now refer by c i,k to the specific action that is at the center of interval µ i,k . Define then the sample value of a set i:

v(i) = Di-1 k=0 γ k r i,k+1 (8)
Each continuous interval µ i,k has a certain length a i,k .

For k ≥ C i , a i,k = 1. For each set, we define its diameter in the semimetric of ( 5):

sup u∞,u ∞ ∈Si |v(u ∞ ) -v(u ∞ )| ≤ δ(i) δ(i) = L ρ Di-1 k=0 a i,k γ k 1-(γL f ) D i -k 1-γL f + γ D i 1-γ (9)
Given the sample value of a set i and its diameter, we have the upper bound:

B(i) = v(i) + δ(i) (10) so that v(u ∞ ) ≤ B(i), ∀u ∞ ∈ S i .
This inequality is shown as follows. By the definition of set S i , for all k < D i we have

d k = σ i,k and c k ∈ µ i,k , so |c k -c i,k | ≤ a i,k
. All the quantities without subscript i refer to the sequence u ∞ . Thus:

∞ k=0 γ k r k+1 = Di-1 k=0 γ k r k+1 + ∞ k=Di γ k r k+1 ≤ v(i) + L ρ Di-1 k=0 a i,k γ k 1 -(γL f ) Di-k 1 -γL f + γ Di 1 -γ = v(i) + δ(i) = B(i) (11)
where the bound on the first summation (in square brackets) follows as in [?], and the bound on the second summation holds because all rewards are at most 1.

Denoting A as the collection of all sets, OPHIS iteratively selects for refinement at each iteration an optimistic set that maximizes the upper bound:

i † = arg max i∈A B(i) (12)
In order to decide whether we have a continuous or discrete split of S i † , and along which dimension, we look at the contribution of each dimension k up to D i † -1 to the diameter (9), as well as at the contribution γ D i † 1-γ of the first unrefined dimension. Whichever contribution is the greatest dictates where we split. Thus:

k † = arg max k∈{0,1,...,D i † } {L ρ a i † ,k γ k 1 -(γL f ) D i † -k 1 -γL f } (13) and if L ρ a i † ,k † γ k † 1-(γL f ) D i † -k † 1-γL f ≤ γ D i † 1-γ ,
we have a discrete split, at dimension D i † . Otherwise, we have a continuous split, along dimension min(k † , C i † ). Recall that for dimensions k between C i † and D i † -1 the size a of the interval is 1. Further, note that by this rule, we always have D ≥ C for any set.

For compactness reasons, we shall denote the contribution of a dimension k in the continuous part of the diameter with

λ k = L ρ a k γ k 1-(γL f ) D-k 1-γL f .
An example is given in Figures 1 and2, for a continuous refinement and a discrete one, respectively. In both cases, we start from a set i, which already had 6 discrete refined dimensions (D i = 6), and 5 continuous discretized dimensions (

C i = 5). The set is S i = {([4/9, 5/9], 1) × ([1/3, 2/3], 0) × ([2/3, 1], 0) × ([0, 1/3], 1) × ([1/3, 2/3], 1) × ([0, 1], 0) × ([0, 1], {0, 1, ..., p}) ∞ }.
We take the center of the interval as the continuous action, and thus we will have the following sequence of actions (1/2, 1), (1/2, 0), (5/6, 0), (1/6, 1), (1/2, 1), (1/2, 0) . Then, f and ρ are called at each step, for this sequence of actions from x 0 , to determine both the sequence of states and rewards for set S i . In this example, there are 2 possible discrete actions, 0 and 1, and M = 3.

Figure 1 shows a continuous split, along dimension 1. In the figure, one can observe in black the parent set, from which M = 3 new sets are formed, having the same continuous intervals for all dimensions other than k = 1, and the same discrete actions. The number of refined continuous dimensions C remains 5 for all children sets and D remains 6. The resulting intervals at k = 1 are shown with different colors and styles. At the top of the figure, we have added, symbolically, a trajectory that represents the states and the rewards. These are the same until step k = 1 (the refined dimension), and differ afterwards,since the continuous actions at step 1 will be 7/18, 1/2 and 11/18, respectively. The middle set will have the same continuous action and trajectories as the parent, so these trajectories can be reused.

Figure 2 shows what a discrete split looks like, starting from the same parent set. A discrete split implies always refining dimension D i , in this case 6. One can see the newly added children with colors and different line styles. They inherit the continuous intervals from the parent, as well as the previous discrete actions. The sample values are the same for the children sets, until dimension k = 6, and then they differ, based on the new discrete action.

Overall, we create a tree structure, where each set is a node in the tree. Whenever a discrete split is done for set i, p + 1 children are added to the node representing set i in the tree. In the case of a continuous split, M children are added instead. An example of such a tree can be seen in Figure 3. Discrete splits are represented by blue, continuous lines and continuous splits by red dotted lines. As we can observe, for the root node, all continuous intervals are [0, 1] and all discrete actions are undefined. Then, by a discrete split, we get two new sets, where the first discrete action is defined as 0 for set S 1 , and 1 for set S 2 , respectively. In the same way, 2 new sets are formed in Figure 2, by defining the discrete action at step D i = 6. Then, set S 2 is chosen for refinement, and a continuous refinement on dimension 0 is done. This adds sets S 3 , S 4 , S 5 , each with a third of interval for the first dimension. Again, this is similar with the continuous split done in Figure 1, where 3 new sets are added, each with the interval for the second dimension being a third of the interval of their parent set.

Note that even though the interval refinements are discrete at each step, asymptotically, after a large number of expansions, the algorithm can reach arbitrarily close to any continuous value; in other words, the set of applicable actions is dense in the unit interval.

We have mentioned before that OPHIS specializes to OPC when there is no discrete action, and to OPD when there is no continuous action. OPC chooses a set to expand based on an upper bound given by b

(i) = v(i) + δ(i) := Ci-1 k=0 γ k r i,k+1 + max( Lρ 1-γL f , 1) ∞ k=0 γ k a i,k .
The maximum is there to cover with a unified formula the contribution of discretized and undiscretized (unsimulated) dimensions [?]. A dimension to refine is chosen in OPC as k † = arg max k=0,...,Ci γ k a i † ,k , without comparing to the discrete-action contribution because there is none. Note that the OPC diameter is differently structured only for convenience reasons, and in fact a tighter diameter can be written:

δ(i) = L ρ Ci-1 k=0 a i,k γ k 1 -(γL f ) Ci-k 1 -γL f + γ Ci 1 -γ (14)
This follows in the same way as the diameter (9) of OPHIS, except that now instead of stopping at D i , we stop at C i because we have not discretized actions further so their rewards are unknown.

For OPD, the optimistic set is chosen for expansion based on an upper bound given as b(i) = v(i) + γ D i 1-γ ; we do not need to take into account any continuous-action deviations from the center sequence.

So that we are able to reuse the center sequence at a continuous split, we impose for the remainder of the paper that M is odd (this will also help in the analysis). Next, a pseudocode of the method is given in Algorithm 4.1. We must pass the model of the system, as well as the initial state. The initial set S 0 represents the root of the tree, with all the continuous intervals being µ = {[0, 1], [0, 1], [0, 1], ...} and the discrete actions not yet defined as certain values, but allowing the sets of all possible values {0, 1, ..., p}. This is because no refinement has been done yet; as such, we also have that D 0 = C 0 = 0. An important parameter is the budget n , which represents the number of simulations of the system dynamics that we are allowed to perform. The algorithm outputs a near-optimal sequence of actions.

SOPHIS

In this section, we discuss an extension of the OPHIS algorithm, in which we have a different set selection rule. To avoid assuming knowledge of the Lipschitz constant for this rule, we will expand all the sets that may be optimistic for any value of this constant. Note however that the Lipschitz constant must still satisfy Assumption 2. Recall that both methods create a tree of sets, and denote by H the depth in this tree. Quantity H represents the sum of continuous and discrete expansions done in order to reach a certain set. Since all sets have the same shape, diameters δ(i) are the same at a given depth, so the maximum-upper-bound set at that depth can only be a set with the largest value v(i). Thus, at each depth Algorithm 1 OPHIS 1: input state x 0 , model f , ρ, split factor M , discrete set {0, 1, ..., p}, budget n, Lipschitz constants L f and L ρ , discount factor γ 2: initialize collection of sets A with S 0 , D 0 = 0, C 0 = 0; 3: while budget still available do 4:

select set i † = arg max i∈A B(i); 5:
select dimension with max contribution for continuous actions k † = arg max k∈{0,1,...,D i † } λ k ;

6:

if λ k † ≤ γ D i † 1-γ then / * split
all children will have D = D i † and C = C i † if k † = C i † , or C = C i † + 1 if k † = C i † ; 16:
end if 17: end while 18: output sequence û of set i * = arg max i∈A v(i)

H that still has nodes unexpanded, we expand the set with the greatest v value among all sets at that depth.

To prevent expansions from continuing indefinitely we also configure a maximum depth H max (n) up to which the expansions are allowed to continue when the budget spent so far is n. If H max grows fast with budget n, the algorithm will favour deep searches. On the contrary, a slower growth with n allows us to do a search focused on breadth. More insight into choosing H max will be provided in the analysis. Of course, expanding several sets per iteration comes at an extra cost. However, as we will show in the analysis, this does not have a great impact on the guarantees and in simulations performance is often improved.

The dimension selection rule for SOPHIS will remain the same as for OPHIS. This means that it will unfortunately still depend on the Lipschitz constant, and there is no way to avoid this.

As previously mentioned, the extension from OPHIS to SOPHIS is similar to the one from OPC to SOPC [?]. Both SOPHIS and SOPC refine several sets per iteration, based on the best value at each depth. Moreover, SOPHIS would specialize to SOPC if there was no discrete action and if SOPC were to use the tighter diameter ( 14 end if 19: end while 20: output sequence û of set i * = arg max i∈A v(i)

select set i † = arg max i∈A v(i) ; 10: select dimension k † = arg max k∈{0,1,...,D i † } λ k ; 11: if λ k † ≤ γ D i † 1-

Discussion of OPHIS and SOPHIS

The following remarks apply for both OPHIS and SOPHIS. We first discuss the algorithms inputs that are selected by the user. The budget should, of course, be taken as large as computationally feasible to get as close as possible to the optimal solution. Moreover, as previously mentioned, we set M to be odd, such that we can reuse the middle sequence when we have a continuous split. This works because we always consider the middle of the interval µ i,k as the continuous action. We suggest taking M to be 3, the smallest feasible odd value.

Regarding now the Lipschitz constants L f and L ρ , while in principle they are given by the problem, in practice they may be difficult to compute (especially L f ), or computing them might give overly conservative values that work poorly across most of the state-action space. Thus we suggest treating them as tuning parameters. Similarly, γ may be fixed by the problem objective; if it is not, it can be treated as a tuning parameter that should not be very far from 1. Larger γ will promote looking for longer-horizon solutions at the expense of refining less the continuous actions, while smaller γ will refine more the continuous actions at the expense of the horizon.

Both algorithms should be applied in receding horizon, by running them at each step, applying the first action from the returned sequence, and then repeating the procedure. Overall, a closed-loop control scheme is obtained. Moreover, we assume a setting in which simulating the system dominates computation, so to obtain a measure of the required computation, we need to look at the number of calls made to the dynamics f and the reward function ρ. For a discrete expansion, we make p + 1 calls to simulate the new discrete step with each of the p + 1 discrete actions and continuous action 0.5. In case of a continuous split of set i † at dimension k † , we need (M -1)(D i † -k † ) calls to simulate the M -1 sequences (except the reused center one) from step k † to step D i † .

Recall that the continuous input c is bounded and in [0, 1]. This unit interval can be obtained by rescaling a different bounded interval, so bound constraints on c are natively supported. Other constraints for the continuous action can also be implemented in the algorithm, by excluding subtrees of solutions that would violate the constraint. However, the analysis that follows in Section 5 will not take such constraints into account.

Analysis

In this section, we prove that the sub-optimality of the OPHIS and SOPHIS algorithms converges to 0 at certain rates with increasing budget n. Recall that the budget represents the number of system dynamics simulations that we are allowed to do. First, we will find an upper bound for the diameter of an arbitrary set i at some depth H in the tree, in Section 5.1. For this, we will look at the way the number of splits per dimension evolves. Then, in Section 5.2, as preliminaries for getting the convergence rates of the algorithms, we define a near-optimal tree that the algorithms focus on expanding, along with its branching factor. Afterwards, using this branching factor, we will establish a relation between the depth H and the budget n for OPHIS and SOPHIS, in Sections 5.3 and 5.4, respectively. Each of these relations, by replacement in the diameter formula, leads to convergence rates: for the OPHIS algorithm in Theorem 11, and for SOPHIS in Theorem 13.

It is important to note that we follow similar steps as [?], but we heavily adapted them to include the fact that we also have discrete splits. The OPC and SOPC algorithms in [?] only deal with continuous splits, and do not have such an intricate space of expansions. We also have fully new questions that did not arise before, such as the difference between the number of discrete split dimensions D and the number of continuous ones C.

Splitting behavior and upper bound for diameter

First, we define the depth of a given set as being the sum of continuous and discrete splits made in order to reach the set. The number of discrete splits for a set will be D, and the total number of continuous splits is denoted by h. Therefore: First, we will derive an upper bound for the diameter. We will focus in the beginning on the number of continuous splits, h. Let us define the continuous splitting function s(k), which is the number of continuous splits per dimension k. An example of this function is given in Figure 4. This function is decreasing with at most 1 at each k, as we will see next. There will be ranges of dimensions for which the splitting function is the same. Quantity q represents the number of ranges at the end, that have a smaller length than the rest.

H = h + D (15)
Assumption 4. M γ > 2.
This assumption is easy to satisfy, since M is at least 3, so as to reuse the center sequence, and generally, we would use values for γ of at least 0.7. Smaller values for γ determine a very short-horizon, almost myopic objective that will likely not work in most problems.

In what follows, we take a better look at the s-ranges.

Lemma 5. a) The first k in a constant s-range is preferred for refinement to later dimensions in the same range. b) s(k) decreases in steps of at most 1.

Apart from the proof of Lemma 10, which was already given in [?], the proofs for all theorems and lemmas, including that of 5, are given in the supplementary material at http://busoniu.net/files/papers/ sophis_suppl.pdf.

Next, we need a better understanding of how the s ranges look. For this, we will denote the lengths of the ranges with τ 0 , τ 1 , ..., τ N , with N being the last range, which is infinitely long and for which s = 0, meaning that there has been no continuous split for those dimensions. We will seek upper and lower bounds on these range lengths.

Lemma 6. For any set, we have:

       τ 0 ≤ τ * τ j ∈ {τ * -1, τ * }, 1 ≤ j < N -q τ j < τ * -1, N -q ≤ j < N τ N = ∞ (16)
where τ = log(M ) log(1/γ) and τ * = τ . This means that other than the first and last q ranges, all s-ranges have either the length τ * or τ * -1.

Next, we look at the difference between C and D, denoted G, which is nearly constant.

Lemma 7. Let G be the smallest positive integer j, for which

L ρ 1-(γL f ) j 1-γL f > γ j 1-γ .
G increases from 0 to G in the beginning, and then it always oscillates between G and G -1. This is important in our analysis, where we will use upper and lower bounds on C, and since G is roughly constant, bounds on C will translate to bounds on D.

Overall, Lemmas 5, 6, and 7 are important in the following way for determining the upper bound on the diameter. Since the contribution to the diameter of a continuous dimension depends on the interval length a k = M -s(k) , we need to find a lower bound on s(k) for a fixed h. Knowing how the s-ranges look like from Lemma 6, we can get a lower bound on s for a fixed C, followed by both upper and lower bounds on C. This will eventually give us an overall lower bound on s, as a function of h and of the lower bound on C. With the help of Lemma 7, this translates to a function of h and the lower bound on D, which can then be utilized in the diameter formula. Detailed proofs are in the supplementary material mentioned above. Overall, this results in the following upper bound for the diameter: Lemma 8. For some positive constants c 1 , c 2 , c 3 and any set i at any depth

H = h + D, δ(i) ≤ c 1 c 2 + 2h(τ * -1) γ √ 2h(τ * -1)τ τ * + c 3 γ √ 2h(τ * -1) .
To understand this result, note that the first term in the summation dominates the second one. Also, due to the fact that γ < 1, the term γ √ 2h(τ * -1)τ τ * asymptotically dominates 2h(τ * -1), causing the convergence to 0 of the diameter. Using the asymptotic notation Õ, we will have:

δ(i) = Õ γ √ 2h(τ * -1)τ τ * (17)
The Õ notation is derived from the Bachmann-Landau notation, O. When the order of complexity of some algorithm is O(f (y)) log(f (y)), we say Õ(f (y)) = O(f (y)) log(f (y)), ignoring the logarithmic term. Here, f (y) is any function; in (17), f is the power of γ and generic argument y is replaced by h.

Overall, we can observe that we have an exponential decrease in √ h, as in [?]. However, the analysis is much more involved, as we must take into account the discrete splits as well.

Near-optimal tree and branching factor

In this subsection, we discuss the tree of near-optimal nodes and define its branching factor m.

Let us define the set of near-optimal nodes at depth H:

T * H = {i at H | v(i) + δ H ≥ v * } (18)
This is a sub-tree of the full tree explained earlier, in Section 4.1, with an example in Figure 3. Both OPHIS and SOPHIS will refine nodes from the near-optimal tree. Now, denote the branching factor as m, as defined below.

Definition 9. The asymptotic branching factor is the smallest m such that ∃C ≥ 1 for which |T * H | ≤ Cm H , ∀H, where | • | represents the cardinality of the set.

Branching factor m represents the complexity of the problem. In case of a smaller branching factor, the problem is simpler. The least complex problems correspond to m = 1, meaning that only the optimal path will be explored. It is important to note that m is not necessarily an integer. It is at least 1, but its maximum value depends on both M , the number of newly formed sets in case of a continuous split, and p + 1, the number of possible discrete actions, representing the number of newly formed sets in a discrete expansion. We define the maximum possible branching factor at any node as Z = max(M, p + 1). Therefore m ∈ [1, Z], different from the branching factors defined in [?] and [?], where there are either only continuous splits, or discrete ones.

Convergence rate of OPHIS

So far, we have discussed the way that the splits look, in the end getting a relation between the diameter and the sum of continuous splits. Now, we want to go further, to get a connection between the depth H and budget n, and link both to the near-optimality. It is important to remember that the depth H is the sum of h and D. Firstly, we provide a bound equal to a diameter on the near-optimality of the sequence returned that is explicitly available for use a posteriori, once OPHIS has run, in Lemma 10. Then, the relation between the depth and the budget, together with G allows us to use (17), which is stated in h. This, in turn, leads to a relation between the diameter and the budget, given in Theorem 11. In the end, taking into consideration also Lemma 10, we get a connection between the sub-optimality and the budget.

The a posteriori bound on the near-optimality of OPHIS is the following [?]: Lemma 10. The sequence i * returned by the algorithm satisfies:

v * -v(i * ) ≤ δ min where δ min is the smallest diameter among all the sets expanded by the algorithm.

Such properties are standard for OP algorithms; both OPD and OPC ensure similar bounds. However, the proofs are different, as it can be seen in the supplementary material. Notably, a required property of increasing set values along continuous splits is not trivial for OPHIS, since we need to handle discrete and continuous splits at the same time. Note that we can also compute lower and upper bounds on the optimal value: v * ∈ [v(i * ), B * ], where B * is the smallest upper bound of any set expanded. Such bounds are popular in hybrid systems, where they are called certification bounds [?].

Next, we give the relation between the budget and suboptimality of the algorithm, an a-priori convergence rate guarantee.

Theorem 11. For large budget n, we have:

a) For m > 1: v * -v(i * ) = Õ γ 2τ 2 (τ * -1) log n τ * 2 log m b) For m = 1: v * -v(i * ) = Õ γ n 1/4 τ τ * 2(τ * -1) ZC
Note that both quantities reach 0 asymptotically. The convergence rate depends on the complexity m. We have a faster decrease for a less complex problem (smaller m) and vice-versa. This is similar to OPC in [?]. For m = 1, we have a convergence to 0 with n 1/4 , as in OPC [?].

Convergence rate of SOPHIS

In this section, we want to find a lower bound on the depth H reached by SOPHIS given a budget n, which ensures we expand a node containing an optimal solution at that depth; and then combine the lower bound with Equation ( 17) to get a convergence rate of this second algorithm. Recall that H also takes into account the discrete expansions, being the sum of all splits needed for a certain set.

Lemma 12. Given the budget n, define H(n) to be the smallest depth for which the following inequality holds:

CZH 2 max (n) H H =0 m H ≥ n (19)
Then, SOPHIS expands a node containing the optimal solution at depth H = min{H(n), H max (h)}, and the sequence returned is δ H -optimal.

Recall that Z = max{M, p + 1}. Reference [?] proves this in the supplementary material. Their proof applies directly here, where instead of M , we have Z, and instead of h, we have H. Next, we find the convergence guarantees of the SOPHIS algorithm.

Theorem 13. Consider the sequence û and the corresponding set i * returned by SOPHIS. For large n: a) for m > 1, we take H max = n , with ∈ (0, 0.5) and we have:

v * -v(i * ) = Õ γˆ τ τ * (τ * -1)(1-2 ) log n log m
b) for m = 1, we take H max = n 1/3 , and we have:

v * -v(i * ) = Õ γˆ n 1/6 τ τ * 2(τ * -1) min{ 1 CZ , 1}
We can observe that for m > 1, the exponent is 1-2 2 times smaller than the one of OPHIS. This means that we lose a bit of convergence speed compared to the first method, but not a lot if is small. Again, we see that the branching factor m is a key factor in the rate, with a smaller m leading to a faster convergence to 0. For m = 1, the rate decreases with n 1/6 , which is a bit of a loss in comparison to OPHIS, but still good, as it is an exponential of a power of the budget n.

One issue that we have to discuss for SOPHIS is the choice of H max . As in practice we do not know the actual branching factor m, we cannot set H max "optimally", i.e. per the rule that for m = 1, n 1/3 should be used, and for m > 1, n with ∈ (0, 0.5) is needed. A practical solution is to first take H max as n 1/3 , in the hope of a simple problem. In case the problem turns out to be more complicated, this leads to a slower convergence rate, but the algorithm remains valid.

This concludes our a-priori convergence rate analysis for the two algorithms. The upper bound for the diameter was a key step in establishing a relation between the budget n and the convergence rate. We were then able to prove that the sub-optimality of both OPHIS and SOPHIS converge to 0 at different rates. This proves the performance of our algorithms. Next, we give some simulation results for both of them.

Simulation results

This section discusses the results of applying the algorithms for a hybrid-input problem. A two-link robot arm example is presented, with simulation results given for both OPHIS and SOPHIS. In a previous paper [?], OPHIS was also applied to a simple, two-tanks problem, where it succeeded in recovering existing results from the literature.

The robot arm has 2 joints, one continuously controlled and one which can either have a brake set or not will either take the value 0 or 1, and the maximum value of τ 1 will be 20, which will be rescaled to 1 for continuous action c. The numerical integration is done using the Euler method, with 5 integration steps per control sampling time, which is ∆ = 0.05s. The goal is to get the state to a desired setpoint x f . The reward function is:

ρ(x k , u k ) = (10 -|x 1 k+1 -x 1 f | -|x 3 k+1 -x 3 f |)/10 (20)
where x k+1 is the state at the next step. We use a nondifferentiable reward to show that the algorithm works in this case.

The following values were used for both algorithms: M = 3, budget n = 5000, L f = 0.8, L ρ = 1.2, γ = 0.8. The starting position is x 0 = [1.2, 0, 0.8, 0]. The desired final state is x f = [π/2, 0, -π/2, 0]. In addition, for SOPHIS, we take H max = n 0.35 . Recall that both algorithms are applied in receding horizon, and the simulation is done over 12 seconds.

The results for OPHIS can be seen in Figure 5, which shows the evolution of the states and actions in time.

The top subplot presents the angles, and the middle one the angular velocities. With blue continuous lines we can see the states corresponding to the first, actuated, joint. The red, dashed lines represent the states of the second joint, which can only be influenced by a holding brake. The last subplot presents the evolution in time of the 2 actions: with blue continuous line c and with red dashed line d. Both angles reach their desired setpoints, which are represented in black dotted lines. Also, the final velocities are oscillating around 0, and the brake is not set in steady state.

Figure 6 presents the results using SOPHIS. with the subplots structured as in Figure 5. As we can see, the angles reach their reference values (in black dotted line) much faster than with OPHIS. This is also the case for the angular velocities, which can be seen stabilizing around 0 much quicker. We notice fewer switches of discrete input d than when using OPHIS. Overall, SOPHIS gives better empirical results. We also compare OPHIS and SOPHIS in terms of discounted return, for several L f values and different budgets. The comparison is presented in Figure 7. For small budgets, OPHIS is better since a focused search on one value of the Lipschitz constant likely makes more sense when computation is limited. When the budgets are larger, SOPHIS starts actually exploiting its potential by spending the extra computation to expand sets for many possible values of the Lipschitz constant, which in effect is similar to automatically finding the best value of this constant (for set selection). So we see that it starts dominating OPHIS. In addition, as presented before, for a large budget, SOPHIS manages to bring the states to the reference value much quicker. Furthermore, we observe flatter curves in Figure 7 for SOPHIS. This means that we have a smaller sensitivity to the value of L f .

In the supplementary material, a comparison is made between the two algorithms when a disturbance is present.

In terms of computational time, for a budget n = 1000, the execution time of one open-loop run is approximately 0.4s. If we increase the budget to n = 2000, the execution time becomes approximately 2.2s, while a budget n = 5000 implies 9.7s. The simulations have been run in Matlab, which is therefore not yet ready for real-time control; however, C would accelerate the algorithm, see [?] for an example of real-time control using SOPC implemented in C and for a larger discussion on the topic.

Conclusions

Two algorithms are proposed for systems with with both continuous and discrete actions: Optimistic Planning for Hybrid-Input Systems (OPHIS) and Simultaneous OPHIS. While OPHIS refines one set per iteration, using computed upper bounds for the returned values, SOPHIS simultaneously refines several sets per iteration. This eliminates the dependency of the Lipschitz constant in the set selection. Theoretical guarantees are given for both methods, where convergence rates are proven, depending on a measure of problem complexity.

In simulations, the algorithms are successful for a twolink robot arm problem, with SOPHIS proving to be the better choice when a larger budget is available.

In future work we aim to study stability properties [?,?] and perhaps even exploit stability in order to achieve tighter bounds.

Fig. 1 .Fig

 1 Fig. 1. Example of continuous split: top -states or rewards trajectories, middle -discrete actions, bottom -continuous intervals

Fig. 4 .

 4 Fig. 4. An example of the splitting function s(k)

Fig. 5 .

 5 Fig. 5. OPHIS: Evolution in time of the states and actions

Fig. 6 .

 6 Fig. 6. SOPHIS: Evolution in time of the states and actions

Fig. 7 .

 7 Fig. 7. Discounted returns of the two algorithms for several L f and n values

  discrete * /

	7:	create p + 1 children sets from i † ;
	8:	children sets inherit continuous intervals and
		discrete actions up to dimension D i † -1;
	9:	create one child set for each d -this action is
		added for dimension D i † ;
	10:	all children will have D = D i † + 1 and
	11:	C = C i † ; else/ * split continuous * /
	12:	expand set i † along k † by creating its M
		children sets;
	13:	children sets inherit continuous intervals and
	14:	discrete actions up to dimension D i † -1; interval at step k † is refined by splitting into
		M equal parts;
	15:

  ) of Section 4.1. A pseudocode of SOPHIS is given in Algorithm 4.2. SOPHIS 1: input state x 0 , model f , ρ, split factor M , discrete set {0, 1, ..., p}, budget n, Lipschitz constants L f and L ρ , discount factor γ, H max (n) 2: initialize collection of sets A with S 0 , D 0 = 0, C 0 = 0; 3: while budget still available do

	Algorithm 2 4: H = smallest depth with unexpanded nodes;
	5:	if H ≥ H max (n) then
	6:	stop and exit the loop;
	7:	else
	8:	while H < H max (n) do
	9:

  γ then

	12:	split discrete (see Algorithm 4.1);
	13:	else
	14:	split continuous (see Algorithm 4.1);
	15:	end if
	16:	H = H + 1;
	17:	end while
	18:

This work was financially supported from EU's H2020 Sea-Clear project No 871295; and by the project 38 PFE in the frame of the programme PDI-PFE-CDI 2021. The work of I.C. Morȃrescu and J. Daafouz was funded by the ANR under grants HANDY ANR-18-CE40-0010 and NICETWEET ANR-20-CE48-0009. Computation resources were provided by the CLOUDUT Project, cofunded by the European Fund of Regional Development through the Competitiveness Operational Programme 2014-2020, contract no. 235/2020.