Second order sliding mode twisting controller tuning based on two-level optimization process

Dominique Monnet 1 , Alexandre Goldsztejn 2 , Franck Plestan 3 and Charles Audet 1

Abstract-State-of-the-art finite time convergence conditions for the sliding mode controllers rely on bounds on perturbation terms. These bounds are often over-approximated, leading to conservative designs, i.e., high gains that amplify undesired behaviors such as chattering. This paper proposes to evaluate precisely the bounds on the perturbation terms to avoid conservative designs by using branch-and-bound algorithms dedicated to nonlinear programming. This leads to non-linear, a priori non-convex, non-differentiable constraints on the controller's gains, which is shown to be solvable using a modern blackbox optimization algorithm. We propose a new methodology employing branch-and-bound and blackbox solvers to generate gains as small as possible ensuring finite time convergence for the twisting algorithm. It is investigated using both a classical and a recently proposed sufficient conditions for finite time convergence. The applicability of the approach is illustrated over a numerical example.

I. INTRODUCTION

One motivation to the development of second order sliding mode (2-SM) techniques has been the reduction of the chattering effect [START_REF] Boiko | Analysis of sliding modes in the frequency domain[END_REF], [START_REF] Levant | Chattering analysis[END_REF], [START_REF] Fridman | Chattering analysis in sliding mode systems with inertial sensors[END_REF] inherent to the classical first order techniques [START_REF] Utkin | Sliding mode control in electromechanical systems[END_REF]. Several 2-SM techniques have been proposed [START_REF] Levant | Principles of 2-sliding mode design[END_REF]. The convergence proofs to the sliding surface rely on bounds on the perturbation terms appearing in the second time derivative of the sliding variable, whether these proofs rely on the majorant curve [START_REF] Shtessel | Sliding mode control and observation[END_REF] or Lyapunov theory [START_REF] Moreno | A lyapunov approach to second-order sliding mode controllers and observers[END_REF], [START_REF] Polyakov | Lyapunov function design for finite-time convergence analysis: "twisting" controller for second-order sliding mode realization[END_REF]. These bounds are generally difficult to obtain, and over-approximating them leads to conservative designs, i.e., large magnitude of the control input that amplifies the chattering effect. As a consequence, the control law is often tuned by simulation rather than derived from the convergence conditions. Although adaptive gain strategies allow to reduce the chattering, investigating sharper theoretical bounds and their usage through numerical optimization is of central importance.

In this paper, the classical 2-SM twisting controller is considered. The aim is to focus on the existing convergence conditions, and to use optimization tools to precisely evaluate the bounds they involve. Using the conditions to tune the twisting controller's gains, the convergence is guaranteed (contrary to a simulation-based tuning), and the convergence conditions are as less conservative as they can be. This approach follows previous works that studied optimization-based tuning for first order sliding mode [START_REF] Rosendo | A global optimization approach for sliding mode tuning and existence maps generation[END_REF], where SIP solvers were used to solve simple instances of first order sliding mode. Such solvers turn out to be unable to solve larger instances of second order SM. We propose a mixed local/global approach in this paper, where a branchand-bound algorithm dedicated to nonlinear optimization is used to solve globally the finite time convergence conditions, and a local blackbox algorithm [START_REF] Audet | Derivative-Free and Blackbox Optimization[END_REF] is used to optimize gains subject to the blackbox constraints.

The paper is organized as follow. Section II states the stateof-the-art convergence conditions of the twisting controller, as well as new less conservative conditions. Section III provides notions on global optimization tools, and how to employ them to compute the bounds involved in the convergence conditions. Section IV illustrates the proposed optimization-based approach on an example, and Section V concludes the paper and proposes directions for future work.

II. TWISTING CONVERGENCE CONDITIONS

Let us consider a dynamical system given by

ẋ = f (p(t), x, u) (1)
where t ≥ 0 is the time variable, p is a perturbation input, x ∈ R n the state, u ∈ R the control input. We assume that the sliding variable σ(p(t), x) depends on perturbations and state variables, as well as its first and second time derivative. We further assume that σ(p(t), x) is affine in the control input and is assumed to be linear in the control input,

σ = h(p(t), x) + g(p(t), x)u. (2)
p is therefore a vector that includes external input signals such as the reference signal, perturbations, etc. and their time derivatives. It is assumed that σ, σ, h and g have known explicit expressions depending on p and x. This assumption is needed to employ optimization tools. Several algorithms have been proposed to ensure the convergence of the system on the sliding surface [START_REF] Levant | Principles of 2-sliding mode design[END_REF]. This paper focuses on the twisting algorithm, whose control structure is given by ( [START_REF] Shtessel | Sliding mode control and observation[END_REF])

u = -r 1 sign(σ) -r 2 sign( σ). (3)
We assume in the following that finite bounds P are known on p(t), p(t) ∈ P, ∀t ≥ 0.

This realistic assumption is required to allow to employ numerical solvers to compute bounds on the right-hand side of Equation [START_REF] Boiko | Analysis of sliding modes in the frequency domain[END_REF]. These bounds are also needed to compute the gains (r 1 , r 2 ) using the convergence conditions stated in Subsection II-A.

Since the control signal is discontinuous and not defined for σ = 0 and σ = 0, the solution of (1) is understood in the Filippov sense [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF].

The following two subsections present conditions that guarantee the convergence of both σ and σ to zero in finite time, hence ensuring the expected sliding motion. The first conditions are classical, the second ones have been recently proposed by the authors in [START_REF] Monnet | Generalization of the twisting convergence conditions to non-affine systems[END_REF], where the conservatism is thoroughly discussed and illustrated. All convergence conditions assume the state is maintained in a given compact X ⊆ R n . This is another condition to employ optimization tools. Note that in the literature, this condition is implicit whenever bounds on g and h are assumed and they are not bounded over the state space.

A. State-of-the-art convergence conditions

The classical convergence conditions for the twisting controller rely on bounds C, k, K ≥ 0 on h(p(t), x) and g(p(t), x) that hold for all t ≥ 0 and x ∈ X. These bounds are defined such that

h(p(t), x) ∈ [-C, C], ∀t ≥ 0, ∀x ∈ X, g(p(t), x) ∈ [k, K], ∀t ≥ 0, ∀x ∈ X. (5)
The convergence conditions based on the bounds ( 5) are given by Theorem 1 [15, pp.148]

Theorem 1: Consider System (2) with u defined by (3). If r 1 > r 2 > 0, C < k(r 1 -r 2 ), (6)
and

(r 1 -r 2 )K + C < (r 1 + r 2 )k -C, (7)
then σ and σ converge to 0 in finite time.

The convergence conditions of Theorem 1 present the advantage to be easy to use in order to tune the gains r 1 and r 2 , since these conditions are expressed as linear inequalities on r 1 and r 2 .

B. Less conservative convergence conditions

Recent work [START_REF] Monnet | Generalization of the twisting convergence conditions to non-affine systems[END_REF] has proposed new convergence conditions for the twisting algorithm that are less conservative than the conditions of Theorem 1. That is, these new conditions define larger sets for the gains (r 1 , r 2 ) ensuring the finite time convergence. Therefore, one is allowed to choose smaller gains limiting chattering effect while ensuring finite time convergence. The new conditions depend on the bounds on σ in the four quadrants Σ i in the space (σ, σ) defined as,

