Orlando Ferrante

Marco Marazza

Alberto Ferrari

Formal Specs Verifier ATG: a Tool for Model-based Generation of High Coverage Test Suites

Keywords: Model-based Automatic Test Generation, High Coverage Test Suites, Formal Methods

In this paper we describe Formal Specs Verifier Automatic Test Generation, a tool generating high coverage test suites for embedded systems. Our tool implements a test case synthesis algorithm using a combination of model checking and optimization techniques starting from a Simulink/Stateflow model of the System Under Test. The main contributions of this paper are the following: we (1) give an extended description of our test generation algorithm, (2) describe the algorithm implementation as part of the Formal Specs Verifier framework, (3) present a concrete application of the tool to a cruise control case study and discuss experimental results comparing our algorithm with a state-of-the art COTS tool.

I. INTRODUCTION

Testing complex hardware-software embedded systems' architectures is one of the most important and costly phases of their entire development life cycle. A significant portion of development time is spent indeed for the verification and validation phases (V&V) during which teams of test engineers aim at discovering requirements' misinterpretation and implementation errors. Early error detection facilitates correction complexity, which in turn favors minimizing the overall system development cost and time. Testing consists in the execution of a set of predefined input vectors on the System Under Test (SUT) and observation of its response for detecting possible deviations from the expected behavior. In order to perform the testing phase, test vectors (or test cases) must be provided to the test execution environment. A test vector consists of a sequence of pairs (inputs, outputs) where inputs is the set of values to be applied to the system and outputs is the set of expected output values. A set of test cases is said to be a test suite. In case the execution of the system produces a set of output values not matching the expected ones the test can be used to highlight an error of the SUT. To quantitatively evaluate the quality of test suites, hence of the entire test execution, coverage metrics are used. Coverage metrics measure the effectiveness of a test suite i.e. how much it covers: (1) the structure of the SUT, as in MC/DC coverage metrics, (2) the set of requirements, as in requirements-based testing, or (3) a meaningful subset of admissible faults of the system, as in fault-injection testing. A good test suite satisfies as much as possible a given coverage metric. As a general rule the more the tests executed, the higher the confidence about the correctness of the system, as long as the executed test cases are of high quality. However, the generation of high quality test suites has a relevant impact on costs, e.g. complexity, setup time and execution time; hence, besides providing test suites maximizing a given coverage metric (optimality), it is of paramount importance to ensure a reduced number of test cases (efficiency). Typically, a test suite represents a tradeoff between optimality and efficiency. Several techniques have been adopted for the test generation process, especially in the area of test generation from system models whose main idea is to represent the SUT using a formal model and use automatic algorithms for the definition of test cases (modelbased test generation). In particular, the use of model checking techniques for test generation has been extensively explored: the test generation procedure is formulated as a problem of reachability, enriching the model with test objectives. Such test objectives must be achieved by the produced test suite in order for a given coverage metric to be satisfied. Once the model is enriched with test objectives, a model checker is used to analyze their reachability. In case a test objective is reachable, a counter-example is produced and stored as a test case, which is said to cover that test objective. Once all test objectives are covered, the generated set of tests is a high coverage test suite. In this paper we describe a novel approach for the automatic test generation from system models combining bounded model checking and optimization to generate efficient test cases and test suites. This paper is an extension of the algorithm presented in [START_REF] Ferrante | Model based generation of high coverage test suites for embedded systems[END_REF]; the main contributions of current paper can be summarized as follows: we (1) provide an extended description of the test generation algorithm, (2) describe the implementation of the algorithm in our Formal Specs Verifier framework and (3) show the application of our tool to a cruise control case study, along with additional experimental results. The proposed use case only covers the MC/DC as an example. The paper is organized as follows: Section II stresses the central role played by model-based test generation and modelbased testing in current industry development processes. Section III summarizes existing approaches for test case generation using exhaustive search techniques. Section IV provides a formal description of the problem our algorithm tries to address and Section VI describes our algorithm in details. Section VII provides an overview of the algorithm implementation as well as an application to a cruise control example and finally Section VIII concludes the paper.

II. MODEL BASED TEST GENERATION

Current practice in industrial systems' design and development process is to create models -also known as virtual prototypes-of the system being developed. A model is the artifact meant to imitate the system of interest. The system can be abstracted (modeled) at different levels of refinement and each level of refinement has its own best-fitting formal models. For example, a mathematical model is a set of formal definitions and mathematical formulas that describe the system under analysis. At the various system development phases, models are gradually refined until these get detailed enough to allow for physical implementation of the system. The system implementation can be achieved by means of either fullyautomated or pseudo-automated synthesis tools. With regards to testing, the main benefits of model-based approach are manifold: (1) tests can be automatically generated from models (model-based automatic test generation), [START_REF] Harman | A theoretical and empirical study of search-based testing: Local, global, and hybrid search[END_REF] testing can be performed at different refinement levels, far before the system gets implemented and, more interestingly, (3) tests generated at early phases of the system design can be refined and reused in later phases of the design to check whether all system properties are still fulfilled (back-to-back testing). Virtual prototyping and related testing are particularly important when the realization of components involves different technology suppliers and manufacturers, and requirements fulfillment must be checked by all the parties at all stages of the design process. This work targets model-based generation of test suites used to bring evidence that functional and non-functional metrics have been achieved by the design implementation.

III. RELATED WORK

Several techniques have been described in literature on the topic of automatic test generation. In this section we provide a brief summary of the available techniques pointing out the differences with the one proposed in this paper. Evolutionary and genetic approaches [START_REF] Harman | A theoretical and empirical study of search-based testing: Local, global, and hybrid search[END_REF] generate tests by randomly exercising the inputs of the SUT and measuring the quality of the test maximizing an objective function. Such function is derived from a structural analysis of the SUT, hence a test that provides better values of the cost function is selected in a set of possible generation process outcomes. In addition, evolutionary testing employs evolutionary algorithm techniques for selection and generation of new tests from a previously generated set of tests. These techniques rely on random search adding structure to the search to avoid generating low coverage tests. However, there are no guarantees that a selected behavior is effectively the best test case with respect to a given coverage metric. Model-checking based techniques are heavily studied for the generation of test cases. In [START_REF] Rayadurgam | Generating mc/dc adequate test sequences through model checking[END_REF] the authors describe the use of a model checking engine for the generation of test cases starting from a modified version of the SUT model enriched with a reset transition (i.e. a variable that, when true, resets the entire state of SUT to the initial state) and test objectives related to the coverage criteria. The proposed approach differs from ours in several ways. At first, the number of satisfied objectives is a non-deterministic result of the execution of the algorithm hence each test case may cover an arbitrary number of test cases, whereas our method allows for selecting the maximum number of satisfiable test objectives at each execution. In addition, the generated tests are produced monolithically starting from the initial model not allowing the application of an incremental test generation methodology; hence, the applicability of the technique to concrete industrial size cases heavily depends on the size of the input model [START_REF] Ferrante | Parallel nusmv: A nusmv extension for the verification of complex embedded systems[END_REF]. In [START_REF] Gargantini | Using model checking to generate tests from requirements specifications[END_REF] a methodology to generate tests from a formalization of requirements based on a tabular representation is described and the counter-example finding capabilities of model checkers is described as a technique to derive such counter-examples. The use of the model checker is straightforward and there is not any guarantee on the quality of generated tests. In [START_REF] Kadono | Using the nusmv model checker for test generation from statecharts[END_REF] the authors describe a method for extracting test cases from Statecharts state machines considering state transitions covering (i.e. covering all states and/or all transitions of the Statecharts). The generation is executed by enriching the model with additional states and searching from counter-examples of a specific CTL formula. The method may produce redundant test cases and after test generation an additional test reduction phase is needed. This approach differs from ours in that it does not guarantee the generation of a minimal set of high coverage test cases and it relies on a specific CTL-based formulation of the problem. As a consequence it cannot be solved using SATbased model checking techniques, that usually perform better than BDD-based techniques when a counter-example search is performed. In [START_REF] Di Guglielmo | Semi-formal functional verification by efsm traversing via nusmv[END_REF] a method exploiting random-simulation and formal verification is described based on a semi-formal method for the traversing of Extended Finite State Machines (EFSMs) that allows for reaching deep test cases exploiting the strength of simulation-based approaches. However, no guarantees are given in terms of efficiency of the generated suite and the use of the model checking engine is limited to the analysis of constraints to guide the main simulation-based test generation algorithm. In [START_REF] Hamon | Generating efficient test sets with a model checker[END_REF] the authors propose a test generation method based on the reachability of given test cases as well as extension of already generated test cases. This approach differs from ours in several ways. First, it does not provide guarantees that every generated test is optimal with respect to a specific coverage metric. The model checker is called in order to satisfy some of the test goals but no guarantees are provided to the number of the goals satisfied. Second, it does not employ an incremental test generation procedure: at every execution no guarantees are given about the length of the found counterexample: there might be cases in which a high coverage test case has a given length, but the same results would be obtained with a shorter test case.

IV. PROBLEM FORMULATION

A. Model And Coverage Criteria Formulation

We focus on the problem of generating high coverage test suites for discrete time models formalized as connections of blocks. Several languages are available to capture embedded control systems using this formalism, e.g. MATLAB Simulink/Stateflow and Esterel SCADE 1 . A model M can be formally represented as a connection of blocks exposing a well-defined interface in terms of input and output ports, each describing its run-time behavior as an Extended Finite State Machine (EFSM). Block's interfaces are formalized as a pair of input and output vectors u In order to be fully specified, an initial value for the state vector should be provided. Connections between blocks act as constraints between the values of inputs and outputs that must match at every time step. The transition relation of a block contains logical and arithmetic expressions for Boolean and bit vectors. For a complete description of the richness of the Simulink/Stateflow languages please refer to [9]. A test suite Π is a set of test traces π 1 , . . . , π n , each representing a finite sequence of admissible input values:

[k] ∈ X, y[k] ∈ Y ,
π k = u[1], . . . , u[m] .
Each test trace π k has a finite length π k = m that corresponds to the maximum time step for which the model is exercised. The objective of test generation is to produce a test suite that exercises the model for maximizing a given coverage criterion. Several coverage criteria have been defined depending on the test generation algorithm input. As an example for finite state machine a relevant coverage criteria is related to the capability of exercising the state machine enabling as much transitions as possible (transition coverage criteria). In model-based approaches several criteria can be defined depending on the input model. In our flow, Simulink and Stateflow models are processed and we use the following coverage criteria (that can be applied to similar languages such as Esterel SCADE). MC/DC coverage: following the definition of MC/DC for software a similar criterion has been followed for model elements associated to Boolean formulas. More precisely, for each block corresponding to a Boolean formula y = f [u 1 , . . . , u N ] to fully satisfy the criterion there should exist a trace such that each input affects the truth value of the output independently of the other inputs [START_REF] Chilenski | Applicability of modified condition/decision coverage to software testing[END_REF]. Relations coverage: for each block that compares two values y = u 1 u 2 (where ∈ {≤, <, =, =, >, ≥}), there should exist tests for which the output value y transitions from FALSE to TRUE and vice-versa. State coverage: for state machines, there should be tests such that the state machine states assume all possible values. Transition coverage: for each state machine, there should be tests such that all the transitions of the state machines are asserted. In our flow, Simulink and Stateflow models are processed and we use the following coverage criteria (that can be applied to similar languages such as Esterel SCADE): MC/DC coverage, relations coverage, state coverage and transition coverage. The first two criteria apply also to guards and actions of the state machines; hence, if a guard is a composite Boolean expression, the test suite should provide tests covering the guard with a maximum MC/DC percentage value.

B. Problem Formulation

The test generation problem we address in this paper can be now formalized as follows. Given a discrete time model M, a coverage metric and a test generation maximum time step m, generate a test suite Π = {π 1 , . . . , π n } of vectors of length π k ≤ m, such that: 1) the achieved coverage is maximal with respect to the coverage metric 2) each test trace contributes to increase the coverage of the test suite

Objective 1) ensures that the generated test suite achieves high test effectiveness, provides high coverage of the input model and captures as much errors as possible, whereas objective 2) ensures that the cost of test execution is well-balanced, meaning that each test trace effectively improves the overall coverage of the model and avoids execution of useless tests. In this section the test generation algorithm presented in [START_REF] Ferrante | Model based generation of high coverage test suites for embedded systems[END_REF] is summarized and extended showing the usage of monitor variables for the automatic synthesis of test cases. In Section VI-A the algorithm is summarized, whereas in Section VI-B the details of the efficient search sub-activity are provided. This section describes the details of the efficient search subactivity of the test generation algorithm we presented in [START_REF] Ferrante | Model based generation of high coverage test suites for embedded systems[END_REF] and differs from a previous implementation [START_REF] Marazza | Automatic generation of failure scenarios for SoC[END_REF] generating Minimal Critical Failure Sets.

A. Test Generation Flow

The overall flow of the test generation algorithm is described in Fig. 2. The initial formal model is elaborated in a model transformation step that enriches the input adding a finite number of test objectives. Each test objective represents a Boolean condition over the discrete time that should be satisfied by at least one generated test case. The formal description of the test objective depends on the coverage criteria selected by the user (e.g. MC/DC, relation coverage, transition coverage, etc.). During the test objectives generation step the input model is instrumented with additional Boolean expressions representing the test objectives the Automatic Test Generation (ATG) step should satisfy. Each test objective is derived by analyzing the input model blocks according to the selected criteria. As an example consider the Simulink model in Fig. 1 containing a Comparator block and a Switch block. The Comparator has an output that is true when input u 1 is greater or equal to zero and false otherwise, while the Switch connects to the output the first input when the control input c is FALSE and the second input (u 2 ) otherwise. In case the user selected the relation coverage criteria, two Test Objectives (TOs) to be satisfied would have been derived for the Comparator block: the first TO is true when the comparator output transitions from FALSE to TRUE and the second is true when it transitions from TRUE to FALSE. Similarly, when the MC/DC criterion is selected, two additional test objectives are derived: one that is true when the control input of the Switch transitions from FALSE to TRUE and another one for the opposite condition. Using this approach we are able to convert different coverage criteria to a common criterion (test objective criterion) and the sub-sequent test generation step will try to produce a high coverage test suite maximizing the number • The test suite coverage percentage always increments with the addition of a new test case;

• Each test case is not redundant: no other test case of the same length covering the same (or a super-set of the) set of satisfied test objectives exists.

The ATG process can be bounded in time or in number of generated test cases by setting a set of parameters controlling its execution.

The algorithm is described in Algorithm 1. The algorithm is described in [START_REF] Ferrante | Model based generation of high coverage test suites for embedded systems[END_REF] and is summarized here for clarity. At start-up the explored depth bound is set to the initial value (it might be 1 or user-defined). The algorithm loops until all test objectives are satisfied or a given resource/time/depth bound is reached and in each iteration a subset of not yet satisfied test objectives (TO) is identified. The subset may be the entire set of unsatisfied TOs or it may be driven by the user needs (i.e. all the TOs related to the coverage of a given sub-function of the SUT, etc.). Once the TOs are selected an inner loop is performed until an explicit exit condition occurs or all the test objectives selected have been covered. In the inner loop, the formal engine is queried to find all the test cases of length equal to l that maximize the number of satisfied TOs among the ones selected at previous step. Each test case satisfies a unique subset of TOs (i.e. there are not two test cases satisfying the same subset of TOs). If new test cases are found, collect all the subsets of satisfied TOs and remove them from the search to avoid search again their satisfaction in next iteration. In case no new test cases are found we can assess that for the given length bound, test cases satisfying the selected test objectives do not exist. This claim is possible because of the exhaustive search (Ψ, tc) = FindTestCaseMaxSat(Ω[l], l)

7: if ({tc} = ∅) then 8: T C[l] ← T C[l -1] ∪ {tc} 9: Ψ[l] ← Ψ[l -1] ∪ Ψ 10: Ω[l] ← Ω[l] -Ψ 11:
bDone = f alse end while 17: end while

18: return Ψ[l], T C[l]
of the underlying formal back-end. Hence the TO selector can be executed again to select a new subset of TOs or exiting from the loop if no more TOs can be selected. If during last iteration no new test cases are found, it is not possible to satisfy any of the current set of selected TOs for the given length, hence the explored depth bound is increased and a new iteration is started if resource limits are not reached (i.e. timeout, memory consumption, etc.).

1) The algorithm loops until all test objectives are satisfied or a given resource/time/depth bound is reached; 2) While this bound has not been reached: a) A subset of not yet satisfied test objectives (TO) is identified. The subset may be the entire set of unsatisfied TOs or it may be driven by the user needs (i.e. all the TOs related to the coverage of a given sub-function of the SUT, etc.) b) The following loop will be executed until new test cases are not found or the selection mechanism has not anymore TOs to select: i) The formal engine is queried to find all the test cases of length equal to the current length bound that maximize the number of satisfied TOs among the ones selected at previous step. Each test case satisfies a unique subset of unsatisfied TOs i.e. there are not two test cases satisfying the same subset of TOs. ii) If new test cases are found, collect all the subsets of satisfied TOs and remove them from the search to avoid search again their satisfaction in next iteration. iii) In case no new test cases are found, we can assess that for the given length bound do not exist test cases satisfying the selected test objectives. This claim is possible because of the exhaustive search of the underlying formal backend. Hence the TO selector can be executed again to select a new subset of TOs or exiting from the loop if no more TOs can be selected. c. If during last iteration no new test cases are found, it is not possible to satisfy any of the current set of selected TOs for the given length. Hence the explored depth bound is increased

The algorithm is guaranteed to terminate provided that the test objective selection performed at line 2 is able to identify the test objectives previously selected even if not satisfied. The simplest admissible selection mechanism picks the entire set of unsatisfied objectives at once and exits from the outer loop at point 3 after the first execution. More efficient test selection mechanisms are admissible and possible but for brevity we will not cover this topic in the report. During the search at line 6 a sub-procedure is called in order to produce queries to the formal engine and storing the counter-examples provided by the model checker as test cases. A Detailed description of the procedure is provided in Algorithm 2. The overall algorithm proceeds in an incremental fashion. Starting from an initial exploration depth bound all the test cases of a given length that maximally satisfies the test objectives are found. After the formal engine proofs that no more test cases of the given length can satisfy additional TOs, the explored bound is incremented and the search is executed again. Due the exhaustive nature of the search, it is guarantee that the set (or a super-set) of the TOs satisfied by a test case of length k cannot be satisfied by test case of length j < k. Hence, the method is efficient in the following sense:

• each new test case covers only new test objectives hence it strictly increments the coverage of the test suite by covering only not yet covered test objectives;

• every test case covers the maximal number of unsatisfied test objectives of a given length

In addition the incremental nature of the algorithm allows for generating test cases by starting from a complex problem in the number of TOs but simpler in the unrolling of the SUT model transition relation and incrementally increasing the complexity due the unrolling of the transition relation but reducing the number of satisfiable test objectives. Our experience with industry sized models shows that this trade off allows for applying ATG on complex models since the satisfied TOs are removed incrementally at every step and the complexity of the formal problem introduced by the increment of the explored bound is partially mitigated by the reduced number of TOs to be satisfied.

B. Maximal Coverage Test Case Generation

The mechanism to find the test case satisfying the maximum number of TOs is an important part of the optimal test suite generation procedure and in this sub-section it will be described in details. The algorithm relies on the concept of monitor variable associated to a test objective TO. A monitor is an integer variable and can be seen as a function m associated to a test objective TO that at every execution step k evolves as follows:

m T O [k] =          1 if m T O [k -1] = 1, or ∃j ≤ k s.t. T O satisfied at step j 0 otherwise and m T O [0] = 1 if T O satisfied at initial step 0 otherwise
Each test objective T k has an associated monitor variable m k . Given a set of unsatisfied test objectives T 1 , T 2 , . . . , T N and their associated monitors m 1 , m 2 , . . . , m N the algorithm that searches for the test case set is described by the pseudo-code shown by Algorithm 2. The algorithm loops until there are counter examples found (lines 1 . . . 9). The search at step 2 searches for a counter example that maximizes the sum of the values of the monitors m 1 , . . . , m N . The counter-example is obtained applying bounded model checking on the input model. The model checker takes into account the dynamics of the SUT and the given coverage metrics (that is used to generate the test objectives and the associated monitors).

The maximization procedure ensures that the found counterexample exercises as much test objectives as possible. If a counter example is found the monitor values' configuration is extracted and a new constraint is added to the model in order to exclude the configuration from future searches (line 4, 5). This step ensures the progress of the iteration loop avoiding the model checker to find a counter-example satisfying the same test objectives over and over (there might be an infinity of them). Then the test case is extracted from the counter example by storing the values of the interesting variables of the system (inputs, outputs and internal variables) and finally it is added to the output set (line 6 and 7). The loop is then executed again looking for additional maximal possible monitor configurations until no more configurations can be found. This guarantees that: 1) every test case produced by the algorithm maximizes the number of test objectives covered and it is not possible that another test case satisfies a super set of the covered ones and 2) for the produced test suite is true by construction that there are not two test cases satisfying the same set of test objectives.

VII. ALGORITHM IMPLEMENTATION IN FSV

Our automatic test generation algorithm has been implemented in the FormalSpecs Verifier (FSV) framework for the verification of embedded systems.

A. FormalSpecs Verifier

The FormalSpecs Verifier is a framework targeting complex embedded systems verification. The core of the tool is based on a translator from Simulink modeling language to NuSMV native language. The transformation process produces a semantically equivalent NuSMV representation of the input model taking into account the non-determinism resolution that may be introduced during the transformation step. In Fig. 4

B. Cruise Control Example

We show the application of our test generation algorithm to a cruise control reference example implemented in MAT-LAB Simulink. We consider a modified version of the model proposed by Aldrich in [START_REF] Aldrich | Coverage analysis for model based design tools[END_REF]. A cruise control is an embedded system that regulates the speed of a vehicle based on a set of commands provided by the driver; the interface of the control algorithm appears to the system as represented in Fig. 5. The ACCEL and BRAKE inputs are Boolean values representing the pressure of the accelerator and brake pedal by the vehicle driver. The CC ON Boolean input is set when the driver wants to engage the cruise control. The KEY ON input represents the presence of the key (Boolean value) and finally the CRUISE SPEED value is an unsigned integer input set by the user representing the desired cruise speed. The outputs of the controller are the ENABLE Boolean value that is true whenever the cruise controller is actively controlling the vehicles speed and the SPEED unsigned integer value representing the reference speed passed to the cascade speed controller that acts directly on the engine throttle based on the reference and current speed values. In Fig. 6 the internal modal logic of the controller is represented as an extended finite state machine (Simulink Stateflow machine). The controller is initially in an OFF state and passing thru an IDLE state can go on a controller active mode (CC MODE) or disengaged mode (CC MODE DISABLED) in which it returns the direct control to the driver. This happens whenever she/he interacts with the car pressing the acceleration or brake pedal. The FormalSpecs Verifier Automatic Test Generation (FSV-ATG) tool is executed to perform the efficient generation of test cases. As a first step the input model is elaborated and automatically enriched with a set of test objectives to enable the subsequent test generation procedure for covering the state machine to obtain MC/DC coverage. The summary of the generation process for the cruise control example is shown in Fig. 7: a total number of 191 test objectives have been added to the formal model to enable the ATG step. The ATG activity is executed and as results generates a set of efficient test case. For the cruise control example the ATG generated 31 test cases given a bound of 100 maximum test cases and a depth bound of 30 steps. To validate the obtained results we simulated the generated test cases on the SUT evaluating the coverage value using the Simulink Verifcation and Validation (V&V) toolbox which provides an independent measure of effectiveness for our approach. The obtained independent measured coverage has been of 100% for decision, condition coverage and MC/DC coverage.

C. Additional experiments

In order to quantitatively compare the proposed algorithm and implementation with respect to state-of-the art tool we compared the test suites generated using the Formal Specs Verifier ATG tool with the Simulink Design Verifier tool (SLDV) 2 that represents an industrial strength tool for the generation of test suites from MATLAB Simulink/Stateflow systems. A set of verification cases has been set up containing blocks ranging over a rich subset of Simulink/Stateflow Libraries. The results for the comparison of the MCDC coverage of the produced suites are presented in Fig. 8. The coverage values are obtained using the MATLAB VnV toolbox 3 . The analysis of the results shows that in some cases the SLDV tool is capable of achieving higher coverage with respect to FSV-ATG whereas in other cases the opposite is true. Cases where the MC/DC coverage percentage is 0% for FSV-ATG indicate that our algorithm was not able to generate any test cases. The analysis of the "losing" cases for the FSV-ATG tool showed that the low coverage is due to an inefficient translation of the Simulink/Stateflow block that does not allow efficient application of our ATG algorithm, as highlighted in Sec. VII-A. As an example, in some cases the FSV-ATG tool produces a structure of comparison/logical blocks that is inefficient for generating a high number of test objectives. The optimization of the translation step for achieving efficient coupling with the test generation engine is one of the activities we have planned as the next development steps of our tool.

VIII. CONCLUSION

Our model-based test generation algorithm produces a test suite starting from a model of the system under test (SUT) that is enriched with a set of test objectives to be satisfied by the test cases derived from a given coverage metric. The produced suite covers the model test objectives in an efficient way such that 1) there are no two test cases satisfying the same set of test objectives and 2) each test case covers the maximum number of test objectives for a given length. Our algorithm relies on the combination of bounded model checking with an optimization-based formulation of the test generation problem. The algorithm is implemented using an incremental execution approach that mitigates the complexity of the problem allowing successful application to complex 2 https://www.mathworks.com/simulinkdv 3 http://www.mathworks.com/products/simverification/ industrial use cases. Several ways of algorithm-improvement are possible, though. From a technical standpoint the algorithm could be improved by employing parallelism at test generation level and not only at the formal back-end level. In addition, the last advancements in pseudo-Boolean SAT solving can be taken into account to explore additional formal back-ends. From a methodological standpoint the use of contract-based design ( [START_REF] Mangeruca | Formalization and completeness of evolving requirements using contracts[END_REF], [START_REF] Benveniste | Contracts for System Design[END_REF]) could allow for exploiting a compositional approach at test generation level. From a methodological standpoint the use of contract-based design could allow for exploiting a compositional approach at test generation level. Finally, the use of formal proofs on the model can be used to ease the process of test generation.

  where X and Y represent the vector domains and k ∈ N the discrete time. Each block behavior is then formalized as a transition function F [u, x, f, k] that at each discrete time k maps the vectors input u[k] and current state x[k] to an output vector y[k] and next-state values x [k], where x [k] = x[k + 1].

Fig. 1 .

 1 Fig. 1. Test Objective example

Algorithm 1 ← 1 ,bDone = f alse 5 :

 115 High Level View of the Test Generator Algorithm Input: Θ = {T O 1 , . . . , T O J } Input: L max explored step, T OU T Algorithm execution timeout Init l Currently Explored Length index T C[0] ← ∅, Test Case Collection @ step 0 Ω[0] ← Θ = set of not yet satisfied TOs @ step 0 Ψ[0] ← ∅ = set of satisfied test objectives @ step 0 1: while ((Ψ = Θ) ∧ (l ≤ L) ∧ (¬T OU T )) do 2: Ω[l] = SelectTestObjectives (Ω[l -1], Θ) while (Ω[l] = ∅) ∧ (bDone = f alse) do 6:

Algorithm 2 3 :

 23 Maximal Coverage Test Generator Algorithm Input: model: The enhanced formal model under analysis Input: m 1 , m 2 , . . . , m N : Set of model variables representing unsatisfied test objectives' monitors Output: testCases : Set of produced test cases 1: repeat 2: find a counter example such that m 1 + m 2 + . . . + m N is maximal if counter example exists then4: define m * = [m * 1 , m * 2 , . . . , m * N ] the found monitor configuration 5: exclude m * from the admissible solutions of the maximization search 6: extract test case values for the found counter-example 7: add extracted test case to the testCases set 8:end if 9: until counter-example has been found

Fig. 3 .

 3 Fig. 3. Test Generation Work Flow

Fig. 4 .

 4 Fig. 4. Formal Specs Verifier Model Transformation Flow

Fig. 5 .

 5 Fig. 5. Controller sub-system

Fig. 6 .

 6 Fig. 6. Controller Extended Finite State Machine

Fig. 8 .

 8 Fig. 8. Comparison of MCDC Coverage Results

http://www.esterel-technologies.com/


Development of an algorithm for energy efficient automated train driving

Keywords: energy efficiency, energy management, automatic train operation, optimization problem, genetic algorithm, smart transportation, real-time systems

Development of an algorithm for energy efficient automated train driving

Artem Ozhigin, Pavel Prunev, Victor Sverdlin, Yulia Vikulina

Introduction

In recent years the question of energy efficiency has become the gravest it has ever been. Railway operation considerably contributes to the overall energy consumption and, therefore, deliberate railway management systems can enhance the global natural resources savings. Significant advances have been made to improve the efficiency of the train motion: minimization of resistant forces, construction improvements, recuperation capabilities [START_REF] Kondo | Recent energy saving technologies on railway traction systems[END_REF]. Nevertheless, there still is opportunity for further savings of energy. In particular, taking into account the track profile, optimization of the driving strategies, measuring of the consumed energy and further analysis to some degree can reduce the cost of trips.

