Pierre Martinon
email: pierre.martinon@inria.fr

Joseph Gergaud
email: joseph.gergaud@enseeiht.fr

Modeling, Optimization, and Control of Dynamic Systems

Keywords: Modeling, Optimization, and optimal control, shooting method, Pontryagin's Principle, control switchings

The shoot2.0 package implements an indirect shooting method for optimal control problems. It is specically designed to handle control discontinuities, with an automatic switching detection that requires no assumptions concerning the number of switchings. Special care is also devoted to the computation of the Jacobian matrix of the shooting function, using the variational system instead of classical nite dierences. The package also features an embedded continuation method and an automatic (parallel) grid shooting in order to reduce the dependency to the initialization.

SHOOT2.0: An indirect grid shooting package for optimal control problems, with switching handling and embedded continuation Résumé : The shoot2.0 package implements an indirect shooting method for optimal control problems. It is specically designed to handle control discontinuities, with an automatic switching detection that requires no assumptions concerning the number of switchings. Special care is also devoted to the computation of the Jacobian matrix of the shooting function, using the variational system instead of classical nite dierences. The package also features an embedded continuation method and an automatic (parallel) grid shooting in order to reduce the dependency to the initialization.

Mots-clés : optimal control, shooting method, Pontryagin's Principle, control switchings

Introduction

The shoot2.0 package implements an indirect shooting method for optimal control problems. It is specically designed to handle control discontinuities, with an automatic switching detection that requires no assumptions concerning the number of switchings. Special care is also devoted to the computation of the Jacobian matrix of the shooting function, using the variational system instead of classical nite dierences. The package also features an embedded continuation method and an automatic (parallel) grid shooting in order to reduce the dependency to the initialization.

1 Algorithmic aspects

Indirect shooting

Direct methods in optimal control rst discretize the problem then solve the resulting nonlinear problem. Indirect methods, on the other hand, rely on the necessary conditions given by Pontryagin's Minimum Principle. These conditions give rise to a boundary value problem that can be solved for instance by shooting methods.

Optimal control problem

As a general framework, we consider an optimal control problem in the Mayer form, in the autonomous case (meaning the dynamics f does not depend explicitly on the time t):

(P )           
M in g(t 0 , x(t 0 ), t f , x(t f )) Objective ẋ = f (x, u) Dynamics u ∈ U Admissible Controls x(t 0 ) = x 0 Initial Conditions c 1 (x(t f )) = 0 T erminal Conditions with x(t) ∈ R n and u(t) ∈ R m .

Remark: a problem with an integral cost t f t0 l(x, u) dt can be expressed in the Mayer form by adding a state variable following the dynamics l.

Optimality necessary conditions: Pontryagin's Principle

We assume in the following that the optimal control is piecewise C 1 , and that we are in the so-called normal case, meaning that the multiplier p 0 associated to the objective is nonzero, and can be therefore set to 1. Let us dene the costate p in the same space as x, and the Hamiltonian H(x, p, u) = (p|f (x, u)).

Pontryagin's Minimum Principle ( [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]): under the assumptions (i) ∃ (x, u) feasible for (P ), with x absolutely continuous and u measurable. (ii) f is continuous with respect to u and C 1 with respect to x.

(iii) g, c 1 are C 1 with respect to x.

Let (x * , u * ) be an optimal pair for (P ), then (i) ∃ p * absolutely continuous such that x * , p * satisfy ẋ * = ∂H ∂p (x, p, u) , ṗ * = -∂H ∂x (x, p, u) (ii) u * minimizes the Hamiltonian ae in [t 0 , t f ].

(iii) ∃ µ 1 such that the transversality condition hold:

p * (t f ) = ∂g ∂x f (t 0 , x(t 0 ), t f , x(t f )) + (µ 1 | ∂c1 ∂x f (x(t f ))) (T C).
Free nal time. If the nal time t f is free, the additional condition

