Nicolas Ducros
email: nicolas.ducros@polimi.it

Cosimo D'andrea

Andrea Bassi

Gianluca Valentini

Simon Arridge

Virtual source pattern method for fluorescence tomography with structured light

Keywords: Fluorescence tomography, optical tomography, dimensionality reduction, image compression

In order to reduce both acquisition and reconstruction times, illumination and detection in fluorescence diffuse optical tomography (FDOT) have recently evolved from a point-based to a pattern-based approach. The use of structured illumination, offering the ability to project any pattern of light onto the object, associated to the compression of the acquired fluorescence images has paved the way to a new generation of fast reconstruction algorithms for FDOT. However, the choice of the most appropriate set of source patterns is still an open problem.

Here, the use of typical source patterns is investigated on experimental data. Reconstructions of similar qualities are obtained for the different types of source patterns. We found that the performances of structured illumination are limited by the required positivity of the source patterns. To alleviate this problem, we introduce a novel method, namely the virtual source patterns method, which allows for considering any kind of patterns, e.g., with negative and complex intensities. This new method provides a significant increase of the contrast of the reconstruction and also a reduction of the reconstruction error, especially when virtual wavelet source patterns are considered.

Introduction

In the last years, the rapid development of specifically designed optical markers and technological improvements of imaging system, has fostered the field of optical molecular imaging [START_REF] Cherry | In vivo molecular and genomic imaging: new challenges for imaging physics[END_REF][START_REF] Ntziachristos | Fluorescence imaging with near-infrared light: new technological advances that enable in vivo molecular imaging[END_REF]. This methodology allows one to follow, using longitudinal studies, functional processes and disease pathways at the molecular level in the living animal. The main advantages of optical techniques are represented by their high sensitivity, low invasiveness and low cost.

The main challenge of optical molecular imaging is to deal with thick tissues for which the images appear blurred due to scattering. Dealing with animal models, this aspect is particularly problematic for mice, whose dimensions and optical properties correspond to a fully diffusive regime, such that transmission imaging suffers from a lack of depth resolution. Fluorescence diffuse optical tomography (FDOT), also referred to as fluorescence molecular tomography (FMT), has been designed to solve these problems allowing to reconstruct the three-dimensional (3-D) fluorescence distribution in small samples [START_REF] Milstein | Fluorescence optical diffusion tomography[END_REF][START_REF] Ntziachristos | Looking and listening to light: the evolution of whole-body photonic imaging[END_REF][START_REF] Herv | Noncontact fluorescence diffuse optical tomography of heterogeneous media[END_REF][START_REF] Montet | Tomographic fluorescence imaging of tumor vascular volume in mice1[END_REF][START_REF] Koenig | In vivo mice lung tumor followup with fluorescence diffuse optical tomography[END_REF][START_REF] Kumar | A time domain fluorescence tomography system for small animal imaging[END_REF][START_REF] Klose | Excitation-resolved fluorescence tomography with simplified spherical harmonics equations[END_REF].

In classical FDOT, part of the sample surface is illuminated by a point source and the fluorescent light reemitted by the tissue is acquired by a detector (CCD cameras or photomultipliers via an optical fibre). Such a measurement is repeated by changing the input source position and/or rotating the sample. To improve the spatial resolution of the fluorescence reconstruction, a common strategy consists in using a large amount of illumination and detection points. However, due to the increase in the size of the data set, this strategy is not always numerically tractable. Moreover it leads to a significant increase of the acquisition time. To remain compatible with in vivo measurements, a trade off between acquisition/reconstruction times and reconstruction performances must be established.

In the last years, different research groups have begun to study the possibility to increase the speed of FDOT by reducing the data-set without sacrificing the reconstruction quality [START_REF] Cuccia | Modulated imaging: quantitative analysis and tomography of turbid media in the spatialfrequency domain[END_REF][START_REF] Bassi | Detection of inhomogeneities in diffusive media using spatially modulated light[END_REF][START_REF] Konecky | Quantitative optical tomography of sub-surface heterogeneities using spatially modulated structured light[END_REF]. To this end, it has been proposed to reduce the dimension of the measurement space by implementing two main concepts, often together. First, a series of works have capitalized on the low-pass filtering of the spatial frequencies of light propagating in highly scattering medium [START_REF] Cuccia | Modulated imaging: quantitative analysis and tomography of turbid media in the spatialfrequency domain[END_REF][START_REF] Bassi | Temporal propagation of spatial information in turbid media[END_REF]. Fluorescent images obtained from biological tissue can be compressed to few coefficients with little degradation using image basis sets such as Fourier [START_REF] D'andrea | Fast 3d optical reconstruction in turbid media using spatially modulated light[END_REF][START_REF] Ripoll | Hybrid fourier-real space method for diffuse optical tomography[END_REF] or wavelets [START_REF] Rudge | Fast image reconstruction in fluoresence optical tomography using data compression[END_REF][START_REF] Ducros | A time-domain wavelet-based approach for fluorescence diffuse optical tomography[END_REF][START_REF] Ducros | Multiple-view fluorescence optical tomography reconstruction using compression of experimental data[END_REF]. Knowing the basis functions corresponding to the retained coefficients, referred to as detection patterns, a compressed version of the inverse problem can be numerically solved. This leads to a significant reduction in the reconstruction time. Analytical expressions for reconstruction with large data sets have also been developed [START_REF] Markel | Inverse problem in optical diffusion tomography. iii. inversion formulas and singular-value decomposition[END_REF][START_REF] Wang | Experimental demonstration of an analytic method for image reconstruction in optical diffusion tomography with large data sets[END_REF] and successfully applied to experimental data [START_REF] Panasyuk | Fluorescent optical tomography with large data sets[END_REF]. However, the latter approach is limited to simple geometries such as slabs, cylinders, or spheres, for which symmetry properties are available [START_REF] Markel | Symmetries, inversion formulas, and image reconstruction for optical tomography[END_REF]). Second, it has been proposed to adopt patterns of illumination instead of point sources [START_REF] Cuccia | Modulated imaging: quantitative analysis and tomography of turbid media in the spatialfrequency domain[END_REF]Joshi et al., 2006a,b;[START_REF] Konecky | Imaging complex structures with diffuse light[END_REF][START_REF] Lukic | Optical tomography with structured illumination[END_REF][START_REF] Konecky | Quantitative optical tomography of sub-surface heterogeneities using spatially modulated structured light[END_REF][START_REF] Dutta | Illumination pattern optimization for fluorescence tomography: theory and simulation studies[END_REF] with a view to reducing the number of illuminations. Indeed, few well-chosen patterns can be used to acquire as much information as a dense raster scanning of a point source [START_REF] Lukic | Optical tomography with structured illumination[END_REF]Ducros et al., 2010a). A collateral advantage is that wide field illuminations allow one to project high light power on the sample without overcoming safety limits based on power density. In this context, we have proposed to use both detection and illumination compression. In this way, fast 3-D reconstructions in arbitrary geometry, has been demonstrated with experimental data acquired from a slab medium (Ducros et al., 2010a) and from a cylinder [START_REF] Ducros | Multiple-view fluorescence optical tomography reconstruction using compression of experimental data[END_REF]. Recently, a similar approach was used to perform tomographic measurement of small animals in time domain in [START_REF] Chen | Time-resolved diffuse optical tomography with patterned-light illumination and detection[END_REF].

An open problem lies in the choice of the most appropriate set of illumination patterns to be projected onto the sample. Up to now different types of source patterns have been proposed. They can be categorized into five groups:

• Uniform or plane wave. It consists in illuminating the object with a spatially uniform pattern (Joshi et al., 2006b). • ad hoc. Different type of patterns have been empirically tested. Among them, moving bars (Joshi et al., 2006b;[START_REF] Bélanger | Real-time diffuse optical tomography based on structured illumination[END_REF], scanning Gaussians and diffractive optics patterns (Joshi et al., 2006a), checkerboard [START_REF] Bélanger | Real-time diffuse optical tomography based on structured illumination[END_REF]. • Sinusoids. The homogeneous infinite slab being a purely low pass filter, the sinusoidal illumination patterns have been proposed to detect the presence of embedded inclusions. This natural candidate has been used for imaging purpose on the slab [START_REF] Cuccia | Modulated imaging: quantitative analysis and tomography of turbid media in the spatialfrequency domain[END_REF], but also for analytical [START_REF] Konecky | Quantitative optical tomography of sub-surface heterogeneities using spatially modulated structured light[END_REF] and numerical [START_REF] Bassi | Detection of inhomogeneities in diffusive media using spatially modulated light[END_REF][START_REF] D'andrea | Fast 3d optical reconstruction in turbid media using spatially modulated light[END_REF][START_REF] Ripoll | Hybrid fourier-real space method for diffuse optical tomography[END_REF] reconstruction in simple geometries. • Wavelets bases. For more complex geometries wavelet bases have been proposed (Ducros et al., 2010a[START_REF] Ducros | Multiple-view fluorescence optical tomography reconstruction using compression of experimental data[END_REF] for their well known compression ability. • Optimal patterns. A fundamental step has been made by [START_REF] Dutta | Illumination pattern optimization for fluorescence tomography: theory and simulation studies[END_REF] who provided optimal source patterns for a slab, a cylinder, and a numerical phantom simulating the optical properties of animal models (Digimouse). However, a comparative study among the different patterns is still missing. In particular, it remains unclear how the choice of the source pattern translates in terms of reconstruction quality when experimental data are considered.

The goal of this paper is twofold. The first goal is to experimentally compare the typically-used source patterns in term of reconstruction quality. In particular, it will be shown that even if different patterns are used the reconstruction quality is similar to that obtained with a uniform illumination. Since physically realisable source distributions are necessarily positive, a second goal of this paper is to introduce the virtual source patterns (VSP) method which is a general method able to determine the diffusive sample response to any negative or complex pattern and therefore to fully take advantage of the structured illumination. In Section 2, we first give the main equations formulating the forward problem of FDOT. Then, a method enabling to consider arbitrary patterns is presented. Last, the four steps of our reconstruction algorithm are described. In Section 3, we present five case studies. The first three case studies deal with actual sources patterns and the last two implement virtual source patterns. In Section 4, we discuss the reconstruction results that show a clear advantage in implementing the VSP method. Finally, our conclusions are reported in Section 5.

R 3 ∂Ω Surface of the medium R 2 ∂Ω e Illumination surface R 2 ∂Ω m Detection surface R 2 Σ e Illumination plane R 2 Σ m Detection plane R 2 φ Photon density C(Ω) d Detection pattern C f Discrete fluorescence distribution R N m measured image R N ccd + m virtually measured image C N ccd s source pattern R N dmd + s virtual source pattern C N dmd B Boundary condition operator C(Ω) → C(Ω) M Measurement operator C(∂Ω m ) → C(∂Σ m ) P Propagation operator C(Ω) → C(Ω) X Illumination operator C(∂Σ e ) → C(∂Ω e )

Method

In this paper, italic letters (a, α or A) indicate functions or constants. Vectors are denoted by lowercase bold italic letters (a or α) and matrices by capital bold italic letters (A or Φ). We denote f, g Ω = r∈Ω f (r)g(r) dr the inner product between two functions f and g, A * the adjoint of operator A, and x the complex conjugate of x. A summary of the notation used throughout the paper is provided in table 1.

Forward Problem

As depicted in figure 1, we consider a device, e.g. a digital micromirror device (DMD), able to illuminate a portion ∂Ω e of the surface ∂Ω of the medium Ω embedding a fluorescence distribution f (r). The fluorescence light emitted from the surface ∂Ω m is recorded by a CCD camera. The forward problem consists in determining the fluorescence image acquired by the CCD knowing the source pattern s loaded on the DMD and the fluorescence distribution f .

Within the medium, the photon density φ (W.m -2 ) obeys the following coupled equations [START_REF] Arridge | Optical tomography: forward and inverse problems[END_REF]:

P e φ e = 0 in Ω
and B e φ e = X s on ∂Ω, (1a)

P f φ f = φ e ημ f a in Ω and B f φ f = 0 on ∂Ω. ( 1b
)
where P is the propagation operator, B is an operator accounting for boundary conditions, and X is a mapping from the DMD plane Σ e to the surface object ∂Ω e . Equation (1a) describes the propagation of light at the excitation wavelength λ e while (1b) describes propagation at the fluorescence wavelength λ f . In (1b), the source term ημ f a φ e (W.m -3 ) represents the part of the photons absorbed at λ e that are remitted at λ f , with η the quantum yield of the fluorescent marker and μ f a its absorption coefficient. In what follows, we denote f = ημ f a the fluorescence distribution, which is to be recovered. In this paper, we consider the diffusion approximation together with partial boundary conditions, which leads to [START_REF] Ishimaru | Theory and application of wave propagation and scattering in random media[END_REF])