Σ 1 = {(σ, σ) | σ > 0 and σ > 0}, Σ 2 = {(σ, σ) | σ > 0 and σ < 0}, Σ 3 = {(σ, σ) | σ < 0 and σ < 0}, Σ 4 = {(σ, σ) | σ < 0 and σ > 0}. ( 8
)
The sub-spaces in X corresponding to these quadrants are denoted X i , for i ∈ 1, 2, 3, 4,

X i (p(t)) = {x ∈ X | (σ(p(t), x), σ(p(t), x)) ∈ Σ i }. (9)
Notice that, since σ and σ depends on p(t), the sub-spaces X i also depend on p(t). The bounds on σ in the four quadrants thus satisfy,

σ(p(t), x) ∈ [m i , M i ], ∀t ≥ 0, ∀x ∈ X i (p(t)). (10)
Theorem 2 [START_REF] Monnet | Generalization of the twisting convergence conditions to non-affine systems[END_REF] states the finite time convergence based on the bounds [START_REF] Moreno | A lyapunov approach to second-order sliding mode controllers and observers[END_REF].

Theorem 2: Consider System (2) with u defined by (3). If the bounds (10) hold,

M 1 , M 2 < 0, (11)
m 3 , m 4 > 0, (12)
and

M 1 m 3 < m 2 M 4 , (13)
then σ and σ converge to zero in finite time.

Remark 1: The proof of Theorem 1 relies on the study of the majorant curve that is obtained from the bounds on σ [15, pp. 148]. With the bounds C, k, K on h and g, and the twisting controller (3), one derives the differential inclusion,

σ ∈ [-C + K(-r 1 -r 2 ), C + k(-r 1 -r 2 )] in Σ 1 , σ ∈ [-C + K(-r 1 + r 2 ), C + k(-r 1 + r 2 )] in Σ 2 , σ ∈ [-C + k(r 1 + r 2 ), C + K(r 1 + r 2 )] in Σ 3 , σ ∈ [-C + k(r 1 -r 2 ), C + K(r 1 -r 2 )] in Σ 4 . (14)
The conditions of Theorem 2 can be retrieved by following the same reasoning than the proof of Theorem 1 by replacing the bounds ( 14) by the generic bounds [START_REF] Moreno | A lyapunov approach to second-order sliding mode controllers and observers[END_REF]. One can remark that under the condition r 1 > r 2 > 0 necessary in Theorem 1, injecting the bounds [START_REF] Sahinidis | Baron: A general purpose global optimization software package[END_REF] in Conditions [START_REF] Polyakov | Lyapunov function design for finite-time convergence analysis: "twisting" controller for second-order sliding mode realization[END_REF] and [START_REF] Rosendo | A global optimization approach for sliding mode tuning and existence maps generation[END_REF] of Theorem 2 leads to Condition (6) of Theorem 1, with the additional constraint C < K(r 1 + r 2 ) which is ignored since it is necessarily satisfied if (6) holds. Condition [START_REF] Rosendo | A global optimization approach for sliding mode tuning and existence maps generation[END_REF] becomes Condition [START_REF] Levant | Principles of 2-sliding mode design[END_REF] in Theorem 1.

Remark 2: It is proved in [START_REF] Monnet | Generalization of the twisting convergence conditions to non-affine systems[END_REF] that the conditions of Theorem 2 are less conservative than the conditions of Theorem 1 provided that the same bounding process is used in both cases. This conservatism is illustrated in Section IV.

In the general case, it is not possible to obtain explicitly the dependancy with respect to r 1 and r 2 for the bounds [START_REF] Moreno | A lyapunov approach to second-order sliding mode controllers and observers[END_REF]. However, since σ has an explicit expression, it is possible to employ numerical solvers to compute such bounds in a guaranteed way.

III. NUMERICAL OPTIMIZATION FOR GAINS

COMPUTATION

This section proposes a two-level optimization process to compute the optimal gains r 1 , r 2 with respect to a given performance criterion, for example minimizing r 1 + r 2 . The lower level of the process consists in computing the bounds [START_REF] Moreno | A lyapunov approach to second-order sliding mode controllers and observers[END_REF] in a guaranteed way. Tightest bounds are aimed to be computed, and since h, g and σ are continuous nonlinear non-convex functions of the states and the perturbations, global solvers must be employed to obtain numerical guaranteed bounds. The upper level of the optimization process consists in exploring the space (r 1 , r 2 ) to find the best gains, that is, the gains that minimize a given performance criterion such that the bounds computed at the lower level respect the convergence conditions. In [START_REF] Rosendo | A global optimization approach for sliding mode tuning and existence maps generation[END_REF], a similar two-level approach for 1-SM tuning is proposed, the upper level being solve with a global solver. However, the underlying optimization problem of the upper level for tuning the twisting controller is more difficult since it involves a greater number of constraints to ensure the convergence, more complex expressions, and more variables. This makes a global approach intractable because of the high complexity of the global optimization (GO) algorithms. That is why we chose here a novel approach by employing an algorithm that converges locally and not globally, but has lower complexity and can be employed for tuning the twisting controller's gains, i.e., solving the upper problem. With

A. Overall optimal gain design

The gains r 1 and r 2 are chosen by optimizing a performance index f (r 1 , r 2 ) subject to the requirement that σ(p(t)) and σ(p(t)) converge to 0 in finite time. Having small gains is a typical requirement, e.g., f (r) = r 1 + r 2 . The finite time convergence requirement is enforced using the sufficient conditions presented in the previous section. Although Theorem 2 is less conservative than Theorem 1, its bounds do not depend linearly on the gains and require additional work 1 : on the one hand, once the bounds C, k, K are computed, Theorem 1 provides linear inequalities on the values of r 1 and r 2 that enforce finite time convergence. The performance index can be optimized inside these linear inequalities. On the other hand, the bounds m i and M i used in Theorem 2 depend non-linearly in r 1 and r 2 , and cannot be expressed explicitly as functions of r 1 and r 2 . It follows that the conditions of Theorem 2 cannot be expressed explicitly, and checking that a pair (r 1 , r 2 ) enforces finitetime convergence requires to employ a global solver to compute m i and M i . Then, it is possible to check if the convergence conditions of Theorem 2 are ensured. This corresponds to solving the lower problem mentioned above.

If one can solve the lower problem, it is possible to discretize the search for r 1 and r 2 , and for each sample check if the conditions of Theorem 2 are satisfied. This is done in Section IV but only to characterize the set of gains satisfying the conditions of Theorem 2 and to illustrate the conservatism of the conditions of Theorem 1 (see Figure 2).

In order to generate gains that are optimal with respect to f (r 1 , r 2 ), we propose to employ a Derivative Free Optimization (DFO) algorithm. The most common optimization algorithms relies on the derivative of f and the constraints given by the convergence conditions to compute (r 1 , r 2 ) that minimizes f . However, in our case, we do not know explicit expressions for m i and M i , as they are the result of the lower level optimization process, and therefore are not differentiable. That is why we need to employ a DFO algorithm that can solve this class of problems. The DFO algorithm chooses iteratively potential candidates in the space (r 1 , r 2 ), provides it to the lower level which computes the bounds m i , M i and returns them to the upper level as illustrated in Figure 1. From the values of f and m i , M i at the provided candidate, the DFO chooses a new potential best candidate at the next iteration and eventually converges to a locally optimal pair of gains.