H(x * (t f ), p * (t f ), u * (t f )) + ∂g ∂t f (t 0 , x(t 0 ), t f , x(t f ) = 0 must hold.
For minimum time problems, this gives the condition H(t f ) = -1.

Boundary Value problem and shooting method

In the following we note y = (x, p) and assume that the Hamiltonian minimization gives the optimal control as a function of y u * (t) = ArgM in w∈U H(x(t), p(t), w) = γ(y(t)).

Then the state-costate dynamics derived from the Hamiltonian system can also be written as a function of y ϕ(y) = ∂H ∂p (x, p, γ(y)), -∂H ∂x (x, p, γ(y)) .

We dene the shooting unknown z = p(t 0 ), such that we have y(t 0 ) = (x 0 , z).

If the nal time is free, the shooting unknown is dened as z = (p(t 0 ), t f ). More generally, the software can solve problems for which the boundary conditions at t 0 for (BV P ) can be written under the form y(t 0 ) = y 0 (z) ∈ R n . At the nal time t f , in most cases the multiplier µ 1 can be eliminated in (T C), leading to a set of equations on y(t f ) = (x(t f ), p(t f )). These equations, in addition to the nal conditions c 1 (x(t f )) = 0, constitute the boundary conditions at t f , that we note B 1 (y(t f )). If the nal time is free, t f is an additional shooting unknown and is part of z. The corresponding equation H(x

* (t f ), p * (t f ), u * (t f )) + ∂g ∂tf (t 0 , x(t 0 ), t f , x(t f )) = 0 is then added to the boundary conditions B 1 .
In the general case, we obtain a Boundary Value Problem on y = (x, p) (P ) Pontryagin's Principle

(BV P )    ẏ = ϕ(y) ae in [t 0 , t f ] y(t 0 ) = y 0 (z) Boundary Conditions at t 0 B 1 (y(t f )) = 0 Boundary Conditions at t f
-------------------→ (BV P ) Shooting method --------------→ S(z) = 0 .
This package uses the HYBRD/HYBRJ routines form Garbow, Hillstrom and More ( [START_REF] Garbow | User Guide for Minpack-1[END_REF]) to solve the nonlinear equation S(z) = 0. Two ODE solvers are available for the evaluation of the shooting function, the classical xed step 4th order Runge Kutta method and the variable step DOPRI5 from Hairer and Wanner ( [START_REF] Hairer | Solving ordinary dierential equations. I[END_REF]). This package is designed to handle two specic diculties: control discontinuities (see 1.2 Switching detection) and initialization (see 1.4 Discrete continuation and 1.5 Grid shooting).

First, the evaluation of S and its Jacobian matrix can be tricky when control discontinuities (or switchings) are present. The key is here to detect the switchings during the solving of (IV P ), and use the so-called variational equation instead of nite dierences to compute the Jacobian.

Then, the matter of nding a suitable initial point is addressed by the means of continuation techniques, coupled to an exploration of a grid of initial points. The latter part is obviously better suited to lower dimension problems, which are the usual situation when using indirect methods as opposed to discretization approaches. This is an alternative to multiple shooting, the latter aiming to increase the convergence radius of the shooting method at the expense of an increased problem dimension.

Switching detection

We consider now the case of a discontinuous optimal control with switchings.

A common situation is when the Hamiltonian is linear in the control and the control is bounded. We assume in the following that the optimal control can take two distinct expressions, depending on the sign of a certain switching function ψ. Namely, the Hamiltonian minimization gives

u = γ 1 (y) if ψ(y) < 0, u = γ 2 (y) if ψ(y) > 0
For instance, for a scalar bounded control, u * is either at the lower or upper bound depending on the sign of ψ = ∂H ∂u . This can be generalized for a compo- nent of a control u ∈ R m , or for a number k > 2 of expressions γ i . The zeros of ψ correspond to the switchings of the optimal control, and we assume a nite number of such switchings. We then obtain y(•, z) as solution of the initial value problem with a discontinuous right hand side

(IV P ) disc    ẏ = ϕ 1 (y) if ψ(y) < 0 ẏ = ϕ 2 (y) if ψ(y) > 0 y(t 0 ) = y 0 (z)
Trying to solve this kind of problem with a xed step integrator (e.g. Runge Kutta 4th order) is in practice nearly impossible except for very simple cases, due to the errors at each switching. On the other hand, a variable step integrator (such as RKF45 or DOPRI5) can usually handle the switchings by itself.

However, this costs many very small steps at the switchings, leading to an increased cpu time, and the remaining small errors at each switching can pile up.

As a result, the shooting method is likely to be slower, and also suers from a loss of precision.

An eective way to solve this problem is to detect the switchings during the integration of (IV P ) disc , as described for instance in [START_REF] Hairer | Solving ordinary dierential equations. I[END_REF][START_REF] Mannshardt | One-step methods of any order for odes with discontinuous rhs[END_REF]. We use here the detection method based on the dense output of the ODE solver, ie a cheap (polynomial) approximation of y on each time interval, noted y dense output . This algorithm can be applied to any xed or variable step integration method with a dense output, and is recalled below.

Main integration loop compute y(t + h) compute switching function ψ(y(t + h)) compare the signs of ψ(y(t + h)) and ψ(y(t))

If sign change Then locate switching τ ∈ [t, t + h] by solving ψ(y dense ouptut (τ )) = 0 switch control pursue integration on [τ, t + h] from y dense ouptut (τ )

End

A simple way to solve ψ(y dense ouptut (τ )) = 0 in order to locate the switching time τ is to use a bisection. This bisection stops either at a xed number of iterations or as soon as the interval length is below a xed tolerance. Switching correction. A small aw of the above algorithm is that the switching point y(τ ) from which we resume the IVP integration is inexact, as given by the dense output. The dense output is typically one order less accurate than the actual ODE method, and these small errors can add up in case of a large number of switchings. We can correct the switching point by solving the equation ψ(y ODE (τ )) = 0, where y ODE (τ ) is given by an actual step of integration from t to τ , instead of the dense output. This equation can be solved by a newton method, using the previously located τ as initial point. In this case, the stopping criterion of the bisection for the detection can be less strict, as this point serves only as an initialization. Remark: It is also possible to perform the whole switching detection while using actual integration instead of the dense output. However, it is faster to use the dense output for an approximate detection and then making the correction.

RR n°7380

This minor modication can improve the convergence of the shooting method, as well as the conservation of numerical invariants of the ODE system (such as the Hamiltonian). Graphs below illustrate the dierence between basic switching detection and switching correction.

Jacobian evaluation and Variational system

The most straightforward method to compute the Jacobian of the shooting function is to use nite dierences. This is done in the HYBRD solver with a step of h j = √ |x j |, where is the error on the shooting function evaluation. However, when using a variable step ODE solver, the time steps can dier at the points for the nite dierences. This can impair the approximation of the Jacobian by nite dierences (see [2,[START_REF] Hairer | Solving ordinary dierential equations. I[END_REF]).

Another option to compute the Jacobian is to use the method described in [START_REF] Hairer | Solving ordinary dierential equations. I[END_REF], which is equivalent to the internal dierentiation 1 from [2].

Smooth (C 1 ) case

In the smooth case, we consider the solution y(•, z) of the ODE system

(IV P ) ẏ(t) = ϕ(y(t)) y(t 0 ) = y 0
Then the derivatives Y = ∂y ∂y0 (y 0 ) are solution of the variational system

(V AR) Ẏ (t) = ∂ϕ ∂y (y(t)) Y (t) Y (t 0 ) = I
Here the initial condition is actually of the form y(t 0 ) = y 0 (z), and we are interested in the derivatives of y(•, z) with respect to z. Assuming that we have 1 while external dierentiation refers to the usual nite dierences RR n°7380 y(t 0 ) = [x 0 ; z] for the sake of simplicity 2 , the variational equation for ∂y ∂zj (t, z)

is (V AR) j Ẏj (t) = ∂ϕ ∂y (y(t))Y j (t) Y j (t 0 ) = 0 • • • 1 • • • 0 T
where the right hand side can be approximated by nite dierences if not available analytically

Ẏ (:, j) ≈ 1 h (ϕ(y + hY (:, j)) -ϕ(y)).
Solving the variational system (V AR) along with (IV P ) provides both y(t f , z)

and ∂y ∂z (t f , z). The Jacobian of the shooting function can be computed from the latter, as S(z) = B 1 (y(t f , z)):

Jac S (z) = ∂B 1 ∂y (y(t f )) ∂y ∂z (t f , z).
With a xed step integration, using nite dierences for the right hand side of (V AR) has the same overall complexity as using nite dierences to compute the Jacobian of the shooting function, namely (n + 1) × N steps evaluations of ϕ. This still holds with a variable step method if we assume that solving both (IV P ) and (V AR) for the Jacobian takes roughly the same number of steps as solving only (IV P ) for the shooting function. This can be enforced in the ODE solver, and is recommended in the discontinuous case.

Discontinuous case

In presence of discontinuities, jumps occur in the variational equation, that must be computed at each control switching. Note that this implies the detection of said control switchings, for instance by the method described in 1.2. Assuming that a switching from ϕ 1 to ϕ 2 occurs at τ ∈]t 0 , t f [, we have to perform the update (see [START_REF] Gergaud | Using switching detection and variational equations for the shooting method[END_REF])

Y ← Y + (ϕ 1 (y τ ) -ϕ 2 (y τ + ))τ (y 0 ).
The switching time is determined by the equation ψ(y(τ, y 0 )) = 0, thus by the implicit function theorem, assuming that ∂ψ ∂y (y τ ) = 0:

τ (y 0 ) = - ∂ψ ∂y (y τ )Y (τ, y 0 ) ∂ψ ∂y (y τ )ϕ 1 (y τ )
.

These updates are done automatically with the switching detection routine.

NB. It should be noted that using the variational system without the switching detection leads to a wrong Jacobian matrix, as the jumps described above are not computed. Therefore, it is not recommended to use the variational system when the switching detection is not available. (IV P ), while its Jacobian matrix needs to solve both (IV P ) and (V AR). This usually leads to dierent time steps when using a variabel step ODE solver. To prevent this, we ignore the variables of (V AR) in the formulas for the initial stepsize and the stepsize control. This has been added in the DOPRI5 code as option ITOL=2.

Discrete continuation method

Continuation techniques are commonly used to nd a suitable initialization for the shooting method, starting with an easier problem and progressively going back to the original problem. Let us dene the continuation parameter λ ∈ [λ 0 , λ f ] and a family of problems (P λ ) such that (P λ f ) = (P ) and we know a solution of (P λ0 ). We dene the homotopy

H : R n × [λ 0 , λ f ] → R n (λ, z) → S λ (z)
where S λ is the shooting function associated to the problem (P λ ). The aim of the continuation is to follow the zero path of H from λ = λ 0 to λ = λ f . Note that the existence of such a path is not guaranteed in general, and of course depends on the family (P λ ). For practical purposes, suitable continuations often involve some physical parameters and / or regularization of the original problem. There are several classes of methods to perform this path following, for instance dierential continuation or simplicial homotopy (see [START_REF] Allgower | Numerical Continuation Methods[END_REF][START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF][START_REF] Gergaud | Low thrust minimum-fuel orbital transfer: an homotopic approach[END_REF][START_REF] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation : résolution par homotopie diérentielle[END_REF][START_REF] Martinon | Numerical resolution of optimal control problems by a Piecewise Linear continuation method[END_REF],

and the hampath3 package). We use here a discrete continuation with a linear prediction, which does not require smoothness of the zero path and can be seen a very simple predictor-corrector method.

End if End while

We note (z n , λ n ) the sequence of zeros of the homotopy computed by the algorithm. The simplest way to compute the initialization z init λn+1 for the next shooting attempt is simply to take the solution of the latest successful shooting, ie z init λn+1 = z n .

A little better is a linear prediction based on the two previous solutions

z init λn+1 = z n + λ -λ n λ n -λ n-1 (z n -z n-1 ).
Numerical experiments indicate that this approach is generally superior to the basic initialization by the previous solution. Higher order (quadratic or cubic)

predictions do not seem to improve the path following.

The maximal step for the continuation parameter λ is typically equal to λ f -λ 0 , so that the path following may be completed in only 1 iteration if the rst shooting attempt succeeds. The default value for the normalized maximal step (set in the .init le see 2.3) is therefore 1. One may want to decrease this value if the path following seems to become unstable, for instance with points branching out from the initial path. This can also be used to force the computation of solutions at certain prescribed values of the continuation parameter.

The minimal step is used to stop the continuation if the stepsize becomes too small, without having to wait until the maximal number of iterations is reached.

As the step for λ can only decrease, it is indeed unlikely that the continuation can succeed when this step becomes very small. Attempts to increase again the step after a certain number of successful shootings seem to make the path following more unstable, and yield no overall benet.

The path following may fail to reach the prescribed value for the continuation parameter in two distinct cases. First, the path stops at a certain limit value for the continuation parameter, and further progress is impossible. This case is usually detected by reaching the minimal stepsize for λ, and the path following should stop at the limit value. Second, the path may go away to innity, meaning that one or more components of the solution tend to innity. This case is usually detected by reaching the maximum number of iterations, and can be checked by plotting the continuation path saved in the .contpath le.

Here is a schematic example of a path following for λ ∈ [0, 1].

Iter Step

Theta Shooting Another way to nd a correct initialization is basically to try several ones, which can be automated in order to explore a grid of initial points. This procedure can be a simple way to obtain convergence, and may also detect local solutions, provided the computational cost remains reasonable.

Two aspects make this idea interesting in our context of solving optimal control problems with an indirect shooting method. First, the dimension of the problem is usually small, typically around 10, as it is more or less the dimension of the state variables.

Then, the shooting method tends to be fast when it converges, but also when it diverges. For instance, the newton method usually fails after only one iteration for a starting point outside the domain of denition of the shooting function, which can be quite small in practice (although not precisely known).

Therefore, the overall cost of making many failed attempts may be quite acceptable, allowing for several thousands of shooting attempts in less than one hour on a standard computer.

This procedure can be combined with the continuation approach described earlier. In this case, for each point of the grid we perform a full path following instead of a simple shooting attempt. At the end of the grid all the obtained solutions are sorted by value of the objective, with a count of the successful shootings for each solution. Besides increasing the chance of nding a suitable initial point for the shooting method, grid shooting also allows for exploring the solution space, and can detect multiples or local solutions. This can be useful for real-life applications, for which a local solution may happen to be more suitable for practical use. Grid shooting is also interesting for benchmarking purposes: due to the high sensitiveness of the shooting method to the initial point, collecting convergence results over a large number of initializations is more reliable than using a single test. For instance, the simulations in [START_REF] Cristiani | Initialization of the shooting method via the hamilton-jacobi-bellman approach[END_REF] study the relevance of using the Hamilton-Jacobi-Bellman approach to initialize the shooting method.

We consider the initialization grid dened by the lower bounds z L ∈ R n , upper bounds z U ∈ R n and range vector r ∈ N n . The grid shooting algorithm will perform an automated sequence of shootings for all initial points

z 0 = z L + Σ n i=1 k i z U i -z L i r i , ∀k i ∈ [0, r i ].
Each initial point is dened by the vector of indices k = (k i ) i=1,n , and the total number of grid points is N = Π n i=1 (r i + 1). In practice, we use a single loop instead of embedded DO loops for each dimension. This loop increments the vector of indexes k from k = (0, . . . , 0) to k = r. k=(0,...,0)

While k = r Do Compute grid point z 0 = z L + Σ n i=1 k i z U i -z L i ri
Perform shooting (or full path following for continuation) with z 0

Increment indexes vector k for next grid point carry = true; i = 1

While (carry & i ≤ n) Do If k i < r i Then k i = k i + 1; carry = false Else k i = 0; carry = true; i=i+1

End if End do End do

It is worth noting that there is no need to generate and store the full grid, as each point can be computed directly from the vector of indexes k = (k i ) i=1,n .

Moreover, this grid shooting can be completely parallelized, as each attempt is independent from the others and the order does not matter. A parallel version using OPENMP is available and can benet from multi-core CPUs. Numerical experiments indicate that the total computation time is divided by the number of cores, as expected.

Software overview

This section presents an overview of the main subroutines of the shoot2.0 package, as well as a description of the user-supplied subroutines required to solve a problem, and the list of the input and output les. The code is written in Fortran90/95, and has been tested with the compilers ifort

Times structure

The shoot2.0 package can handle problems with elaborate times structures, comprising xed or free nal time as well as multi-phase dynamics. Control switchings are handled specically by the switching detection method, and therefore do not appear explicitly in the times structure.

Final time

The nal time can be either xed of free, in which case it is part of the shooting unknown, as the last component of z. A free t f implies the additional optimality condition H(x(t f ), p(t f ), u(t f )) + ∂g ∂tf (t 0 , x(t 0 ), t f , x(t f )) = 0, which is added as the last component of the shooting function S(z). The type of nal time is set in the input le (see 2.3): 1 stands for a xed nal time, and -1 for a free nal time.

Multi-phase problems

We refer here as a multi-phase problem a system whose dynamics is governed by several sets of equations, depending on the time. For instance, in the case of trajectory optimization for a multi-stage space launcher, the ight is divided into several phases corresponding to the each propulsion system. Each time interval with a certain set of equation is called a phase, and the times delimiting the phases can be either xed or free. In the latter case, the optimal value for this phase time is part of the shooting unknowns, and the corresponding equation must be dened in the subroutine PhaseCondition.

The current phase is stored in the global variable ipar(19), and can be used in the user supplied subroutines, see 2.4. The phase is typically required to compute the control and/or dynamics, as well as the phase conditions in the free time case. Phase times are specied in the time structure in the input le (see 2.3): a 2 indicates a xed phase time, whose value must be specied; -2 indicates a free phase time, whose value is part of the shooting unknown z.

For instance, a problem with free initial costate, two free phase times t 1 , t 2 and free nal time would give z = [p(0); t 1 ; t 2 ; t f ].

Input and output les

Continuation

Initial, nal value and max normalized step for continuation parameter Initial and nal values λ 0 , λ 1 for the continuation parameter λ, and maximal normalized stepsize for λ during the path following (see 1.4).

Max iterations and iterations output frequency maximal number of allowed iterations for the continuation. output frequency (< 1: no output, n > 0: display every n iteration(s)).

Prediction type prediction mode for the continuation (0: constant, 1: linear).

Output sol and path generate solution le .sol (1: enabled, other: disabled). generate path following le .contpath (1: enabled, other: disabled).

Grid shooting

Lower bounds for grid shooting (see 1.5). Upper bounds for grid shooting (see 1.5). Range for grid shooting (see 1.5).

Problem specic parameters

Parameters size and value of problem specic parameters.

User supplied subroutines

The following subroutines are specic to each problem and must be provided:

-InitialConditions: computes the initialization y(t 0 ) = y 0 (z).

-FinalConditions: computes the boundary conditions B 1 (y(t f )).

-Control: computes the optimal control minimizing the Hamiltonian.

-Dynamics: computes the dynamics for the state and costate variables.

These optional subroutines can be left empty 9 if unused:

-InitPar: specic initializations, such as table interpolation.

-Switch: computes the switching function used for switching detection.

-BangBangControl: control according to the switching function.

-PhaseCondition: equations to be satised by the free phase times.

General notations:

y(nxp) denotes the state-costate pair (x, p).

u(m) denotes the control. -nfree0 is the dimension of the shooting unknown z.

ipar(lipar),rpar(lrpar) are global variables arrays described in 2.5.

InitialConditions

Subroutine InitialConditions(z,dim,y,lipar,ipar,lrpar,rpar) implicit none integer, intent(in) :: dim, lipar, lrpar integer, intent(inout) :: ipar(lipar) real(kind=8), intent(in) :: z(nz) real(kind=8), intent(inout) :: rpar(lrpar) real(kind=8), intent(out) :: y(dim) ... This subroutine denes the initial conditions for the state and costate at initial time, namely y(t 0 ) = (x(t 0 ), p(t 0 )), in the rst nxp components of y. Typically, the state x is set according to x(t 0 ) = x 0 , while the free initial costate is given by the shooting unknown as p(t 0 ) = z. More generally, y(t 0 ) is dened by the initial and transversality conditions at t 0 , with the missing parts coming from the shooting unknown z.

If the variational system is used for the Jacobian, then the derivatives of y(t 0 ) with respect to z must also be lled (the ∂y(t0) ∂z are stored by lines after the rst nxp components of y).

FinalConditions

Subroutine FinalConditions(y,s,dsdy,lipar,ipar,lrpar,rpar) implicit none real(kind=8), intent(in) :: y(nxp) integer, intent(in) :: lipar, lrpar integer, intent(inout) :: ipar(lipar) real(kind=8), intent(inout) :: rpar(lrpar) real(kind=8), intent(inout) :: s(nfree0), dsdy(nfree0,nxp) ... This subroutine computes the value of the boundary conditions at t f for the shooting function. On entry the state-costate pair y = (x, p) at the nal time is provided in y. In the free nal time case, the last component of s contains H(t f ) on entry. On exit the shooting function value s must be lled with the value of B 1 (y(t f )). If the variational system is used to compute the Jacobian matrix of the shooting function, then the derivatives of theses conditions must also be lled in the matrix dsdy. This subroutine is used to compute the optimal control when using the automatic detection of switchings. Switching detection is enabled by setting the corresponding ag to 0 or 1 in the input le .init (see 2.3), and requires the subroutine Switch above that denes the switching function. The user has to provide the two expressions of the optimal control depending on the sign of the switching function. This sign at the current time step is automatically stored in the global variable ipar(17).

Global variables

The two arrays ipar and rpar, of size lipar and lrpar, contain the integer and real global variables for the shoot2.0 package. Most settings and parameters for the algorithms are stored in these arrays, while other parts are for internal use. Below is a short description of the values that can be of use in the usersupplied subroutines.

Integer global variables: IPAR

ipar(2): mode for embedded continuation (0: disabled, 1: enabled) -ipar(3): Jacobian mode (0: Finite dierence, 1: Variational System) RR n°7380

Orbital transfer

We show here an orbital transfer problem studied in [START_REF] Gergaud | Using switching detection and variational equations for the shooting method[END_REF], from an elliptic transfer orbit to the geostationary orbit. We consider a satellite with a low thrust electroionic propulsion, with thrusts ranging from 10 Newtons to 0.1 Newton. The forces applied to the satellite are the Earth attraction and the engine thrust, giving the dynamics

r = -µ r |r| 3 + T m .
The objective is to maximize the payload, i.e., to minimize the fuel consumption during the transfer. Unlike the minimum-time transfer, the optimal trajectories present a bang-bang control, with either full thrust (at apogees and perigees)

or no thrust (see gure). An interesting continuation approach is to go from an energy type criterion to the mass criterion as

M in t f t0 λ|u(t)| + (1 -λ)|u(t)| 2 dt.
The graph below shows the smooth control for λ = 0 and the discontinuous control for λ = 1 (for a 10N thrust).

Orbital transfer -smooth and discontinuous control depending on objective

The solution from the energy criterion is usually sucient to initialize the shooting method and solve the problem for the mass criterion. This continuation is quite eective, but solving the problem for the energy criterion becomes dicult for low thrusts. We test the grid shooting with embedded continuation on the 1N transfer, with dopri5 (tolerances 10 -8 ) as ODE solver. With the simple grid {-0.1, 0.1} n for the initial costate p(0) (ie for each component of p(0) we only try the two values ±0.1, for a total of 128 shooting attempts), we are able to solve the problem for a maximal thrust as low as 0.1N . For thrusts of 10N , 1N and 0.1N , the optimal trajectories present 14, 119 and 1195 control switchings respectively. The graph below shows the optimal trajectory for a 1N thrust, with the thrust arcs (in red) located at apogees and perigees of the orbits, and arcs with no thrust (in green).

RR n°7380

Orbital transfer -Trajectory for a 1N thrust (119 switchings)

We detail now the tests on the 1N transfer to compare several options. For the Jacobian, we use either nite dierences (FD) or variational system (VAR).

For the switchings, we use either no detection, basic detection, or detection with correction. The table below sums up the convergence results for these six congurations. For each grid shooting we indicate the success ratio (percentage of successful shootings over the grid) and the best convergence obtained (ie lowest norm of the shooting function).

Jacobian

No detection Detection Correction FD 25% 1.25 10 -5 25% 4.83 10 -8 25% 3.06 10 -9 VAR 87% 6.51 10 -5 89% 2.80 10 -8 89% 2.19 10 -10

Orbital transfer (1N) -grid shooting CV results

Here the interest of switching detection is clear, as well as the eectiveness of the grid shooting with embedded continuation. Overall, we once again observe that the switching detection and correction improve the precision of the shooting method. The variational system, on the other hand, gives better chances of success than the basic nite dierences. We also notice that using the variational system without the switching detection is much slower, which is probably due to the wrong Jacobian (see 1.

For a given value

  of the shooting unknown z, we note y(•, z) the solution of the Initial Value Problem (IV P ) ẏ = ϕ(y) y(t 0 ) = y 0 (z) and then dene the shooting function S that maps z to the value of the boundary condition at the nal time S : R p → R p z → B 1 (y(t f , z)) with p = n if the nal time is xed and p = n + 1 for a free nal time. Finding a zero of S gives a trajectory that satises the optimality conditions from Pontryagin's Principle. The indirect shooting method thus consists in solving the equation S(z) = 0, as summarized below:

  Time steps. Computing the shooting function and its Jacobian matrix should use the same sequence of time steps, at least to ensure that the same switchings are detected. The shooting function only requires the solving of 2 ie z is actually the initial costate p(0) RR n°7380

  continuation can be used to explore a range of values for a certain parameters, by setting the lower, upper bounds and maximal stepsize for the continuation parameter accordingly.

  All input and output les for a given problem will use the same prex, with dierent extensions. The algorithm requires a single initialization le, say Problem.init, and typically generates a solution le Problem.sol. Discrete continuation can produce an additional path le Problem.contpath storing the path following. Grid shooting also outputs a Problem.gridcv le which summarizes the convergence results over the grid, as well as Problem-cvxx.sol les for each solution found. The Matlab scripts sol.m and contpath.m can visualize the .sol and .contpath les. The .init le contains the following entries (see appendix A for an example). Shooting method Shooting, continuation and verbose mode Shooting: -1: only compute S(z), 0: basic shooting, 1: grid shooting. Continuation: 0: no continuation, 1: use continuation. Verbose: from -1 (no output) to 2 (most info). State, Control, Switch dimensions dimension n, m of state and control variables, and size of the vector returned by the subroutine Switch (should be set to 1 if unused). Number of arcs and times structure number of arcs (at least 1, see below for arc denition). time structure (at least initial and nal times): 0: (xed) initial time 1: xed nal time, -1: free nal time 2: xed phase time, -2: free phase time. values for xed and free times 7 . Shooting unknown dimension and initial guess size of the shooting unknown 8 , size of initial unknown values. starting point for shooting unknown. Integrator choice, xed steps, relative and absolute tolerances 4: Runge Kutta 4th order, 5: dopri5. number of steps for RK4, tolerances for DOPRI5. Switching detection mode -1: disabled, 0: use switching detection, 1: use switching correction Jacobian mode 0: usual nite dierences, 1: use variational system. Target convergence tolerance (on shooting function norm) for a successful shooting.

1 A. 2

 12 3.93 1106 / 281 3.94 1152 / 293 3.93 VAR 5479 / 1384 3.96 1025 / 261 3.93 1017 / 259 3.93 Orbital transfer (1N) -Total CPU and clock times (s), cpu / clock time ratio All tests were run on a quad-core 3GHz Xeon processor, using the parallel (OPENMP) version of the grid shooting. The times indicate that the actual computation time (clock time) is roughly 4 times smaller than the total CPU time (summed on all cores), as expected on a quad-core CPU. RR n°7380 A Input le and user supplied subroutines A.1 Input le sample.init *** Shooting method *** Shooting, continuation and verbose mode 1 choice, fixed steps, relative and absolute tolerances 5 100 1d-8 1d-8 Switching detection mode (-1: disabled 0: dense output) -1 Jacobian mode (0:FD 1:VAR) 0 Target convergence 1d-4 *** Embedded continuation *** Initial, final value and max normalized step for homotopic parameter 0d0 1d0 1d0 Max iterations and iterations output frequency 100 Problem specific parameters *** 1 Subroutines in le sample.f90 Subroutine InitialConditions(y,lipar,ipar,lrpar,rpar) implicit none integer, intent(in) :: lipar, lrpar integer, intent(inout) :: ipar(lipar) real(kind=8), intent(inout) :: rpar(lrpar) real(kind=8), intent(inout) :: y(nxp) if (contmode == 0) then lambda = rpar(roffset1+1)

  subroutine provides the optimal control that minimizes the Hamiltonian, according to Pontryagin's Minimum Principle. On entry the state-costate pair y = (x, p) at current time is provided in y. On exit the control value must be provided in u. subroutine provides the dynamics for the state and costate variables.On entry the value of the state, costate and control at the current time are given in y,u. If mode=0, the full dynamics ẏ = ( ẋ, ṗ) is required in f. If mode=1, only the state dynamics ẋ is required in the rst half of f. optional subroutine provides the switching function ψ used for the switching detection. On entry the value of the state and costate at the current time are given in y. On exit the switching function must be provided in the rst component of psi. Additional components may be used to store other values that the user would like to compute at each time step for visualization purposes.

	real(kind=8), intent(in) :: y(nxp)
	integer, intent(in) :: lipar, lrpar
	integer, intent(inout) :: ipar(lipar)
	real(kind=8), intent(inout) :: rpar(lrpar) 2.4.3 Control real(kind=8), intent(out) :: psi(npsi)
	Subroutine Control(y,u,lipar,ipar,lrpar,rpar) ...
	implicit none
	integer, intent(in) :: lipar, lrpar
	integer, intent(inout) :: ipar(lipar)
	real(kind=8), intent(inout) :: rpar(lrpar)
	real(kind=8), intent(in) :: y(nxp)
	real(kind=8), intent(out) :: u(m)
	...
	2.4.7 BangBangControl
	Subroutine BangBangControl(y,u,lipar,ipar,lrpar,rpar)
	implicit none
	integer, intent(in) :: lipar, lrpar
	integer, intent(inout) :: ipar(lipar)
	real(kind=8), intent(inout) :: rpar(lrpar) 2.4.4 Dynamics real(kind=8), intent(in) :: y(nxp)
	Subroutine Dynamics(y,u,f,mode,lipar,ipar,lrpar,rpar) real(kind=8), intent(out) :: u(m)
	implicit none select case (ipar(17)) integer, intent(in) :: lipar, lrpar case (-1) integer, intent(inout) :: ipar(lipar) !switching function is negative real(kind=8), intent(inout) :: rpar(lrpar) u(1) = ... integer, intent(in) :: mode case (1) real(kind=8), intent(in) :: y(nxp), u(m) !switching function is positive real(kind=8), intent(out) :: f(nxp) u(1) = ... ... end select
	end Subroutine BangBangControl
	2.4.5 InitPar
	Subroutine InitPar(lipar,ipar,lrpar,rpar)
	implicit none
	integer, intent(in) :: lipar, lrpar
	integer, intent(inout) :: ipar(lipar)
	real(kind=8), intent(inout) :: rpar(lrpar)
	...
	This optional subroutine is called at the beginning of each shooting attempt
	(and not for every shooting function call, unlike InitialConditions). It may
	be used for one-time initializations, for instance related to discrete continuation
	or table interpolations (splines computations).
	2.4.6 Switch
	Subroutine Switch(y,psi,lipar,ipar,lrpar,rpar)
	implicit none

This

This This

RR n°7380

http://apo.enseeiht.fr/hampath/ RR n°7380

, gfortran

and g95

. 4 http://software.intel.com/en-us/intel-compilers/ 5 http://gcc.gnu.org/fortran/ 6 http://www.g95.org/ RR n°7380

http://www.inria.fr ISSN 0249-6399

ipar(19): current phase number (starts at phase 1) -ipar(30): total number of phases, computed from time structure 2.5.2 Real global variables: RPAR rpar(1): objective value (from last call to ivp) -rpar(2): norm of the shooting function (idem) -rpar(3): λ (continuation parameter) -rpar(4): Hamiltonian value (from last call to rhs) -rpar [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]: initial value λ 0 for continuation parameter rpar(13): nal value λ 1 for continuation parameter rpar(roffset1+1:roffset1+npar): optional parameters dened at the end of the .init le, see 2. [START_REF] Cristiani | Initialization of the shooting method via the hamilton-jacobi-bellman approach[END_REF] 3 Illustration problems

Step by step example

We now illustrate the method on a very simple optimal control problem:

The Hamiltonian is dened by

and the costate p is solution of the adjoint equation ṗ1 = 0, ṗ2 = -p 1 .

The optimal control is discontinuous

The unknown for the shooting method is the initial costate and the shooting function is dened by

The user-supplied subroutines and the input le corresponding to this simple example are given in appendix A. The objective M in

a third component for the state, whose corresponding costate is equal to 1.

We test the shooting method on a 50 × 50 grid over [-10, 10] 2 for the initial costate. We compare the two approaches for the Jacobian, nite dierences and variational system, with or without switching detection. We note that depending on the initial costate p(0), we have 9 possible control structures with 0, 1 or 2 switchings, and that convergence is related to these regions. The table summarizes the convergence results for the grid shootings, with tolerances of 10 -8 for the ODE solver. For each run we indicate the percentage of successful shootings over the grid, and the best convergence (norm of the shooting function) obtained. We observe that the switching detection and correction improve the precision of the shooting method, with a better norm. Using the variational system gives better chances of success than the basic nite dierences.

Jacobian

No detection Detection Correction FD 53%


C Fiter

T Floquet

J Rudolph
email: j.rudolph@lsr.uni-saarland.de

Sensorless Control of a Stepper Motor Based on Higher Order Sliding Modes

Keywords: stepping motor, observer, control, stability, robustness, sliding mode

A robust control for a stepper motor with no position nor velocity sensors and only needing current and voltage measurements is designed. Second order sliding mode based observers are realized to estimate both rotor angular position and velocity. Moreover, a robust control law, which is also based on second order sliding modes and which uses the estimates of the observer, is designed. The stability of the observer based control loop is discussed. The results obtained in simulations indicate the usefulness and the robustness of the method.

INTRODUCTION

Stepper motors are widely used in industry and in consumer goods. In order to design a robust control law for such a motor, the knowledge of the angular position and velocity of the motor is required. However, in many applications, the implementation of mechanical sensors measuring these signals is not desirable. This may be due to the space they require, the wires, their fragility, and also because of non negligible additional cost. For all these reasons, during the last years there has been an increasing interest in designing observers and control laws for the stepper motor which rely on measurements of electrical variables only.

One can find work in the literature which deals with the design of a control supposing that the rotor position is known and using a speed observer [START_REF] Chiasson | Nonlinear speed observer for the PM stepper motor[END_REF], [START_REF] Defoort | A third-order sliding-mode controller for a stepper motor[END_REF], [START_REF] Nollet | Observerbased second order sliding mode control laws for stepper motors[END_REF]). Recently, some other work deals with "sensorless" control based on different velocity and position estimation methods, as in [START_REF] Chen | New adaptive sliding observers for position and velocity-sensorless controls of brushless DC motors[END_REF], [START_REF] Poulain | An observer for permanent magnet synchronous motors with current and voltages as only measurements[END_REF], and [START_REF] Xu | An adaptive sliding stator flux observer for a direct-torque-controlled IPM synchronous motor drive[END_REF].

The subject of the present work is resumed in Fig. 1. An observer based on second order sliding modes is presented, which provides accurate estimates of the velocity and the position of the stepper motor and uses electric variables measurements only. Then, a robust control law using the estimated variables, based on second order sliding modes and on the notion of flatness (see Fliess et al. [1992]) is designed in order to achieve trajectory tracking. The performance of this control scheme is illustrated by simulation results showing the robustness of the method, even in the presence of perturbation and noise. The paper is organized as follows. In Section 2, the different models of the stepper motor that will be used for the observation and the control are displayed. The notion of higher order sliding modes is recalled in Section 3. In Section 4, observers are presented, and in Section 5, a control law based on those observers is analyzed. The results obtained are illustrated by simulations.

STEPPER MOTOR MODELS

Model in the fixed frame (αβ)

The following model will be used to design the observer:

di α dt = 1 L (v α -Ri α + KΩ sin N θ) di β dt = 1 L (v β -Ri β -KΩ cos N θ) dΩ dt = 1 J (K(i β cos N θ -i α sin N θ) -f v Ω -τ L ) dθ dt = Ω (1)
in which i α and i β are the coil currents, v α and v β are the voltage inputs, θ is the rotor position, and Ω its velocity. The parameter R represents the resistance, L the inductance, N the rotor number of teeth, K the motor torque constant, J the rotor moment of inertia, f v the coefficient of friction, and τ L the torque load, which will be assumed to be an unknown perturbation.

Only the voltage inputs and the stator currents will be assumed to be measured. Note that the stepper motor without mechanical sensors is observable except for zero velocities.

Model in the moving frame (dq)

From the model presented above, one can derive a model in a moving frame, linked to the rotor, thanks to a simple rotation:

[i d , i q ] T = M p [i α , i β ] T [v d , v q ] T = M p [v α , v β ] T with M p = cos N θ sin N θ -sin N θ cos N θ .
Using this change of coordinates, the system (1) becomes

di d dt = 1 L (v d -Ri d + N LΩi q ) di q dt = 1 L (v q -Ri q -N LΩi d -KΩ) dΩ dt = 1 J (Ki q -f v Ω -τ L ) dθ dt = Ω.
(2)

This model will be used to design the control law in Section 5.

HIGHER ORDER SLIDING MODE

Sliding modes (see [START_REF] Perruquetti | Sliding mode control in engineering[END_REF]) are often used to design robust nonlinear observers or control laws. The aim of sliding mode control is, by means of a discontinuous control, to robustly constrain the system to reach and stay, after a finite time, on a sliding surface where the resulting behavior has some prescribed dynamics. The sliding surface is defined by the vanishing of a corresponding sliding variable s and its successive time derivatives up to a certain order, i.e. the r th order sliding set

S r = x ∈ R n : s(x) = ṡ(x) = • • • = s (r-1) (x) = 0 .
A control law leading to such a behavior is called a r th order ideal sliding mode algorithm with respect to s. Higher order sliding modes, that are characterized by a discontinuous control acting on the r th , r > 1, time derivatives of the sliding vector (instead of the first time derivative in classical sliding mode, r = 1), can reduce the chattering phenomenon while preserving the robustness properties.

A wide range of applications using sliding modes for either observation or control, as in mechanics, robotics, or electric machines, can be found in the literature (see for example [START_REF] Bartolini | A survey of applications of second-order sliding mode control to mechanical systems[END_REF], [START_REF] Butt | Estimation of gasoline-engine parameters using higher order sliding mode[END_REF], [START_REF] Canale | Vehicle yaw control via second-order sliding-mode technique[END_REF], [START_REF] Defoort | Sliding mode formation control for cooperative autonomous mobile robots[END_REF], [START_REF] Drakunov | Stabilization and tracking control for an extended Heisenberg system with a drift[END_REF][START_REF] Floquet | Higher order sliding mode stabilization for a class of nonholonomic perturbed systems[END_REF], [START_REF] Martinez | Hybrid sliding mode-based control of underactuated systems with dry friction[END_REF], [START_REF] Pisano | Cascade control of PM DC drives via second order sliding mode technique[END_REF], [START_REF] Riachy | Second order sliding mode control of underactuated mechanical systems I: local stabilization with application to an inverted pendulum[END_REF]). In what follows, higher order sliding mode observers will be designed because they allow one to obtain estimates in finite time without introducing a low-pass filter (which is not the case with a first order sliding mode, for which it is necessary to filter the high frequency oscillations). Especially, the super twisting algorithm (see [START_REF] Levant | Universal SISO sliding-mode controllers with finite-time convergence[END_REF]) will be used, which is a second order sliding mode defined by

u st (s) = u 1 (s) + u 2 (s) (3) with u1 (s) = -αsign(s) u 2 (s) = -λ|s| 1 2 sign(s) ( 4
)
where s is the sliding variable, and where α and λ are two parameters which have to satisfy sufficient conditions (see [START_REF] Levant | Universal SISO sliding-mode controllers with finite-time convergence[END_REF]) depending on the system and the chosen sliding variable in order to ensure the finite time convergence towards the second order sliding set

S 2 = {x ∈ R n : s(x) = ṡ(x) = 0}.
(5)

POSITION AND VELOCITY OBSERVER BASED ON SECOND ORDER SLIDING MODES

In order to design the observer, consider the two electric equations of model ( 1), and replace the unknown terms by an output injection based on the second order sliding mode algorithm described above. The observer dynamics is described by

d îα dt = 1 L [v α -Ri α ] -χ 1 d îβ dt = 1 L [v β -Ri β ] -χ 2 (6)
The sliding variables are chosen as the current estimation errors:

ε 1 = i α -îα ε 2 = i β -îβ .
and the output injections χ 1 and χ 2 are functions of ε 1 and ε 2 given by the super twisting algorithm (3):

χ 1 = u st (ε 1 ) χ 2 = u st (ε 2 ). ( 7
)
The observation error dynamics is obtained using ( 1) and ( 6):

ε1 = K L Ω sin N θ + χ 1 ε2 = - K L Ω cos N θ + χ 2 . (8)
In order to satisfy the conditions given in [START_REF] Levant | Universal SISO sliding-mode controllers with finite-time convergence[END_REF] that insure the finite time convergence of the estimation errors onto the sliding set defined in (5), one needs to set the parameters α and λ from both super twisting algorithms described in (3) and (4):

α > C 0 and λ 2 > 4C 0 α + C 0 α -C 0 , with C 0 = K L (| Ω| max + N |Ω 2 | max ).
With these conditions, one has ε 1 = ε1 = 0 and ε 2 = ε2 = 0 after a finite time. This yields

χ 1 = - K L Ω sin N θ χ 2 = K L Ω cos N θ. (9)
From the two latter equations, estimates of the rotor position θ and velocity Ω are given by

Ω = ± L K χ 2 1 + χ 2 2 θ =            θ0 if Ω = 0, - 1 N arctan χ 1 χ 2 if Ωχ 2 > 0, 1 N arctan χ 1 χ 2 otherwise ( 10
)
where θ0 is the latest calculated value of θ before Ω = 0.

Two problems caused by the non-uniqueness of the solutions of the electrical system of (1) can be pointed out.

The first is that the sign of the velocity is undetermined, since if (θ, Ω) is a solution of ( 9), (θ ± π N , -Ω) is a solution too. This also means that if the estimate of the sign of the velocity is not correct, the position estimate will be shifted by ± π N . The second problem is that the value of the position estimate is only obtained modulo 2π N because θ is calculated according to (10) within the interval ] -π N , π N [. Assuming the sign of the velocity is known, the observer can achieve good performance even in presence of perturbations. However, the two problems mentioned above are still to be solved.

The idea for solving the sign issue is to study the evolution of the position. To this end, the estimate θ of the position is used and another super twisting algorithm based observer is introduced:

d θ dt = χ 3 ( θ -θ) ( 11
)
where χ 3 is given by the super twisting algorithm.

Writing ε 3 = θθ as the observation error, the error dynamics is:

ε3 = d θ dt -χ 3 (ε 3 ).
Setting the parameters for this other super twisting algorithm N if the estimate of the velocity sign is wrong. In this case, the time derivative of θ is equal to Ω. Then, the output of the third observer Ω converges towards the actual value of the velocity Ω. As a consequence, the sign of the velocity estimate Ω is fully determined, except locally around possible discontinuity points of θ, for which the sliding mode could be destroyed. Those discontinuity points can create estimation errors for this observer. In order to solve this problem, one can use an algorithm to detect these discontinuity points on the sliding variable ε 3 = θθ (see Fig. 2).

α > C ′ 0 and λ 2 > 4C ′ 0 α + C ′ 0 α -C ′ 0 , with
C ′ 0 = d Ω

Fig. 2. Discontinuity point detection and position reconstruction algorithm

Such an algorithm to detect the discontinuity points is by construction sensitive to measurement noise: a "jump" caused by the noise could be interpreted as a passage to another interval by the algorithm. However, the super twisting algorithm used for the position observation in the observer (11) plays the role of a filter and provides a continuous position estimate θ. The results obtained in simulation for a system with large perturbations (τ L = τ h 2 sin(6πt), with τ h the motor holding torque) and measurement noise (white noise -up to 5% of the nominal current) highlight the interest of the method (see Fig. 3).

SENSORLESS CONTROL

In this section, a robust control law without any mechanical sensor is designed using the estimate obtained with the observers defined above.

The control objective is to accurately track a given trajectory. First, note that the system (2) is flat (see Fliess et al. [1992]) with z = [θ, i d ] T as a flat output. This implies that all the system states and inputs can be expressed as functions of a flat output z and a finite number of its time derivatives. Thanks to this property, one can easily define a reference trajectory (denoted Γ r ), satisfying the nominal system dynamics without load torque. The dynamics of the tracking error

e = [i d -i dr , i q -i qr , Ω -Ω r , θ -θ r ] T = [e 1 ,
= 1 L (v d -Re 1 + N L(e 3 e 2 + e 3 i qr + e 2 Ω r )) ė2 = 1 L (v q -Re 2 -N L(e 3 e 1 + e 3 i dr + e 1 Ω r ) -Ke 3 ) ė3 = 1 J (Ke 2 -f v e 3 -τ L )
ė4 = e 3 (12) with vd = v dv dr and vq = v qv qr the new system inputs.

State feedback control

Before designing the sensorless control law, let us recall a second order sliding mode control described in [START_REF] Nollet | Observerbased second order sliding mode control laws for stepper motors[END_REF] that requires the knowledge of the full state.

Two sliding variables are defined as

S θ = ke 4 + ė4 and S i d = e 1 , for k > 0.
It can be seen that the sliding motion S θ = S i d = 0 implies the vanishing of the position tracking error e 4 and the direct current tracking error e 1 .

Because the system has relative degree 2 with respect to S θ , one has to compute the second time derivative of S θ in order to derive the control law used for the tracking of the position:

Sθ = K JL vq + µ θ (e, Γ r ) - k J - f v J 2 τ L - 1 J dτ L dt where µ θ (e, Γr) = k J (Ke 2 -fve 3 ) - fv J 2 (Ke 2 -fve 3 ) - K JL (Re 2 + N L(e 3 e 1 + e 3 i dr + e 1 Ωr) + Ke 3 ).
Then, a control law based on the sampled twisting algorithm (see [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]) is considered (it is well suited to systems with relative degree 2):

K JL vq = -µ θ (e, Γ r ) + u te (S θ , ∆S θ ) (13) u te (S, ∆S) = -λ m sign(S) if S∆S ≤ 0 -λ M sign(S) if S∆S > 0 ( 14
)
where ∆S is the variation of S a sampling period.

To guarantee the finite time convergence in finite time on the surface S θ = 0, λ m and λ M have to verify

λ m > k J - f v J 2 τ L - 1 J dτ L dt max λ M > λ m + 2 k J - f v J 2 τ L - 1 J dτ L dt max (15)
Because the system has relative degree 1 with respect to S i d , the first time derivative of S i d is calculated in order to derive the control law used for the tracking of the direct current:

Ṡi d = ė1 = 1 L vd + µ i d (e, Γ r ) ( 16
)
where

µ i d (e, Γ r ) = 1 L (-Re 1 + N L(e 3 e 2 + e 3 i qr + e 2 Ω r )).
The super twisting algorithm (well suited to systems with relative degree 1) is then used to design the control law 1

L vd = -µ i d (e, Γ r ) + u st (S i d ). (17)
Setting the parameters α > 0 and λ > 0 of this super twisting algorithm ensures the finite time convergence towards the surface S i d = 0.

Therefore, the designed control law guarentees the exponential convergence of the tracking error for a system in which the whole state is measured. It is also shown (see [START_REF] Nollet | Observerbased second order sliding mode control laws for stepper motors[END_REF]) that the system is robust with respect to perturbations and parametric uncertainties.

Observer based control

Suppose now that only the electrical variables are measured and derive a sensorless control law based on the state feedback one described in Section 5.1. Denote the observation error vector as

ε = [i α -îα , i β -îβ , θ -θ] T = [ε 1 , ε 2 , ε 3 ] T .
Then, the dynamics of the whole system state, composed of the tracking error and the observation errors

Ξ = [e 1 , e 2 , e 3 , e 4 , ε 1 , ε 2 , ε 3 ] T is given by Ξ =              - R L e 1 + N (e 3 e 2 + e 3 iq r + e 2 Ωr) - Re 2 + Ke 3 L -N (e 3 e 1 + e 3 i dr + e 1 Ωr) 1 J (Ke 2 -fve 3 -τ L ) e 3 K L
(e 3 + Ωr) sin(N (e 4 + θr))

-K L (e 3 + Ωr) cos(N (e 4 + θr))

e 3 + Ωr              +          1 L vd 1 L vq 0 0 χ 1 (ε 1 ) χ 2 (ε 2 ) χ 3 (ε 3 )          . ( 18
)
Consider the control given in ( 13) and ( 17). The values of the control inputs vq and vd have to be adapted to the sensorless case. For this, the following error between the estimated values of the state and the reference values is introduced:

ξ = [ îd -i dr , îq -i qr , Ω -Ω r , θ -θ r ] T = [ξ 1 , ξ 2 , ξ 3 , ξ 4 ] T . Then, one has K JL vq = -μ θ (ξ, Γ r ) + u te ( Ŝθ , ∆ Ŝθ ) 1 L vd = -μ i d (ξ, Γ r ) + u st ( Ŝi d ) (19)
where Ŝθ and Ŝi d are the new sliding variables (similar to S θ and

S i d ) Ŝθ = kξ 4 + ξ4 Ŝi d = ξ 1 , k > 0
and where μθ (ξ, Γ r ) and μi d (ξ, Γ r ) are given by

μθ (ξ, Γ r ) = k J (Kξ 2 -f v ξ 3 ) - f v J 2 (Kξ 2 -f v ξ 3 ) - K JL (Rξ 2 + N L(ξ 3 ξ 1 + ξ 3 i dr + ξ 1 Ω r ) + Kξ 3 ) μi d (ξ, Γ r ) = 1 L (-Rξ 1 + N L(ξ 3 ξ 2 + ξ 3 i qr + ξ 2 Ω r )).
The superposition principle cannot be used to prove the closed-loop system stability because the system is nonlinear. Therefore, the convergence of both systems (observer/control), considered separately, does not imply the convergence of the whole system.

In order to prove the exponential stability of the closedloop system, it will first be shown that the state is bounded in finite time. Then, since the designed observers converge in finite time, after a finite time the system will behave exactly like the system controlled with the state feedback described in the subsection 5.1 and will hence be exponentially stable.

Let us prove that the system (18) is bounded in finite time.

To do this, write ξ as a function of the system state Ξ, reference trajectory Γ r , and the control inputs χ 1 (ε 1 ), χ 2 (ε 2 ), χ 3 (ε 3 ), u te ( Ŝθ , ∆ Ŝθ ) and u st ( Ŝi d ). One has

ξ 1 ξ 2 = îd îq - i dr i qr .
Since one also has The variable ξ 3 can also be expressed as

îd îq = cos(N θ) sin(N θ) -sin(N θ) cos(N θ) îα îβ with θ = - 1 N arctan χ 1 χ 2 , îα = i α -ε 1 îβ = i β -ε 2 as well as i α i β = cos(N (e 4 + θ r )) -sin(N (e 4 + θ r )) sin(N (e 4 + θ r )) cos(N (e 4 + θ r )) i d i q and i d i q = e 1 + i dr e 2 + i qr then one obtains ξ 1 ξ 2 = Φ(χ 1 , χ 2 ) Ψ(e 4 + θ r ) e 1 + i dr e 2 + iq r - ε 1 ε 2 - i dr iq r (20) with Φ(χ 1 , χ 2 ) =          1 1 + χ1 χ2 2 - χ1 χ2 1 + χ1 χ2 2 χ1 χ2 1 + χ1 χ2 2 1 1 + χ1 χ2 2         
ξ 3 = Ω -Ω r = χ 3 -Ω r . (21)
After substituting in the system (18) the equations of the controls vd and vq ( 19) with the expressions ( 20) and ( 21) of the errors between the estimated values and the reference values, one gets a model depending only on the state and the control inputs:

Ξ = F (Ξ, Γ r , χ 1 , χ 2 , χ 3 , u te ( Ŝθ , ∆ Ŝθ ), u st ( Ŝi d )).
Without loss of generality, it can be assumed that the chosen reference trajectory Γ r is bounded. Then one has, using ( 20) and ( 21)

ξ 1 = O(e 1 + e 2 + ε 1 ) ξ 2 = O(e 1 + e 2 + ε 2 ) ξ 3 = O(ε 3 ). (22)
Since the currents i d and i q are saturated, as well as their estimated values îd and îq , the tracking errors of the currents, e 1 and e 2 , as well as the estimation errors ε 1 and ε 2 are bounded. Hence, ξ 1 and ξ 2 are also bounded according to (22).

It can be seen that the sampled twisting control input ( 14) is also bounded, as well as the super twisting algorithms (7), ( 11) and ( 17) since they are continuous functions with values in a compact set.

As a consequence, the following inequality can be obtained:

Ξ ≤ Q Ξ + ḡ ( 23
)
where Q and ḡ are positive constants. where C is a positive constant.

Thus, the trajectories of the whole system state Ξ are bounded in finite time, which concludes the proof of exponential stability of the closed-loop system (18)-( 19).

In Fig. 4 are reported simulation results obtained with large perturbations (τ L = 2 3 τ h sin(6πt), with τ h the motor holding torque), measurement noise (white noise -up to 5% of the nominal current) and parametric errors (5% on parameters R, L, K, J and f v ). It can be seen in Fig. 4 that good tracking of both position and velocity is achieved, even in the presence of perturbations, noise and parameter uncertainties.

CONCLUSION

Higher order sliding mode observers have been designed to estimate the mechanical state of a stepper motor. Then, a robust sensorless control law using those observers has been presented. The results obtained in simulation for the observation scheme as well as the control are convincing and show the potential of the method.

Fig. 1 .

 1 Fig. 1. Sensorless control scheme for the stepper motor

  time convergence of the observer towards the second order sliding set {ε 3 = ε3 = 0} and thus: θ = θ and Ω = χ 3 (ε 3 ) = d θ dt . It has been seen that one has θ = θ ± π

Fig. 3 .

 3 Fig. 3. Real (in red) and estimated (in blue) values of θ and Ω ė1 = 1 L (v d -Re 1 + N L(e 3 e 2 + e 3 i qr + e 2 Ω r ))

  and Ψ(e 4 + θ r ) = cos(N (e 4 + θ r ))sin(N (e 4 + θ r )) sin(N (e 4 + θ r )) cos(N (e 4 + θ r )) .

Then

  Fig. 4. Reference (in red) and real (in blue) values of θ and Ω with measurement noise and parametric errors

  e 2 , e 3 , e 4 ] T

	Real/estimated position (rad)	-1 0 1 2	0	0.2	0.4	0.6	0.8	1
						Time t (s)
	Real/estimated velocity (rad/s)	-10 0 10 20	0	0.2	0.4	0.6	0.8	1
						Time t (s)
									is then given by