P = -∇[D(r)∇] + μ a (r), and (2a)
B = I -AD(r) n • ∇, ( 2b
)
where D is the diffusion constant of the medium and μ a its absorption coefficient. In (2b), I stands for the identity operator, A is a constant that depends on the index mismatch, and n indicates the outward-pointing normal to the surface ∂Ω.

The fluorescence image acquired by the CCD is m = M[φ f ] where M is a mapping from the object surface ∂Ω m to the CCD plane Σ m . In our approach we compress the acquired image

m m = K k=1 mk d k (3)
where the d k are the basis of the image compression method; in analogy to the source compression basis used herein, we consider the image compression basis vectors as detection patterns.

Given a fluorescence distribution f (r), the components mk can be obtained by means of an integral operator with kernel w. Mathematically:

mk = w(s, d k ), f Ω , ( 4
)
where the integration kernel w only depends on the source and detection pattern s and d k . It can be shown that (see appendix A for proof)

w(s, d k ) = G e [u] G f * [v k ], (5)
where G e (resp. G f ) is the Green's operator for (1a) (resp. (1b)), u = X [s] is the source pattern mapped onto surface ∂Ω e , and v k is the adjoint detection pattern defined on surface ∂Ω m . The adjoint detection pattern is defined by

v k = M * [d k ] ( 6 )
In a general problem, the object is rotated and multiple-illumination acquisitions are performed. We denote I the number of rotation angles and J the number of source patterns that are projected onto the object for each view. Such a process results in the acquisition of a set of IJ images. Each of them is compressed to K components, which leads to a total number of K tot = IJK image components. Discretizing the volume Ω in N voxels, the forward problem given in (4) simplifies to the discrete product m = W f.

(

) 7
where m ∈ C Ktot is the measurement vector containing the retained image components, W ∈ C Ktot×N is a matrix referred to as the weight matrix, and f ∈ R N is the discrete concentration vector. We denote by w i,j,k the discretised kernel of ( 5) for the k th detection pattern of the j th source and i th rotation, with k the fastest varying index, followed by j, and with i the slowest varying index. Equation ( 7) can be written:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ m1,1,1 . . . m1,1,k . . . m1,j,k . . . mi,j,k . . . mI,J,K ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ w 1,1,1 1 . . . w 1,1,1 n . . . w 1,1,1 N . . . . . . w 1,1,k 1 . . . w 1,1,k n . . . w 1,1,k N . . . . . . w 1,j,k 1 . . . w 1,j,k n . . . w 1,j,k N . . . . . . w i,j,k 1 . . . w i,j,k n . . . w i,j,k N . . . . . . w I,J,K 1 . . . w I,J,K n . . . w I,J,K N ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ • ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ f 1 . . . f n . . . f N ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (8)

Virtual source patterns method

Here we describe a simple approach, namely the virtual source pattern (VSP) method, that offers the possibility to exploit the response of the diffusive sample to any negative or complex pattern, although only patterns with positive intensities can be projected on the sample. The method capitalize on the linearity of the forward problem.

Let S = [s 1 . . . s J ] represent the J (desired) virtual source patterns. Our approach consist in finding J actual source patterns S = [s 1 . . . s J ] and a J × J (possibly complex) matrix T such that ST = T S T .

(9)

The matrix T is referred to as transformation matrix. As a consequence of linearity, the same transformation can be applied to the actual measurements M = [m 1 . . . m J ] so as to get the virtual measurements M = [ m1 . . . m J ]. The virtual images would have been measured if the virtual patterns S were projected onto the medium. Formally:

M T = T M T . ( 10
)

Virtual wavelet patterns

In this section, we describe how to generate wavelet patterns by means of the virtual source method. The reader can refer to the work of [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] for a general discussion on wavelet bases. The wavelet functions (mother wavelets) can be derived from the so-called scaling functions (father wavelets). In dimension 2, at a given resolution, the i)M vertical wavelets {q V }, ii)M horizontal wavelets {q H }, and iii)M diagonal wavelets {q D } can be constructed from the 4M scaling functions {p} taken at a higher resolution. Hence, to obtain virtual wavelet patterns S in ( 9), the actual source patterns S can be chosen as scaling functions and the transformation matrix T can be constructed accordingly. We choose:

S = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ q V 1 q H 1 q D 1 . . . q V M q H M q D M ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ T = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ t V 1,1 . . . t V 1,4M t H 1,1 . . . t H 1,4M t D 1,1 . . . t D 1,4M . . . . . . t V M,1 . . . t V M,4M t H M,1 . . . t H M,4M t D M,1 . . . t D M,4M ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ S = ⎛ ⎜ ⎝ p 1 . . . p 4M ⎞ ⎟ ⎠ (11)
As detailed in Appendix A.2, each row of T contains the coefficients of a 2D vertical, horizontal or diagonal high-pass filter, reshaped into a 4M × 1 vector. Hence, any DC component added to the scaling functions S cancels out after multiplication with T . This guarantees that any wavelet pattern can be generated even if the corresponding scaling functions contain negative value. The technical details to computing T are reported in appendix Appendix A.2.

FDOT reconstruction algorithm

The sketch of our algorithm is presented in figure 2. The different steps are detailed below. We denote m i,j the fluorescence image measured at the ith view after illumination of the object with the jth source pattern s j .

Step 1 The I J virtual fluorescence images mi,j are calculated from the IJ fluorescence images m i,j according to (10). The same way, the J virtual source patterns sj are calculated from the J actual source patterns s j according to (9). The virtual source patterns are mapped from the DMD plane Σ e onto the object surface ∂Ω e i and the resulting ũi,j = X i sj are concatenated to form the virtual source matrix

Ũ = [ũ 1,1 . . . ũi,j . . . ũI, J ] ∈ C N ×I J .
Step 2 Each of the I J virtual fluorescence image mi,j is wavelet transformed and retaining the K i,j (absolute value) largest components leads to the compressed image mi,j =

Ki,j k=1 mi,j,k d k . ( 12
)
The retained components mi,j,k are stacked into the vector m ∈ C Ktot . In this paper, the detection patterns d k are chosen from a Daubechie 4 basis.

Next, the adjoint detection patterns v i,j,k , defined on the object surface ∂Ω m i ,

Weight Matrix computation Step 3 The computation of the weight matrix W ∈ C Ktot×N requires the calculation of ϕ i,j ∈ C N and ψ i,k ∈ C N that are the vectorised versions of G e [u] and G f * [v k ] in (5), respectively. The source density ϕ i,j is needed for all couples (i, j), while the adjoint detection density ψ i,k must be obtained for any couple (i, k). This may be done in one shot inverting the matrix equations

P e Φ = Ũ and P f Ψ = V . ( 13
)
The matrix Φ = [ϕ 1,1 . . . ϕ i,j . . . ϕ I, J ] ∈ C N ×I J contains the source densities ϕ i,j and the matrix Ψ = [ ψ1 . . . ψ . . . ψL ] ∈ C N ×L contains the non-redundant adjoint detection densities ψ . Any vector ψ i,k can be obtained by picking up the appropriate ψ in Ψ . In the latter equation, P e and P f are two N × N sparse matrices that result from the discretization of the propagation operator P e and P f and account for the boundary conditions. Here, P e and P f have been computed thanks to the Matlab TOAST package [START_REF] Schweiger | TOAST package[END_REF].

The weight matrix W is finally built up row by row, according to (5). Precisely, the row of W corresponding to the triplet (i, j, k) is given by

w i,j,k = ϕ i,j • ψ i,k , ( 14
)
where • stands for the entrywise (Hadamard) product.

Step 4 The inverse problem in FDOT is well-known to be ill-posed and requires regularization to ensure the stability of the solution f in the presence of noise.