The two-level approach we propose generates a pair (r 1 , r 2 ) which is guaranteed to ensure finite time convergence of σ and σ to zero. The pair is locally optimal with respect to f in the set of gains ensuring finite time convergence.

B. Lower level: computing optimal bounds with global optimization

The bounds ( 5) and ( 10) have the form b(t, x) ∈ [l, u], where b is alternatively g, h, or σ. The bounds l and u are to be computed so that b(p(t), x) ∈ [l, u] holds for all t ≥ 0 and all x ∈ X(p(t)). When b is either g or h in (5) then the set X(p(t)) does not depend on t, but we discuss here this more general case only. One can remark that, the tighter are the computed bounds the larger is set of gains proved to enforce finite time convergence by theorem. Therefore, we aim at computing the optimal bounds l * = min b(p(t), x).

(15)

Local minima or maxima cannot be used because inequalities would not be guaranteed to hold. Therefore, deterministic global optimization such as branch-and-bound algorithms are required to maintain a global lower bound (respectively an upper bound) that converges to the minimum (respectively to the maximum).

The dependency with respect to time in the optimization problems [START_REF] Shtessel | Sliding mode control and observation[END_REF] causes difficulties. The first one is that explicit expressions with respect to time of perturbations are unavailable, but only bounds on perturbations are known. The second one is that t ∈ [0, ∞) is unbounded while branch-andbound algorithms require bounded feasible search spaces. As explained in Section I, we assume that all time dependencies arise through some unknown function p(t) bounded in P holds for all t ≥ 0. Theorem 3 allows solving problems involving a variable p ∈ P instead of t ≥ 0.

Theorem 3: Define

l + = min p∈P x∈X(p) b(p, x), u + = max p∈P x∈X(p) b(p, x). (16)
Then

[l + , u + ] ⊇ [l * , u * ].
Proof: Let t * ≥ 0 x * ∈ X(t * ) be minimizers of ( 15), i.e., b(p(t * ), x * ) ≤ 0. Then p * ∈ P , and x * ∈ X(p(t * )) = X(p * ) therefore p * and x * are feasible for the problem [START_REF] Utkin | Sliding mode control in electromechanical systems[END_REF], which implies l + ≤ l * . A symmetric reasoning shows that u + ≥ u * .

Note that the bounds [l + , u + ] are conservative with respect to [l * , u * ] because the impact of time correlation on b(p(t)) is lost when only considering the information p(t) ∈ P . One cannot expect better if no accurate knowledge on p(t) is known.

Off-the-shelf global solvers are available [START_REF] Sahinidis | Baron: A general purpose global optimization software package[END_REF], [START_REF] Chabert | Contractor programming[END_REF], and discussing in detail how they are implemented is beyond the scope of this paper.

In the end, global solvers enable to compute guaranteed lower (resp. upper) bounds on l + (resp. u + ) as defined in Theorem 3 within a given relative precision, and by extension guaranteed bounds on the optimal bounds l * and u * as stated by Theorem 3. One can obtain the bounds C, k, K, and from the conditions of Theorem 1 define a set of gains (r 1 ,r 2 ) by linear inequalities. Furthermore, provided numerical value for (r 1 , r 2 ), one can compute guaranteed bounds m i , M i on σ, and thus concluding if the conditions of Theorem 2 are satisfied.

C. Upper level: optimizing the gains

The upper problem consist in finding r 1 , r 2 that minimize a performance criterion f (r 1 , r 2 ), such that (r 1 , r 2 ) satisfies the convergence constraints of Theorem 2. We recall that the convergence constraints are nonlinear, non convex, and non-differentiable and do not have explicit expressions. Such problems can be solved by DFO algorithms, that do not need to have access to the derivative of f or the constraints to converge to a local solution. In this paper, we propose to use the NOMAD solver [START_REF] Digabel | Algorithm 909: Nomad: Nonlinear optimization with the mads algorithm[END_REF] to solve the upper level problem. NOMAD has been developed to solve blackbox (BB) problems, meaning that the evaluation of the performance criterion f and the constraints is the result of some BB time-consuming process, for example simulations [?]. In our case, the BB process is the lower problem, which takes as input a pair (r 1 , r 2 ) and returns the value of the bounds by employing a GO solver as explained in the previous subsection. From an initial pair (r 0 1 , r 0 2 ), NOMAD explore the space of gains by repeatedly choosing points in that space and calling the GO solver to ensure that this points respect the convergence constraint, and eventually converges to a local optimum (r * 1 , r * 2 ). Figure 1 schematizes the two-levels optimization process.

NOMAD implements a Mesh Adaptive Direct Search (MADS) algorithm. Devising on the principle of this kind of algorithm is beyond the scope of this paper, but the interested reader can refer to [?] for the latest development in MADS.

IV. NUMERICAL EXAMPLE

Consider the pendulum example adapted from [START_REF] Levant | Principles of 2-sliding mode design[END_REF]. The system is a rigid pendulum loaded with a mass m tracking a reference signal θ c (t), with a varying arm length R(t) to model a perturbation. The states of the pendulum are x = (θ, θ), the angle of the arm and its angular velocity. The state θ = 0 corresponds to the vertical position (down). We suppose that an actuator provides a torque u considered as the control input. The system dynamic is given by the Input: (r 0 1 , r 0 2 )

Output:

(r * 1 , r * 2 )
Upper Level solver: NOMAD min r1,r2

f (r 1 , r 2 ) M 1 (r 1 , r 2 ) < 0 M 2 (r 1 , r 2 ) < 0 m 3 (r 1 , r 2 ) > 0 m 4 (r 1 , r 2 ) > 0 M 1 m 3 -M 4 m 2 < 0
Lower Level (Black Box) slover: Branch-and-Bound

m i = min p∈P x∈Xi(p) σ(r 1 , r 2 , p, x) M i = max p∈P x∈Xi(p) σ(r 1 , r 2 , p, x) (r 1 , r 2 ) f (r 1 , r 2 ) M i , m i Fig. 1.
Two-levels optimization process for optimizing the twisting controller gains. differential equation

θ = -2 Ṙ R θ -g 1 R sin(θ) + 1 mR 2 u. (17)
We take the mass m = 1 and g = 9.81.

The reference signal is given by θ c (t) = 0.5 sin(0.5t) + 0.5 cos(0.5t),

and the perturbation is modeled by R(t) = 0.8 + 0.1 sin(2t) + 0.2 cos(t).

The sliding variable is

σ = θ -θ c . (20)
It follows that σ = θθc

and σ = h(t, x) + g(t, x)u

with

g(p, x) = 1 mR 2 , h(p, x) = -2 Ṙ R θ -9.81 1 R sin(θ) -θc . ( 23
)
The vector of perturbations is

p(t) = (R(t), Ṙ(t), θ c (t), θc (t), θc (t)). ( 24
)
From the expressions of R(t) and θ c (t), we derive the bounds

p(t) ∈ P = [0.5, 1.1] × [-0.4, 0.4] × [-1, 1] ×[-0.5, 0.5] × [-0.25, 0.25]. (25)
We propose now to generate the set of controllers satisfying the convergence conditions of Theorem 1 and 2. We limit the study to x ∈ [-π, π] × [-4, 4].

In the following, all the problems are solved with the global solver implemented in the C++ library IBEX [START_REF] Chabert | Contractor programming[END_REF] with a relative precision ϵ = 0.01. We obtain the bounds C = 26.65, k = 0.82 and K = 4.04. The resulting set of feasible values of (r 1 , r 2 ) with respect to the linear constraints of Theorem 1 is displayed by the gray set in Figure 2.