Innumerable studies have contributed to the energy efficiency in train control. A detailed review of works on optimal train control can be found in [START_REF] Yang | A review on intelligent control for energy-efficient train operation[END_REF]. Diverse methods and approaches can be found in literature concerning optimal train control, optimal schedule construction and railway system management. Beginning with the pioneer works of Milroy [START_REF] Milroy | Aspects of automatic train control[END_REF], Howlett [START_REF] Howlett | The Optimal Control of a Train[END_REF], Asnis et al. [START_REF] Asnis | Solution of the problem of the energetically optimal control of the motion of a train by the maximum principle[END_REF], numerous scientists developed and investigated train control optimization problem, giving birth to great number of specific methods and approaches. Two different models can be highlighted among them: a model with continuous control [START_REF] Milroy | Aspects of automatic train control[END_REF], [START_REF] Howlett | The Optimal Control of a Train[END_REF], [START_REF] Asnis | Solution of the problem of the energetically optimal control of the motion of a train by the maximum principle[END_REF]; and a model with discrete control [START_REF] Benjamin | Energy-efficient operation of long-haul trains[END_REF]. The first model considers applied acceleration as a control variable whereas the second one considers switches of traction and brake for this purpose. The second model has received much attention in recent decades, due to the opportunity to model the real train control directly and to obtain as a solution a sequence of a driver's actions for optimal train control.

One of the most effective approaches to energy efficient driving is speed regulation through coasting control, minimization of brakes or a combination of optimal strategies [START_REF] Howlett | Energy-efficient train control[END_REF]. In [START_REF] Albrecht | Energy-efficient train control: from local convexity to global optimization and uniqueness[END_REF] driving strategy 'power-speedhold-coastbrake' is considered and the algorithm of optimal switching points finding is proposed. The "speed regulation without braking" strategy is implemented for efficient manual driving on long distance lines in [START_REF] Sicre | Modelling and optimising energy efficient manual driving on high speed lines[END_REF].

Among the methods of solution of the train control problem one can point out two main groups widely used by the researchers: exact methods and meta-heuristical.

The first approach was utilized, for instance, by Howlett in [START_REF] Howlett | The Optimal Control of a Train[END_REF]. He has elaborated the algorithm to minimize the energy consumption using Pontryagin's principle and Kuhn-Tucker equations. Assuming that the trip shall be completed with a fixed time he finds local optimal points of control switching for each part of the track with a particular slope and obtains a globally optimal strategy. In the paper [START_REF] Matsuura | Optimal train speed profiles by dynamic programming with parallel computing and the fine-tuning of mesh[END_REF] Matsuura investigated an algorithm optimizing a train speed profile by the Bellman's Dynamic Programming (DP), taking into account track profile incline, regenerative brake and so on. The parallel computing technique is introduced to deal with the computation complexity issue. Comparison of Dynamic Programming (DP), gradient method and Sequential Quadratic Programming (SQP) is introduced in [START_REF] Miyatake | Optimization of Train Speed Profile for Minimum Energy Consumption[END_REF] in application to the complete optimization of multiple trains' speed profiles and energy storage devices with constraints. However, the increase of control input dimension causes explosion of computation time and use of memory space and the introduced methods cannot be applied to real-time control of trains.

As long as the train motion is a multicriterial process, constructing an analytic optimization algorithm, taking into account all the constraints, is still very difficult. Using Pontryagin's principle of maximum or the Bellman equation requires significant simplifications of the problem considered. Furthermore, these methods reduce the problem to the system of differential or even partial differential equations which need to be solved by means of some numerical methods, e.g. Runge -Kutta methods or finite element method respectively. Therefore, the final solution inevitably contains inaccuracies.