Here, the inverse problem is efficiently solved considering the underdetermined Tikhonov solution

f = W T (W W T + αI) -1 m ( 15
)
where α is the regularization parameter that enforces the stability of the solution.

The procedure for choosing its value is detailed in Section 2.6.

Sample and Experimental set-up

We consider a cylindrical phantom mimicking the dimensions and optical properties of the mouse. The phantom is an epoxy resin cylinder of radius 20 mm and height 45 mm. Its optical properties at λ e = 630 nm have been measured with a time-resolved spectroscopy system [START_REF] Bassi | Portable, large-bandwidth time-resolved system for diffuse optical spectroscopy[END_REF] and the following values have been obtained: μ a = 0.022 mm -1 μ s = 1.35 mm -1 . The geometry of the fluorescent inclusions are represented at figure 3b.

The experimental set-up has been reported in detail in D'Andrea et al. ( 2010). It is based on three main components: illumination device, detection system and rotational stage sample holder. Light is emitted by a He:Ne laser source (05LHP991, Melles Griot) at wavelength of 633 nm with optical power of about 10 mW coupled to an optical fiber guide (LLG0538-4, Thorlabs). Then, the light is focused by a lens and, through a total internal reflection prism, is delivered to a Digital Micromirror Device (DMD) (Discovery 1100 -ALP1, Vialux). DMD consists of a matrix of titlable mirrors (720 × 521), whose reflecting elements can be independently controlled. By uploading the image of the desidered spatial pattern, the arrangement of the tiltable mirrors and the light beam reflected by the DMD change accordingly. Through the control of the duty cycle, patterns with a dynamic range of 8 bits can be generated by the DMD. Finally, the pattern is magnified and projected on the cylindrical phantom by means of an a projection lens. The illumination area ∂Ω e is a rectangle of width x = 13 mm and height y = 26 mm, 2.5 mm apart from the top of the cylinder. The sample holder is mounted on a motorized rotational stage in order to perform multiple view measurements.

The light exiting the opposite side of the cylindrical phantom is imaged by an objective lens (f = 50 mm, f #= 2.8, Nikon Co.) on a low noise high sensitivity 16 bit CCD camera (Versarray 512, Princeton Instruments,) cooled to -40˚C. The detection area ∂Ω m is a rectangle of width 16 mm and height 32 mm. The illumination and detection area are depicted in figure 3a. In order to select the fluorescence signal, a combination of an interference filter (XF 3076, Omega) and a high pass filter RG-695, Schott) is placed in front of the objective lens. The signal at the excitation wavelength can be acquired by removing the filters. The system is placed in a lightproof cabinet in order to prevent stray light. Particular attention has been devoted to the automatization of the measurements system, which is fully computer controlled, in order to reduce the acquisition time.

Comparison of the reconstructions

In order to fairly compare the reconstructions obtained from different types and numbers of illumination patterns, the choice of the regularization parameter is crucial. The proposed approach is inspired by [START_REF] Ducros | A comprehensive study of the use of temporal moments in time-resolved diffuse optical tomography: part ii. three-dimensional reconstructions[END_REF], where the regularization parameter is set by imposing the reconstruction variability. Here, the regularization parameter is set by maximizing the contrast-to-noise ratio (CNR) of the reconstructed fluorescence distribution f . The CNR is a figure of merit that indicates how well structures of interest are rendered in a reconstruction [START_REF] Song | Automated region detection based on the contrast-to-noise ratio in near-infrared tomography[END_REF][START_REF] Baritaux | An efficient numerical method for general l p regularization in fluorescence molecular tomography[END_REF]. It is defined by

CN R = μ roi -μ back (w roi σ 2 roi + w back σ 2 back ) 1/2 , ( 16
)
where μ roi (μ back ) and σ roi (σ back ) are respectively the average and standard deviation of f in the ROI (background). The weights w roi and w back are the volume ratios

w roi = V roi /V tot and w back = V back /V tot .
Here, the region of interest (ROI) refers to the volume occupied by the fluorescence inclusions and the background refers to the rest of the volume.

In addition to the CNR, we consider other two figures of merit, namely the contrast resolution and the reconstruction error. The contrast resolution C is defined by

C = μ roi -μ back μ roi + μ back . ( 17
)
It measures the ability to tell the difference between the ROI and the background in terms of mean values. The higher the contrast resolution, the better. Finally, we consider the reconstruction error E r . Expressed in decibels, the following definition is retained:

E r (dB) = 20 log f μ true roi μ roi -f true -20 log f true . ( 18
)
where μ true roi is the average of f true in the ROI. The smaller the reconstruction error, the better.

Results

We present five experiments for five different types of source patterns. The first three experiments implement real source patterns while the last two implement virtual source patterns. The profiles of some source patterns are displayed for each experiment in figure 4. The five experiments have been performed in the same conditions. Every source pattern has been projected onto the the cylindrical phantom for Δt = 35 s. Sixteen rotation angles from 0˚to 337.5˚have been considered. For objective comparison among reconstruction qualities, the three performance metrics are calculated for each experiment and reported in table 2. Some reconstruction profiles for the five experiments are provided at figure 10. Another important aspect is the amount of time required to perform FDOT. In table 3 both acquisition and reconstruction times are reported for each experiment. Central column: optimal regularization parameter αopt as defined in Section 2.6; Left-hand column: under-regularized situation for which α = α under < αopt, ; Right-hand column: over-regularized situation for which α = αover > αopt. The maximum and minimum value of each colorbar have been set to 1.5μ roi and -0.5μ roi , respectively. b) Graph of the CNR of the reconstruction with respect to the regularization parameter α. The blue asterisk indicates α under , the red dot αopt, and the black triangle αover.

Case 1 -Uniform illumination

The first experiment is performed loading onto the DMD a constant image whose intensity level is set to the maximum value (255). This experiment is considered as the reference experiment. First, because it enables for measuring the highest signal in transmission. Indeed, the zero spatial frequency is the less attenuated one (see for example [START_REF] Bassi | Portable, large-bandwidth time-resolved system for diffuse optical spectroscopy[END_REF]). Second, since a DMD intrinsically cuts light off, this experiment benefits from the highest amount of light available. The fluorescence images experimentally acquired are depicted in figure 5. The choice of the regularization parameter follows the rule given in Section 2.6. In figure 6a we present some reconstruction slices corresponding to three different regularization parameters: a small α, the optimal α, and a large α. The graph of the CNR with respect to the regularization parameter is given in figure 6b. The whole set of reconstruction slices for the optimal reconstruction parameter is depicted in figure 7a.

Case 2 -Sinusoids

In this experiment, we considered sinusoidal patterns of the form cos 2π(n x x/ x + n y y/ y )+1, where n x and n y are the numbers of cycles along x and y, respectively. In diffuse optics, the considered media behave as low-pass filters. Though source patterns of increasing spatial frequencies could increase the spatial resolution of reconstructions, such frequencies are increasingly attenuated and cannot be experimentally measured. Thus, a common strategy consists in considering few low frequency patterns. In this case, patterns up to 1 cycle along x (the smallest dimension) and up to 3 cycles along y (the largest dimension) are regarded. Precisely, a set of 8 bi-dimensional frequencies is considered: n x ∈ {0, 1} and n y ∈ {0, 1, 2, 3}. Note that this experiment includes the uniform pattern of case 1, i.e., (n x , n y ) = (0, 0). Note also that, apart from this pattern, all the other patterns lead to an illumination of half the amount of total light. The resulting reconstructed fluorescence distribution is presented in figure 7b.

Case 3 -Haar scaling functions

In this experiment, we consider the 2-D scaling functions of a wavelet basis. It is well known that such functions can provide a multiresolution approximation of images [START_REF] Unser | Wavelet theory demystified[END_REF]. Since there is generally no closed form for these functions, they are more easily obtained from the inverse wavelet transform. Let As denotes the wavelet transform of s at a given resolution. The matrix A may be written A = [A a , A d ] where A a s provides an approximation of s, while A d s give the details of s. Henceforth, the scaling patterns may be simply obtained as s = A -1 a ŝ where ŝ is a unit vector chosen from the natural basis {e i }.

In this experiment, we consider the Haar wavelet. At the resolution we consider, the patterns simplify to height black images containing a white square of side x /2 = y /4. The reconstructed fluorescence distribution obtained from eight Haar scaling patterns is presented at figure 7c.

Case 4 -Phasor methods

The phasor method has been extensively used for demodulation in communications systems [START_REF] Carlson | Communication systems[END_REF] where it is known under square-law detection. Referred to as phase shifting technique, it has been a method of choice in profilometry [START_REF] Srinivasan | Automated phase-measuring profilometry of 3-d diffuse objects[END_REF] and was successfully applied in confocal microscopy [START_REF] Neil | Method of obtaining optical sectioning by using structured light in a conventional microscope[END_REF]. In diffuse optics, it was first introduced by [START_REF] Cuccia | Modulated imaging: quantitative analysis and tomography of turbid media in the spatialfrequency domain[END_REF] for filtering out the DC component in the measured images. The phasor method consists in i)projecting onto the medium a set of P shifted source patterns s[p] = cos 2π(n x x/ x + n y y/ y + p/P ) + 1, with p ∈ {0, P -1}, p/P being the phase shift and ii) post processing the resulting P images m to obtain a demodulated image m = p m[p] exp(i2πp/P ). From the VSP point of view, the phasor method can be interpreted as the use of some actual source patterns s to gain the response to a virtual source pattern s = exp[i2π(n x x/ x + n y y/ y )]. Note the virtual source pattern no longer contain additive DC component. For a given pair of frequencies (n x , n y ) and choosing P = 3, the transformation matrix defined in Section 2.2 simplifies to:

T = 1 -1 2 -1 2 0 1 -1 . ( 19
)
The patterns used in experiment 2, with n x ∈ {0, 1} and n y ∈ {0, 1, 2, 3}, are used together with three phase shifts p ∈ {0, 1, 2}, which results in a set of twenty-four actual source patterns. It is worth noting that after processing the measurements not only the real part, but also the imaginary part of the virtual fluorescence image are available.

The reconstruction of the fluorescence distribution obtained from the phasor method is available in figure 9a.