We propose to generate a set of feasible gains from the conditions of Theorem 2 for illustrative comparison with the conditions of Theorem 1. To do so, we discretize the search space (r 1 , r 2 ), and for every point we compute the optimal bounds m i , M i , as explained in Subsection III-B and verify if the conditions of Theorem 2 are satisfied. The discretization points are r 1 = {1, 3, 5, ...199}, r 2 = {1, 3, 5, ...199} and generate a total number of 10,000 points. The points that satisfy the convergence conditions of Theorem 2 are represented by the dots in Figure 2. Processing these 10,000 points takes approximately 70 minutes on a standard laptop. In Figure 2, the set of feasible gains described by the conditions of Theorem 2 is larger than the one described by the conditions of the classical conditions provided by Theorem 1. However, a grid search procedure is very inefficient to locate the point that minimizes r 1 + r 2 .

We now focus on reducing the chattering effect by employing the two-level optimization process described in Section III to compute a pair of gains (r 1 , r 2 ) that minimizes the function f (r 1 , r 2 ) = r 1 + r 2 , corresponding to the maximum amplitude of the control input. We select as initial point (r 0 1 , r 0 2 ) for the NOMAD solver the point that minimizes f and enforces the convergence condition of Theorem 1. This point (r 0 1 , r 0 2 ) = (129, 96.5) corresponds to the leftmost point of the gray triangle in Figure 2. We run NOMAD with the default parameters and set the maximum number of blackbox evaluations, i.e., solving the lower level problem, to 500. In 8 minutes, NOMAD produces the point (47.1, 36.6) represented by the star in Figure 2.

The simulations of (17) with (3) tuned with these pairs of gains are performed using an explicit Euler method with 10 -5 s time step. The control input is updated every 10 -2 s to simulate the sampling time of a real system. The simulation results are displayed in Figure 3. The simulations show that tuning the controller's gains using the new conditions decrease the chattering effect. Between 2 and 10 seconds (sliding motion set), the L2 norm of the error (x c -x) is decreases 59.99 % using the gains obtained with Theorem 2.

To conclude, using the new conditions, although computationally expensive, provides a larger set of gains than the classical conditions. Having a larger set provides more freedom in choosing the controller's gains, that can help reducing the chattering effect while guaranteeing the convergence to the sliding surface.

V. CONCLUSION This paper showed how numerical optimization can be employed for tuning the twisting controller's gains. Using a global solver enabled to compute tight bounds on which are based the convergence conditions, and therefore helped in reducing the conservatism due to over-approximations of the bounds. Employing global optimization was necessary since the bounds are found by solving programs involving non-linear functions. Blackbox optimization has been used to solve the higher optimization problem consisting in finding optimal gains subject to finite time convergence. The study of a numerical example showed that using such optimizationbased methods with the new convergence conditions of Theorem 2, despite being computationally expensive, enables more freedom in choosing controller's gain while preserving the finite-time convergence.

Future works include finding even less conservative conditions than the ones of Theorem 2 for the twisting algorithm, and extending this approach to other 2-SM controllers.

Fig. 2 .

 2 Fig. 2.Sets of gains (r 1 , r 2 ) satisfying the convergence conditions of Theorem 1 (classical conditions) and the ones of Theorem 2 (less conservative) tested over a grid of 10,000 points. The star represents the solution obtained with the proposed Black-Box optimization process.

This additional work is worthwhile in an optimal design process.

ACKNOWLEDGMENT This work is supported by Région Pays de la Loire, Centrale Nantes, D-ICE engineering, West Atlantic Marine Energy Community (WEAMEC), Atlanstic 2020, LS2N and CNRS as part of the O2GRACE project. Part of this work is supported by the NSERC Alliance grant 544900-19 in collaboration with Huawei-Canada. The authors thank Sébastien Le Digabel for his help with the NOMAD solver.


Ai-Hua Fan
email: ai-hua.fan@u-picaride.fr

AND Ling-Min Liao
email: lingmin.liao@u-picardie.fr

Ji-Hua Ma
email: jhma@whu.edu.cn

ON THE FREQUENCY OF PARTIAL QUOTIENTS OF REGULAR CONTINUED FRACTIONS

We consider sets of real numbers in [0, 1) with prescribed frequencies of partial quotients in their regular continued fraction expansions. It is shown that the Hausdorff dimensions of these sets, always bounded from below by 1/2, are given by a modified variational principle.

Introduction

Let Q c denote the set of irrational number. It is well-known that each x ∈ [0, 1) ∩ Q c possesses a unique continued fraction expansion of the form

x = 1 a 1 (x) + 1 a 2 (x) + 1 a 3 (x) + . . . , (1.1)
where a k (x) ∈ N := {1, 2, 3, • • • } is the k-th partial quotient of x. This expansion is usually denoted by x = [a 1 (x), a 2 (x), • • • ]. For each j ∈ N, define the frequency of the digit j in the continued fraction expansion of x by τ j (x) := lim n→∞ τ j (x, n) n , whenever the limit exists, where

τ j (x, n) := Card{1 ≤ k ≤ n : a k (x) = j}.
This paper is concerned with sets of real numbers with prescribed digit frequencies in their continued fraction expansions. To be precise, let p = (p 1 , p 2 , . . . ) be a probability vector with p j ≥ 0 for all j ∈ N and ∞ j=1 p j = 1, which will be called a frequency vector in the sequel. Our purpose is to determine the Hausdorff dimension of the set

E p := {x ∈ [0, 1) ∩ Q c : τ j (x) = p j ∀j ≥ 1}.
Let us first recall some notation. For any a 1 , a 2 , • • • , a n ∈ N, we call

I(a 1 , a 2 , • • • , a n ) := {x ∈ [0, 1) : a 1 (x) = a 1 , a 2 (x) = a 2 , • • • , a n (x) = a n }
a rank-n basic interval. Let T : [0, 1) → [0, 1) be the Gauss transformation defined by T (0) = 0, T (x) = 1/x (mod 1) for x ∈ (0, 1). For a given frequency vector p = (p 1 , p 2 , . . . ), we denote by N ( p) the set of Tinvariant ergodic probability measures µ such that | log x|dµ < ∞ and µ(I(j)) = p j for all j ≥ 1.

(1.2)

Let h µ stand for the measure-theoretical entropy of µ, and dim H A for the Hausdorff dimension of a set A. The main result of this paper can be stated as follows.

Theorem 1.1. For any frequency vector p, one has

dim H (E p ) = max 1 2 , sup µ∈N ( p) h µ 2 | log x|dµ ,
where the "sup" is set to be zero if N ( p) = ∅. By virtue of log |T ′ (x)| = 2| log x|, we see that 2 | log x|dµ is the Liapunov exponent of the measure µ. Therefore, the "sup" in the above is a variational formula which relates the Hausdorff dimension to the entropy and Liapunov exponent of measures.

Theorem 1.1 provides a complete solution to the long standing problem that requests an exact formula for dim H (E p ). Let us recall some partial results in the literature. In 1966, under the condition that ∞ j=1 p j log j < ∞, Kinney and Pitcher [START_REF] Kinney | The dimension of some sets defined in terms of f -expansions[END_REF] showed that

dim H (E p ) ≥ -∞ j=1 p j log p j 2 | log x|dµ p ,
where µ p is the Bernoulli measure on [0, 1] defined by