The aforesaid is the reason of the ubiquity of the second group of methods involving diverse metaheuristics: genetic algorithm (GA), simulated annealing, tabu search, etc. Although these approaches do not guarantee to find the global optimum point, in most cases it is sufficient for real-time applications to know the quasi-global optimum. During the last decade a lot of papers were published where that kind of algorithm was used for solving the efficient train control problems ( [START_REF] Sicre | Modelling and optimising energy efficient manual driving on high speed lines[END_REF], [START_REF] Sicre | Real time regulation of efficient driving of high speed trains based on a genetic algorithm and a fuzzy model of manual driving[END_REF], [START_REF] Boschetti | Optimizing the Energy Efficiency of Electric Transportation Systems Operation Using a Genetic Algorithm[END_REF], [START_REF] Besinovic | A simulation-based optimization approach for the calibration of dynamic train speed profiles[END_REF], [START_REF] Martinis | The evaluation of energy efficient solutions in train operation: a simulation-based approach[END_REF], [START_REF] Wei | Energy Saving Train Control for Urban Railway Train with Multi-population Genetic Algorithm[END_REF], [START_REF] Bigharaz | Simultaneous Optimization of Energy Consumption and Train Performances in Electric Railway Systems[END_REF], [START_REF] Tuyttens | Simulation-Based Genetic Algorithm towards an Energy-Efficient Railway Traffic Control[END_REF], [START_REF] Cucala | Fuzzy optimal schedule of high speed train operation to minimize energy consumption with uncertain delays and driver's behavioural response[END_REF], [START_REF] Carvajal-Carreno | Efficient driving algorithms for non-distributed and distributed trains with the CBTC signalling system[END_REF], etc.).

The genetic algorithm, for example, was utilized in works of Sicre et al. [START_REF] Sicre | Modelling and optimising energy efficient manual driving on high speed lines[END_REF], [START_REF] Sicre | Real time regulation of efficient driving of high speed trains based on a genetic algorithm and a fuzzy model of manual driving[END_REF], where the problem of constructing an on-line advisory system for efficient manual train driving is solved. The manual driving is modelled by means of fuzzy parameters as long as the exact moment of control switching depends on the reaction speed of the driver.

The GA-based approach is proposed in [START_REF] Boschetti | Optimizing the Energy Efficiency of Electric Transportation Systems Operation Using a Genetic Algorithm[END_REF] to search for a control strategy under different operating conditions and the realistic electromechanical simulation model is utilized for energy consumption calculation. GA is also utilized in the works of Besinovic et al. [START_REF] Besinovic | A simulation-based optimization approach for the calibration of dynamic train speed profiles[END_REF] for calibration of the train input parameters using observed data; De Martinis et al. [START_REF] Martinis | The evaluation of energy efficient solutions in train operation: a simulation-based approach[END_REF] proposed a framework operating on timetable, rolling stock characteristics, signalling system and infrastructure characteristics and an optimization subroutine composed by the genetic algorithm.

In most works the control sequence is represented by a pre-defined array of positions or as a mesh with fixed size. This evidently causes the loss of generality. In this work we propose an approach where the control switching points are computed during the GA routine, allowing to select not only a type of control, but also a particular switch moment minimizing the consumed energy.

Among the works utilizing other meta-heuristical algorithms we should mention the PhD thesis of Kim [START_REF] Kim | Optimal train control on various track alignments considering speed and schedule adherence constraints[END_REF] with an approach to optimize train motion strategies for various track alignments (single or mixed) and maximum operating speed (constant or variable) using simulated annealing. Modified and hybrid meta-heuristic methods applied to optimal train control were considered in works of Wei et al. [START_REF] Wei | Energy Saving Train Control for Urban Railway Train with Multi-population Genetic Algorithm[END_REF], Bigharaz et al. [START_REF] Bigharaz | Simultaneous Optimization of Energy Consumption and Train Performances in Electric Railway Systems[END_REF].

Considering the uncertain disturbances arising from the weather, route and locomotive rolling conditions, Li et al. [START_REF] Li | Existence of an optimal strategy for stochastic train energy-efficient operation problem[END_REF] proved the existence of an optimal operation strategy for stochastic train energy-efficient operation.

The main advantages of meta-heuristic methods are relatively simple logic, the capability to expand the problem and an opportunity to take into account any new circumstances and conditions without affecting the optimization algorithm. These methods do not have requirements of differentiability or continuity on the functions involved. The possibility of mutation in the genetic algorithm allows to not be stuck in the local optimum. Moreover, as a result of computations a set of solutions is obtained and if the best one is not suitable for some reason, another candidate can be taken right away. Meanwhile, analytical methods provide only one solution and in order to satisfy new conditions, the necessity of the recalculation of the whole problem arises. In some cases, reconstruction of the method itself may be needed. This can be impossible or difficult to carry out.

Until recently, the question of efficient train control was a responsibility of the driver, who used his own experience, known heuristics and driving patterns as a key to optimality and efficiency. Autopilot (automated train driving) algorithm providing efficient driving will take the responsibility of routine actions and will allow the driver to concentrate on more important trip factors.

Despite the fact that there are a lot of works devoted to energy efficient train control (see [START_REF] Sicre | Modelling and optimising energy efficient manual driving on high speed lines[END_REF], [START_REF] Martinis | The evaluation of energy efficient solutions in train operation: a simulation-based approach[END_REF], etc.), none known to us deals with completely automated control. Apparently, this is a consequence of the local regulations prescribing that the train must be always under human control as a guarantee of safety. An optimization algorithm can be implemented only as an advisory system. On the contrary, in this work we are trying to develop an algorithm which will be able to implement automated real-time train control capable to realize the whole trip without human interference. The genetic algorithm is applied to the problem in order to find an optimal sequence of controls for energy efficient train driving.

Safety is a crucial matter in transportation systems. As mentioned before, the automated train driving function is surrounded by train control systems of a high safety integrity level and it does not control the motion of the train directly. The commands generated by the system under consideration are sent to train control systems and treated by them as recommendations which will only be implemented if safety is not violated. This way of integration isolates the automated train driving system and prevents any impact on train safety. Due to the isolation there are no real safety requirements imposed on the autopilot. On the other hand, in order to be in actual control of motion and perform designated functions the autopilot shall not violate criteria implied by safety systems for its requests to be implemented. Inappropriate operation of autopilot can affect on passengers' comfort (e.g. frequent switching between traction and braking modes is perceived as not comfortable) and even cause engagement of emergency modes of safety systems (e.g. when autopilot operation causes violation of speed limit). This must be taken in serious consideration in the development of an automated driving algorithm.

The objective of our work is to develop an automatic real-time train control system with a view to reducing of energy consumption. For that purpose we need to find an equilibrium between three antithetic objectives: energy efficiency, fulfilment of schedule constraints and speed restrictions and passengers' comfort.

Problem Statement

The aim of this section is to state an optimization problem for energy efficient train control. Firstly, let us consider the mathematical model of the train which is one of the most important parts of a system.

Train Model Description

A train model is designed to define the velocity and position of a train for each moment of a trip basing on the following parameters: traction intensity (%), brake intensity (%), track profile, adhesion coefficient, weather conditions such as wind velocity. Generalized traction and brake curves were determined from the train documentation.

From the Newton's law we obtain the following dynamic equation of train motion:

M d 2 S dt 2 = (u tr u br ) F tr F br -F f r -F h , ( 1
)
where M is train mass, S is train coordinate, F tr is traction effort, F br is brake effort, F f r is friction force, F h is projection of the gravity force, u tr and u br are traction and brake control correspondingly. Traction force values are determined by the train specification. Brake force F br comprises two components: electrodynamic and pneumatic brake effort. Train control system is responsible for actual value of the traction effort and for dynamic distribution of braking efforts between pneumatic and electrodynamic brakes. Control system is trying to fulfil the requests from Autopilot on required speed and acceleration values while simultaneously preventing the violation of safety requirements.

The next term in the eq. ( 1) is friction force F f r modelled using "Davis equation" [START_REF] Rochard | A review of methods to measure and calculate train resistances[END_REF]:

F f r = A + Bv + Cv 2 ,
where v = dS dt is velocity and the coefficients A, B and C are determined by aerodynamic characteristics and rolling resistance. Their values were derived from the train specification. Aerodynamic characteristics include influences of the wind velocity and tunnels. It is possible to take into account other weather conditions, e.g. precipitation, as easy as adding one more term in right-hand side of the eq. ( 1).

The more sophisticated train models, including non-linear ones, can also be considered and the more realistic electromechanical simulation can be embedded. Complexity of the simulation will not require any changes in the optimization algorithm. Thus, in this work we use the standard simplification in order to pay all the attention to the energy efficiency and finding an optimal speed profile.

Solution of the optimization problem U (S) is represented as a sequence of commands u i and corresponding switching point s i .

U (S) = [s i u i ] =      s 1 u 1 s 2 u 2 . . . . . . s n u n      (2)
Control u i can be equal to {0, 1, 2, 3} where 0 designates the brake, 1 -coasting mode, 2 -cruising mode and 3 -traction. Coasting mode is implemented if traction and brake forces are equal to zero and only friction and gravitational forces act on a train. Cruising mode means the keeping by all means of the current speed using traction or brake force depending on the track profile and other factors. Each mode defines the value of control variables u tr and u br : both equal to zero in coasting mode, one of them has some value from an interval [0, 1] in cruising mode, u tr = 1 in traction mode, u br = 1 in braking mode.

The train model provides information about energy consumption and recuperated energy. Energy consumed from the beginning of a trip up to the current point, is computed according to the relation: Recuperated energy depends only on the electrodynamic brake force which is a component of the whole brake effort as was described above:

E = s 0 F tr dS,
E r = η s 0 F edbr dS,
where E r is recuperated energy, η is the maximum percentage of the electrodynamic brake energy which can be returned to the grid, F edbr is electrodynamic brake force.

The regenerative brake opportunity enables to significantly reduce the energy consumption, although, there are some difficulties in efficient utilization of recuperated energy, the first of which is a necessity to absorb the energy by accelerating train or by energy storage device. If the energy is not absorbed, energy consumption is larger because of regenerative failure [START_REF] Miyatake | Optimization of Train Speed Profile for Minimum Energy Consumption[END_REF]. However, due to the lack of necessary infrastructure facilities in Russian railways the recuperated energy can not be spent efficiently. Therefore, in this work the regenerative brake energy is not a part of fitness function and is considered as loss.

Optimization Problem Statement

To complete an optimization problem statement let us define the constraints and restrictions imposed on it. The information defining the train run parameters is described in a digital schedule uploaded before the trip into the automatic train driving system. This schedule contains set of route points with prescribed arrival and departure time, values of speed restriction on different route segments, altitude of the route points. The intention of the system is to make decisions on the driving mode to be implemented at every moment along the route based on current coordinate and schedule data.

Firstly, the optimal train control should satisfy the schedule constraints. We can describe the time condition using the following inequality:

T schedule -∆ -≥ T ≥ T schedule + ∆ + ,
where ∆ -and ∆ + characterise the time tolerance interval, T is trip time for the particular solution, T schedule is time demanded by the schedule. The boundary conditions define departure and destination points:

S(0) = S depart , S(T ) = S dest ,
where S depart and S dest are respectively the departure and destination stations coordinates. Finally, the velocity shall satisfy the speed restrictions:

V ≤ V max (S),
where V max (S) is maximum permitted speed for the coordinate S. Here it is important to note that the length of a train should be taken into account because that constraint should not be violated by any part of a train.

Objective Function

We understand the term optimal first of all as energy efficient. Nevertheless, the time of the trip should be also included into consideration by an optimization algorithm. Otherwise, the best chosen solution will always Figure 2: Feasibility check and analytical improvement implement a trip with longest duration. Therefore, we include in the objective function not only an energy related term, but also a time related one [START_REF] Bocharnikov | Optimal driving strategy for traction energy saving on DC suburban railways[END_REF]:

O = O time + O energy → min, O time = (T advise -T ) 2 (∆ -+ ∆ + ) 2 , O energy = E E f lat-out + B B f lat-out , ( 3
)
where O is objective function, O time and O energy are terms, estimating correspondingly the time and energy optimality.

It is evident that T = T advise provides the minimum for O time , see fig. 1. By default T advise = T schedule -∆ - 2 , but it is to be adjusted by the algorithm of time buffers optimization. This special form of O time gravitates the solution to T advise , not to T schedule + ∆ + as it would be if only the energy related term in the objective function would be considered.

The energy optimality component consists of two terms. The first is consumed energy E divided by the E f lat-outenergy spent during the flat-out trip, i.e. intense mode with minimal trip time. The second one is estimation of brake energy B divided by the one spent during the flat-out mode B f lat-out . Adding of that component to the objective function allows us to minimize the number and duration of brake modes.

Energy spent in flat-out mode is considered as maximum possible. Although, technically, one can exceed it by using an enormously inefficient solution. Nevertheless, flat-out mode gives us a perfect tool for simple normalization of the energy component of the objective function. The typical schedule for the trip contains several fixed time points. Some points demand the stop, others -only passing the station within particular time interval defined by tolerances. In order to control the passage time for each point, the optimization process should be separately implemented between all of them instead of considering only points with stops.

Optimization Algorithm

For the case when the schedule is unbalanced (too large difference between the route segments' time buffers) we need some preliminary process of estimation and rearranging of the time buffers. For that purpose we move the parameter T advise forward and backward within the tolerance interval for each trip part between two fixed time points.

At the first stage of the optimization process several typical solutions are analysed:

• flat-out solution,

• flat-out solution with cruising modes replaced by coasting modes,

• solution with maximum possible coasting mode duration (maximum time consumption estimation).

If the time buffer after time tolerances optimization is negative, an optimization will not be started and the flat-out mode will be taken as the solution for current trip interval. Otherwise, the time consumption for the next typical solution is checked. If the time buffer is negative, the second typical solution will be taken as a result. If there is a time for improvements, the last check is implemented and if the time buffer for the third typical solution is still positive, then additional brakes will be added to the control sequence during the optimization. Then, the optimization process based on the genetic algorithm (GA) is started, see its main stages fig. 3. The first step is filling of the initial population which consists of the randomly generated solutions. During that process not only controls u i are generated randomly, but also the moment of switching the command, under assumption that the minimal time between switching the commands equals to 4 sec. This allows to analyse the whole search space of the solutions. Most of the randomly generated solutions are infeasible. But simply discarding the unsuitable solution leads to enormous computation time due to the low value of the relation "feasible solutions / all generated". Therefore a special procedure is necessary to generate feasible random solutions more frequently.

For that purpose, a feasibility check process corrects the solutions if the velocity constraints V max (fig. 2,cases A,B) or V min (fig. 2, case C) are violated. For the V max violation the coasting or brake is added to the control sequence. In case if V min is reached the traction mode is switched on. Secondly, the feasibility check process tracks the time consumed by a randomly generated solution. If the resulting time misses the tolerance interval [T schedule -∆ -; T schedule + ∆ + ], then some random controls are added to the control sequence till the solution becomes feasible or the allowed number of attempts is over. During that process the flat-out solution with minimal time consumption allows to have an estimation of the minimal time needed to get to the destination with current velocity constraints. Finally, an algorithm of analytical improvement attempts to optimize new solutions using several heuristics, see an example in the fig. 2, case D. Starting from the shortest switches, it iteratively changes the control sequence parts with frequent re-switching traction-brake to coasting mode while the solution is feasible and only if the fitness function has not gotten worse due to that improvement.

As a result, the constructed population is a number of feasible and analytically enhanced solutions represented by the sequence of commands for the train. Then, using a train model described in Section 3.1. we estimate the energy and time consumption for each solution in population. The objective function (or fitness function in terms of the genetic algorithm) is computed according to eq. ( 3).

The next step is sorting of the solutions by fitness and selecting of the best or elite members in the population to save them till the next iteration while the worst of them are to be deleted. After that the freed place in the population is filled by the solutions obtained by one-point crossover, see fig. 4. Here a random mutation can happen, which is in fact an addition of a random control: cruise, traction or coasting or replacing one of the previously defined controls, see fig. 6. If the time buffer was too big, together with these commands, an additional brake can also be added during the mutation. Mutation by replacing a control New "children" are checked and corrected by the feasibility check algorithm. From this point the algorithm repeats the iteration until the necessary number of steps is reached.

Further, we continue the optimization for the next trip part and repeat these steps.

At the end of a process we have a set of the solutions which provide the minimum for the objective function. In other words, we have energy efficient train control sequence.

Results

In order to implement the optimization algorithm described above a C++ program was developed. Estimated time of computations for a 320 km trip is about 40 seconds with Intel Core i5 @ 2.77 GHz, 8 Gb RAM.

An example of the solution for the first 15 minutes of a trip with a real schedule is depicted in the Fig. 7,8. Despite the fact that the schedule Figure 7: Velocity profile and control modes is tough, this solution uses the sequences "traction mode → long coasting mode" one after another. Brakes are used only for satisfying the velocity constraints, which is very reasonable. The train is late at the first station (the cross is to the right from the box), nevertheless, it arrives early at the next one (the cross is to the left from the box) without use of the flat-out mode. Cruise mode is not used here at all. Track profile advantages are utilized to hold the speed during the coasting mode: see the interval between 60 and 210 seconds. Here the algorithm does not choose to use brakes or traction, despite the significant change of the speed. The time interval between 420 and 480 seconds also worth mentioning. An algorithm avoids use of a short traction modes here, replacing that with one traction turned on at the end of the 100 km/h constraint and turned off at the middle of the next constraint -160 km/h. During this speed up the train perfectly fits all the requirements, taking into account its own length (see the green line).

In the optimization process 20 iterations were used. The greater number of iterations would be inefficient because the subsequent steps will lead to diversity reducing in the population whereas no significant improvement of the objective function will be earned. The population size was chosen as a compromise between the computation time and accuracy of the results, i.e. it shall provide the same final solution after multiple experiments. The 200 members population satisfies both requirements.

Conclusion

The developed algorithm of automated train control produces an intelligent control sequence which provides an energy efficient trip. The track profile inclines are taken into account, consumed and recuperated energy is computed. There is a capability of consideration of weather conditions and any other additional constraints and circumstances.

We still consider only two possible positions of the traction and brake lever -0 and 100% intensity. This constraint seems to be reasonable according to the hypothesis that the energy efficient driving means that traction is used for short periods of time but with 100% intensity, where it is possible. However, consideration of several intermediate levels is one of the directions of improvement.

The design of the autopilot with a feature of energy optimality is under discussion and consideration for now. Rough theoretical estimations have shown savings of consumed energy up to 14.5% as opposed to manual driving. Actual values will be computed after the first field trials which are planned for the near future. The authors are looking forward to any reviews and comments.

Figure 1 :

 1 Figure 1: Time component of the objective function O time

Figure 3 :

 3 Figure 3: Genetic algorithm chart

Figure 4 :

 4 Figure 4: Crossover of two solutions

Figure 5 :

 5 Figure 5: Mutation by adding a control

Figure 6 :

 6 Figure 6: Mutation by replacing a control

Figure 8 :

 8 Figure 8: Forces acting on a train


Tom Kristensen

Ludovic Pricoupenko

One-dimensional ultracold atomic gases: Impact of the effective range on integrability

Keywords: numbers: 05.30.Fk, 05.30.Jp, 03.75.-b, 67.85.-d

Three identical bosons or fermions are considered in the limit of zero-range interactions and finite effective range. By using a two channel model, we show that these systems are not integrable and that the wave function verifies specific continuity conditions at the contact of three particles. This last feature permits us to solve a contradiction brought by the contact model which can lead to an opposite result concerning the integrability issue. For fermions, the vicinity of integrability is characterized by large deviations with respect to the predictions of the Bethe ansatz.

I. INTRODUCTION

The one-dimensional (1D) Bose and Fermi gases with zero-range interactions are celebrated examples of exactly solvable many-body problems [1][2][3][START_REF] Sutherland | Beautiful Models, 70 Years of Exactly Solved Quantum Many-Body Problems[END_REF]. Ultracold atoms offer the stupendous possibility to achieve these systems in the degenerate regime by using highly elongated cigar traps [START_REF] Bloch | [END_REF][6][START_REF]Proceedings of the school -Quantum Gases in Low Di-mensions J. Phys. IV[END_REF]. Moreover, using magnetic Feshbach resonances and/or tuning the trap parameters make it possible to study 1D systems in strongly correlated regimes. This way, the Tonks-Girardeau and the super Tonks-Girardeau phases have been achieved [START_REF] Girardeau | [END_REF][9][10][11][12]. In addition, the existence of confinement induced resonances and resonances shifts have been confirmed [13,14]. Analogously to the three-dimensional (3D) case [15][16][17], it has been shown recently that considering a large 1D effective range parameter permits one to enrich the phase diagram of the Bose and Fermi 1D degenerate gases [18,19]. For bosons, this regime is achieved for narrow resonances i.e., in the limit of small Feshbach coupling between atoms and diatomic molecules [20]. For fermions, the p-wave scattering resonance is intrinsically narrow and this regime can be easily reached [21]. As in 3D systems, for small energy processes, the regime of large effective range can be studied in the limit where the actual radius of the interacting potentials is formally zero whereas the effective range and the scattering length of the model are finite. Using a contact model (CM), it was shown for one-component fermions in Ref. [18] and for identical bosons in Ref. [19] that the eigenstates of these systems are given by the Bethe ansatz (BA) and are thus integrable in the limit of large effective range [22]. This result is in strong contradiction with the McGuire-Yang-Baxter criterion which when applied in this regime shows diffractive effects in multiple scattering [2,[START_REF] Yurovsky | Advances in Atomic Molecular and Optical Physics[END_REF]. Hence integrability is inherently not possible.

In this paper, we consider the three-body problem, which is intimately related to the integrability issue [START_REF] Sutherland | Beautiful Models, 70 Years of Exactly Solved Quantum Many-Body Problems[END_REF]. To this end we use a Hamiltonian two-channel model (HTCM), which encapsulates the Feshbach mechanism. Whereas the CM and the HTCM are strictly equivalent at the two-body level, in the three-body problem the HTCM gives large deviations with respect to the predictions based on the BA. We show that in the limit of the contact of three particles, all the solutions of the HTCM have the same type of singularity not satisfied by the BA. The behavior of the wave function in the limit where the three particles fall one on top of the other appears then as a key ingredient in the violation of the integrability. We show that equivalence of the CM and of the HTCM can be achieved at the three-body level by imposing continuity conditions on the wave function.

II. TWO-BODY PROBLEM AND MCGUIRE-YANG-BAXTER CRITERION

Our modeling of the system is based on a parameterization of the two-body 1D asymptotic scattering states including the effective range term. For an incoming wave of relative wave number k 0 and relative coordinate z, we write it as

z|ψ k0 = e ik0z + f 0 (k 0 ) + sgn(z)f 1 (k 0 ) e ik0|z| . ( 1
)
In Eq. (1), sgn(z) is the sign function and f 0 (f 1 ) is the scattering amplitude in the even (odd) sector, parameterized as

f η (k 0 ) = -(ik 0 a η ) η 1 + ik 0 a η + b η (ik 0 ) 3-η (-a η ) η . (2)
For ultracold atoms in a 1D waveguide, the scattering lengths a η and the effective range parameters b η in Eq. ( 2) can be expressed as a function of 3D scattering parameters in the homogeneous space [START_REF] Olshanii | [END_REF][26][27][28][29][30][31][32]. In what follows, we consider only positive values of the effective range parameter b η > 0, an assumption justified in the limit of narrow resonances [20,21]. From the analyticity of the scattering amplitude, one finds a single bound state, i.e., a dimer |φ η of energy -2 κ 2 η /m in the even sector for all values of a 0 and in the odd sector only for positive values of a 1

z|φ η = [1 -2ηθ(z)] × κ η 1 + 2b η κ 3-2η η × e -κη|z| . (3)
In Eq. (??) θ(z) is the Heaviside function and in the odd sector (η = 1), one recognizes the sign function sgn(z) = 2θ(z) -1. The dimer binding wave number κ η in Eq. ( 3) is the positive root of

1 -a η κ η -b η (a η ) η κ 3-η η = 0. ( 4
)
We now come to the integrability issue for a system of N one-component bosons (fermions) where the two-body scattering occurs only in the even (in the odd) sector [33]. Integrability means that the eigenstates are given by the BA and there is thus no diffractive scattering i.e., the wave numbers of the particles are globally conserved after multiple collisions in the system [START_REF] Sutherland | Beautiful Models, 70 Years of Exactly Solved Quantum Many-Body Problems[END_REF]. The expressions of the transmission (t η ) and reflection (r η ) coefficients in the scattering of two identical particles are thus particularly relevant. They are defined by an alternative expression of the asymptotic scattering state in Eq. ( 1) where the interaction occurs only in one of the sectors η = 0 or η = 1:

z|ψ k0 = e ik0z + r η (k 0 ) × e -ik0z for z < 0 t η (k 0 ) × e ik0z for z > 0 (5)
From Eq. ( 1), considering a pair of particles (i, j) of wave numbers (k i , k j ), the transmission and reflection coefficients are related to the exchange of momentum between the scattering particles with:

t ij η = 1 + f η k i -k j 2 (6)
r ij η = (-1) η f η k i -k j 2 . (7)
A necessary condition for integrability is given by the McGuire-Yang-Baxter criterion, which follows from the absence of diffractive scattering in the three-body integrable problem [2,[START_REF] Yurovsky | Advances in Atomic Molecular and Optical Physics[END_REF]:

r 12 η r 13 η t 23 η + t 12 η r 13 η r 23 η = r 23 η t 13 η r 12 η . (8)
For η = 0 this last equality is verified if and only if b 0 = 0 (Lieb Liniger model) and for η = 1 if and only if a 1 = 0 or |a 1 | = ∞ i.e., in the Fermi Tonks-Girardeau (FTG) regime [36,37]. This is in strong contradiction with the results of Refs. [18,19] where the BA was used as an eigenstate of contact models in regimes where Eq. ( 8) is not satisfied. To understand this discrepancy, in the rest of this paper we focus on the three-body problem which has the advantage of the simplicity while being a cornerstone of the integrability.

III. CONTACT MODEL A. Contact conditions and pseudo-potentials

We first use a CM which includes the effective range as a straightforward generalization of the Lieb Liniger model and introduced in Refs. [20,21]. It is analogous to the one used in the context of narrow Feshbach resonances for atoms in the three dimensional space [15,38,39]. For convenience we introduce the shorthand notations (z) N for the N coordinates of the system and (z ij ) for the relative coordinate of the pair of particles (i, j):

(z) N ≡ (z 1 , z 2 , . . . z N ), z ij = z i -z j . (9)
The center of mass of the pair (i, j) is denoted (Z ij ), and the relative distance between the pair and the particle (k) is denoted (Z ij,k ):

Z ij = z i + z j 2 ; Z ij,k = z i + z j 2 -z k ( 10
)
The CM is defined as follows: firstly, for all the configurations where ∀i = j, z i = z j , the wave function (z) N |Ψ verifies the Schrödinger equation without any interaction between particles; secondly, for each pair of interacting particles (i, j) the wave function verifies the contact conditions

lim zij →0 + 1 + a η ∂ zij + (-a η ) η b η ∂ 3-η zij (z) N | Πij η |Ψ = 0 ( 11
)
where for η = 0 (for η = 1) the operator Πij η symmetrizes (antisymmetrizes) the state |Ψ in the exchange of the particles i and j:

(z) N | Πij η |Ψ = 1 2 [ z 1 , . . . z i , . . . z j , . . . z N |Ψ +(-1) η z 1 , . . . z j , . . . z i , . . . z N |Ψ ] . (12)
In Eq. ( 11), the positions Z ij and z k (where k = i, j) are kept fixed [40]. One can verify that the exact expressions of the scattering amplitudes in Eq. ( 2) are deduced from the contact conditions of Eq. ( 11) by using the wavefunction of Eq. ( 1).

Another equivalent way to implement the contact model is to include directly the contact condition in the Schrödinger equation by using the Λ potentials V ij η (Λ) for each pair of interacting particles (i, j). For a pair of particles of reduced mass µ, in the even sector of the interaction:

(z) N |V ij 0 (Λ)|ψ = - 2 µ f 0 (iΛ)δ(z ij ) × lim zij →0 + Λ + (1 -Λ 3 b 0 )∂ zij + Λb 0 ∂ 3 zij (z) N | Πij 0 |ψ (13)
and in the odd sector of the interaction:

(z) N |V ij 1 (Λ)|ψ = - 2 µ f 1 (iΛ) Λ δ (z ij ) × lim zij →0 + Λ + Λ 2 b 1 + ∂ zij -b 1 ∂ 2 zij (z) N | Πij 1 |ψ . ( 14
)
In Eqs. ( 13) and ( 14), Λ is an arbitrary parameter, i.e., the action of the pseudo-potential on exact eigenstates do not depend on the value of Λ [21,41,42].

B. Puzzling result for the trimers

In the regime where a dimer exists and if the system is integrable, then the ground state for three identical particles of mass m is a trimer of energy -4 2 κ 2 /m given by the BA [2]

z 1 , z 2 , z 3 |ψ p,BA η = e -κ i<j |zij | i<j [1 -2ηθ(z ij )] . (15)
Following the standard method in Refs [1,43], one considers the contact condition in Eq. (11) for each pair (i, j) in configurations where the third particle k is distinct from the center of mass Z ij (i.e., z k = Z ij ). For instance, in the case where z 2 < z 1 < z 3 :

z 1 , z 2 , z 3 |ψ p,BA η = e -κ(z12-2Z12,3) . ( 16
)
Applying the contact condition for the pair (1, 2) on Eq. ( 16) gives κ = κ η . The same reasoning for the other configurations give the same result. Moreover, the wave function in Eq. ( 15) is a solution of the free Schrödinger equation almost every where excepted at the contact of two or three particles. Thus surprisingly, the BA for the trimer appears as an eigenstate of the CM with the binding wavenumber q t,BA η = 2κ η , in deep contradiction with the McGuire-Yang-Baxter criterion.

IV. THE CONSISTENCY OBTAINED FROM THE TWO-CHANNEL MODEL A. Two-channel Hamiltonian

The consistency of the CM is thus puzzling and to go further we now use a HTCM which is a more conventional approach. In this model, the scattering process between two particles is only due to the coherent coupling between the pair of particles and a molecular state of mass 2m. For a plane wave of wave number k, we choose the convention z|k = exp(ikz) and we denote the creation operator in the open channel â † η,k , where η = 0 for bosons and η = 1 for fermions. The creation operator for molecules in the closed channel is denoted by b † η,k , where the index η permits one to distinguish the composite boson (i.e., the molecule) made of two fermions, from the molecule made of two bosons. We consider only pure systems with identical particles and for each system (η = 0 or η = 1), the Hamiltonian is

Ĥη = dk 2π k â † η,k âη,k + k 2 + E m η b † η,k bη,k + 2 λ η m dKdk (2π) 2 k|δ η â † η, K 2 +k â † η, K 2 -k bη,K + h.c. ( 17
)
In Eq. (17

) k = 2 k 2
2m is the single particle kinetic energy, λ η is the strength of the coherent coupling between the two channels and E m η is the internal energy of the molecular state. The function k|δ η in the second line of Eq. ( 17) is a cut-off for the inter-channel coupling

k|δ η = (ik) η e -k 2 2 /4 . (18)
Physically, the short-range parameter represents the length scale below which the collisional properties have a 3D character. For atoms moving in the monomode regime of a 1D harmonic waveguide of atomic frequency ω ⊥ , it is typically of the order of the transverse length a ⊥ = /(mω ⊥ ). At this scale the 1D effective model of Eq. ( 17) is no more relevant. This explains the fundamental interest of considering the zero-range limit ( → 0) which permits one to capture the universal 1D properties for energies much smaller than the level spacing in the waveguide i.e., 2 /(ma 2 ⊥ ). In the zero-range limit, the scattering lengths and the effective range parameters of the HTCM are given by

a 0 = mE m 0 2 |λ 0 | 2 ; a 1 = →0 1 2 π 1 - mE m 1 2 |λ1| 2 ; b η = 1 |λ η | 2 . ( 19
)
The molecular energy in the odd sector (E m 1 ) is a bare parameter which diverges in the zero-range limit in such a way that a 1 keeps a desired finite value, whereas the parameters E m 0 and λ η stay finite in this limit.

B. Integral equation of the three-body problem

In the HTCM, a three-body state is the coherent superposition of a particle state (denoted by |ψ p η ) in the open channel and of a mixed channel state (denoted by √ 3!|ψ m η /λ η ). In the center of mass frame, it can be written as

|Ψ η = dkdK (2π) 2 k, K|ψ p η √ 3! â † η, K 2 +k â † η, K 2 -k â † η,-K |0 + dK (2π) K|ψ m η √ 3! λ η b † η,K â † η,-K |0 . ( 20
)
For a positive energy (E > 0), |Ψ η is a scattering state and we denote the three-particle incoming state by |ψ 0 η . In Eq. ( 20) k, K|ψ p η is symmetric (for η = 0) or antisymmetric (for η = 1) in the exchange of two particles i.e., in the transformation (k → -k) or in the transformations (k → k ± ; K → K ± ) where

k + = -3K 4 + k 2 K + = -K 2 -k k -= -3K 4 -k 2 K -= -K 2 + k. ( 21
)
The projection of the Schrödinger equation at energy E onto the open channel gives

(E - 3 K 2 -2 k ) k, K|ψ p η = 2 2 m k|δ η K|ψ m η + -k + |δ η K + |ψ m η + k -|δ η K -|ψ m η . (22)
Combining Eq. ( 22) with the projection of the Schrödinger equation onto the one atom plus one molecule space, one obtains in the zero-range limit a 1D Skornyakov Ter-Martirosyan equation [44] (

ik rel K ) 2η-1 K|ψ m η f η (k rel K ) + dk 2π M η (K, k, E) k|ψ m η = dk 2π (-ik) η k, K|ψ 0 η . ( 23
)
In Eq. ( 23) we have introduced the kernel

M η (K, k, E) = 4(k + K/2) η (K + k/2) η -m 2 (E + i0 + ) + K 2 + kK + k 2 (24)
and the relative momentum

k rel K = mE 2 -3 4 K 2
, where for a negative argument of the square root one uses the standard analytic continuation in scattering theory i.e., -q 2 = -i|q|. For a state of negative energy (E < 0), there is no incoming three-particle state (|ψ 0 η = 0) and the prescription E → E + i0 + in Eq. ( 24) can be omitted. For b η = 0, one can deduce from Eq. ( 23) the large momentum behavior (|K| → ∞) of the mixed channel wave function solution of the problem as a function of its value at the contact of the three particles:

K|ψ m η ∼ 8 3(-2) η Z = 0 + |ψ m η + Z = 0 -|ψ m η b η K 4-2η . (25)

C. Mapping with the contact model

In the zero-range limit, the cut-off function in the configuration space z|δ 0 [ z|δ 1 ] tends to the δ distribution denoted by δ 0 (z) [to its first derivative δ (z) denoted by δ 1 (z)]. Hence Eq. ( 22) gives the singular behavior of the particle wave function at the contact of two particles. For a pair (i, j) located at the distance Z ij,k = Z ij -z k from the third particle k, one has

∂ 2 zij (z) 3 |ψ p η = 2δ η (z ij ) Z ij,k |ψ m η + "non δ terms". ( 26
)
To achieve the mapping with the contact model, one first imposes that the three-body wave function coincides with the particle wave function of the two-channel model. The Schrödinger equation obtained with the Λ potential must also coincide exactly with the zero range limit of Eq. ( 22). Thus, the action of the Λ potential on a three-particle state verifies

(z) 3 | V ij η (Λ)|ψ p η = 2 2 m δ η (z ij ) Z ij,k |ψ m η ∀Λ (27)
By using the particular expressions of the Λ potentials of Eq. ( 13) and Eq. ( 14), in the limit Λ → 0 for η = 0 and in the limit Λ → ∞ for η = 1, one also finds: a result which does not depend on the value of Λ for states which are solutions of the Schrödinger equation in the contact model. In the three-body problem and for the contact model, the mixed channel wave function is thus obtained from the three-particle state with the following relation:

(z) N | V ij η (Λ)|ψ = 2 2 m δ η (z ij ) lim zij →0 + ∂ 1-η zij (z) N | Πij η |ψ , (28)
Z ij,k |ψ m η = lim zij →0 + ∂ 1-η zij z 1 , z 2 , z 3 | Πij η |ψ . (29)
From Eq. ( 29), in the case of the BA, we denote the mixed channel state |ψ m,BA η and using Eqs. ( 29) and ( 15), one obtains

Z|ψ m,BA η = -κ 1-η η [1 -2ηθ(Z)] 2 e -2κη|Z| , (30)
a result which can be obtained also directectly from Eq. ( 26) [START_REF]We use the standard identities ∂z|z| = sgn(z) and ∂zsgn(z) = 2δ(z)[END_REF]. Applying a Fourier transform on Eq. ( 30), one finds the momentum representation of the mixedchannel state

K|ψ m,BA η = - 4κ 2-η η 4κ 2 η + K 2 . (31)

D. Trimers

We are now ready to compare the trimers obtained from Eq. ( 23) with the BA. It is clear that for b η = 0, the BA in Eqs. (30) and (31) does not fulfill the correct asymptotic behavior in Eq. ( 25) which confirms the non-integrability. For bosons, this result was found in a model including also the direct particle-particle interaction [START_REF] Yurovsky | [END_REF]. The fact that for fermions, the three-body problem is ill-defined when both b 1 = 0 and the numerator of Eq. ( 25) is not zero, shows also that the BA in Eq. ( 31) can never be an exact solution of Eq. (23) [47]. This can be shown as follows: firstly, for η = 1 and b 1 = 0 at large momentum, Eq. ( 23) is scale invariant and the mixed-channel wave function can be searched as a power law: K|ψ m η ∝ K s ; secondly, the integral in the first line of Eq. ( 23) is definite at least if s < -1; thirdly, implementing the limit of large momentum in Eq. ( 23) one finds K|ψ m η ∝ 1/K unless the numerator of Eq. ( 25) equals zero, which completes the proof. Similarly to the integrable case, we have found numerically that whenever a dimer exists, there exists also one and only one trimer characterized by an even symmetry (i.e., K|ψ m η = -K|ψ m η ). We denote the trimer energy by E t η = -( q t η ) 2 /m. In Fig.

(1) the wave number q t η is plotted as a function of the dimer wave number κ η .

In the limit of large scattering length (a η b

1/(3-2η) η
), the binding wave number of the trimer tends to the integrable limit q t η ∼ 2κ η ∼ 2/a η . For bosons, the convergence is fast and one can verify straightforwardly that for b 0 = 0, Eq. ( 31) is the trimer solution of Eq. ( 23). For fermions, the convergence toward the integrable limit is very slow: one finds the approximate law q t 1 ∼ 2κ 1 [1 + 1.86/ ln(0.466 × b 1 /a 1 )]. The shape of the mixed-channel wave function converges also slowly toward the BA of Eq. ( 31) even for very large values of the ratio a 1 /b 1 . In Fig. 2, the wave function in configuration space is plotted for the ratio a 1 /b 1 = 10 3 and a 1 /b 1 = 10 12 . Even for this last value, there is a clear distinction between the solution of the HTCM and the BA (dotted line) in the vicinity of the three-body contact. In this last region, the deviation with respect to the integrable solution is large due to the discontinuity of the BA at Z = 0 for η = 1 in Eq. (30) [sgn(0) = 0 and thus Z = 0|ψ m,BA 1 = 0]. Nevertheless, our numerical solutions of the HTCM show that the discontinuity is asymptotically recovered in the integrable limit. In the opposite limit of a large dimer wave number, q t η ∼ κ η and the mixedchannel wave function tends to the expected results for a shallow two-body (i.e., atom-dimer) bound state Z|ψ m,BA η ∼ exp(-δq t η |Z|) and δq t η = 2 (q t η 2 -κ 2 η )/3. We find numerically δq t 0 ∼ -2.66/a 0 for a 0 → -∞ (where κ 0 ∼ -a 0 /b 0 ) and δq t 1 ∼ 0.835 × a 1 /b 3 1 for a 1 → 0 + (where κ 1 ∼ 1/ √ a 1 b 1 ).

E. Atom-dimer scattering lengths

We have also solved the dimer-particle scattering problem, for an incoming wave of momentum k 0 and a negative energy:

E = 3 2 4m k 2 0 - 2 κ 2 η m < 0. ( 32
)
The mixed-channel state |ψ m,scat η can be written for an arbitrary overall normalization as

K|ψ m,scat η = 2πδ(K -k 0 ) + 2ik 0 g η (K, k 0 ) k 2 0 + i0 + -K 2 (33)
FIG. 2: Mixed-channel wave function for the fermionic trimer (η = 1). Solution of Eq. ( 23) plotted in the configuration space for (a1 = 10 3 b1) (dashed line) and (a1 = 10 12 b1) (continuous line). Dotted line: BA mixed-channel state in Eq. ( 30) with a normalization chosen for convenience. Inset: same plot where Z/a1 is in logarithmic scale. The cross indicates that the solutions of the HTCM do not vanish at the three-body contact.

and in the configuration space one has for Z = 0

Z|ψ m,scat η = e ik0Z + g η (sgn(Z)k 0 , k 0 )e ik0|Z| . (34)
The comparison of Eq. ( 34) with the generic definition of the 1D scattering amplitudes in Eq. ( 1) and Eq. ( 2) gives the relation between the function g η and the even (a dp η,0 ) or odd (a dp η,1 ) dimer-particle scattering lengths in the bosonic (η = 0) and fermionic (η = 1) cases:

a dp η,η = lim k→0 1 -η ik + g η (k, k) + (1 -2η )g η (-k, k) 2ik (35)
Using Eq. ( 23), the function g η (k, k 0 ) is obtained from the integral equation

(q rel K ) 2η-1 g η (K, k 0 ) (k 2 0 -K 2 )f η (iq rel K ) -P ∞ -∞ dk 2π M η (K, k, E) k 2 0 -k 2 g η (k, k 0 ) -M η (K, k 0 , E) g η (k 0 , k 0 ) 4ik 0 -M η (K, -k 0 , E) g η (-k 0 , k 0 ) 4ik 0 = M η (K, k 0 , E) 2ik 0 , (36)
where P denotes the Cauchy principal value and

q rel K = κ 2 η + 3 4 (K 2 -k 2 0 ). ( 37
)
The four types of particle-dimer scattering lengths are plotted in Fig. (3). In the integrable limit (small values of κ η ), the even scattering lengths (η = 0) diverge, whereas the odd scattering lengths (η = 1) tend toward -2/κ η : a result expected from Ref. [34]. In the opposite limit (large values of κ η ), the trimer is shallow and the even scattering length is thus given by 1/δq t η [see Fig.

V. THE DOMAIN OF THE CONTACT MODEL

After this study of the three-body problem with the HTCM, we point out that making the assumption that for b η = 0:

∂ 1-η Z Z|ψ m η is continuous at Z = 0, (38)
filters out unphysical solutions of the CM. For bosons, the continuity of the derivative ∂ Z Z|ψ m η=0 at the threebody contact is a necessary condition to avoid the spurious K -2 behavior at large momentum of K|ψ m η=0 . For fermions and in the even sector, the condition in Eq. ( 38) is necessary to recover Z|ψ m η=1 from K|ψ m η=1 by the inverse Fourier transform for all values of Z, including at the three-body contact Z = 0. In the odd sector, a first-order discontinuity at Z = 0 is incompatible with the K -2 behavior at large momentum of K|ψ m η=1 in Eq. [START_REF] Olshanii | [END_REF].

The continuity condition in Eq. ( 38) excludes the BA from the set of eigenstates of the CM, because it does not belong to the correct domain, whereas it makes it possible to derive Eq. ( 23) from the CM as follows. Using the standard method, one uses in the Schrödinger equation the δ source terms which are related to the two-body singularities of the wave function [44] and correspond to the expressions of the Λ potential in Eqs. (28) and (29). In our case, one obtains Eq. ( 22) in the exact zero-range limit ( = 0). The particle wave function can be then expressed in terms of the mixed-channel wave function with

E + i0 + -3 K 2 -2 k ( 39
)
where k ± , K ± are defined in Eq. ( 21). The Skornyakov Ter-Martirosyan equation ( 23) follows from the application of the contact condition in Eq. ( 11) on this last expression.

VI. CONCLUSIONS

To conclude, we have used the HTCM to study the trimers and the particle-dimer scattering lengths, for three identical particles in one dimension when one takes into account the effective range term. Considering the zero-range limit of the HTCM has been fruitful to define the correct domain of the corresponding 1D CM including the effective range parameter. This way, we conciliate these two different approaches in accordance with the McGuire-Yang-Baxter criterion about the integrability issue. Exploring the phase diagram of the 1D atomic gas from small to large effective ranges is an open issue both experimentally and theoretically. Current experimental techniques make it possible to explore few-and many-body properties in regimes of large effective ranges [12-14, 48, 49]. One expects large deviations from the integrable dynamical properties, observable in the thermalization or in the response functions [50][51][52].

FIG. 1 :

 1 FIG.1: Spectrum of the trimers as a function of the wave number of the dimer. Continuous line: bosons (η = 0), dashed line: fermions (η = 1). Inset: plot of δq t η = 2 (q t η 2 -κ 2 η )/3 in the region where q t η ∼ κη.

  (1)] .

FIG. 3 :

 3 FIG. 3: Dimer-particle scattering lengths defined in Eq. (35) plotted as a function of the dimer binding wave number. Continuous line: bosonic system (η = 0); dashed line: fermionic system (η = 1). (a): Odd sector (η = 1); (b): even sector for bosons; (c): even sector for fermions.

(z) 3 2 e ikz12 e iKZ12, 3 ×

 33 |ψ η = (z) 3 |ψ 0 (ik) η K|ψ m η + (-ik + ) η K + |ψ m η + (ik -) η K -|ψ m η