Case 5 -Virtual Haar wavelets

While scaling functions are used in case 3, wavelet functions are desired in this experiment. Wavelet functions possess a higher frequency content than scaling functions and provide details (at different resolution level) in image compression. Implementing such patterns requires the VSP method since wavelets contain negative values. Interestingly, wavelet functions at a given resolution level are easily obtained from scaling functions at a finer resolution level (see Section 2.3). Capitalizing on this property, the measurements of case 3 can serve in this experiment to get the response to Haar wavelet patterns at a lower resolution. Considering the two filters

h = 1 √ 2 [1, 1] T and g = 1 √ 2 [1, -1]
T and following the four step described in Section Appendix A.2 leads to the following transformation matrix:

T = 1 2 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 1 0 0 -1 -1 0 0 1 -1 0 0 1 -1 0 0 1 -1 0 0 -1 1 0 0 0 0 1 1 0 0 -1 -1 0 0 1 -1 0 0 1 -1 0 0 1 -1 0 0 -1 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 20
)
The virtual source patterns, i.e., the columns of S = T S are displayed in figure 8 together with the associated virtual measurements.

The reconstruction of the fluorescence distribution from the virtual Haar wavelets is available in figure 9b.

Discussion

Using only a uniform illumination, the overall quality of the reconstruction is quite good since the contours of the three inclusions can be recovered. Inclusion 2 is reconstructed with a correct localization in the xy plane as well as along z. However, the reconstruction of inclusions 1 and 3 is of lower quality, especially in top slices down to z = 9.5 mm, where inclusion 1 joins inclusion 3. Reconstructed fluorescence is observed to vanish in the first top slices (from z = 19.5 mm to z = 17.5 mm, in particular). The latter effect is a typical problem that can be attributed to the vicinity of boundary. It must be outlined that the top slices also suffer from receiving a limited amount of light since the illumination area starts 3 mm underneath the cylinder summit.

The reconstructions obtained from cases 1, 2, and 3 indicate the performance of the real source patterns. First of all, note that only one source pattern is considered in case 1 while height source patterns are considered in both cases 2 and 3. Hence, the reconstruction quality is expected to be higher in case studies 2 and 3 than in case 1 (at an extra cost in terms of acquisition and computation times). By visual inspection of the reconstruction slices displayed in figure 7, the benefit of using 8 source patterns appears limited. To get a better picture of the situation, it is informative to refer to the figures of merit reported in table 2 and to the reconstruction profiles of figure 10. An improvement is found with a diminution of the reconstruction error for both case studies 2 and 3. While case 3 also exhibits a slight improvement in terms of CNR and contrast, case 2 stands at a CNR of 0.77 and decreases in terms of contrast.

Regarding now the two cases implementing the VSP method, it can be seen that the reconstruction quality is greatly improved with respect to the other three cases. Indeed, the best reconstruction obtained with actual source patterns, namely case 3, is outperformed by both VSP reconstructions of case studies 4 and 5 (compare figure 7a to figure 9a andb and the profiles at figure 10). In particular, inclusion 1 is nicely individuated. Moreover, in spite of the lack of illumination light, the first three slices show little blurring. These visual observations are confirmed by the figures of merits reported in table 2. Among the two VSP implementations, the wavelet approach offers the best performance, whether in terms of CNR, contrast, or reconstruction error. Very similar reconstruction qualities have been obtained from different types of actual source patterns. However, the VSP method results in a better reconstruction quality. For a given set of fluorescence images, the interest of combining real images to form virtual ones prior to reconstructing is demonstrated (compare case 2 to case 4 and also case 3 to case 5). A key feature of the VSP method lies in its ability to remove the DC components in the fluorescence images. This advantage is shown comparing case studies 1, 2, and 4. First, note that each pattern of case 2 contains a DC component plus an AC component. Hence, the similarity of reconstruction 1 and 2 can indicate that the DC components in case 2 overwhelm the potentially useful AC components. This is confirmed by the results of case 4 that show a significant improvement of the reconstruction quality in the absence of the DC component. It must be outlined that removing the DC component is likely to be further profitable in experimental scenarios since it also permits to cancel out undesirable light. To capitalize on this important feature, a satisfying transformation matrix must highpass filter the actual fluorescence images. Hence, the coefficient of the columns of T must satisfy

J j=1 t j,j = 0, ∀j ∈ {1, .., J }, ( 21
)
where the t j,j 's are the matrix elements of T . The superiority of the VSP method is also clearly shown comparing cases 3 and 5, that exploit the same data set. Note that the transformation matrix of case 5 (see (20)) is chosen to satisfy the previous condition (21). To further test the importance of removing DC components, we have also built the virtual wavelet patterns of case 5 choosing a transformation matrix that does not verify condition (21). This results in a poor reconstruction quality (results not shown), which support our interpretation. Although it exploits fewer patterns, the virtual wavelet approach outperforms the phasor method. An interpretation is that the phasor method, which works in the Fourier space, only exploits the frequency resolution of the source patterns and it is known that high frequencies are increasingly damped. On the other hand, the wavelet approach offers a trade-off between the spatial and frequency resolution of the source patterns. In this sense, the wavelet approach provide a trade-off between the pointsource approach and the phasor-method approach. Another important aspect is the noise amplification. The virtual fluorescence images are by construction noisier than the actual fluorescence images they results from. Despite this, the reconstruction quality is improved using virtual patterns, which indicates the applicability of the method on real data.

While satisfying condition ( 21), there are many ways to choose T and S so as to get the desired S . An important aspect of this choice is the cost of the VSP method in terms of extra measurements, which can be indicated through the ratio J/J. The most straightforward approach for obtaining a negative pattern s is to separate its positive part s + from its negative part s -and to consider s = s +s -, which leads to J/J = 1/2. The phasor method provides a more efficient approach for sinusoids with a ratio J/J of 2/3 > 1/2. In Section 2.3, a method is given that enables for building 3M wavelets from 4M scaling functions. With this approach, we have: J/J = 3/4 > 2/3, which is a better ratio than for phasor method. Note that even if the scaling functions contain negative values, the approach of Section 2.3 can still be used provided an appropriate background is added to the scaling functions.

The speed of our FDOT algorithm can be evaluated by means of table 3. We limit the discussion to the case 4 and 5 for which the gain of reconstruction quality is substantial. Since the acquisition time is proportional to the number of real source patterns J, the phasor method requires an acquisition 3 times longer than the virtual wavelet approach. As far as the reconstruction time is concerned, the reduction of the set of detection patterns at step 3 allows for a great speed-up of the algorithm and the use of different source patterns for a given view comes at a limited extra cost. Anyway, the wavelet approach is about three times faster than the phasor method.

Conclusions

This work addresses the problem of the choice of the source patterns in fluorescence optical tomography considering experimental measurements. We have shown that different types of actual source patterns, which are limited to positive intensities, provides very similar reconstruction qualities. A positive pattern can be considered as a uniform illumination (DC component) summed to a non-uniform pattern (with positive and negative intensities). The presence of the DC component limits the advantage of using structured light. Due to the inherent tissue response, the DC component is much less attenuated than the higher spatial frequencies of the pattern. Therefore, the fluorescence images obtained with high spatial frequencies illumination are very similar to those obtained with uniform illuminations and the reconstruction quality is slightly improved by the use of patterned illumination. To alleviate this problem, the so-called virtual source pattern (VSP) method, which allows for considering any source pattern, has been introduced. Implementing this method, virtual source patterns with no DC components are shown to improve significantly the reconstruction. In particular, patterns belonging to 2-D wavelet bases offer the best performance. It is worth outlining that the VSP method has the ability to reject the undesirable sources of light. This is a key feature in experimental scenarios for which the illumination of the object is unavoidably corrupted by some light background. Even though the results of this paper have been verified in FDOT experiments, we believe that many aspects of the results and discussion can be translated to optical system based on photon migration, such as diffuse optical tomography and lifetime tomography. In future work, we will investigate the benefit of the proposed method in small animal imaging.

(ii) Computing 2-D filters. The vertical, horizontal and diagonal 2-D filters, denoted G V , G H , and G D , respectively, are obtained by:

G V = h V (g H ) T , G H = g V
T = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
vect(G V 0,0 ) T vect(G H 0,0 ) T vect(G D 0,0 ) T vect(G V 0,1 ) T vect(G H 0,1 ) T vect(G D 0,1 ) T . . .

vect(G V

MV-1,MH-1 ) T vect(G H MV-1,MH-1 ) T vect(G D MV-1,MH-1 )

T ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (A.11)
Example In the Haar case, we have the following low-pass and hight-pass filters: .12) Choosing M V = 2 and M H = 1, we build (step i):

h = 1 √ 2 [1, 1] T and g = 1 √ 2 [1, -1] T . (A
h V = 1 √ 2 [1, 1, 0, 0] T , g V = 1 √ 2 [1, -1, 0, 0] T , h H = 1 √ 2 [1, 1] T , g H = 1 √ 2 [1, -1] T . (A.13)
From the previous four 1-D filters, we build the following three 2-D filters (step ii):

G V = ⎛ ⎜ ⎜ ⎝ 1 -1 1 -1 0 0 0 0 ⎞ ⎟ ⎟ ⎠ , G H = ⎛ ⎜ ⎜ ⎝ 1 1 -1 -1 0 0 0 0 ⎞ ⎟ ⎟ ⎠ , G D = ⎛ ⎜ ⎜ ⎝ 1 -1 -1 1 0 0 0 0 ⎞ ⎟ ⎟ ⎠ . (A.14)
Next, we have to shift the 2-D filters (step iii):

G V 0,0 = ⎛ ⎜ ⎜ ⎝ 1 -1 1 -1 0 0 0 0 ⎞ ⎟ ⎟ ⎠ , G H 0,0 = ⎛ ⎜ ⎜ ⎝ 1 1 -1 -1 0 0 0 0 ⎞ ⎟ ⎟ ⎠ , G D 0,0 = ⎛ ⎜ ⎜ ⎝ 1 -1 -1 1 0 0 0 0 ⎞ ⎟ ⎟ ⎠ , G V 1,0 = ⎛ ⎜ ⎜ ⎝ 0 0 0 0 1 -1 1 -1 ⎞ ⎟ ⎟ ⎠ , G H 1,0 = ⎛ ⎜ ⎜ ⎝ 0 0 0 0 1 1 -1 -1 ⎞ ⎟ ⎟ ⎠ , G D 1,0 = ⎛ ⎜ ⎜ ⎝ 0 0 0 0 1 -1 -1 1 ⎞ ⎟ ⎟ ⎠ . (A.15)
Vectorising the previous six matrices, we obtain the matrix T as given in (20) (step iv).

Figure 1 .

 1 Figure 1. Principle of acquisition and notations

Figure 2 .

 2 Figure 2. Sketch of the FDOT reconstruction algorithm.

Figure 3 .

 3 Figure 3. a) Source surface ∂Ω e (13 × 26 mm 2 , red line) and detection surface ∂Ω m (16 × 32 mm 2 , yellow line) superimposed on the cylinder. b) Cylindrical phantom (height 45 mm, diameter 20 mm). The three inclusions are cylinders with radius 2 mm. Thirty slices are displayed from z = -9.5 mm to z = 19.5 mm. The fluorophore concentration is about 10 -4 M of Nile Blue in Ethylene Glycol. The concentration has been chosen in order to achieve a peak signal of about 5×10 4 electrons per pixel in the camera sensor.