µ p (I(a 1 , a 2 , • • • , a n )) = n j=1 p aj .
The above lower bound is just the Hausdorff dimension of the Bernoulli measure µ p . However, by a result of Kifer, Peres and Weiss in 2001, this is not an optimal lower bound. Indeed, it is shown in [START_REF] Kifer | A dimension gap for continued fractions with independent digits[END_REF] that, for any Bernoulli measure µ p ,

dim H µ p ≤ 1 -10 -7
This surprising fact indicates that the collection of Bernoulli measures are insufficient for providing the correct lower bound for dim H (E p ).

In 1975, under the same condition ∞ j=1 p j log j < ∞, Billingsley and Henningsen [START_REF] Billingsley | Hausdorff dimension of some continued-fraction sets[END_REF] obtained an improved lower bound

dim H (E p ) ≥ sup µ∈N ( p) h µ 2 | log x|dµ .
Moreover, they proved that, for any fixed N ∈ N, this lower bound is the exact Hausdorff dimension of the set

{x ∈ E p : a n (x) ≤ N for all n ≥ 1}
provided that p j = 0 for all j > N . It is therefore quite natural to guess that this lower bound is the right value for dim H E p in general. However, as will be shown in Theorem 1.1, this is not the case. Actually, the lower bound due to Billingsley and Henningsen is only a half of the correct lower bound. The other half of the lower bound, namely, dim H (E p ) ≥ 1/2, can be proved basing on Lemma 2.4 in [START_REF] Liao | Dimension of some non-normal continued fraction sets[END_REF]. However, we will give a direct proof in this paper.

The upper bound estimate is more difficult. In its proof, we will use techniques from [START_REF] Liao | Dimension of some non-normal continued fraction sets[END_REF] and [START_REF] Billingsley | Hausdorff dimension of some continued-fraction sets[END_REF] to estimate the lengths of basic intervals. Not incidentally, an entropy-involved combinatorial lemma borrowed from [START_REF] Glasner | On the interplay between measurable and topological dynamics[END_REF] (see Lemma 2.7) will play an important role.

The paper is organized as follows. In Section 2, we give some preliminaries on the basic intervals and on the entropy of finite words. In Section 3, we establish the upper bound in Theorem 1.1. In Section 4, we prove that dim H (E p ) ≥ 1/2 and show that we can drop the condition ∞ j=1 p j log j < ∞ in Billingsley and Henningsen's theorem and then obtain the lower bound in Theorem 1.1. The last section serves as a remark.

Preliminary

Let x = [a 1 (x), a 2 (x), • • • ] ∈ [0, 1) ∩ Q c .
The n-th convergent in the continued fraction expansion of x is defined by

p n q n := p n (a 1 (x), • • • , a n (x)) q n (a 1 (x), • • • , a n (x)) = 1 a 1 (x) + 1 a 2 (x) + . . . + 1 a n (x)
.

For ease of notation, we shall drop the argument x in what follows. It is known (see [START_REF] Ya | Continued Fractions[END_REF] p.9) that p n , q n can be obtained by the recursive relations:

p -1 = 1, p 0 = 0, p n = a n p n-1 + p n-2 (n ≥ 2), q -1 = 0, q 0 = 1, q n = a n q n-1 + q n-2 (n ≥ 2).
By the above recursion relations, we have the following results.

Lemma 2.1 ([6]). Let q n = q n (a 1 , • • • , a n ) and p n = p n (a 1 , • • • , a n ), we have (i) p n-1 q n -p n q n-1 = (-1) n ; (ii) q n ≥ 2 n-1 2 , n k=1 a k ≤ q n ≤ n k=1 (a k + 1). Lemma 2.2 ([14]). For any a 1 , a 2 , • • • , a n , b ∈ N, b + 1 2 ≤ q n+1 (a 1 , • • • , a j , b, a j+1 , • • • , a n ) q n (a 1 , • • • , a j , a j+1 , • • • , a n ) ≤ b + 1 (∀1 ≤ j < n).
Recall that for any a 1 , a 2 , • • • , a n ∈ N, the set

I(a 1 , a 2 , • • • , a n ) = {x ∈ [0, 1) : a 1 (x) = a 1 , a 2 (x) = a 2 , • • • , a n (x) = a n }
is a rank-n basic interval. We write |I| for the length of an interval I.

Lemma 2.3 ([10] p.18). The basic interval I(a 1 , a 2 , • • • , a n ) is an interval with endpoints p n /q n and (p n + p n-1 )/(q n + q n-1 ). Consequently, one has

I(a 1 , • • • , a n ) = 1 q n (q n + q n-1 ) , (2.1
)

and 1 2q 2 n ≤ I(a 1 , • • • , a n ) ≤ 1 q 2 n . (2.2)
Lemma 2.4. We have

|I(x 1 , • • • , j, • • • , x n )| ≤ 8 (j + 1) 2 I(x 1 , • • • , j, • • • , x n ) ,
where the notation j means "deleting the digit j".

Proof. By Lemma 2.2 and (2.2), we have

|I(x 1 , • • • , j, • • • , x n )| ≤ q -2 n (x 1 , • • • , j, • • • , x n ) ≤ 2 j + 1 2 q -2 n-1 (x 1 , • • • , ĵ, • • • , x n ) ≤ 8 (j + 1) 2 I(x 1 , • • • , j, • • • , x n ) .
We will simply denote by I n (x) the rank n basic interval containing x. Suppose that a n := a n (x) ≥ 2 and consider

I ′ n (x) = I(a 1 , • • • , a n-1 , a n -1) and I ′′ n (x) = I(a 1 , • • • , a n-1 , a n + 1
) which are two rank n basic intervals adjacent to I n (x). By the recursive equation of q n and (2.1), one has the following lemma.

Lemma 2.5. Suppose that a n := a n (x) ≥ 2. Then the lengths of the adjacent intervals I ′ n (x) and I ′′ n (x) are bounded by 1 3 I n (x) from below and by 3 I n (x) from above.

For any x ∈ [0, 1) \ Q and any word i

1 • • • i k ∈ N k , (k ≥ 1), denote by τ i1•••i k (x, n) the number of j, 1 ≤ j ≤ n, for which a j (x) • • • a j+k-1 (x) = i 1 • • • i k .
For N ∈ N, define Σ N := {1, . . . , N }. We shall use the following estimate in [START_REF] Billingsley | Hausdorff dimension of some continued-fraction sets[END_REF].

Lemma 2.6 ([2]). Let N ≥ 1 and n ≥ 1. For any x = [x 1 , x 2 , • • • ] ∈ [0, 1] ∩ Q c with x j ∈ Σ N for 1 ≤ j ≤ n. Then for any k ≥ 1, we have log |I n (x)| ≤ 2 i1•••i k ∈Σ k N τ i1•••i k (x, n) log p k (i 1 , • • • , i k ) q k (i 1 , • • • , i k ) + 8 + 8n 2 k . (2.3)
Now we turn to the key combinatorial lemma which will be used in the upper bound estimation. Let φ : [0, 1] → R denote the function φ(0) = 0, and φ(t) = -t log t for 0 < t ≤ 1.

For every word ω ∈ Σ n N of length n and every word u ∈ Σ k N of length k, denote by p(u|ω) the frequency of appearances of u in ω, i.e.,

