Juan Mario Restrepo

Gary K Leaf

INNER PRODUCT COMPUTATIONS USING PERIODIZED DAUBECHIES WAVELETS

Keywords: Galerkin, differential operators, connection coefficients]

Inner products of wavelets and their derivatives are presently known as connection coe cients. The numerical calculation of inner products of periodized Daubechies wavelets and their derivatives is reviewed, with the aim at providing potential users of the publicly-available numerical scheme, details of its operation. The numerical scheme for the calculation of connection coe cients is evaluated in the context of approximating di erential operators, information which is useful in the solution of partial di erential equations using wavelet-Galerkin techniques. Speciÿc details of the periodization of inner products in the solution di erential equations are included in this presentation. [wavelets;

INTRODUCTION

Wavelets have found a well-deserved niche in such areas of applied mathematics and engineering as approximation theory, signal analysis, and projection techniques for the solution of di erential equations. While wavelets are not conceptually new, [START_REF] Chui | An Introduction to Wavelets[END_REF][START_REF] Daubechies | Ten Lectures on Wavelets[END_REF][START_REF] Meyer | Ondelettes et Operateurs I[END_REF] the past ÿfteen years have produced much of the theoretical underpinnings for the concept, as well as the generation of new wavelet families and the exploration of their potential in various areas of applied science (see References 4 and 5 for wavelets in signal processing, References 6 and 7 for wavelet applications in approximation theory, References 8 and 9 for wavelet applications in the solution of di erential equations, and Reference 10 for other types of applications).

Wavelets have several advantages: (1) they have compact support or exponentially decaying support; (2) their continuity properties may easily be increased, albeit at the expense of a larger domain of support; (3) for a given spline order, an complete basis may easily be generated by simple recurrence relations; (4) in the context of projection techniques, their convergence properties are as good as or better than Fourier methods, and they permit the analysis of extremely local functional behaviour without the need for windowing and with little or no bias from global behaviour; and (5) the manner in which the space is broken down into a family of multiplyenclosed subspaces enables spatial or temporal function multiresolution analysis.

A multiresolution analysis is a nested sequence

V 0 ⊂ V 1 ⊂ • • • ⊂ L 2 (R)
satisfying the following properties:

j∈Z V j = 0 (1) clos L 2 j∈Z V j = L 2 (R) (2)
f(x) ∈ V j ⇔ f(2x) ∈ V j+1 (3)
and there is a function ' ∈ V 0 such that{' 0; k (x) = '(x -k)} k∈Z forms a Riesz basis for V 0 (4)

The term ' is called the mother scaling function since, from (3), there exist {h k } ∈ l 2 such that

'(x) = √ 2 k∈Z h k '(2x -k)
This relation, called the scaling relation, will also hold for '(2x) and, by induction, for '(2 j x):

In accordance with the notation in (4), we denote the translates and dilations of ' by ' j; k (x) = 2 j=2 '(2 j x -k)

The set {' j; k } forms a Riesz basis for V j . We deÿne W j to be the orthogonal complement of V j with respect to V j+1 : Just as V j is spanned by dilations and translations of the mother scaling function, so are the W j 's spanned by translations and dilations of the mother wavelet. The mother wavelet is deÿned by

(x) = √ 2 k (-1) k-1 h -k+1+2M ' 1; k (x)
with M a particular integer. Daubechies 2 constructed compactly supported wavelets and scaling functions using a ÿnite set of non-zero {h k } N -1 k=0 scaling parameters with N =2M and

N -1 k=0 h k = √ 2.
M is the order of the wavelet function. With these scaling parameters, the recursion formulas generate the desired orthogonal wavelets and scaling functions with supp(') = [0; N -1]. Henceforth, these will be the wavelets we shall use, which we refer to as Daubechies wavelets of genus N .

Values for ' are calculated using the scaling relation as indicated in the following procedure. First, the values it takes are determined at integer points. Then at the dyadic rationals at level 1 (the dyadic rationals at level j are D j = {k=2 j } k ∈ Z ); using that information, we calculate the values at D j+1 and so on, until ' is deÿned at the dyadic rationals at all levels. Since the dyadics are dense in the reals, we simply extend ' continuously to R. This procedure creates a function that is continuous but not di erentiable for N = 4, di erentiable but not twice di erentiable for N = 6, and with increasing regularity for increasing N . [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF] The nature of the scaling relation guarantees that ' will be discontinuous in some derivative. [START_REF] Press | Numerical recipes in Fortran[END_REF] This procedure is easily accomplished computationally.

Another pleasing feature of these wavelets is their compact support. Whereas Fourier methods return global results, with compactly supported wavelets one can easily analyse short-lived events or pulses. Wavelet projection methods avoid distortion that might result from a local analysis with a windowed Fourier transform. As will be shown in the next section, compact support also makes the periodization of these wavelets an elegant process.

In addition to items (1)-( 4) mentioned above, wavelets have a number of other interesting properties. These will be given without proof. For further details, see Daubechies 2 or Chui: 1 {' j; k } j¿0; k∈Z is an orthonormal basis for L 2 (R)

(5)

V j+1 = V j ⊕ W j (6)
L 2 (R) = clos L 2 V 0 ∞ j=0 W j (7) {' 0; k ; j; k } j¿0; k ∈ Z is an orthonormal basis for L 2 (R) (8) ∞ -∞ '(x) dx = 1 (9) k∈Z ' 0; k = 1 (10) ∞ -∞ (x)x k dx = 0 : k = 0; : : : ; M -1 ( 11
)
{x k } M -1 k=0 ∈ V N (12)
Item ( 6) is really at the heart of multiresolution analysis and provides wavelet-based analysis with a distinctly di erent resolving quality in contrast to spectral methods: to go to a higher resolution of spatial scale, one simply adds on the next wavelet level (the next W j ) as implied by item (6). At some given level (say, with a representation in V j ), the multiresolution property guarantees all spatial scale information at all coarser levels. In contrast, with Fourier methods, information about one frequency gives no information about other frequencies. This study is devoted to the numerical calculation of connection coe cients involving periodized Daubechies wavelets. Connection coe cients are matrix structures that result from the evaluation of inner products

j; k ≡ ' (d0) j; k0 ' (d1) j; k1 • • • ' (dn) j; kn dx
where d i is the number of di erentiations with respect to x of the scaling function ' = '(x). Inner products arise naturally in the context of the Galerkin solution of di erential equations. 9; 12 The name for these inner products was coined by Latto et al., [START_REF] Latto | The evaluation of connection coe cients of compactly supported wavelets[END_REF] who developed the computational method presented herein. It replaces a quadrature problem by a linear algebra problem. Since the method is simple, fast, and general, it was considered worthy of implementing numerically. The technique developed by Latto et al. is by no means the only alternative available for the calculation of inner products involving wavelets. Of note is Dahmen and Micchelli's 14 algorithm which is similar to the approach favoured in this study in that the quadrature problem is recast as a linear algebra problem. The algorithm is more simply explained in Reference 4 (p. 396 -401). This latter algorithm was numerically implemented very recently by Kunoth and is also publiclyavailable software. [START_REF] Kunoth | Computing reÿnable integrals: documentation of the program, version 1.1[END_REF] The code is capable of performing the inner products of up to four wavelets and their derivatives. Sections 2 and 3 of this study describe how periodized Daubechies wavelets on the unit interval are assembled and present estimates of their approximating characteristics. Furthermore, they will be shown to inherit most of the properties of their more standard Daubechies wavelets cousins, and capable of serving as basis set for L 2 functions on the unit interval which are periodic. Readers interested in the construction of a basis sets deÿned over the semi-inÿnite interval, or with di erent edge conditions over a ÿnite interval, may ÿnd Cohen et al. [START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF] a useful point of reference. Section 4 describes the connection coe cient algorithm. This section also covers several useful strategies for the periodization of general inner product matrices which arise in Galerkin solutions of di erential equations, as well as some properties of the connection coe cients, which are useful in testing numerical codes and results, and in the manipulation of these coe cients.

In Section 5 the connection coe cient algorithm is evaluated by examining how wavelet-Galerkin methods approximate several commonly used di erential operators. Finally, several appedices are provided with information on how to obtain the connection coe cients numerical code as well as further details on the properties and approximating characteristics of the periodized Daubechies wavelets.

PERIODIZED WAVELETS

Daubechies wavelets are deÿned on the whole real line. For many applications, however, wavelets deÿned on a periodic domain are needed. Interestingly, the wavelets deÿned above can be periodized with a Poisson summation technique to give periodic wavelets 2 that possess many of the same properties of their non-periodic kin. Moreover, the periodized ' j; k and j; k are identical to their non-periodic forms except for wrapping around the edges of the domain; for large enough j, this too can be reduced to the non-periodic case for most calculations. Hence, many of the abovementioned properties are preserved in the periodic case as a result of the construction by the "scaling" property of the non-periodic functions and their compact support.

The wavelets are periodized as follows:

' j; k (x) ≡ l∈Z ' j; k (x -l) and j; k (x) = l∈Z j; k (x -l)
By construction ' and are periodic and are well deÿned on [0; 1], since ' and have compact support. Note that ' j; k = ' j; k if k ≡ k mod(2 j ). Thus we shall restrict our attention to 06k ¡2 j . The same holds for the 's. In what follows, the properties of the periodic wavelets will be investigated in some detail. Periodized wavelet bases are not generated in quite the same way as the non-periodic versions. In the non-periodic case, bases are generated by repeated translation and dilation of the mother functions; but this approach is not possible in the periodic case because periodization does not commute with dilation. Therefore, the wavelet must be ÿrst dilated, then periodized. Although proof can be shown for the general case, Appendix II includes a proof that shows instead that the elements in V j for j 6 0 are all constant functions. If dilation commuted with periodization, this would not be true.

Periodized Daubechies wavelets form an orthonormal basis that still has a scaling relation, as shown in Appendix II. This means that in comparison with the non-periodic case, we have a chain of spaces

V 0 ⊂ V 1 ⊂ • • • ⊂ L 2 [0; 1] with the following properties: j¿0 V j = L 2 [0; 1] with V j = span{ ' j; k } 2 j -1 k=0 (13) j∈Z V j = {constant functions} (14) f(x) ∈ V j ⇔ f(2x) ∈ V j+1 (15)
By deÿning W j = span{ j; k } 2 j -1 k=0 we see that W j is the orthogonal

complement of V j in V j+1 . So then clos V 0 ∞ j=0 W j = L 2 [0; 1] (16)
Clearly, di erences exist between the properties of the periodic case and the non-periodic case. While they are both multiresolution spaces, the basis functions in the non-periodic case are all formed by translations and dilations of the mother scaling function, ', while in the periodic case it is often impossible to derive ' j+1 from ' j (for example, consider ' 1 and ' 0 ; the latter is a constant function and thus unable to represent the former). It turns out, however, that there is no relation between ' j+1 and ' j for very small j only. For j suitably large, the periodic case actually looks exactly the same as the non-periodic case. Thus, for large enough j, the periodization will a ect the functions only by 'wrapping them around' the edges of the domain. This is also a strong argument for using the scaling functions as trial-and-test functions when using periodized wavelets. By choosing j well, calculations can be performed as in the non-periodic case; if desired, a multiresolution can then be performed easily. Calculations will not be so simple if V 0 j k=1 W k is used as a test space, because, for low k, the basis for W k is not equivalent to the basis for W k .

APPROXIMATION RESULTS

A function f ∈ L 2 [0; 1] may be projected into the wavelet basis and expressed as f = k a kj ' j; k , where a kj = f; ' j; k : The calculation of a j k is usually performed numerically. In Reference 7, a method was developed which employs Taylor series expansions to approximate f. The method requires use of the moment equations to make O(h n ) approximations for f ∈ C n : Unfortunately, in applications such as signal processing or any area where only a ÿnite number of samples of f are provided, this method may produce unacceptable results. As an alternative, a function f given as samples on D j ∩ [0; 1] may be approximated by a function f ∈ V j .

Deÿne the samples of

f as → f ∈ R, with the kth component of → f = f((k -1)=2 j ): Construct → ' j; k
from ' j; k as shown previously. This yields, for j ¿ log 2 (N -2), a linearly independent spanning set for R 2 j . Further, since supp( ' j;0 ) = [0; 2 -j (N -1), → ' j; k takes only N -2 values on D j , and → ' j; k is just the kth forward cyclic permutation of the elements of → ' j; 0 . The problem is thus reduced to ÿnding the unique representation of → f in terms of { → ' j; k } 0 6 k ¡ 2 j , which is simply the solution of

A → v = → f
where A is the very sparse transformation matrix from the standard basis to the → ' basis. Deÿne

f(x) = k → v (k) ' j; k (x):
By construction, f| Dj = f. In summary, this method involves solving the inverse problem with a sparse matrix, and results in a function f with f| Dj = f: In essence, f will be equal to f at all the sampled values.

We next show the manner in which the periodized wavelets may be used in the context of functional approximation. Since the wavelets form an orthonormal basis, the orthogonal projection operators onto V j and W j are deÿned respectively as

P j (f) = 2 j -1 k=0 f; ' j; k ' j; k ; Q j (f) = 2 j -1 k=0 f; j; k j; k
As we have already seen, periodized wavelets provide a basis for L 2 [0; 1] so we have f -P j f 2 → 0, as j → ∞. This is a property of any orthonormal basis of L 2 , but this particular periodized basis has some additional properties. If f is a continuous function on the torus, then [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF] The proofs for these two results appear in Appendix III of this study. These two last statements are strong, in contrast to the convergence properties of Fourier functional approximations. In fact it has been shown that the continuous functions whose Fourier series do not uniformly converge are dense in C(T). [START_REF] Rudin | Real and Complex Analysis[END_REF] In this sense, wavelets provide a much more general basis than Fourier bases and hence have potentially broader applications. These results suggest that wavelets should do a better job at pointwise approximation, especially for continuous functions.

f -P j f ∞ → 0 as j → ∞. If f ∈ L 1 [0; 1], then f -P j f 1 → 0 as j → ∞.
To measure the error to which a truncated projection will approximate a desired function, we shall estimate its convergence. The natural choice of norms with which to measure convergence is the Sobolev norms. The sth Sobolev norm of a function f is deÿned as

f H s = (1 + 2 ) s f2 ( ) d 1=2
where H s consists of those functions whose s Sobolev norm exists and is ÿnite. Daubechies 2 states that the norm for H s is equivalent to

f H s [0;1] =    j ¿ 0 0 6 k ¡ 2 j (1 + 2 2js ) f; ˆ j; k 2    1=2
Using this result, we easily ÿnd a bound for f

-P p f 2 . For f ∈ H s [0; 1], f -P p f 2 = j ¿ p0 6 k ¡ 2 j f; ˆ j; k 2 = j; k f; ˆ j; k 2 1=2 = j; k 2 2js 2 2js f; ˆ j; k 2 1=2 6 j; k 2 2js 2 2ps f; ˆ j; k 2 1=2 =2 -ps j; k 2 2js f; ˆ j; k 2 1=2 62 -ps f H s [0; 1]
The same technique may be used to find error bounds in the H l norm for l6p:

f -P p f H l [0; 1] =    j ¿ p 0 6 k ¡ 2 j (1 + 2 2jl ) f; ˆ j; k 2    1=2 = j; k (1 + 2 2jl ) (1 + 2 2j(s-l) ) (1 + 2 2j(s-l) ) f; ˆ j; k 2 1=2 6(1 + 2 2p(s-l) ) -1=2 j; k (1 + 2 2sj ) f; ˆ j; k 2 1=2 62 -p(s-l) f H s [0; 1]
These bounds are similar to those on Fourier series. [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]

METHOD FOR COMPUTING CONNECTION COEFFICIENTS

In wavelet applications, one often must represent operators in terms of wavelets. [START_REF] Beylkin | On the Representation of operators in bases of compactly supported wavelets[END_REF] An example of such an application is the Galerkin solution of di erential equations. The formulation of solutions will require integrations of the form d1; d2; : : : ;dn k1; k2; : : : ;kn = ' (d0) j; k0 ' (d1) j; k1 : : : ' (dn) j; kn = 1 0 ' (d0) j; k0 ' (d1) j; k1 : : : ' (dn) j; kn d x where ' (d) = d d '=dx d . This expression is an n-term connection coe cient. Since ' cannot be represented in closed form for N ¿2 and, by construction, has limited regularity, analytic calculation of the integral is impossible, and numerical quadrature is often inaccurate as a result of the wildly oscillating nature of the resulting kernels. An alternative approach developed by Latto et al. [START_REF] Latto | The evaluation of connection coe cients of compactly supported wavelets[END_REF] circumvents some of the di culty by exploiting the scaling relation and the moment condition to reduce the calculation to an eigenvector problem. Their method is designed for non-periodic compactly supported wavelets. However, by invoking an extension of the earlier result regarding the equivalence of periodized and non-periodized wavelets in Section 2, one may infer that for j ¿ log 2 ((N -1)n), the periodized case yields the same result as the non-periodized case. For illustration and for completeness in what follows we adopt closely the general procedure given in detail in Reference 13 to the 2-tuple case. Several tabulated connection coe cients using periodized wavelets are included in an appendix of this study. First, integration by parts is performed repeatedly on the above integral to obtain From these relations it is clear that any 2-tuple can be represented by a d k .

To construct the eigenvector problem, ÿx d, then solve for { d k } 0 6 k ¡ 2 j , creating a system of 2 j homogeneous relations in d k and enough inhomogeneous equations to reduce the dimension of the associated eigenspace to 1. Although we are using the connection-coe cient method for the non-periodized case, we are computing them for the periodic case (by equivalence), which is where the bounds on k come into play.

Connection coe cients algorithm

To generate the homogeneous equations ÿx d; j ∈ N, such that ' (d) j is well deÿned. To simplify notation, denote ' (d) j; k ≡ d k . In Reference 13 it is suggested, without proof, that this method also holds for the ÿrst d for which d is discontinuous. For low-order di erential equations, however, N = 6 or N = 8 wavelets should provide su cient regularity. Since for every 0 6k ¡2 j ,

d k = 0 (x) d k (x) dx = N -1 m=0 h m m (2x) N -1 l=0 h l d l+2k (2x) 2 d d(2x) = 2 d m l h m h l m (2x) d l+2k (2x) d(2x) = 2 d m l h m h l 0 ( ) d l+2k-m ( ) d ; d k = 2 d N-1 m=0 N -1 l=0 h m h l d l+2k-m
In the above discussion, the integration is over the real line. This linear homogeneous system can be represented as

A → d = 2 -d → d where → d = { d k } 0 6 k ¡ 2 j :
It is worth noting here that if one needs to compute an n-tuple connection coe cient for j ¡ log 2 ((N -1)n), then the periodic scaling relation can be used to resolve each ' j; k into the sum of ' j ; k 's with j ¿ log 2 ((N -1)n): Thus, the reduction to the non-periodic case is a universally applicable method.

To generate the inhomogeneous equations, we must ÿrst assume d 6M -1. The moment condition then guarantees that with

M d l = 2 -j(2d+1)=2 M d l
The linear system formed by the 2 j homogeneous equations and the above inhomogeneous equations has eigenspace dimension equal to 1. Thus, all that remains to complete the system is to calculate M for all l.

The linear system is now complete and ÿxes the values of d k . Note that while the scaling relation, which is used to generate the homogeneous relations, exists for periodized wavelets, currently nothing is analogous to the moment condition that may be used to generate the necessary inhomogeneous equations. One possible approach is to use the scaling equation linked with the fact that V j = {constant functions} for j 60 to ÿnd additional inhomogeneous relations. Problems arise with relating ' j+1; 0 dx to ' j; 0 dx; however, since dilation does not commute with periodization. While this method could probably be worked out, the periodic case can always be reduced to an equivalent non-periodic case for which the method is already well deÿned. Thus, to compute an n-term connection coe cient for periodized wavelets, one need only resolve the terms into V j , for some j ¿ log 2 (N -1)n and apply the above method. This approach takes full advantage of the equivalence of periodic and non-periodic scaling functions and circumvents the need for a connection coe cient method particular to periodized wavelets.

Periodization of inner product matrices

The two most common situations that arise in the wavelet Galerkin solution of di erential equations involve 2-tuples and 3-tuples. Periodic conditions require wraparound of the entries in in the upper right-hand corner and the lower left-hand corner, assuming that the matrix row index increases from top to bottom corresponding to the dyadic points in [0; 1). The wraparound will produce at most N (N + 1) additional entries in the matrix, where [0; N ] is the support of the scaling functions. In this section we present explicit schemes that enable the construction of these matrices. The 2-tuple case is easily generated, since is a circulant matrix of the form = circ( 0 ; 1 : : : N ; : : : ; 0; 0; : : : -N ; -N+1 ; : : : -1 ) and = ± T , depending on whether the operator is symmetric or skew-symmetric (see Section 4.3).

The 3-tuple is somewhat more involved. The j ; k have the structure illustrated schematically in Figure 1(a). We call this structure the 'index pad'. Let K = 2 p -1, where the superscript p corresponds to the resolution of the scaling functions. A typical situation in a Galerkin discretization when 3-tuples are involved would be the calculation of

s k = K l=0 a l ! k; l where ! k; l = K j=0 h j k; j; l
where a l and h j are the wavelet coe cients of the projection of real quantities A(x; •) and H (x; •) into V p , and is the inner product d1 k d2 j d3 l . For a given k there is an index pad in ( j; l) over which these indices range when forming ! k; l and s k . Moreover, associated with each index ( j; l) in the pad there is a j ; k , where j = j -k and l = l -k, which gives the value of k; j; l . It is easier to think of the index pad as having a rectangular structure in which the entries lying outside of the hexagon of the true index pad are zero.

First consider the case in which there is no wraparound, that is, for N 6 k¡k , where k = K -(N -1), the largest index k for which the support of the scaling function lies entirely in [0; 1]. Associated with each (j ; l ) in the index pad, there is a connection coe cient j ; l . As k is varied, the pad centred at k, moves along the line l = j in the ( j; l) plane. This is the 'regular case' and is shown in Figure 1(a).

Next, consider the case k ¿ k . Let k = k + , 06 6N -1. The situation is illustrated in Figure 1(b). The index pad is beyond (j; l) = (K; K). For convenience identify the following subregions.

A : k -N 6 j 6K k -N 6l6K B : j = K + 1 + q; 06q 6 l = K -t; 06t 6 N -1 C : j = K -t; 06t 6N -1 l = K + 1 + q; 06q 6
D : j = K + 1 + q; 06q 6 l = K + 1 + q; 06q 6
Consider a given k, the corresponding index pad is subdivided into the four subsets A; B; C; D deÿned above. When index pairs (j; l) lie in subregions B; C or D, periodization comes into play since the corresponding basis functions have been periodically extended. We now consider the consequences of this periodization on the s k in detail. Note that s k will be the sum of three partial sums, summed over subsets of the index pad. Next we observe that the range of l is k -N 6l 6K +1+ . We divide this range into two parts: k -N 6l 6 K and K +16l6K +1+ . For the range k -N 6 l 6 K, the index l is not a ected by the periodization since l 6K. The index j, on the other hand, produces a pair (j; l) that ranges over the set A and B. Let

! k; l = ! r k; l + ! p k; l . Then ! r k; l = K j=k-N h j j ; l (the regular part) ! p k; l = K+1+ j=K+1
h j-K-1 j ; l (the periodic adjustment) and thus

s (1) k = K l=0 a l ! k; l
where, throughout this section, we set j = j -k and l = l -k.

Next we deal with the range K + 16l6K + 1 + . In this case the pair ( j; l) ∈ C ∪ D. When ( j; l) ∈ C,

! k; l = K j=k-N h j j ; l
Now the index l is a ected by periodization, so that

s (2) k = K+1+ l=K+1 a l-K-1 ! k; l
When (j; l) ∈ D, both indices are a ected by periodization, thus

! k; l = K+1+ j=K+1 h j-K-1 j ; l
Then,

s (3) k = K+1+ l=K+1 a l-K-1 ! k; l
So, for k = k + , 06 6N -1, combining, we have

s k = s (1) k + s (2) k + s (3)
k

We now consider the case when the index pad encounters the left boundary, j = 0, l = 0, which occurs for 0 6k 6N -1. Schematically, the situation is portrayed in Figure 1(c). Denote the following subregions:

A : 06j; l6 k + N B : 06l6N -1 --16j6 -1 C : 06j6N -1 --16l6 -1 D : --16j6 -1 --16l6 -1 where 0 6 6N -1.
Again we divide the entire l range --16l 6k + N into two parts: 0 6l 6k + N and --1 6l 6 -1. Consider the range 0 6l6 k + N for which the index pair ( j; l) ranges over B and A. When (j; l) ∈ B, the index j is a ected by periodization so that ! k; l has a regular part and a periodic part. Let ! k; l = ! r k; l + ! p k; l . Then ! r k; l = k+N j=0 h j j ; l (the regular part)

! p k; l = -1 j=--1 h K+1+j j ; l thus s (1) k = k+N l=0 a l ! k; l
Next we deal with the range --16 l 6 -1. In this case (j; l) ∈ C ∪ D. For (j; l) ∈ C, we have

! k; l = N -1 j=0 h j j ; l Then s (2) k = -1 l=--1 a K+1+l ! k; l
When (j; l) ∈ D we have

! k; l = -1 j=--1 h K+1+j j ; l and
s 3 k = -1 l=--1 a K+1+l ! k; l
Combining, for k = N --1, 06 6N -1, we have

s k = s (1) k + s (2) k + s (3) k
The periodization is then complete. The result is the vector {s k }. This procedure generalizes to the n-tuple case in a straightforward manner.

Useful connection coe cient relations

We list a number of identities and relations that are useful in the manipulation of connection coe cients. Many of these relations appear in References 13 and 20. Except for the operator inversion relationship, most of the inner product relations can easily be derived by invoking integration by parts, translation, and change of variables. d1; d2 l; m = d2; d1 m; l (17) d1; d2 l = -d1-1; d2+1 l (18) 1; d2; d3 = -0; d2+1; d3 -0; d2; d3+1 (19) d1; d2; d3 l; m = d2; d1; d3

-l; m-l (20)
d1; d2; d3 l; m = d3; d2; d1 l-m; -m (21) d1; d2; d3 l; m = (-1)

d1 d1 i=0 d 1 i 0; d2+i; d3+d3+d1-i l; m (22)
Another relation that is very useful in checking the construction and accuracy of 3-tuple matrices is the check-sum procedure: the column sum of a 3-tuple matrix must equal a corresponding 2-tuple vector component for component. For example, K l=0 1; 0; 0 l; m = 1; 0 m . Lastly, we mention an e cient procedure for the inversion of a ÿrst-order operator, which exploits the symmetric or skew-symmetric nature of the matrix, = ± T , so that the matrix may be diagonalized, = D T . A concrete example of operator inversion appears in Reference 20.

APPROXIMATION OF DIFFERENTIAL OPERATORS

The spectrum { j } of the continuous di erential operator L with periodic boundary conditions on 0 and 1 is discrete. A discrete approximation L of the operator may be found by projecting the operator onto the space spanned by the periodized scaling functions. The discrete and continuous operators can be compared by looking at their spectra. By the Galerkin procedure outlined in Section 4.1 one can form the matrix , an approximate representation of L, in the subspace of periodized scaling functions of resolution p and of genus N . In this section we examine the usefulness and properties of wavelet projection techniques for di erential operators. A comparison of the spectrum of the di erential operator L and its approximation using the connection coecients algorithm reveals information important within the context of solving ordinary and partial di erential equations using wavelet-Galerkin techniques.

The projection of the di erential operator with periodic boundary conditions leads to a 2-tuple . Eigenvalues of this matrix may be found by noting [START_REF] Davis | Circulant Matrices[END_REF] that the circulant matrix , has eigenvalues j given by the expression

j = K k=0 a k e i2 kj=K
where a k is an entry in the ÿrst row in the circulant matrix, namely, an element of the set 0 ; 1 : : : N ; 0; : : : ; 0; -N ; -N+1 ; : : : -1 and j = 1; 2; : : : ; K + 1. The 's were obtained by the procedure outlined in the preceding section. The overdetermined system was solved using a QR algorithm from LAPACK [START_REF] Anderson | LAPACK User's Guide[END_REF] called dgess. For 0; 1 , the residual and the condition number as a function of the resolution p and the genus N appears in Table I.

Some of the eigenvalues of 0; 1 are given in Tables II andIII. Regarding these tables, one can make three general observations. First, the imaginary part of the eigenvalues grows only as O(2 p ), just as the Fourier-basis computation would, and this is in contrast to spectral approximations based on Chebychev or Legendre bases. Second, we observe that the real parts of the eigenvalues are comparable to the size of machine precision zero, since the real part of each eigenvalue is proportional to the 0 , again in sharp contrast to Chebychev and Legendre bases investigated by Trefethen et al. [START_REF] Trefethen | An instability phenomenon in spectral methods[END_REF] Third, little evidence of contamination due to round-o error was seen even for high resolution (in our experiments, we tried values of p as high as 11), in sharp contrast to the experience of Trefethen et al. with the above-mentioned bases functions. The spectrum of : =35•8. We examined the dependence of the eigenvalue convergence on both the resolution p and the genus N . We did this by examining the imaginary part of several eigenvalues. We found that the rate of convergence was almost exactly 2 -p , while the dependence of the rate of convergence on N was essentially quadratic. Tables II and III illustrate the above-mentioned rates of convergence. Qualitatively, we ÿnd the same type of behaviour in higher odd-ordered di erential operators. Table IV shows the approximation to the spectrum of the operator d 3 =dx 3 with periodic boundary conditions. We found the rates of convergence very similar to those of the ÿrst derivative, with the proviso that the basis functions had to be chosen with su cient smoothness.

Finally we consider approximation to even-ordered di erential operators. For example, Table V shows the approximation of the spectrum of d 2 =dx 2 for periodic boundary conditions for the same values of p. This case corresponds to the eigenvalues of 0; 2 . Each eigenvalue has a complex conjugate, for which the imaginary part is no greater than machine precision zero. Again, we found the same rates of convergence as in the odd-ordered cases.

The condition number (cn) of the vector calculation was found to have the following correspondence: for p = 5 it was cn = 52 778, p = 6 it was cn = 298 557, p = 7 it was cn = 1 688 892. In general it was found that the condition number grew as the size of p, N , and on the order of the derivatives was increased. The source of the high condition number is the ill-posed nature of the over-determined matrix. The high condition number is especially troublesome in the calculation of the 3-tuples. There are few options available to improve the ill-conditioned nature of the matrix. One option could have been to use Chebychev polynomials in connection with the moment conditions. However, the ill-posedness is in the linear homogeneous equations, its exact nature being that several of the rows will have nearly the same entries for some values of p, N , and the order of the derivatives. Hence, little or no improvement would result in the use of the Chebychev polynomials as a way to mix the moment conditions.

Using the eigenvalue data corresponding to the three di erential operators, we have also been able to assess numerically the size of the largest eigenvalues. This is a useful estimate, for example, in the determination of time stability in numerical schemes for the solution of di erential equations. The estimate is that the largest eigenvalue is comparable in size to Fourier spectral approximations; that is, the largest eigenvalue =(2 ) n = O(2 pn ), for d n =dx n , where n = 1; 2; 3. It was also possible to conÿrm that the eigenvalues of the discretized matrix of odd-ordered di erential operators have real parts of magnitude no greater than machine precision zero. The same can be said of the imaginary part of discretizations of even-ordered operators.

SUMMARY

The primary aim of this study has been to complement the report of Latto et al. [START_REF] Latto | The evaluation of connection coe cients of compactly supported wavelets[END_REF] in which they propose a simple and useful technique for the calculation of inner products involving Daubechies wavelets. It is hoped that this study will be a helpful guide to a publicly available numerical code which implements their scheme. Details on how to obtain the code are given in Appendix I.

The periodization of matrices involving inner products of two and three wavelets and their derivatives was presented in detail. The three wavelet inner product matrix is not straightforward to implement, especially in periodized problems. In this study we have provided several useful algorithms to construct such matrices and have listed several properties of the inner products that are useful in checking the actual implementation of the wavelet-Galerkin solution of di erential equations.

Lastly, we have shown numerically that common di erential operators may be well approximated using wavelets. In light of other studies that compare several spectral approximations to the above operators, [START_REF] Trefethen | An instability phenomenon in spectral methods[END_REF] the wavelet approximations compare favourably. Information on the approximation of di erential operators is also useful in the solution of partial di erential equations using wavelet-Galerkin techniques, for example, in making choices of time integrators in the solution of evolution equations, and in estimating the size of the time stepping in order to achieve time-stability in the calculation. It was found that the eigenvalues, which approximate commonly used di erential operators, where fairly close in value and magnitude to those calculated by means of discrete Fourier projection techniques.

APPENDIX I

Tables of connection coe cients

Connection coe cients: 2-tuples for p = 0, N = 6. Computed in double precision on a Sun Sparc1 workstation using LAPACK solver routines.

  where the periodicity of the wavelets has been invoked. By changing variables, we further reduce the equation to 0; d k1; k2 = 0; d 0; k2-k1 ≡ d k2-k1 where d = d 1 + d 2

  d ; ' 0;l : Setting x = 2 j and deÿning M d l = x d ; l , we haveM d l = 2 dj 2 j=2 d ; l = 2 dj 2 d = d! ' j; 0 ( ) d = d!2 -j=2 Thus l M d l d l = d!2 -j=2 :The sum over l is actually over |l|6 N -2 since the ''s are compactly supported. Thus, by changing the indices of summation by m = l + 1 + (N -2)

  x d '(x -l) dx= (y + n) k '(y) a property of the ', the above relation can be used to evaluate recursively M d l

Figure 1 .

 1 Figure 1. Structure of j ; k index pad for 3-tuples in the ( j; l) plane: (a) the index pad corresponds to N = 4 in the region where no wraparound arises; (b) corresponding to N = 4 for k = k + 2; (c) corresponding to N = 4 for k = N -3, = 2

0; 1 - 4 = 1 - 3 =

 1413 -3•4246575342471D -04 0; -1•4611872146119D -02

Table I .

 I Residual and condition number for the leastsquares calculation of 's used in 0; 1 , as a function of p and N

				Condition
	N	p	Residual	number
	4	4	-2•1926904736347D -15	137
	4	5	-8•4099394115356D -15	387
	4	6	-8•3266726846887D -15	1096
	4	7	-2•6922908347160D -14	3100
	4	8	-2•1510571102112D -16	8770
	6	4	4•3905851176973D -15	189
	6	5	4•3905851176973D -15	531
	6	6	2•2332830029725D -14	1503
	6	7	3•2980562725271D -14	4252
	6	8	-1•6563139748627D -14	12 029
	8	4	-1•6766260279716D -11	229
	8	5	-3•3533817109376D -11	642
	8	6	6•7080697513231D -11	1816
	8	7	1•3417001243474D -10	5138
	8	8	2•6830899981444D -10	14 532

Table II .

 II Imaginary part of eigenvalues divided by 2 for increasing resolution 2 p calculated using N = 6 periodized Daubechies wavelets symmetric operator d=dx is purely imaginary and equal to 2 k, where k ∈ Z. Periodicity will lead to the eigenvalue 0 having multiplicity 2. TableIIshows the magnitude of the imaginary part of the ÿrst few eigenvalues, divided by 2 , as a function of p. It is clear from the table that approximations to the eigenvalue improve as p is increased and the number of reasonably correct eigenvalues grows as well with p. Parenthetically, we remark that with p = 6 we can go as far k = 11 for which 11 : = 10•8 while retaining at least one digit of accuracy. In the case p = 8 we can go as far as k = 36 for which 36

	p = 4	p = 6	p = 8
	0•0D0	0•0D0	0•0D0
	0•99997222310553	0•99999999300555	0•99999999999689
	1•9967784015701	1•9999991102504	1•9999999997795
	2•0423307943545	2•9999849535965	2•9999999962590
	2•9532542889943	3•9998888924219	3•9999999720194
	3•4758292699748	4•9994799193486	4•9999998667457
	3•7208826878105	5•9981781285634	5•9999995232107
	4•0025895152829	6•9947814463569	6•9999985996877
	the skew-

Table III .

 III Imaginary part of eigenvalues of d=dx as a function of N for p = 6

	N = 4	N = 6	N = 8
	0•99999690699928	0•99999999300555	0•99999999998407
	1•9999013641150	1•9999991102504	1•9999999916456
	2•9992552692402	2•9999849535965	2•9999996834005
	Table IV. Spectrum approximation of d 3 =dx 3 with N = 6. Imaginary
		part of eigenvalues
	p = 4	p = 6	p = 8
	0•0D0	0•0D0	0•0D0
	39•411980740054	39•478149261094	39•478416554374
	308•17967443293	315•7932221828	315•82720643905
	956•10277711077	1065•3407907411	1065•9149806290
	1670•0988161963	2522•3667673633	2526•6015527133
	1877•4681030848	4914•9287548852	4934•7204114901
	2625•9035313813	8457•8432634425	8 527•0455871351
	2641•9768303508	13 342•439258378	13 540•237953821

Table V .

 V Real part of eigenvalues for d 2 =dx 2 divided by 4 2 for increasing resolution p. N = 6

	p = 4	p = 6	p = 8
	1•1518884962509D-14	2•7645323910020D-13	0•0D0
	1•0033162940157	1•0000135268157	1•0000000529814
	4•1847329915128	4•0008582644240	4•0000033890207
	10•666776606377	9•0096362703930	9•0000385683104
	22•726823972803	16•053060704248	16•000216429062
	41•582358732169	25•197220008678	25•000824274322

ACKNOWLEDGEMENTS Special thanks to George Schlossnagle for his help at various stages of this work. The authors thank John Weiss for making his preprints available to us. Gilbert Strang provided extremely helpful comments on a draft of this report. This research was supported in part by an appointment to the Distinguished Postdoctoral Research Program sponsored by the U.S. Department of Energy, O ce of University and Science Education Programs, and administered by the Oak Ridge Institute for Science and Education. Further support was provided by the Atmospheric and Climate Research Division and the O ce of Scientiÿc Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38.

for the least-squares solution of the overdetermined system.

1; 1

A code for the computation of connection coe cients is available. The URL is http: //www.math.ucla.edu/∼restrepo.

APPENDIX II

Proofs related to periodizing Daubechies wavelets are presented in this appendix.

Proposition. To prove that dilation and periodization of Daubechies wavelets do not commute it su ces to show that for j60; ' j; k = 2 -j=2 :

Letting y = 2 j x and summing over l, we obtain

Periodized Daubechies wavelets form an orthonormal basis for L 2 [0; 1]:

Proposition. { ' j; k : j¿0; 06k ¡ 2 j } forms an orthonormal basis for L 2 [0; 1].

Deÿnition. f; g = f(x)g(x) dx, the standard L 2 inner product.

Proof. We begin by showing that ' j; k ; 'j ; k = 0:

' j; k+2 j r ; ' j ; k = jj kk with l -l = r. Thus ' j; k and ' j ; k are orthonormal. Next we show that they form a basis for

f ∈ L 2 (R) and {' j; k } form an orthonormal basis for L 2 (R); so we have

which, when periodized, becomes

This ÿnal result is actually a ÿnite sum since, for k¿2 j and k61 -N , supp('

Thus f has a representation in the periodized wavelets.

The proof that { ' j; k ; j ; k : j ¿j¿0; 06k ¡ 2 j } also form an orthonormal basis is nearly identical. Since '(x) = N -1 k=0 h k '(2x -k), we can periodize both sides to get

The connection coe cients code may be used for both periodized as well as the non-periodized Daubechies wavelets inner products. This result is formalized as follows.

Proposition. For j¿ log 2 (N -1); ' j; 0 = ' j; 0 ; where ' is extended to R by setting it to 0 away from the unit interval.

Proof. supp(' j; 0 ) = [0; 2 -j (N -1)], so for j¿ log 2 (N -1); supp(' j; 0 ) = [0; ÿ]; ÿ61. Thus ' j; 0 (x) = l∈Z ' j; 0 (x -l) = ' j; 0 : APPENDIX III Several approximating characteristics in Section 3 are proved below. As we have already seen, periodized wavelets provide a basis for L 2 [0; 1] so we have f -P j f 2 → 0, as j → ∞. This is a property of any orthonormal basis of L 2 , but this particular periodized basis has some additional properties.

Theorem. If f is a continuous function on the torus; then f -P j f ∞ → 0 as j → ∞.

Proof. We begin the proof by showing that our projection operator is bounded. P j is an integral operator of the form P j f(x) = 1 0 2 j -1 k=0 ' j; k (y) ' j; k (x)f(y) dy. Thus,

' ∞ since for j¿ log 2 (N -1) there are at most (N -1) k's such that for a given x, {x} ∩ supp( ' j; k ) = ∅. Hence, we have

If we take f ∈ j∈N V j , then ∃J such that ∀j¿J; Q j f = 0. Thus, P j f = f for j¿J . j∈N V j is dense in L 2 [0; 1] which is dense in C(T), continuous functions of period 1 on the unit interval. Finally, by the boundedness of P j , the theorem follows.

Theorem (Daubechies). If f ∈ L 1 [0; 1], then f -P j f 1 → 0 as j → ∞ [10].

Proof. Since L 1 [0; 1] ⊂(C[0; 1]) ? , we have P j f 1 = sup{| P j f; g |; g continuous; g ∞ 61} = sup{| f; P j g |; g continuous; g ∞ 61} 6 f 1 P j g ∞ P j ∞ is bounded by the previous theorem. Since j∈N V j is dense in L 2 [0; 1], which is also dense in L 1 [0; 1], the uniform bound on P j is su cient to prove our result.


Plamen Petrov
email: ppetrov@tu-sofia.bg

Fawzi Nashashibi
email: fawzi.nashashibi@inria.fr

Automatic Vehicle Perpendicular Parking Design Using Saturated Control

Keywords: automatic parking, perpendicular parking, path planning, saturated control

This paper considers the perpendicular reverse parking problem of front wheel steering vehicles. Relationships between the widths of the parking aisle and the parking place, as well as the parameters and the starting position of the vehicle for planning a collision-free reverse perpendicular parking in one maneuver are presented. A nonlinear saturated (tanh-type) feedback steering controller for straight-line tracking is proposed and evaluated. It is demonstrated that the saturated controller, which is continuous, achieves quick steering can be successfully used in solving parking problems. Simulation results and first experimental tests confirm the effectiveness of the proposed control scheme.

I. INTRODUCTION

The perpendicular parking is efficient and economical because it accommodates the most vehicles per linear meter [1]. However, due to the special constraint environments, much attention and driving experience is needed to control the vehicle, and this parking maneuver may be a difficult task. For this reason, automated operation attracts significant attention from research view point, as well, and from the automobile industry. One of the difficulties in achieving automatic parking is the narrow operating place for collision-free motion of the vehicle, and planning of optimal trajectories is often used in the applications. In [START_REF] Gutjahr | Automatic collision avoiding during parking maneuver -an optimal control approach[END_REF], an optimal stopping algorithm was designed for parking using an approach combining an occupancy grid with planning optimal trajectories for collision avoidance. The geometry of the perfect parallel parking maneuver is presented in [START_REF] Blackburn | The geometry of perfect parking[END_REF]. In [START_REF] Pradalier | Path planning for parking assistance system : Implementation and experimantation[END_REF], a practical reverse parking maneuver planner is given. A trajectory planning method based on forward path generation and backward tracking algorithm, especially suitable for backward parking situations is reported in [START_REF] Moon | A trajectory planning method based on forward path generation and backward tracking algorithm for automatic parking systems[END_REF]. A car parking control using trajectory tracking controller is presented in [START_REF] Lee | Car parking control using a trajectory tracking controller[END_REF]. In [START_REF] Petrov | Saturated feedback control for an automated parallel parking assist system[END_REF], a saturated feedback control for an automated parallel parking assist system is reported. In recent years, automatic parking systems have been also developed by several automobile manufacturers [9,10].

In this paper, we deal with the problems of geometric collision-free path planning and feedback steering control for automatic perpendicular reverse parking. An important condition to achieve effective collision-free parking maneuver is to determine a convenient start position of the vehicle depending on the widths of the parking aisle and the parking place, as well as the dimensions and the steering capabilities of the car. Geometric path planning based on admissible circular arcs within the available spot is presented in order to steer the vehicle in the direction of the parking place in one maneuver. A saturated (tanh-type) feedback control is proposed which achieves quick steering in tiny spots and stabilizes (practical stabilization) the vehicle in the parking place. The parking problem is seen as an extension of the tracking problem, since the vehicle has to track a straight line which passes through the goal position in the parking space and is parallel to the parking place with varying velocity depending of the position of the vehicle with respect to the goal position. The rest of the paper is organized as follows: In Section II, geometric considerations for planning perpendicular reverse parking in one maneuver are presented. In Section III, a saturated feedback steering controller is described. Simulation results and first experimental tests are reported in Section IV. Section V concludes the paper.

II. COLLISION-FREE PATH PLANNING WITH CONSTANT TURNING RADIUS

A. Vehicle Model

In this paper, front-wheel steering vehicles are considered and a rectangular model of the vehicle is assumed. To illustrate the results presented in this paper, we conduct simulation and experiments using a CyCab vehicle developed at INRIA, which is shown in Fig. 1. The vehicle parameters which affect the parking maneuver, as well as the parameter values used in the simulations, are presented in Table I.

B. Geometric Considerations for Collision-Free Perpendicular Parking in One Maneuver

The geometry of the reverse perpendicular parking in one collision-free maneuver is shown in Fig. 2.

Fig. 2. Geometry of the collision-free perpendicular parking maneuver

Two inertial coordinate systems are assigned to the parking place: The frame Axy with center A which is placed in the middle between the borders of the parking space and which has its Fy-axis aligned with the boundary c 2 of parking lot L2; the frame Gxy with center G placed at the goal position of the parking place which has its Gx-axis aligned with the Ax-axis of Axy, as shown in Fig. 2.

1) Feasibility conditions for collision-free perpendicular parking in one maneuver

In the perpendicular parking scenario, the vehicle starts to move backward from an initial position 1 in the parking aisle, with constant steering angle α c which is smaller or equal to the maximum steering angle α max and has to enter in the parking place without colliding with the boundary c 1 of parking lot L1 and boundaries c 2 , and c 3 of parking lot L2. When a parallel orientation with respect to the parking place is attained, the vehicle continues to move backward in a straight line into the parking place until it reaches the final position 2 (Fig. 2).

Point P placed in the middle of the rear wheel axle is assigned as the reference point of the vehicle (Fig. 2). Assuming a vehicle circular motion with center O and turning radius ρ c (ρ c = OP) calculated from the formula

c c l α ρ tan = .
(1)

The boundaries of the turning path during perpendicular parking are determined by the dimensions of the traces (circular arcs) formed by the left corner of the front bumper B 2 with radius r B2 , the left corner of the rear bumper B 4 with radius r B4, and the end of the rear wheel axle C 1 , as shown in Fig. 2. Since the vehicle motion is rotation in the plane, the trajectories of these points form arcs of concentric circles.

From the ΔOC 2 B 2 , applying the Pythagorean theorem, we obtain an expression for the radius r B2 of the circular arc traced by the left corner of the front bumper B 2 in terms of the vehicle parameters l, l 1 , b (Table I), and the constant turning radius ρ c , as follows ( )

2 2 1 2 2 2 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + + = = b l l OB r c B ρ . ( 2
)
Similarly, from the ΔOC 2 B 4 , we determine the radius r B4 , of the circular arc traced by the left corner of the rear bumper B 4

2 2 2 4 4 2 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + = = b l OB r c B ρ . ( 3
)
Let O denotes the center of rotation (the Instantaneous Center of Rotation (ICR)) of the vehicle when it starts the parking maneuver with constant front-wheel steering angle α c . Depending on the sign of the x-coordinate of ICR O with respect to coordinate frame Axy, i.e., the offset s (Fig. 2), different formulas can be derived in order to determine the required width h p of the parking place and the width of the parking aisle (the corridor) h c , in order to ensure collision-free perpendicular parking in one maneuver. Due to a limited space, in this paper, we consider the case when s is negative and belongs to the interval

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - ∈ 0
C 2 B 4 • P • • • 1 • h p • F
are not too large. The lower value corresponds to the case when the right side of the vehicle B 1 B 3 (Fig. 2) lies on the boundary line c 2 of parking lot L2.

In order to avoid collision between the left corner B 2 of the front bumper with the boundary c 1 of L1 (Fig. 2), using (2), we obtain an expression for the width of the parking aisle h c , as follows ( )

s b l l s r h c B c - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + + = - = 2 2 1 2 2 ρ . ( 4
)
The function h c = f(s) is linear in s, positive and monotonically increasing in the above-mentioned closed interval for s. Therefore, it takes its minimum and maximum values at the ends of this interval.

To avoid collision between the right point C 1 of the rear axle with vertex A of the obstacle L2, from the ΔOAD, applying the Pythagorean Theorem, the distance OD is determined as follows

2 2 2 s b OD c - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = ρ . ( 5
)
In order to avoid collision between the left corner B 4 of the rear bumper with the edge c 3 of the parking place, using (3) and ( 5), the following expression for the width of the parking space h p , is obtained

2 2 2 2 2 4 2 2 s b b l OD r h c c B p - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + = - = ρ ρ . (6)
Consider the function h p = f(s), which is continuous on the closed interval of s mentioned above. The function is

differentiable on the open interval ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - ∈ 0 , 2 b s c ρ , and its derivative is 0 2 2 2 < - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = ∂ ∂ s b s s h c p ρ . ( 7
)
Therefore h p = f(s) is strictly decreasing on the closed interval [-(ρ c -b/2), 0]. The maximum and minimum values of h p can be found by replacing the limit values s = -(ρ c -b/2) and s = 0 of the interval, respectively, in [START_REF] Lee | Car parking control using a trajectory tracking controller[END_REF].

2) Collision-free parking with assigned values for the widths of the parking aisle and the parking place

From a practical point of view, it is important to determine the starting positions of the vehicle for automatic parking without collision in one maneuver, in the case, when the widths h c and h p of the parking aisle and the parking space, are specified in advanced. Suppose that the widths of the parking aisle and the parking place are set as h c = h cd and h p = h pd , respectively, and also that h cd < r B2 . In this case, from (2) and (4), it follows that

2 max B cd r h s - = - . ( 8
)
From ( 3) and ( 6), we obtain a formula for the minimum value of s as follows

( ) 2 4 2 min 2 pd B c h r b s - - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - = - ρ . ( 9
)
It should be noted that, when s = -|s| min the vehicle will park without collision, but asymmetrically (not centered in the middle) in the parking space with respect to its boundaries, i.e., h pl ≠ h pr (Fig. 2). When the vehicle is parallel to the parking space, the distances between the car and the parking space boundaries h pl and h pr are determined as follows

- ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = ρ ρ , ( 10
) pr pd pl h b h h - - = .
(11)

3) Symmetrical collision-free parking with assigned values for the widths of the parking aisle and the parking place

From a practical view point, it is better to park the car symmetrically with respect to the boundaries of the parking place, since it is not very wide. For this end, we calculate the minimum value of the offset s = s m , in order to park symmetrically the vehicle in the center of the parking space (Fig. 2). We set

2 : b h h h h pd pl pr ps - = = = . ( 12
)
From the ΔOAD (Fig. 3), the distance OD is determined as

2 2 2 m c s b OD - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = ρ . ( 13
)
Since the turning radius can be expressed as

OD h b pr c + + = 2 ρ , ( 14
)
and substituting h pr from (10) and OD from (13) into (14), we arrive to an expression for s m , as follows

2 2 2 2 | | ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - = - pd c c m h b s ρ ρ . ( 15
)
The new offset -|s| m is bigger than minimum offset -|s| min given by ( 9), (-|s| m > -|s| min ). In general, it must be checked whether the new offset -|s| m is smaller than -|s| max given by [START_REF] Petrov | A hybrid control for automatic docking of electric vehicles for recharging[END_REF]. If it is the case, the car can park symmetrically without collision in reverse with a minimum offset s = -|s| m . In this case, however, the boundary c 3 of the parking place will not be tangent to the arc of circle traced by point B 4 of the left corner of the rear bumper; nevertheless, point A (vertex A of obstacle L2) will lie again on the arc of circle traced by point C 1 of the rear vehicle axle. Therefore, given specified values h c = h cd and h p = h pd , the offset s can take values in the closed interval -|s|∈ [-|s| m , -|s| max ], where -|s| m and -|s| max are determined by formulas ( 15) and ( 8), respectively.

Hence, in order to perform reverse perpendicular parking without collision in one maneuver and to place the vehicle symmetrically into the parking place with assigned constant steering angle α c which is smaller or equal to the maximum steering angle α max , the starting position, i.e., the reference point P of the vehicle has to be on any one of the arcs of circles with radius ρ of center O(x O , y O ) with respect to the coordinate system Fxy attached to the parking place (Fig. 2), where x O belongs to the interval [-|s| m , -|s| max ] and y O = -ρ c . The vehicle initial orientation has to be tangent to this arc. In particular, if the orientation of the vehicle is (almost) perpendicular to the parking place, the initial coordinates of the reference point P with respect to the frame Fxy have to be x P (0) = ρ -|s|, where -|s|∈[-|s| m , -|s| max ] and y P (0) = -ρ c . The reference path of the parking maneuver consists of two parts. The first one is a circular arc with center O connecting the staring position of the vehicle and the Fx-axis of Fxy, which is tangent to the circular arc. At the tangent point T with coordinates (-|s|, 0), the car will be parallel to the parking place. The second part of the reference path is a straight line with distance TG along the Fyaxis of the coordinate frame Fxy between the tangent point T and the goal position G of the parking place, since point G lies on the Fx-axis of Fxy.

III. STEERING CONTROL

For a low speed motion, which is the case of the parking maneuver, we assume that the wheels of the vehicle roll without sliding, and the velocity vectors are in the direction of the orientation of the wheels. We consider a simplified (bicycle model) of the vehicle, where the front and rear wheels are replaced by two virtual wheels, placed at the longitudinal axis of the vehicle. The inertial coordinate system Gxy is attached to the parking place in the goal position G (Fig. 2). The coordinates of the reference point P in Gxy are denoted by (x P , y P ). The orientation of the vehicle θ is defined as an angle between the x-axis of Gxy and the longitudinal vehicle base. Information from the encoders mounted on the wheels was used for estimation the position of the vehicle in the parking area. The steering angle of the front virtual wheel is denoted by α. The equations of motion of the vehicle in the plane have the form [START_REF] Petrov | Saturated feedback control for an automated parallel parking assist system[END_REF] where v P is the velocity of point P. We consider a practical stabilization of the vehicle in the parking place. Our approach is based on controlling the motion of the vehicle along a straight line (the Gx-axis of Gxy) passing through the goal point G in the parking place and aligned with the orientation of the place. The velocity profile is dependent of the distance between the vehicle and the goal position [START_REF] Petrov | Saturated feedback control for an automated parallel parking assist system[END_REF]. Since the reference path for the first part of the parking maneuver is a circular arc, a bang-bang feedback controller [START_REF] Petrov | A hybrid control for automatic docking of electric vehicles for recharging[END_REF], where the front wheel steering angle is constrained by magnitude and takes only two constant values may be a solution, since the vehicle trajectories will represent circular arcs. However, in practice, due to the discontinuity of the bang-bang control law, an undesirable behavior of the system (chattering) will occur. In order to avoid the chattering, in this paper, a saturated control based on hyperbolic tangent function is proposed, which is also constrained by magnitude, but the control function is continuous.

A. Saturated Control

In this paper, we propose a differentiable saturation in the form of hyperbolic tangent (tanh(.)) constraint. This function is bounded by ± 1. Also tanh(x) ≥ 0 if x ≥ 0, and tanh(x) < 0 if x < 0. Tanh(x) is close to the signum function, when in tanh(K t x) the gain K t is large, as shown in Fig. 3.

] ) tanh( tan v K lu a t = α , (17)
where

) ( 0 P y a K v - = θ , ( 18
) l u c α tan = , ( 19
)
and K t, K and a 0 are positive constants. The proposed controller achieves asymptotic stability of the closed loop system composed of equations ( 16)-( 19). The proof is similar to those presented in [START_REF] Petrov | Saturated feedback control for an automated parallel parking assist system[END_REF] and is omitted here due to space limitation.

It should be noted that the constant value α c for the front wheel steering angle can be smaller or equal of its maximum value α max . If the starting position of the vehicle does not satisfy the conditions mentioned in Section II for parking in one maneuver, multiple maneuvers have to be performed in order to reach the goal position in the parking place.

IV. SIMULATION RESULTS

Simulation results using MATALAB are presented to illustrate the effectiveness of the geometrical path planning procedure for collision-free perpendicular parking in one maneuver, as well as the performance of the proposed saturated steering controller for practical stabilization of the vehicle into the parking place. The parameters of the test vehicle (CyCab) are given in Table I.

A. Path Planning Simulation

In order to demonstrate the proposed path planning procedure for collision-free perpendicular parking in one maneuver presented in Section II, in the first simulation, we present the relationships between the widths h c and h p of the parking aisle and the parking space, as functions of the offset s which belongs to the interval [-(ρ c -b/2), 0] where s ≤ 0 using parameters of the test vehicle with α c = α max , ( ρ c = ρ min ). The assigned values of h c and h p were chosen to be: h cd = 3m and h pd = 2m. As seen from Fig. 4, the function h p = f(s) (the solid blue line) decreases in the interval and converges to b=1.2m (the red dotted line), which is exactly the length of the wheel base of the vehicle. Meanwhile, the graph intersects the horizontal line for the assigned value of h pd = 2m (the blue dotted line) at s = -|s| min = -1.3016m, which is the minimum value of s obtained from (9) for collision-free parking. In order to park the vehicle in one maneuver for s = -|s| min = -1.3016m, from (4), the required minimum width h c of the parking aisle has to be h c = 1.793m which is less than the specified value of h cd = 3m.

The function h c = f(s) (the green solid line) increases linearly in this interval and the graph intersects the horizontal line for the assigned value of h cd = 3m (the green dotted line) at s = -|s| max = -0.0946m, which is the maximum value of s, obtained from [START_REF] Petrov | A hybrid control for automatic docking of electric vehicles for recharging[END_REF]. For s = -|s| max = -0.0946m, from (6), the required minimum width h p of the parking place has to be h p = 1.2258m, which is less than the specified value of h pd = 2m.

Therefore, given specified values h c = h cd = 3m and h p = h pd = 2m for the parking aisle and the parking space, respectively, for collision-free parking, the offset s can take values in the interval [-|s| min , -|s| max ] = [-1.3016, -0.0946], where the boundary values are determined by ( 9) and ( 8), respectively.

For s = -|s| min , the vehicle will park without collision but not centered in the middle of the parking place. In this case, Using (10) and (11), simulation results indicate that the distances between the car and boundaries of the parking place are h pr = 0.7772 m and h pl = 0.0228m for h pd = 2m. In order to park symmetrically in the parking place (h pr = h pl = 0.4m), by using ( 15

B. Simulation of Perpendicular Parking using Saturated Control

Simulation results are presented to illustrate the performance of the proposed saturated steering controller for perpendicular reverse parking. Using (1), for α c = α max , for the minimum turning radius ρ c =ρ min is obtained ρ min =2.0785m. The dimensions of the parking aisle and the parking place were the same used before: h cd = 3m and h pd = 2m.

1) Simulations for perpendicular parking in one maneuver

The initial coordinates of the vehicle reference point P with respect to the inertial frame Gxy attached to the goal position were chosen to be (x P (0), y P (0)) = (3m, -ρ) =(3m, -2.07m) where y P (0) = -ρ min . In this case the offset -|s| = -1.0785m ∈ [-|s| m , -|s| max ] = [-1.3016, -0.0946]. The initial orientation of the vehicle is chosen to be θ(0) = -π/2rad. Therefore, the vehicle is able to park in one maneuver. The maximum value of the timevariant vehicle velocity was |v P | = 0.3m/s. The values of the saturated tanh-type controller were K t = 8, K = 5.85, a 0 = 0.17.

The planar path of the vehicle (a) and an animation of the perpendicular reverse parking in one maneuver are presented in Fig. 5.Evolution in time of the vehicle velocity (a) and the front-wheel steering angle (b) are presented in Fig. 6. Since the saturation control is continuous, there is not chattering when the position of the vehicle is in the vicinity of the tracking line and in addition, the orientation error is small too. As seen from the simulation results (Fig. 6,b), the saturated control can be used instead of bang-bang control in order to steer the vehicle into the parking place without collision, since it is close to the signum function when the gain K t is large enough.

2) Simulations for perpendicular parking in multiple maneuvers

In this case, the initial coordinates of the vehicle reference point P with respect to the inertial frame Gxy attached to the goal position were chosen to be (x P (0), y P (0)) = (3m, -(ρ+0.5)) = (3m, -1.57m). The offset -|s| = -1.0785m belongs to the interval [-|s| min , -|s| max ] = [-1.3016, -0.0946], but y P (0) ≠ -ρ min and the vehicle cannot park in one maneuver. The initial orientation of the vehicle was chosen to be θ(0) = -π/2rad. In this case, the vehicle cannot park in one maneuver and multiple maneuvers are necessary. The values of the controller gains were like in the first simulation for automatic parking in one maneuver. The planar path of the vehicle (a) and an animation of the perpendicular reverse parking in three maneuvers are presented in Fig. 7. The overall parking maneuver was achieved by three consecutive maneuvers (backward-forwardbackward). The saturated controller has been implemented on an experimental automatic electric vehicle CyCab and initial tests of perpendicular reverse parking in one maneuver have been initialized (Fig. 1). In the first tests, only information from the encoders mounted on the wheels were used for estimation of the the position and localisation of the vehicle in the parking area. First experiments confirm the effectiveness of the proposed controller. V. CONCLUSION

In this paper, the problem of perpendicular collision-free reverse parking of front wheel steering vehicles has been considered. Geometric considerations for collision-free perpendicular parking in one reverse maneuver have been first presented, where the shape of the vehicle and the parking environment have been expressed as polygons. Relationships between the widths of the parking aisle and parking place, as well as the parameters and the initial position of the vehicle have been given, in order to plan a collision-free maneuver, in the case, when the car has to be centered into the parking place. A tanh-type steering controller for straight-line tracking has been proposed and evaluated. It was demonstrated that this controller which is continuous, was able to achieve quick steering avoiding chattering and can be successfully used in solving parking problems. Simulation results, as well as the first experimental tests with an automatic CyCab vehicle confirm the effectiveness of the proposed control scheme.

Fig. 1 .

 1 Fig. 1. The automatic CyCab vehicle developed at INRIA equiped with odometers, DGPS, and laser range finder

  the case when the dimension of the parking aisles and the parking place b

Fig. 3 .

 3 Fig. 3. The function tanh(Ktx) for Kt = 1, 3, and 10.The proposed bounded controller has the form

[

Fig. 4 .

 4 Fig. 4. Collision-free interval for s

  ), the minimum offset -|s m | was obtained to be -|s m | = -1.0113m, and this value is bigger than -|s| min = -1.3016, (-|s| m > -|s| min ), obtained from (9), but smaller than -|s| max = -0.0946m. Therefore, given specified values h cd = 3m and h pd = 2m, for centered parking into the parking place, the offset s can take values in the closed interval -|s|∈ [-|s| m , -|s| max ] = [-1.0113, -0.0946].

  parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane

Fig. 5 .Fig. 6 .

 56 Fig. 5. Perpendicular parking in one maneuver: Planar paths of the vehicle (a) and animation of the parking maneuver (b)

  parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane Perpendicular parking -Vehicle motion in the plane

Fig. 7 .

 7 Fig. 7. Perpendicular parking in three maneuvers: Planar paths of the vehicle (a) and animation of the parking maneuver (b) Evolution in time of the vehicle velocity (a) and the frontwheel steering angle (b) are presented in Fig. 8.

Fig. 8 .

 8 Fig. 8. Perpendicular parking in three maneuvers: Evolution in time of the vehicle velocity (a) and the front-wheel steering angle (b)

TABLE I .

 I VEHICLE PARAMETERS

	Vehicle parameters	Notation	Value
	Longitudinal vehicle base	l	1.2m
	Wheel base	b	1.2m
	Distance between the front axle and the front bumper	l1	0.35m
	Distance between the rear axle and the rear bumper	l2	0.35m
	Maximum steering angle	αmax	π/6rad