Figure 4 .

 4 Figure 4. Some of the source patterns considered in the five cases.

Figure 5 .

 5 Figure 5. Fluorescence images acquired at θ = 0˚, 90˚, 180˚, 270˚, from right to left. A uniform illumination pattern is considered. Original images are represented in the top row while the compressed images are represented in the bottom row. The colorbar indicates the number of counts.

  Figure 6. a) Reconstructions obtained for three different regularization parameters. Six slices are displayed from z = 17.5 mm to z = -7.5 mm. Central column: optimal regularization parameter αopt as defined in Section 2.6; Left-hand column: under-regularized situation for which α = α under < αopt, ; Right-hand column: over-regularized situation for which α = αover > αopt. The maximum and minimum value of each colorbar have been set to 1.5μ roi and -0.5μ roi , respectively. b) Graph of the CNR of the reconstruction with respect to the regularization parameter α. The blue asterisk indicates α under , the red dot αopt, and the black triangle αover.

Figure 7 .

 7 Figure 7. Reconstructions from actual illumination patterns: a) Plane wave (case 1); b) Sinusoids (case 2); c) Wavelet scaling functions (case 3). For each reconstruction, thirty slices are displayed from z = 19.5 mm to z = -9.5 mm. Reconstructions are given in arbitrary units. For fair comparison between the different reconstructions, the maximum and minimum value of each colorbar have been set to 1.5μ roi and -0.5μ roi , respectively.

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Acquisitions using the virtual illumination patterns of case 5. First row, virtual illumination patterns s on the DMD. The colorbar indicates gray levels. Second row, virtual fluorescence images m on the CCD. Third row, virtual fluorescence images after compression m. The colorbar indicates the numbers of counts.

  (h H ) T , and G D = g V (g H ) T . (A.9) (iii) Circularly shifting 2-D filters. The 2-D filters obtained at the previous step must be circularly shifted by multiples of 2. We denote G i,j the 2-D filter resulting from a circular shift of G by i elements along its first dimension and by j elements along its second dimension. The following shifted version of G V , G H , and G D are required:G V 2i,2j , G H 2i,2j , and G D 2i,2j for (i, j) ∈ {0, ..., M V -1} × {0, ..., M H -1}. (A.10) (iv) Vectorising. Each of the shifted 2-D filter is vectorised to provide a row of T . Noting vect(G) the vectorised version of G, we have:

Table 1 .

 1 Mathematical notations

	Variable	Meaning	Dimension
	I	Number of views	1
	J	Number of actual source patterns	1
	J	Number of virtual source patterns 1
	K	Number of retained components	1
	N	Number of voxels	1
	Ω	Volume of the medium

Table 2 .

 2 Performance metrics for the different illumination patterns

	Case CN R	C	E r (dB)
	1	2.10	0.77	-17.2
	2	2.09	0.77	-18.3
	3	2.13	0.78	-18.1
	4	2.49	0.83	-18.6
	5	2.76	0.87	-19.0

Table 3 .

 3 Acquisition and reconstruction times (in sec.) for the different case studies

	Case J Acquisition time (s)	J	L	Reconstruction time (s)
						Step 2 Step 3 Step 4 Total
	1	1	35	1 1024	6.1	19	2	28
	2	8	280	8 1189	6.5	44	54	110
	3	8	280	8 1352	6.7	44	54	112
	4	24	840	16 1705	8.5	93	206	329
	5	8	280	6 2220	9.7	54	30	97

Acknowledgment

This work was partly supported by Royal Society International Joint Project 2009/R2, Cariplo Foundation (Grant 2009(Grant -2626) ) and MIUR under the project Futuro in Ricerca (prot. RBFR08XH0H 002).

Appendix A.1. Integration kernel of the forward operator

The retained component m k can be expressed as the projection of the measured image m onto the detection pattern d k . Generally speaking we have:

Substituting m = Mφ f in the previous equation lead to:

By definition of the adjoint operator, we have:

Isolating f that is the quantity of interest, we obtain

and finally

. Derivation of Matrix T

Any wavelet basis is associated to a pair of high-pass and low pass filters, denoted h and g, respectively. We denote L the length of the filters, M V the number of desired wavelets along the first direction (vertical), and M H the number of desired wavelets along the second direction (horizontal). The matrix T , which is of size 3M V M H × 4M V M H , can be derived from the pair of filters following the steps:

(i) Resizing 1-D filters. We first build two filters of size M V × 1 denoted h V and g V .

If L < M V , the filters are obtained by padding h and g with zeros so as to get the desired length. If L > M V , the filters are 'periodised' with a period M V , i.e. h V [i] = k=0 h[i + kM V ] and g V [i] = k=0 g[i + kM V ] are considered.

The same way, we build h H and g H that are two filters of size M H × 1.


Pierre Riedinger

A switched LQ regulator design in continuous time

In this paper, the design of a LQ regulator for linear switched systems in continuous time is investigated. From a relaxation of the optimal control problem, a Lyapunov based switching law is provided. Even if the subsystems are all unstable, the state feedback switching law can be applied subject to a positiveness condition. In any cases, the real cost is always upper bounded by the Lyapunov function value. The optimality of the switching law is also discussed and we prove that the switching conditions are optimal in some generic cases. This point explains why the obtained results over examples approach finely the optimal solutions. Finally, a design strategy is also given that extends the results to the cases where the subsystems are controlled linear systems.

I. Introduction

In the last decade, many contributions have addressed the design of stabilizing switching laws for switched systems both in continuous and discrete time (see for examples [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF], [START_REF] Lincoln | Relaxing dynamic programming[END_REF], [START_REF] Seatzu | Optimal control of continuous-time switched affine systems[END_REF], [START_REF] Zhang | On the value functions of the discrete-time switched lqr problem[END_REF] for dynamic programming approaches, [START_REF] Bengea | Optimal control of switching systems[END_REF], [START_REF] Axelsson | Gradient descent approach to optimal mode scheduling in hybrid dynamical systems[END_REF], [START_REF] Shaikh | On the hybrid optimal control problem: Theory and algorithms[END_REF] for variational approaches, and [START_REF] Deaecto | Switched state feedback control for continuous time-varying polytopic systems[END_REF], [START_REF] Colaneri | Stabilization of continuous-time switched nonlinear systems[END_REF], [START_REF] Geromel | Stability and stabilization of discrete time switched systems[END_REF] for Lyapunov based approaches). This problem cannot be considered trivial even numerically [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF], [START_REF] Shaikh | On the hybrid optimal control problem: Theory and algorithms[END_REF] and the goal to design closed loop control, based on the optimization of a criteria, is a challenging task.

LQ regulators play a central role in the control theory of linear systems due to their simple design to meet performance requirements and their robustness properties. Up to now, the exact solution of a switched LQ problem is not available. Only approximations of the state feedback switching law have been proposed in the literature using the dynamic programming [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF], [START_REF] Lincoln | Relaxing dynamic programming[END_REF], [START_REF] Seatzu | Optimal control of continuous-time switched affine systems[END_REF], [START_REF] Zhang | On the value functions of the discrete-time switched lqr problem[END_REF]. The practical obstacle to the application of these methods relies on the difficulties to compute numerically a good approximation of the solution not only for small dimensional problems. Alternatively, open loop control can be achieved using direct *The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement n • 257462: HYCON2 Network of Excellence "Highly-Complex and Networked Control Systems". 1 Université de Lorraine, CNRS-CRAN UMR 7039, 2, avenue de la forêt de Haye, 54516 Vandoeuvre-lès-Nancy Cedex, France.

or indirect optimization methods [START_REF] Shaikh | On the hybrid optimal control problem: Theory and algorithms[END_REF], [START_REF] Axelsson | Gradient descent approach to optimal mode scheduling in hybrid dynamical systems[END_REF] but singular solutions [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF], [START_REF] Bengea | Optimal control of switching systems[END_REF] entail numerical difficulties [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF].

In [START_REF] Geromel | Stability and stabilization of discrete time switched systems[END_REF], Lyapunov based approaches are developed for a mode independent quadratic cost and the authors show that the design conditions are not convex and the problem cannot be solved by Linear Matrix Inequality (LMI) tools. A switched quadratic Lyapunov approach is also proposed in [START_REF] Deaecto | Switched state feedback control for continuous time-varying polytopic systems[END_REF], based on Lyapunov-Metzler inequalities. If some LMI conditions can be guaranteed, the method allows to define a state feedback switching law and can be applied even when subsystems are unstable. An upper bound on the cost is also provided but the optimality of the proposed stabilizing switched law is not evaluated.

In this paper, a variational approach has been used to analyze the optimal necessary conditions occurring in the switched LQ problem. A first result yields the right expression of the optimal cost with respect to the state and the co-state. From this analysis, a control Lyapunov function as well as a state feedback switching law is determined.

An important point which is discussed, is that the obtained switching conditions are actually optimal in some generic cases when the optimal control is constant. This explains why the state feedback switching law leads to some appealing results close to the optimal. In any case, the value of Lyapunov function defines an upper bound on the real cost and gives a guarantee on performances.