p(u|ω) = τ u (ω) n -k + 1 ,
where τ u (ω) denote the number of j, 1 ≤ j ≤ nk + 1, for which

ω j • • • ω j+k-1 = u. Define H k (ω) := u∈Σ k N φ(p(u|ω)).
We have the following counting lemma.

Lemma 2.7 ([5]).

For any h > 0, ǫ > 0, k ∈ N, and for any n ∈ N large enough, we have

Card{ω ∈ Σ n N : H k (ω) ≤ kh} ≤ exp(n(h + ǫ)).
3. Upper bound

3.1. Some Lemmas. Let (p(i 1 , • • • , i k )) i1•••i k ∈N k be a probability vector indexed by N k .
As usual, we denote by q k (a 1 , • • • , a k ) the denominator of the k-th convergent of a real number with leading continued fraction digits a 1 , . . . , a k .

Lemma 3.1. For each k ∈ N and each probability vector

(p(i 1 , • • • , i k )) i1•••i k ∈N k , i1•••i k ∈N k -p(i 1 , • • • , i k ) log p(i 1 , • • • , i k ) ≤ i1•••i k ∈N k -p(i 1 , • • • , i k ) log |I(i 1 , • • • , i k )|.
Proof. Applying Jesen's inequality to the concave function log, we have

i1•••i k ∈N k p(i 1 , • • • , i k ) log |I(i 1 , • • • , i k )| p(i 1 , • • • , i k ) ≤ log i1•••i k ∈N k |I(i 1 , • • • , i k )| = 0.
Lemma 3.2. Let p = (p 1 , p 2 , . . . ) be a probability vector and q = (q 1 , q 2 , . . . ) a positive vector. Suppose -∞ j=1 p j log q j = ∞ and ∞ j=1 q s j < ∞ for some positive number s. Then

lim sup n→∞ - n j=1 p j log p j -n j=1 p j log q j ≤ s.
Proof. This is a consequence of the following inequality (see [START_REF] Walters | An introduction to ergodic theory[END_REF], p.217): for nonnegative numbers s j (1 ≤ j ≤ m) such that m j=1 s j = 1 and any real numbers

t j (1 ≤ j ≤ m), we have m j=1 s j (t j -log s j ) ≤ log( m j=1 e tj ). (3.1)
Fix n ≥ 1. Let s j = p j for 1 ≤ j ≤ n and s n+1 = ∞ j=n+1 p j . Let t j = s log q j for 1 ≤ j ≤ n and t n+1 = 0. Applying the above inequality (3.1) with m = n + 1, we get

s n j=1 p j log q j - n j=1 p j log p j -( ∞ j=n+1 p j ) log( ∞ j=n+1 p j ) ≤ log(1 + n j=1 q s j ).
Consequently, -

n j=1 p j log p j -n j=1 p j log q j ≤ s + ( ∞ j=n+1 p j ) log( ∞ j=n+1 p j ) -n j=1 p j log q j + log(1 + n j=1 q s j ) -n j=1 p j log q j .
Using the facts -∞ j=1 p j log q j = ∞ and ∞ j=1 q s j < ∞, we finish the proof by letting n → ∞.

Lemma 3.2 implies the following lemma. Recall that Σ k N = {1, . . . , N } k . Lemma 3.3. Let k ≥ 1. If (p(i 1 , • • • , i k )) i1•••i k ∈N k is a probability vector such that i1•••i k ∈N k p(i 1 , • • • , i k ) log q k (i 1 , • • • , i k ) = ∞, then we have lim sup N →∞ -i1•••i k ∈Σ k N p(i 1 , • • • , i k ) log p(i 1 , • • • , i k ) 2 i1•••i k ∈Σ k N p(i 1 , • • • , i k ) log q k (i 1 , • • • , i k ) ≤ 1 2 .
Proof. By Lemma 2.1, for any k ∈ N and any s > 1/2, we have

i1,••• ,i k q k (i 1 , • • • , i k ) -2s ≤ i1,••• ,i k (i 1 • • • i k ) -2s = ( ∞ j=1 j -2s ) k < ∞.
Thus we get the result by Lemma 3.2.

3.2. Proof of the upper bound. To prove the upper bound, we shall make use of multi-step Markov measures. Let k ≥ 1, by a (k -1)-step Markov measure, we mean a T -invariant probability measure P on [0, 1) satisfying the Markov property

P (I(a 1 , • • • , a n )) P (I(a 1 , • • • , a n-1 )) = P (I(a n-k , • • • , a n )) P (I(a n-k , • • • , a n-1 )) (3.2)
for all n ≥ 1 and a 1 , • • • , a n ∈ N (see [START_REF] Cajar | Billingsley dimension in probability spaces[END_REF], p.9). We may regard a Bernoulli measure as a 0-step Markov measure.

For each N ≥ 2, we denote by P k N = P k N ( p) the collection of (k -1)-step Markov measures satisfying the condition

P (I(j)) = p j for 1 ≤ j ≤ N -1 and P (I(N )) = 1 - N -1 j=1 p j . (3.3)
These Markov measures are supported by the set of continued fractions for which the partial quotients are bounded from above by N .

For each

i 1 • • • i k ∈ {1, 2, • • • , N } k , write p(i 1 , . . . , i k ) = P (I(i 1 , . . . , i k )). Put α N,k := sup P ∈P k N -1 k p(i 1 , • • • , i k ) log p(i 1 , • • • , i k ) -2 p(i 1 , • • • , i k ) log(p k (i 1 , • • • , i k )/q k (i 1 , • • • , i k )) . (3.4)
The argument in [START_REF] Billingsley | Hausdorff dimension of some continued-fraction sets[END_REF] (pp.171-172) shows that the following limit

α N := lim k→∞ α N,k ,
exists and coincides with each of the following three limits:

α ′ N := lim k→∞ sup P ∈P k N -p(i 1 , • • • , i k ) log p(i 1 , • • • , i k ) 2 p(i 1 , • • • , i k ) log q k (i 1 , • • • , i k ) , α ′′ N := lim k→∞ sup P ∈P k N -p(i 1 , • • • , i k ) log p(i 1 , • • • , i k ) -p(i 1 , • • • , i k ) log |I(i 1 , • • • , i k )| , and
α ′′′ N := lim k→∞ sup P ∈P k N h P 2 | log x|dP . Let α := lim sup N →∞ α N = lim sup N →∞ α ′ N = lim sup N →∞ α ′′ N = lim sup N →∞ α ′′′ N . (3.5)
To prove the upper bound, we need only to prove the following two propositions.

Proposition 3.4. For any N ∈ N large enough, we have

dim H (E p ) ≤ max 1 2 , α N .
Proposition 3.5. We have

α ≤ max 1 2 , sup µ∈N ( p) h µ 2 | log x|dµ .
Remark that we will finally establish the formula in Theorem 1.1, thus by Proposition 3.4 and Proposition 3.5 we have max 1 2 , lim inf

N →∞ α N = max 1 2 , lim sup N →∞ α N = max 1 2 , sup µ∈N ( p) h µ 2 | log x|dµ ,
and if lim inf N →∞ α N ≥ 1/2, then the limit of α N exists and equals to sup

µ∈N ( p) h µ 2 | log x|dµ .
Proof of Proposition 3.5. By virtue of (3.5), we set Without loss of generality, we may suppose that there is a sequence of Markov measures P N,k ∈ P k N converging to a measure µ ∈ N ( p) in the weak * -topology such that lim sup

α = lim sup N →∞ lim k→∞ sup P ∈P k N -p(i 1 , • • • , i k ) log p(i 1 , • • • , i k ) 2 p(i 1 , • • • , i k ) log q k (i 1 , • • • , i k ) . Denote D := lim sup N →∞ lim k→∞ sup P ∈P k N 2 p(i 1 , • • • , i k ) log q k (i 1 , • • • , i k ). If D = ∞,
N →∞ lim k→∞ h P N,k 2 | log x|dP N,k = α.
Then by the upper semi-continuity of the entropy function and the weak convergence, we have

α ≤ sup µ∈N ( p) h µ 2 | log x|dµ .
Proof of Proposition 3.4. For any fixed integer N which is large enough, and any ǫ > 0, we have

E p ⊂ ∞ ℓ=1 ∞ n=ℓ H n (ǫ, N ),
where

H n (ǫ, N ) := x ∈ [0, 1) \ Q : τ j (x, n) n -p j < ǫ, 1 ≤ j ≤ N .
For any γ > max {1/2, α N } and for any integer k ∈ N, we have the γ-Hausdorff measure (see [START_REF] Falconer | Fractal Geometry, Mathematical Foundations and Application[END_REF], for the definition of H γ )

H γ ∞ n=ℓ H n (ǫ, N ) ≤ | τ j (x,n) n -pj |<ǫ,1≤j≤N |I n (x)| γ (∀n ≥ ℓ) = n(pj -ǫ)<mj<n(pj +ǫ),1≤j≤N x1•••xn∈A |I(x 1 , • • • , x n )| γ ,
where

A := {x 1 • • • x n ∈ Σ n N : τ j (x 1 • • • x n ) = m j , 1 ≤ j ≤ N } . (We recall that τ j (x 1 • • • x n ) denotes the times of appearances of j in x 1 • • • x n .)
Let ñ := N j=1 m j . By Lemma 2.4, we have the following estimate by deleting the digits j > N in the first n partial quotients x 1 , . . . , x n of x ∈ I n (x):

x1•••xn∈A |I(x 1 , • • • , x n )| γ ≤   ∞ j=N +1 8 (j + 1) 2γ   n-ñ x1•••xñ∈ à |I(x 1 , • • • , x ñ)| γ , where à := x 1 • • • x ñ ∈ Σ ñ N : τ j (x 1 • • • x ñ) = m j , 1 ≤ j ≤ N . Since γ > 1/2, for N large enough ∞ j=N +1 8 (j + 1) 2γ < 1. (3.6)
By applying Lemma 2.6, and noticing that

τ i1•••i k (x, ñ) ≤ τ i1•••i k (x 1 • • • x ñ) + k, we have |I(x 1 , • • • , x ñ)| = exp{log |I(x 1 , • • • , x ñ)|} ≤ exp    2 i1•••i k ∈Σ k N (τ i1•••i k (x 1 • • • x ñ) + k) log p k (i 1 , • • • , i k ) q k (i 1 , • • • , i k ) + 8 + 8ñ 2 k    . Thus x1•••xñ∈ Ã |I(x 1 , • • • , x ñ)| γ ≤ mi 1 •••i k x1•••xñ∈B exp    2γ i1•••i k ∈Σ k N (m i1•••i k + k) log p k q k + 8γ + 8ñγ 2 k    ,
where

B := x 1 • • • x ñ ∈ Ã : τ i1•••i k (x 1 • • • x ñ) = m i1•••i k ∀i 1 • • • i k ∈ Σ k N . Take h = 1 k i1•••i k ∈Σ k N φ m i1•••i k ñ -k + 1 (3.7)
in Lemma 2.7. We have for any δ > 0 and for ñ large enough

x1•••xñ∈B exp    2γ i1•••i k ∈Σ k N (m i1•••i k + k) log p k q k + 8γ + 8ñγ 2 k    ≤ exp    ñ (h + δ) + 2γ i1•••i k ∈Σ k N (m i1•••i k + k) log p k q k + 8γ + 8ñγ 2 k    .
Rewrite the right side of the above inequality as

exp {ñ (L(γ, k, m i1•••i k ))} ,
where

L(γ, k, m i1•••i k ) := h + 2γ i1•••i k ∈Σ k N m i1•••i k + k ñ log p k q k + 8γ ñ + 8γ 2 k + δ.
Since there are at most (ñ

-k + 1) N k possible words of i 1 • • • i k in Σ ñ N , we have x1•••xñ∈ Ã |I(x 1 , • • • , x ñ)| γ ≤ (ñ -k + 1) N k exp ñ sup mi 1 •••i k L(γ, k, m i1•••i k ) .
Notice that by the definition of à and B, the possible values of m i1•••i k are restricted to satisfy the condition that the frequency of digit j in

x 1 • • • x ñ is about p j . Then when ñ → ∞, for i 1 , • • • , i k ∈ Σ k N m i1•••i k ñ -k + 1 → p(i 1 , • • • , i k ), (3.8)
and

{p(i 1 , • • • , i k ) : i 1 , • • • , i k ∈ Σ k N }
defines a probability measure P in P k N . Now take δ > 0 small enough and k large enough such that

8γ 2 k < δ, (3.9) and γ
> -1 k p(i 1 , • • • , i k ) log p(i 1 , • • • , i k ) + 5δ -2 p(i 1 , • • • , i k ) log(p k (i 1 , • • • , i k )/q k (i 1 , • • • , i k )) . (3.10)
The last inequality comes from the definition of α N , (3.4) and the assumption γ > α N . By (3.8), for sufficiently large ñ, we have 8γ ñ < δ and

1 k i1•••i k ∈Σ k N φ m i1•••i k ñ -k -1 - 1 k i1•••i k ∈Σ k N φ (p(i 1 • • • i k )) < δ, i1•••i k ∈Σ k N m i1•••i k + k ñ log p k q k - i1•••i k ∈Σ k N p(i 1 • • • i k ) log p k q k < δ. (3.11)
By (3.7) and (3.11),

L(γ, k, m i1•••i k ) < 1 k i1•••i k ∈Σ k N φ (p(i 1 • • • i k )) + 2γ i1•••i k ∈Σ k N p(i 1 • • • i k ) log p k q k + 5δ. Thus (3.10) implies L(γ, k, m i1•••i k ) < 0.
Hence finally, by (3.6), we can obtain for any γ > max {1/2, α N },

H γ ∞ n=k H n (ǫ, N ) < ∞.
This implies that dim H (E p ) ≤ max {1/2, α N } as desired.

Lower bound

In this section, we first prove dim H (E p ) ≥ 1/2. Then we examine what happens if the condition ∞ j=1 p j log j < ∞ in Billingsley and Henningsen's theorem is violated. We will see that if the condition is not satisfied, then dim H (E p ) = 1/2 and

N (p) = ∅.
The following is the key lemma for proving that dim H (E p ) ≥ 1/2.

Lemma 4.1. For any given sequence of positive integers {c n } n≥1 tending to the infinity, there exists a point z = (z 1 , z 2 , . . . ) ∈ E p such that z n ≤ c n for all n ≥ 1.

Proof. For any n ≥ 1, we construct a probability vector (p

k , . . . ) such that p For each digit k ≥ 1, consider the random variables of x ∈ N N , X i (x) = 1 {k} (x i ), (i ≥ 1). By Kolmogorov's strong law of large numbers (see [START_REF] Shiryaev | Probability, Second Edition[END_REF] p.388), we have for each digit k,