These results are established in a relaxed framework for which the singular arcs [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF], [START_REF] Bengea | Optimal control of switching systems[END_REF] that appear in the optimal solution, are properly taken into account. This is a key point since for this class of systems, the optimal solutions are frequently singular as indicating by the big quantity of randomly tested examples.

The paper is organized as follows. Section II is dedicated to the problem formulation and to the necessary conditions that must be fulfilled to solve the optimal control problem.

A preliminary result establishes the right expression of the optimal cost when the optimal switching law is said regular. In Section III, the main result of this paper is given. A control Lyapunov function and an ad hoc switching law are determined. In Section IV, we discuss the optimality of the proposed switching law. In Section V. several examples illustrate that the proposed switching law is effectively optimal. We propose also a design strategy in the case where the subsystems are, eventually all non stabilizable, controlled linear systems. [START_REF] Seatzu | Optimal control of continuous-time switched affine systems[END_REF] DRAFT

II. Problem Formulation and preliminary results

We consider the class of linear switched systems in continuous time:

ẋ(t) = A σ(t) x(t) x(0) = x 0 (1)
where σ : [0, +∞) → S = {1, • • • , s} denotes the switching law that selects the active mode at time t by choosing among a finite collection of matrices A i ∈ R n×n , i ∈ S . Our aim is to design a state feedback switching law (i.e. x → σ(x)) for system (1) that approaches the optimal solution of the following optimization problem:

Problem 1: Minimize the switched quadratic criterion:

min σ(•) 1 2 ∞ 0 x T (t)Q σ(t) x(t)dt ( 2
)
where [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF], [START_REF] Bengea | Optimal control of switching systems[END_REF] to solve optimal control problem for switched systems ( ẋ = f i (x), i ∈ S ) is to solve its relaxed version, replacing the vector field set ( f i (x)) by its convex hull ( ẋ = co{ f i (x)}). At least, three reasons justify the convexification of the problem: (i) the solutions are well defined (Fillipov; [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF]); (ii) the density of the switched system trajectories into the trajectories of its relaxed version [START_REF] Ingalls | An infinite-time relaxation theorem for differential inclusions[END_REF]; (iii) the existence of singular optimal solutions are taking into account [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF], [START_REF] Bengea | Optimal control of switching systems[END_REF].

Q i = Q T i > 0, i ∈ S subject to ẋ(t) = A σ(t) x(t), x(0) = x 0 . A usual framework
The relaxed version of Problem 1 is then given as a pure continuous time optimal control problem consisting of replacing the matrices Q i and A i respectively by their convex combination i.e.

Q(λ(t)) = s i=1 λ i (t)Q i and A(λ(t)) = s i=1 λ i (t)A i where λ(t) ∈ Λ = {λ ∈ R s : s i=1 λ i = 1, λ i ≥ 0}
H(x, λ, p) = s i=1 λ i H i (x, p) (3)
with H i (x, p) = p T A i x + 1 2 x T Q i x and where p defines the co-state. This leads to the following classical necessary conditions for optimality [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]:

Theorem 1: Suppose that λ * is optimal with the corresponding state x * . Then, there exists an absolutely continuous function p * , named co-state, such that:

1) p * 0, 2) ṗ * = s i=1 λ * i (t)(-A T i p * -Q i x * ) for almost all t ∈ R + , 3) λ * (t) = arg min λ∈Λ H(x * (t), λ, p * (t)), 4) H(x * (t), λ * (t), p * (t)) = 0.
From Item 3 and equation (3), it is clear that if there exists i ∈ S at time t such that H i (x * (t), p * (t)) < H j (x * (t), p * (t)), ∀ j ∈ S \ {i}, then the optimal control has to satisfy λ * i (t) = 1 and λ * j (t) = 0, ∀ j ∈ S \ {i}. A switching instant can occur at time t if there exists at least a pair (i, j) ∈ S 2 such that H i = H j = 0. At this time, the value of λ cannot be determined directly. Actually, if we suppose that 0 = H i = H j < H k , ∀k ∈ S \ {i, j} then the values that satisfy the relation, λ i +λ j = 1, are potential candidate for optimality. Moreover, a so called singular control λ can exist, for which 0 = H i = H j < H k , ∀k ∈ S \ {i, j} on a non empty time interval (a, b). This is a well known situation in the literature [START_REF] Robbins | A generalized legendre-clebsch-condition for the singular cases of optimal control[END_REF], [START_REF] Bryson | Applied optimal control[END_REF], [START_REF] Bonnard | Singular Trajectories and their Role in Control Theory[END_REF] and second order necessary conditions given by the Generalized Legendre-Clebsh Condition [START_REF] Michel | Singular optimal control -the state of the art[END_REF], [START_REF] Krener | The higher order maximal principle and its application to singular extremals[END_REF] are generally necessary to solve the optimal control problem. Definition 1: A singular control λ(.) is such that there exist at least two indices i, j, for which H i = H j = 0 on a non zero measure time interval (a, b) and satisfying λ(t) ∈ Λ,

λ k (t) 1 ∀k ∈ S , ∀t ∈ (a, b).
The corresponding part of the trajectory is named a singular arc.

Remark 1: A singular control defines a Fillipov solution [START_REF] Cortes | Discontinuous dynamical systems[END_REF] for the original switched system [START_REF] Xu | Optimal control of switched systems based on parameterization of the switching instants[END_REF]. Hence, it allows to extend properly the notion of optimal solution for switched systems. Roughly speaking when an optimal solution of the relaxed problem possesses singular arcs, these arcs define sliding surfaces for the switched system (1) which lead to chattering if the surface is attractive. It is noteworthy that only suboptimal solutions can be achieved for the switched systems due to the limited switching frequency; see for example [START_REF] Bengea | Optimal control of switching systems[END_REF].

Theorem 2: If p * (0) is an optimal value for the co-state at a given x 0 , then ηp * (0) is an optimal value for the initial state ηx 0 , for all η ∈ R. Moreover, the optimal cost V * is homogenous of degree 2 :

V * (ηx 0 ) = η 2 V * (x 0 ).
Proof: First, observe by linearity that applying the same control λ from two initial positions x 0 and ηx 0 that the resulting trajectories are homogeneous of degree one i.e.

x(t, ηx

0 ) = ηx(t, x 0 ). It follows that the associated cost V(x 0 , λ) := 1 2 ∞ 0 x T (t)Q λ(t) x(t)dt
are also given by : V(ηx 0 , λ(.)) = η 2 V(x 0 , λ(.)). In particular, it is obvious that if λ * is an optimal control for initial state x 0 , then it is also an optimal one for the initial position ηx 0 (easy to show by contradiction). It remains to show that if p * (0) is an optimal value for the co-state for a given initial position x 0 then ηp * (0) is for ηx 0 . This is achieved if the two corresponding control laws are the same. Due to the homogeneity of (x, p), the switching condition λ * (t) = arg min λ∈Λ H(x * (t), λ, p * (t)) is not modified when we integrate the Hamiltonian system from (x 0 , p(0)) or from (ηx 0 , ηp(0)).

Remark 2: If the value function V * is homogeneous with degree two then its gradient is also homogeneous with degree one (at points where it can be defined).

Theorem 3: Assuming that Problem 1 admits a solution for each x 0 , then the value function

V * (x 0 ) is continuous.
Proof: It can be achieved using the fact that for given control λ ∈ Λ and a given R > 0, A(λ)x is clearly Lipschitz on any ball B(0, R) as well as x T Q(λ)x. Then, an adapted proof and inspired by [chap. 8 pp. 188-189, [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF]] can be applied to get the result.

Remark 3: Unfortunately, we fail to prove that V * is also Lipschitz and then differentiable almost everywhere (which is the case for the optimal control problem states in finite time). Now, consider the set of Lyapunov equations corresponding to each subsystem (A i ), i ∈ S for which a symmetric solution P i exists and is unique:

A T i P i + P i A i + Q i = 0. (4)
Let us define the variables p i from the co-state p as:

p = p i + P i x, ∀i ∈ S , (5)
then, the Hamiltonian function H can be simplified to

H = s i=1 λ i H i with H i = p T i A i x. It is straightforward to verify that when subsystem i is active on a time interval (a, b) i.e. ∀t ∈ (a, b), λ i (t) = 1 and λ j (t) = 0, ∀ j ∈ S \ {i}, the dynamic of p i is given by: ṗi = -A T i p i , (6)
while the evolutions of the others p j for all j ∈ S \ {i} follow from (5):

p j = p i + (P i -P j )x. (7)
Item 4 in above Theorem 1 yields:

H i = p T i A i x ≡ 0, on (a, b) (8)
and a switch can only occur from subsystem (A i ) to subsystem (A j ) at time t = b, if

H i = H j ⇔ p T i A i x = p T j A j x = 0, (9)
and before this switch, the following conditions are fulfilled:

H j = p T j A j x > 0, for all j ∈ S \ {i}. ( 10
)
Theorem 4: Assume there exists a non singular stabilizing control law λ and a co-state p such that Theorem 1 is satisfied. Denote by {t 0 , t 1 , t 2 , • • • } with t 0 = 0, the associated switching time sequence (possibly infinite) and by {i 0 , i 1 , i 2 , • • • } the corresponding mode sequence. Then, the cost function is determined by:

V(x 0 , λ) = 1 2 (x T 0 P i 0 x 0 + p i 0 (0) T x 0 -L). ( 11
)
where L = lim k→N x T k p i k (t k ) with x k := x(t k ) and N ∈ N {+∞} the number of switchings. If the trajectory is optimal, L = 0 and the co-state p i 0 (t) for t 0 ≤ t ≤ t 1 can be identified with:

p i 0 (t) = N k=1 Π T (t, t k )(P i k -P i k-1 )x k . (12)
where Π(t, t k ) is the transition function from time t to t k .

Proof: As it is assumed that the switching sequence is non singular, the cost function can be written:

V(x 0 , λ) = 1 2 t 1 0 x T Q i 0 xdt + 1 2 t 2 t 1 x T Q i 1 xdt + • • •
Using (4) and integrating, we get:

V(x 0 , λ) = 1 2 x T 0 P i 0 x 0 + 1 2 N k=1 x T k (P i k -P i k-1 )x k ( 13
)
where

x k := x(t k ) is the state at switching times t k , k = {0, ..., N}.
Define G := N k=1 x T k (P i k -P i k-1 )x k . Using [START_REF] Shaikh | On the hybrid optimal control problem: Theory and algorithms[END_REF], G can be written as:

G = N k=1 x T k (p i k-1 -p i k )(t k ) = x T 1 p i 0 (t 1 ) -x T 1 p i 1 (t 1 ) + x T 2 p i 1 (t 2 ) + • • • -x T N-1 p i N-1 (t N-1 ) + x T N p i N-1 (t N ) -x T N p i N (t N ) = x T 1 p i 0 (t 1 ) -x T N p i N (t N )
The last equality is due to the fact that the product x T (t)p i k (t) is constant on time interval (t k , t k+1 ) which implies that all the differences disappear. Denoting L by L = x T N p i N (t N ) establishes the cost expression [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF]. Now, we assume that the switching sequence is optimal. If N is finite then it means a last mode is used. Consequently, as for a classical LQ problem, the co-state after the last switch equals to p = P N x and it follows that p i N = 0 and L = 0.

If N is infinite, the term G defines a telescoping series which is convergent if and only if x T N p i N (t N ) converges to a limit L. The value of the series is then x T 1 p i 0 (t 1 ) -L. As the control problem is homogenous (By homogenous, we means that if p 0 is the optimal co-state value corresponding to an initial condition x 0 then αp 0 is for αx 0 ), the only possible value for L is zero when the extremal is optimal. Else it implies that p goes to infinity when

x is steered to the origin which contradicts the homogeneity of the problem. This can be also viewed as follows: Considering the optimal control problem in finite time t f , then the necessary condition p(t f ) = 0 implies from (5) that p i(t f ) = -P i(t f ) x(t f ). Then, taking the limit

lim t f →+∞ p i(t f ) = lim t f →+∞ -P i(t f ) x(t f ) = 0 since lim t f →+∞ x(t f ) = 0 and i(t f ) ∈ S is finite.
Using one more time the fact that the product x T (t)p i 0 (t) is constant between (t 0 , t 1 ), it can be concluded that the value function equals to:

V * (x 0 ) = 1 2 (x T 0 P i 0 x 0 + p i 0 (t 0 ) T x 0 ). ( 14
)
Moreover, denote by Π(a, b) the transition function from time t = a to time t = b, then the cost function can be explicitly defined in terms of x 0 and of the switching time sequence T = (t 0 , t 1 , t 2 , ...) by:

V(x 0 , T ) = 1 2 (x T 0 P i 0 x 0 + x T 0 ( N k=1 Π T (t 0 , t k )(P i k -P i k-1 )Π(t 0 , t k ))x 0 ), (15)
and p i 0 (t) in ( 11) could be identified as in [START_REF] Bonnard | Singular Trajectories and their Role in Control Theory[END_REF].

Remark 4: Equation [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF] shows that the cost function not only depends on x 0 as in the standard LTI LQ Problem. Unfortunately, as we can see in equation [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF], to obtain the value of the cost, the switching states x k (or the switching instants t k ) must be available.

Theorem 5: The variable p i k defined by [START_REF] Bengea | Optimal control of switching systems[END_REF] in the time interval (t k , t k + 1) is given by

p i k (t) = Pi k (t)x(t)
where the time variant symmetric matrix Pi k (•) follows the differential Lyapunov equation:

-Ṗi k (t) = A T i k Pi k (t) + Pi k (t)A i k (16) Pi k (t k+1 ) = Pi k+1 (t k+1 ) + P i k+1 -P i k . ( 17
)
Proof: The proof is obvious in view of equation ( 15) and from the fact that for

t k ≤ t ≤ t k+1 , Pi k (t) = e A T i k (t k+1 -t) Pi k (t k+1 )e A i k (t k+1 -t) .
In view of these preliminary results, it appears that the state space is partitioned following the mode i ∈ S which is used. Along an optimal trajectory, when mode i k ∈ S is used, the optimal cost is of the form

V * (x) = x T P i k x+c i k with c i k = x T Pi k x a constant. At switching time
i.e. on the boundary of the partition, the next constant is given by

c i k+1 = c i k + x T (P i k -P i k+1 )x.
Unfortunately, as the given expression of the optimal cost is related to a trajectory (particularly the constants c i k ), the right expression of the value function is not available on the whole state space and it is not possible to use directly Hamilton-Jacobi-Bellman equation to determine the state space partition. In the next section, a control Lyapunov function is exhibited that allows to approach finely the optimal solution in some generic cases as the section IV and V will show.

III. Lyapunov based switching law

The aim of this part is to define a Lyapunov function as a tight upper bound on the value function. We means tight in the sense that the two functions may coincide at some points.

Consider the set of Lyapunov equations corresponding to each subsystems (A(λ)), λ ∈ Λ for which a symmetric solution P λ exists and is unique if the eigenvalues,

α k , k = 1, • • • , n
of A(λ), satisfy α i + α j 0 for all pair (i, j):

A(λ) T P λ + P λ A(λ) + Q(λ) = 0. ( 18
)
Definition 2: For a fixed M > 0 (arbitrary large), let us define the set:

Λ + (M) = {λ ∈ Λ : P λ > 0, max(eig(P λ )) ≤ M}, ( 19
)
where eig(P) stands for the set of eigenvalues of the matrix P.

The following proposition holds:

Proposition 1: If Q i > 0,
for all i ∈ S then Λ + (M) is a compact set and there exists a α min = min λ∈Λ + (M) min(eig(P λ )) > 0.

Proof: Λ + (M) is obviously bounded. Let us show it is a closed set. Suppose that there exists a sequence λ k ∈ Λ + (M) such that min(eig(P λ k )) → 0. As Λ is a compact set, there exists a subsequence that converges to an element λ ∈ Λ and such that min(eig(P λ )) = 0. Denote by v 0 the associated eigenvector. It follows the contradiction from 2v

T 0 P λ A(λ)v 0 + v T 0 Q(λ)v 0 = v T 0 Q(λ)v 0 = 0 which is absurd since Q(λ) = i∈S λ i Q i > 0.
Then, Λ + (M) is a closed set and there exists α min > 0 such that ∀λ ∈ Λ + (M), min(eig(P λ )) ≥ α min .

It can now be introduced the following min-type function and proposition:

Proposition 2: The function defined by

x → V min (x) := min λ∈Λ + (M) x T P λ x, ( 20
)
is locally Lipschitz continuous and proper.

Proof: The solution of ( 18) are always well defined on Λ + (M) and vary continuously with λ. Then, by compacity and for any fixed x, the minimum exists and inf λ∈Λ + (M) x T P λ x = min λ∈Λ + (M) x T P λ x.

Moreover, the familly of function indexed by λ ∈ Λ + (M) and defined by x → x T P λ x is uniformly locally Liptchitz on any ball B(0,R) with a Lipschitz constant equals to K := 2MR. This is obvious since for all λ ∈ Λ + (M),

α min x 2 ≤ inf λ∈Λ + (M) x T P λ x ≤ x T P λ x ≤ M x 2 . (21)
So, as any x T P λ x have a Lipschitz constant equals to K on a ball of radius R, we can conclude that it is also the case for the minimum. Indeed, for any pair (x, y) (∈ B(0, R)), there exists a pair (λ 1 , λ 2 ) ∈ Λ + (M) 2 such that:

|V min (x) -V min (y)| = |x T P λ 1 x -y T P λ 2 y|. ( 22
)
Thus, by continuity, there exists z on the line segment [x y] such that z T P λ 1 z = z T P λ 2 z and it follows:

|V min (x) -V min (y)| ≤ |x T P λ 1 x -z T P λ 1 z| + |z T P λ 2 z -y T P λ 2 y| (23) ≤ K x -z + K z -y = K x -y . ( 24
)
The min type function V min is also proper in view of the left part of inequality [START_REF] Cortes | Discontinuous dynamical systems[END_REF]. is achieved and I(x) as the set of indices where the minimum over i ∈ S of x T (P λ (x)

-P i )A i x, λ ∈ L(x)
, is reached:

I(x) := arg min i∈S ,λ∈L(x)
x T (P λ -

P i )A i x. ( 25
)
The following lemma is useful latter in order to define an appropriate descent direction for

V min :
Lemma 1: For every λ ∈ Λ, such that equation ( 18) admits a solution, min i∈S x T (P λ -P i )A i x is non positive.

Proof: For a given λ ∈ Λ, left-multiplying by x T and right-multiplying by x the equation ( 18) and using Q(λ) = i∈S λ i Q i and A T i P i + P i A i + Q i = 0, we get:

i∈S λ i x T (P λ -P i )A i x = 0. ( 26
)
So, we cannot have x T (P λ -P i )A i x > 0 for every i ∈ S because in this case the left-hand member of equality (26) would be positive.

We can now state one of the main result of this paper:

Theorem 6: Assume there exists at least a λ ∈ Λ such that A(λ) is Hurwitz. Then, the state feedback switching rule defined by:

i * (x) ∈ I(x) = arg min i∈S ,λ∈L(x)
x T (P λ -P i

)A i x, ( 27
)
stabilizes the switching system (1) with a cost smaller than 1 2 V min (x 0 ). Proof: Let λ ∈ Λ such that A(λ) is Hurwitz, then the Lyapunov equation ( 18) admits a (unique) solution P λ > 0. The matrix A(λ 0 ) is still Hurwitz if λ 0 belongs to a sufficiently small neighborhood of λ and so equation ( 18) admits a unique solution for every λ 0 in some neighborhood of λ. Thus, the interior of Λ + is non empty and V min is well defined.

Let us now consider the directional derivative of V min (x(t)) with respect to ẋ = A(ν)x, ν ∈ Λ as:

V ′ min (x; ν) = lim h→0;h>0 V min (x + hA(ν)x) -V min (x) h . ( 28
)
Denote by g λ (x) := x T P λ x the continuous function defined on Λ + (M) × R n . Then, as the following conditions are met: 1. the set Λ + (M) is compact, 2. for each λ ∈ Λ + (M), the directional derivative g ′ λ (x, •) exists, 3. ∀w ∈ R n , g λ (x+hw)-g λ (x) h → g ′ (x; w) uniformly in λ as h → 0, it can be concluded (from Theorem 6.1, p. 350-353 in [START_REF] Furukawa | Optimality conditions in nondifferentiable programming and their applications to best approximations[END_REF]),