lim n→∞ 1 n n i=1 1 {k} (x i ) - n i=1 E(1 {k} (x i )) = 0 P -a.s., which implies lim n→∞ 1 n n i=1 1 {k} (x i ) = lim n→∞ 1 n n i=1 p (i) k = p k P -a.s.. (4.2)
That is to say, for P almost every point in the space ∞ n=1 {1, . . . , c n }, the digit k has the frequency p k . Considering each point in N N as a continued fraction expansion of a number in [0, 1], we complete the proof.

Proof of dim H (E p ) ≥ 1/2. Take c n = n in Lemma 4.1, we find a point z ∈ E p , such that

z n = a n (z) ≤ n (∀n ≥ 1). (4.3)
For a positive number b > 1, set

F z (b) := {x ∈ [0, 1) : a k 2 (x) ∈ (b k 2 , 2b k 2 ]; a k (x) = a k (z) if k is nonsquare}.
It is clear that F z (b) ⊂ E p for all b > 1. We define a measure µ on F z (b). For n 2 ≤ m < (n + 1) 2 , set

µ(I m (x)) = n k=1 1 b k 2 . (4.4)
Denote by B(x, r) the ball centered at x with radius r. We will show that for any θ > 0, there exists b > 1, such that for all x ∈ F z (b), lim inf

r→0 log µ(B(x, r)) log r ≥ 1 2 -θ. (4.5)
In fact, for any positive number r, there exist integers m and n such that

|I m+1 (x)| < 3r ≤ |I m (x)| and n 2 ≤ m < (n + 1) 2 . (4.6)
By the construction of

F z (b), a n 2 (x) > b n 2 > 1. Let x = [x 1 , x 2 , • • • ]
. By Lemma 2.5, B(x, r) is covered by the union of three adjacent rank n 2 basic intervals, i.e.,

B(x, r) ⊂ I(x 1 , x 2 , • • • , x n 2 -1) ∪ I(x 1 , x 2 , • • • , x n 2 ) ∪ I(x 1 , x 2 , • • • , x n 2 + 1).
By the definition of µ, the above three intervals admit the same measure. Hence by (4.6), we have

log µ(B(x, r)) log r ≥ log 3µ(I(x 1 , x 2 , • • • , x n 2 )) log 1 3 |I(x 1 , x 2 , • • • , x m+1 )| . (4.7)
However, on the one hand, by (4.4) Combining (4.8) and (4.9), for any θ > 0, take b > 1 to be large enough, we have lim inf

-log µ(I(x 1 , x 2 , • • • , x n 2 )) = -log n k=1 1 b k 2 = n k=1 k 2 log b. ( 4
n→∞ log µ(I(x 1 , x 2 , • • • , x n 2 )) log |I(x 1 , x 2 , • • • , x m+1 )| ≥ 1 2 -θ, ∀x ∈ F z (b).
Hence by (4.7), we obtain (4.5). Since θ can be arbitrary small, by Billingsley Theorem ([1]), we have

dim H (E p ) ≥ 1 2 .
Thus Theorem 1.1 is already proved under the condition ∞ j=1 p j log j < ∞. Now assume ∞ j=1 p j log j = ∞. Then for any invariant and ergodic measure µ such that µ(I(j)) = p j for all j ≥ 1, we have

| log x|dµ(x) ≥ ∞ j=1 µ(I(j)) log j = ∞ j=1 p j log j = ∞. (4.10)
This implies D = ∞ (see the proof of Proposition 3.5 for the definition of D).

Then by Proposition 3.4, we have dim

H (E p ) ≤ 1 2 . Since we have already proved dim H (E p ) ≥ 1 2 , we get dim H (E p ) = 1 2 if ∞ j=1 p j log j = ∞.
This is in accordance with the formula of Theorem 1.1 under the convention that sup ∅ = 0, because (4.10) implies N ( p) = ∅. Finally, we remark that N ( p) = ∅ if and only if ∞ j=1 p j log j = ∞. We have seen the "if" part. For the other part, assume which implies that N ( p) = ∅.

A remark

As suggested by the referee, we add a remark on a problematic argument appearing in the literature. To obtain an upper bound of the Hausdorff dimension of a set, one usually applies the Billingsley's theorem by constructing a finite measure P on the set such that |U | s ≤ P (U ) (see [START_REF] Falconer | Fractal Geometry, Mathematical Foundations and Application[END_REF], p.67). For the set E p , where p j = 0 for some j, the Markov measure P satisfying (3.3) does not match because the cylinders starting with j do not charge the measure and the above inequality is obviously not true. This appeared unnoticed for long (see the remarks of Kifer [START_REF] Kiefer | Fractal dimensions and random transformations[END_REF], p. 2012).

This problem did exist in the proof of Theorem 2 in [START_REF] Billingsley | Hausdorff dimension of some continued-fraction sets[END_REF]. Let us briefly indicate how to get around the problem in the proof of Theorem 2 of [START_REF] Billingsley | Hausdorff dimension of some continued-fraction sets[END_REF] when p j = 0 for some j's. The basic idea is similar to that of Cajar (see [START_REF] Cajar | Billingsley dimension in probability spaces[END_REF], p.67) and that of Kifer [START_REF] Kiefer | Fractal dimensions and random transformations[END_REF]. Recall that P is a (k -1)-step Markov measure supported on the set {x ∈ [0, 1) ∩ Q c : a n (x) ≤ N for all n ≥ 1} .

It is uniquely determined by its values on the k-cylinders, namely, p(i 1 , . . . , i k ) = P ([i 1 , . . . , i k ]), (i

1 • • • i k ) ∈ {1, 2, • • • , N } k .
Let 0 < ǫ < 1 and P ǫ be the perturbed (k -1)-Markov measure determined by

p ǫ (i 1 , . . . , i k ) = (1 -ǫ)P ([i 1 , . . . , i k ]) + ǫ N k , (i 1 • • • i k ) ∈ {1, 2, • • • , N } k .
Now, we can apply the Billingsley's theorem with P ǫ to find an upper bound, and then get the desired result by letting ǫ → 0.

In the present paper, we have intentionally avoided using the Billingsley's theorem. The proof of upper bound consists of Propositions 3.4 and 3.5. Proposition 3.5 concerns some calculations for which the zero frequency of some digits will not cause any trouble. In the proof of Proposition 3.4, we have used a "covering argument" depending on the estimate (2.3) instead of using the Billingsley's theorem. This enables us to get the upper bound of the Hausdorff dimension.

  then by Lemma 3.3, we have α ≤ 1/2. Now suppose that D < ∞, which is equivalent to lim sup

k

  > 0 for all 1 ≤ k ≤ c n and cn k=1 p (n) k = 1, and that for any k ≥ 1, lim n→∞ p (n) k = p k . (4.1) Consider a product Bernoulli probability P supported by ∞ n=1 {1, . . . , c n }.

∞

  j=1 p j log j < ∞. Then the Bernoulli measure µ such that µ(I(j)) = p j satisfies| log x|dµ(x) ≤ ∞ j=1 µ(I(j)) log(j + 1) = ∞ j=1 p j log(j + 1) < ∞.

Acknowledgments :

The authors are grateful to the referee for a number of valuable comments. This work was partially supported by NSFC10728104 (A. H. Fan) and NSFC10771164 (J. H. Ma).