V ′ min (x; ν) = min λ∈L(x) 2x T P λ A(ν)x. ( 29
)
Taking ν = i * (x) ∈ I(x), it follows using Lemma 1 that:

V ′ min (x, i * (x)) = min λ∈L(x) 2x T P λ A i * (x) x ≤ -x T Q i * (x) x < 0, x 0.
Therefore, for any initial condition x 0 ,

V min (x(t)) + t 0 x(τ) T Q i * (x(τ)) x(τ)dτ ≤ V min (x 0 ), ∀t ≥ 0.
As Q i > 0, ∀i ∈ S and as V min (•) is proper, it follows that: V min (x(t)) → 0 and x(t) → 0 when t → +∞.

IV. Discussion concerning the switching law and its optimality

In this part we want to discuss the degree of optimality of the provided switching law.

Observe first that in the case where all matrices A i are Hurwitz, then the matrices P i are definite positive for each mode i and the Lyapunov function V min satisfies always the following inequality: V min (x) ≤ min i∈S x T P i x. One can also observe that for a given initial state x, the value 1 2 V min (x) is the best cost related to every constant convex combination that stabilizes the relaxed system. As we will see, this is an important point if a singular control occurs in the solution. Now, in the general case, when can we say that 1 2 V min (x) is optimal? The answer is :'Along the part of trajectories where the optimal control λ * is constant to reach the origin".

In order to prove this last sentence, consider an optimal solution and assume that the optimal control is piecewise smooth with a countable number of discontinuities occurring at time t 1 , t 2 , .... Then, on time interval (t k , t k+1 ), the optimal control λ * is smooth and we can denote by Ā(t) = A(λ * (t)), Q(t) = Q(λ * (t)), the time varying matrices. The co-state p(t) can be identified as for a time varying linear system to P(t)x(t) (as suggested by Theorem 5) October 17, 2013 DRAFT where the time varying matrix P(t) satisfies the Lyapunov equation:

-Ṗ = ĀT P + P Ā + Q. ( 30
)
As it is mandatory that the hamiltonian H = p T Āx + 1 2 x T Qx remains equal to zero, it implies that along the optimal trajectory: x T Ṗx = 0. Thus, H can be rewritten as in the case of Equation (26) as follow:

H = i∈S λ * i x T ( P -P i )A i x = 0, (31)
where the optimal control λ * satisfies the complementarity constraint:

0 ≤ λ * i ⊥ x T ( P -P i )A i x ≥ 0, i ∈ S . (x ⊥ y means xy = 0) (32)
Theorem 7: If after a given time instant t, the optimal control λ * stays constant on (t, +∞)

and if P λ * > 0 then the switching law provided by Theorem 6 coincides with the optimal one.

Proof: This is obvious since as in the classical LQ problem and after the last switch, the matrix P is given by the solution P λ * of the Lyapunov equation ( 18) where λ * stands for the constant optimal control (eventually singular). Thus, as it is assumed that P λ * > 0, by optimality, x T P λ * x = min λ∈Λ:P λ >0 x T P λ x and it can be concluded that λ * ∈ L(x). Along the optimal trajectory, the optimal switching condition is now based on the minimization of

x T (P λ * -P i )A i x, i ∈ S which is exactly the switching condition provided by the switching rule of Theorem 6. The resulting trajectories are then the same if λ * is admissible for the switched system. Else a chattering motion occurs around the sliding surface defined by conditions (32) and the trajectory can be considered to be similar of the optimal one provided that the switching frequency is sufficiently high [START_REF] Bengea | Optimal control of switching systems[END_REF].

Remark 5: Theorem 7 applies at least for two important generic cases: when the switching sequence is finite, and, in dimension 2, along the singular arcs which are necessarily defined by lines passing through the origin and by a constant control [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF].

Formally, we can justified the design of V min as follow. Assuming known the value function, one can write for any T > 0, V * (x 0 ) = min σ

1 2 T 0 x T (t)Q σ(t) x(t)dt +V * (x(T )).
The transversality condition of PMP implies at time T , p(T ) = ∂V(x(T )) ∂x (if exists). Now suppose that V * (x(T )) is approximated by V min . Then, an approximate of p(T ) is given by p(x(T )) ≈ P λ x(T ) with λ ∈ L(x). Thus, it is easy to check that the minimization of the Hamiltonian at time T leads October 17, 2013 DRAFT to the switching law (27). As the problem is homogenous and if the approximation is "good", one can infer that p(x) ≈ P λ(x) x with λ(x) ∈ L(x) on the whole state space.

V. Illustrative examples

Before presenting some examples, it is important to mention that it is not necessary to ensure a stabilizing switched law to determine all the possible values of the set Λ + (M). Only one value is sufficient to guarantee the stability. So, a reasonable finite number of values ensures performances. A possibility to get this finite number is to discretize the set Λ.

A. Example 1 : a regular case

In [START_REF] Riedinger | Linear quadratic optimization for hybrid systems[END_REF], we have proposed a periodic solution that approximated the optimal one for a switched system defined by two modes. The matrices corresponding to this example are:

A 1 =               
2.15 0.66 0.66 1.33

         .
For this example, there is no singular arc and the solution is always regular (and, thus, admissible for the switched system). Figure 1 shows the state space trajectories for the switching law given by ( 27) and the optimal one. The later is obtained by NL programming in a suitable formulation taking into account singular arcs [START_REF] Riedinger | A numerical framework for optimal control of switched input affine nonlinear systems subject to path constraint[END_REF]. If not, numerical difficulties in the control determination are often encountered. We can see that the two solutions match well together. Thus, even if no argument has been advanced in this case, it seems that the provided switching conditions are close to the optimal. Figure 2 compares the optimal cost with the costs obtained by using the switching law, only mode 1 and only mode 2, respectively. This comparison is made for initial states taken on the unit ball, the x-axis represent the angle θ. It can be observed that the cost associated to the switching law coincides the cost of the optimal numerical solution. Of course, the essential difference is that the numerical solution is an open loop control while the switching law defines a closed loop control.

B. Example 2 : a singular case

Let us take an example (chosen randomly) for which singular arcs occur in the solutions. The matrices for this example are: A

(! ("#$ ("#% ("#& ("#' " "#' "#& "#% "#$ ! (! ("#$ ("#% ("#& ("#' " "#' "#& "#% "#$ ! ()*)+,(-*.+
1 =          -0.96 -9.93 0.68 -5.08          , A 2 =          -1.10 0.95 1.88 -2.76          Q 1 =         
1.17 0.90 0.90 5.60

         , Q 2 =          2.11 -1.07 -1.07 2.04          .
As shown in Figures 3 and4, the result is clearly near optimal. By comparison, a min switching strategy defined by: i(x) = arg min i∈S {x T P i x} leads to a cost given by 1 2 min i∈S {x T P i x} which is clearly not optimal as the two examples show.

C. Example 3 : extension to controlled linear switched systems

Consider the controlled linear switched systems in continuous time:

ẋ(t) = A σ(t) x(t) + B σ(t) u σ(t) (t), x(0) = x 0 , (33)
where σ : [0, +∞) → S = {1, • • • , s} denotes the switching law that selects the active mode at time t by choosing among a finite collection of controlled linear systems. For i ∈ S , A i ∈ R n×n , B i ∈ R n×m i , x(∈ R n ) is the state and u i (t) ∈ R m i (m i ≤ n) is the control associated to mode i.

We assume that the switched system is stabilizable but eventually all the pair (A i , B i ), i ∈ S cannot be stabilized.

         .
Clearly the controllable subspace is generated by the canonical basis vector {e 1 = (1, 0) T } for the first system while it is {e 2 = (0, 1) T } for the second. Thus, following the first step, an LQ design is used on each controllable subspace to determine two static gains K 1 and K 2 such that ∞ 0 x T Q i x + u T i R i u i is minimized by u i = -K i x. We have chosen R i = Q i = 1, i = 1, 2 then K 1 = K 2 = 0.2361. Thus, the switched linear system is now defined by: ẋ = (A i -B i K i )x, i = 1, 2.

The second step consists in applying the switching law with weight matrices Qi . For this example, we have chosen : Qi = Id + K T i R i K i , i = 1, 2. The Figures [START_REF] Bengea | Optimal control of switching systems[END_REF] and [START_REF] Axelsson | Gradient descent approach to optimal mode scheduling in hybrid dynamical systems[END_REF] show once again that the obtained switching law is optimal by comparison with the numerical solution.

It can be noticed that the solutions of (18) for mode 1 and 2 are not positive definite as expected, see Figure [START_REF] Axelsson | Gradient descent approach to optimal mode scheduling in hybrid dynamical systems[END_REF].

VI. Conclusion

A state feedback switching law based on control Lyapunov function for switched LQ regulator problems in continuous time has been proposed. The stabilizing feedback can be applied even if the subsystems are all unstable. The only condition that is required, is the existence of at least a stable convex combination of the subsystems. Even if the exact optimal solution is not determined, we have shown that the switching conditions involved by the switching law can be optimal. More precisely, this generic situation occurs along arcs (singular or not) ending to the origin with a constant optimal control. As in dimension two, the singular controls are constant, the switching law is really optimal in most of encountered 2-dimensional examples. In any case, a guarantee on the cost is provided by an upper bound given by the value of the Lyapunov function.

Definition 3 :

 3 Let us define L(x) as the set of values λ ∈ Λ + (M) where the minimum V min (x)

Fig. 1 . 1 Fig. 2 .

 112 Fig. 1. Example1: State space trajectories: (red) optimal solution (NLP); (blue) switching law

Fig. 3 .Fig. 4 .

 34 Fig. 3. Example 2 : State space trajectories: (red) optimal (NLP); (blue) switching law

  of the proposed switching law consists to solve:

Fig. 5 .Fig. 6 .

 56 Fig. 5. Example 3 : State space trajectories: (red) optimal (NLP); (blue) switching law

[START_REF] Seatzu | Optimal control of continuous-time switched affine systems[END_REF] DRAFT