Y Moudden

J Bobin

J.-L Starck

J Fadili

DICTIONARY LEARNING WITH SPATIO-SPECTRAL SPARSITY CONSTRAINTS

Devising efficient sparse decomposition algorithms in large redundant dictionaries has attracted much attention recently. However, choosing the right dictionary for a given data set remains an issue. An interesting approach is to learn the best dictionary from the data itself. The purpose of this contribution is to describe a new dictionary learning algorithm for multichannel data analysis purposes under specific assumptions. We assume a large number of contiguous channels as in so-called hyperspectral data. In this case it makes sense to consider a priori that the collected data exhibits sparse spectral signatures and sparse spatial morphologies in specified dictionaries of spectral and spatial waveforms. Building on GMCA, the proposed algorithm gives a practical way to enforce the additional a priori spectral sparsity constraint on the dictionary space. Numerical experiments with synthetic and real hyperspectral data illustrate the efficiency of the proposed algorithm.

INTRODUCTION

Generalized Morphological Component Analysis (GMCA) is a recent algorithm for multichannel data analysis described in [START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF][START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF]. It was derived as a multichannel extension to MCA [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA)[END_REF] and as such it is a powerful algorithm for the sparse decomposition of multichannel data over multichannel dictionaries. A major feature of the GMCA algorithm, in comparison to other sparse decomposition algorithms, is that it alternates sparse decomposition steps with dictionary learning steps, in a fast iterative thresholding loop with a progressively decaying threshold leading to a very robust salient to fine estimation process. In fact, learning dictionaries for the sparse representation of given data sets is now a growing field of interest [START_REF] Gribonval | Dictionary Identifiability from Few Training Samples[END_REF][START_REF] Rubinstein | Learning sparse dictionaries for sparse signal representations[END_REF] although it is worth noting that this problem has a long history [START_REF] Olshausen | Sparse coding of natural images produces localized, oriented, bandpass receptive fields[END_REF][START_REF] Cardoso | Some experiments on independent component analysis of non-gaussian processes[END_REF][START_REF] Hyvärinen | Sparse code shrinkage: denoising by nonlinear maximum likelihood estimation[END_REF]. Obviously, it is required to set a priori some constraints on the dictionary space in which the optimal dictionary is to be looked for. In the case of GMCA, it is assumed that the size of the multichannel dictionary Ω, i.e. the number of atoms n × t ′ , is specified beforehand and that the multichannel atoms are all rank one matrices, product of a spectral signature a k ∈ R m,1 and a spatial density profile φ k ∈ R 1,t . The column vectors a k are grouped into a matrix noted A ∈ R m,n and the line vectors φ k are taken from matrix Φ ∈ R t ′ ,t . GMCA models the data X ∈ R m,t as follows :

X = AνΦ + N = k k ′ ν k ′ k a k φ k ′ + N (1)
where the entries of the sparse matrix of coefficients ν representing X in the multichannel dictionary Ω = A ⊗ Φ are noted ν k ′ k and N ∈ R m,t is included to account for Gaussian instrumental noise or modeling errors. GMCA further assumes that the dictionary of spatial waveforms Φ is known beforhand while the spectral components A, also called the mixing matrix in blind source separation (BSS) applications, is learned from the data. The image from the p th channel is represented here as the p th row of X, x p . The successful use of GMCA in a variety of multichannel data processing applications such as BSS [START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF], color image restoration and inpainting [START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF] motivated research to extend its applicability. In particular, there are instances where one is urged by additional prior knowledge to further constrain the dictionary space. For instance, one may want to enforce equality constraints on some atoms, or the positivity or the sparsity of the learned dictionary atoms. Builiding on GMCA, the purpose of this contribution is to describe a new dictionary learning algorithm for so-called hyperspectral data processing. Hyperspectral imaging systems collect data in a large number (up to several hundreds) of contiguous regions of the spectrum so that it makes sense to consider for instance that some physical property will show some regularity from one channel to the next. In fact, the proposed algorithm, referred to as hypGMCA, assumes that the multichannel atoms to be learned from the collected data exhibit diversely sparse spatial morphologies as well as diversely sparse spectral signatures in specified dictionaries Φ ∈ R t,t ′ and Ψ ∈ R m,m ′ of respectively spatial and spectral waveforms. The proposed algorithm is used to learn from the data rank one multichannel atoms which are diversely sparse [START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF] in a given larger multichannel dictio-

nary.

In what follows, regardless of other models living in other scientific communities, the term hyperspectral denotes multichannel data following model [START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF] with the above two specific properties i.e. that the number of channels is large and that these achieve a regular sampling of some additional and meaningful physical index (e.g. wavelength, space, time) which we refer to as the spectral dimension. We describe next the proposed modified GMCA which we devised to account for the a priori sparsity of columns a k in Ψ, a given dictionary of spectral waveforms. Accounting for this prior requires a modified objective function, discussed in section 2. The resulting hypGMCA algorithm is given in section 3. Finally, numerical experiments in section 4 demonstrate the efficiency of the proposed method.

OBJECTIVE FUNCTION

With the above assumptions, equation ( 1) is rewritten as follows :

X = k X k + N = k Ψγ k ν k Φ + N = ΨαΦ (2)
where X k = a k s k are rank one matrices sparse in Ω = Ψ⊗Φ such that a k has a sparse representation γ k in Ψ while s k has a sparse representation ν k in Φ. In a BSS context, the rows s k of S are commonly referred to as sources. For the sake of simplicity, we assume that Ψ and Φ are orthonormal bases ant that the noise N is uncorrelated inter-and intrachannel with variance σ 2 . Extensions to more general cases are readily derived. Denote α k = γ k ν k the rank one matrix of coefficients representing X k in Ω . Looking at the above matrix equation column-wise, each column of Ψα = XΦ T is represented as a sparse linear combination of columns in A with coefficients in the corresponding column of ν = SΦ T . Indeed, the assumption that the lines of ν are diversely sparse (e.g. sparse and independent, or sparse with disjoint support, etc.) results in the columns of ν also being sparse. Usual sparse decomposition algorithms work on each column of XΦ T separately to seek its sparsest representation in A. GMCA comes into play when the dictionary is not fully known beforehand and A needs to be learned from the data. GMCA constrains the dictionary space assuming that matrix A is unknown yet of specified size m × n and that Φ is given. This leads to a joint sparse coding and dictionary learning objective which can be expressed as a minimization problem in augmented Lagrangian form using an ℓ 1 sparsity measure :

min A,ν K k=1 λ k ν k 1 + X - K k=1 a k ν k Φ 2 2 (3)
which is clearly a difficult non-convex optimization problem. Nonetheless, the GMCA algorithm is able to provide a practical approximate solution as reported in [START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF]. Problem (3) is readily interpreted as a MAP estimation of the model parameters A and ν where the ℓ 1 penalty terms imposing sparsity come from a Laplacian prior on the sparse coefficient matrix ν.

Building on GMCA, we want here to learn spectral dictionary A so that again on average the columns of XΦ T are simultaneously sparse in A and, what is more, the columns of A are sparse in Ψ.

A well known property of the linear mixture model ( 2) is its scale and permutation invariance. A consequence is that unless a priori specified otherwise, information on the separate scales of γ k and ν k is lost, and only a joint scale parameter for γ k , ν k can be estimated. This needs to be translated into a practical prior on γ k ν k . We propose here that the following p π is a good and practical candidate joint sparse prior for γ k and ν k after the loss of information induced by multiplication :

p π (γ k , ν k ) ∝ exp(-λ k γ k ν k 1 ) ∝ exp(-λ k i,j |γ k i ν j k |) (4)
where γ k i is the i th entry in γ k and ν j k is the j th entry in ν k . The proposed distribution has the nice property, for subsequent derivations, that the conditional distributions of γ k given ν k and of ν k given γ k are Laplacian distributions which are commonly and conveniently used to model sparse distributions. Finally, inserting the latter prior distribution in a Bayesian MAP estimator leads to the following minimization problem, expressed in coefficient space :

min {γ k ,ν k } 1 2σ 2 X - k Ψγ k ν k Φ 2 + k λ k γ k ν k 1 (5)
Let us first note that the above can be expressed slightly differently as follows :

min {α k } 1 2σ 2 X -k X k 2 + k λ k α k 1 with X k = Φα k Ψ and ∀k, rank(X k ) ≤ 1 (6)
which uncovers a nice interpretation of our problem as that of approximating the data X by a sum of rank one matrices X k which are sparse in the specified dictionary of rank one matrices. This is the usual ℓ 1 minimization problem [START_REF] Donoho | Optimally sparse representation in general (non-orthogonal) dictionaries via ℓ 1 minimization[END_REF] but with the additional constraint that the X k are all rank one at most. The latter constraint is enforced here mechanically through a proper parametric representation of X k = a k s k or α k = γ k ν k . A similar problem was previously investigated by [START_REF] Zhang | Low-rank approximations with sparse factors I: Basic algorithms and error analysis[END_REF] with a very different approach.

We also note that rescaling the columns of ν ← ρν while applying the proper inverse scaling to the rows of S ← 1/ρS, leaves both the quadratic measure of fit and the ℓ 1 sparsity measure in equation ( 5) unaltered. Although renormalizing is still worthwhile numerically, it is no longer dictated by the lack of scale invariance of the objective function and the need to stay away from trivial solutions, as in GMCA.

Finally, adopting a BSS point of view, the objective function ( 5) is fully symmetric in its treatment of the mixing matrix A and the source processes S. The great majority of BSS methods invoke a uniform improper prior distribution for the spectral parameters A. Truly, A and S often have different roles in the model and very different sizes. However, dealing with so-called hyperspectral data, such an asymmetry is questionable. There have been previous reports of a symmetric treatment of A and S for BSS [START_REF] Stone | Spatiotemporal independent component analysis of event-related fmri data using skewed probability density functions[END_REF][START_REF] Hyvärinen | Imposing sparsity Fig. 7. Estimated spatial concentration maps of H 2 O (left) and CO 2 (right) ice. on the mixing matrix in independent component analysis[END_REF] however in the noiseless case. We also note that very recently, the objective function ( 5) was proposed in [START_REF] Rubinstein | Learning sparse dictionaries for sparse signal representations[END_REF]. However the algorithm used in [START_REF] Rubinstein | Learning sparse dictionaries for sparse signal representations[END_REF] is very different from the method proposed here which benefits from all the good properties of GMCA, notably its speed and robustness which come along the iterative thresholding with a decreasing threshold.

ALGORITHM

Unfortunately, there is no obvious closed form solutions to optimization problem (5) which is again clearly nonconvex. Similarly to the GMCA algorithm, we propose here a numerical approach by means of a block-coordinate relaxation iterative algorithm, alternately minimizing with respect to γ and ν. Indeed, thanks to the chosen prior, for fixed γ (resp. ν), the marginal minimization problem over ν (resp. γ) is convex and is readily solved using a variety of methods. Inspired by the iterative thresholding methods described in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Hale | A fixedpoint continuation method for l1 -regularized minimization with applications to compressed sensing[END_REF][START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF], akin to Projected Landweber algorithms, we obtain the following system of update rules :

ν (+) = ∆ η ν (-) + R ν α -γν (-) γ (+) = ∆ ζ γ (-) + α -γ (-) ν R γ (7)
where R ν and R γ are appropriate relaxation matrices for the iterations to be non-expansive. Assume left invertibil-ity of A and right invertibility of S. Then, taking R ν = γ T γ -1 γ T and R γ = ν T νν T -1 , the above are rewritten as follows :

ν (+) = ∆ η γ T γ -1 γ T α (8) γ (+) = ∆ ζ αν T νν T -1 (9)
where vector η has length n and entries η 8) and ( 9) rules are easily interpreted as thresholded alternate least squares solutions. Finally, in the spirit of the fast GMCA algorithm [START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF][START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF], it is proposed that a solution to problem ( 5) can be approached efficiently using the following symmetric iterative thresholding scheme with a progressively decreasing threshold, which we refer to as hypGMCA :

[k] = λ k γ k 1 / γ k 2 2 ,
1. Set the number of iterations I max and initial thresholds λ

(0) k 2. Transform the data X into α 3. While λ (h) k ≥ λ min
k , -Update ν assuming γ is fixed using eq. ( 8).

-Update γ assuming ν is fixed using eq. ( 9) . -Decrease the thresholds λ The salient to fine estimation process is again the core of hypGMCA. With the threshold successively decaying towards zero along iterations, the current sparse approximations for γ and ν are progressively refined by including finer structures spatially and spectrally, alternatingly. The final threshold should vanish in the noiseless case or it may be set to a multiple of the noise standard deviation as in common detection or denoising methods. Soft thresholding results from the use of an ℓ 1 sparsity measure, which comes as an approximation to the ℓ 0 pseudo-norm. Applying a hard threshold instead towards the end of the iterative process, may lead to better results as was noted experimentally in [START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF][START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF]. When non-unitary or redundant transforms are used, the above is no longer strictly valid. Nevertheless, simple shrinkage still gives satisfactory results in practice as studied in [START_REF] Elad | Why simple shrinkage is still relevant for redundant representations?[END_REF]. In the end, implementing the proposed update rules requires only a slight modification of the GMCA algorithm given in [START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF][START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF]. Where a simple least squares linear regression was used in the GMCA update for a k , the proposed update rule applies a thresholding operator to the least squares solution thus enforcing sparsity on the estimated spectral signatures as a priori desired.

NUMERICAL EXPERIMENTS

Toy experiment with synthetic data

In this section, we compare the ability of hypGMCA and GMCA to recover and learn the dictionary atoms that were used to synthesize a given data set, generated according to model [START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF]. Using the language of BSS, we consider synthetic 2D data consisting of m = 128 mixtures of n = 5 image sources. The sources were drawn at random from a set of 128×128 structured images shown on Figure 1. These images provide us with 2D structured processes which are sparse enough in the curvelet domain [START_REF] Starck | The curvelet transform for image denoising[END_REF]. The spectra were generated as sparse processes in some orthogonal wavelet domain given a priori. The wavelet coefficients of the spectra were sampled from a Laplacian probability density with scale parameter µ = 1. Finally, white Gaussian noise with variance σ 2 was added to the pixels of the synthetic mixture data in the different channels. The graph on figure 3 traces the evolution of C A = I n -P à † A 1 , which we use to assess the recovery of the spectral dictionary A, as a function of the SNR which was varied from 0 to 40dB. Matrix P serves to reduce the scale and permutation indeterminacy inherent in model ( 2) and à † is the pseudo-inverse of the estimated spectral dictionary. In simulation, the true source and spectral matrices are known and so that P can be computed easily. Criterion C A is then strictly positive and null only if matrix A is correctly estimated up to scale and permutation. Finally, as we expected since it benefits from the added a priori spectral sparsity constraint it enforces, the proposed hypGMCA is clearly more robust to noise. A visual inspection of figure 4 a further qualitative assessment of the improved dictionary and source recovery provided by correctly accounting for a priori spatial as well as spectral sparsity. The top images were obtained with GMCA while the bottom images, were obtained with hypGMCA. In all cases, both methods were run in the curvelet domain [START_REF] Starck | The curvelet transform for image denoising[END_REF] with the same number of iterations.

Application to real data

We applied the proposed dictionary learning algorithm to hyperspectral data from the 128 channels of spectrometer OMEGA on Mars Express (www.esa.int/marsexpress), at wavelengths ranging from 0.93µm to 2.73µm. Example maps collected in four different channels are shown on figure 5. Model 2 is clearly too simple to describe this hyperspectral reflectance data set. Non-linear instrumental and atmospheric effects are likely to contribute to the true generative process. In any case, we use hypGMCA to learn a dictionary of spectral atoms, sparse in some orthogonal wavelet dictionary, that well represent the data. The spectral dictionary was learned assuming it is sparse in an orthogonal wavelet representation. Preliminary results are shown on figure 6. Interestingly, exactly two atoms were found to be respectively strongly correlated with reference spectra for H 2 O and CO 2 ice. Given the simplicity of the model, the close match between the learned spectra and the reference spectra is is remarkable. Figure 7 shows the corresponding spatial maps which were assumed sparse in an orthogonal wavelet basis. We note that the CO 2 ice appears spatially concentrated at the poles which is in close agreement with the results presented in [START_REF] Moussaoui | On the decomposition of mars hyperspectral data by ica and bayesian positive source separation[END_REF].

CONCLUSION

We described a new dictionary learning algorithm, hypGMCA, for sparse signal representation in the case where it is known a priori that the spatial and spectral features in the data have sparse representations in known dictionaries of template waveforms. The proposed method relies on an iterative thresholding procedure with a progressively decreasing threshold. This alone gives the method true robustness to noise. As expected, taking into account the additional prior knowledge of spectral sparsity leads to enhanced performance. Numerical experiments focus on a comparison between GMCA and hypGMCA. GMCA is compared to state-of-the-art dictionary learning and BSS algorithms in [START_REF] Bobin | Blind Source Separation: the Sparsity Revolution[END_REF][START_REF] Bobin | Morphological diversity and sparsity for multichannel data restoration[END_REF].

Fig. 1 .

 1 Fig. 1. Image data set used in the experiments.

2 .

 2 while ζ has length m and entries ζ[k] = λ k ν k 1 / ν k 2 The multichannel soft-thresholding operator ∆ η acts on each row k of ν with threshold η[k] and ∆ ζ acts on each column k of γ with threshold ζ[k]. Equations (

5 .

 5 Transform back γ and ν to estimate A and S.

  Figure 2 displays four typical noisy simulated mixture data with SNR = 20dB.

Fig. 2 .

 2 Fig. 2. Four 128 × 128 mixtures out of the 128 channels. The SNR is equal to 20dB.

  allows

Fig. 3 .

 3 Fig. 3. Evolution of the mixing matrix criterion C A as a function of the SNR in dB. Solid line : recovery results with GMCA. • : recovery results with hypGMCA.

Fig. 4 .

 4 Fig. 4. Left column : Estimated sources using the original GMCA algorithm. Right column : Estimated sources using the new hypGMCA.

Fig. 6 .

 6 Fig. 6. Left picture : Reference (solid line) and estimated (•) spectra for H 2 O ice. Right picture : Reference (solid line) and estimated (•) spectra for CO 2 ice.

Acknowledgments :

The authors are grateful to Olivier Forni for providing the hyperspectral data from Omega on Fig. 5. From left to right : Mars Express observations at wavelengths = 1.38 -1.75 -1.94 and 2.41µm. Mars Express.


Antoine Cizeron
email: antoine.cizeron@geeps.centralesupelec.fr

Hugo Milan

Javier Ojeda •

Olivier B Éthoux

A Cizeron •

H Milan •

J Ojeda

PWM-Induced Current Modelling in Stator Slots with Multiple Stacked Coils

This paper deals with the PWM-induced current and losses in a specific segmented winding structure. The proposed segmentation process enables to split a winding into several coils. These latter are supplied independently by H-bridge converters and are wound around the same magnetic circuit. This process leads to a deeper segmentation of electric drives for enhanced modularity and reduced voltage rating. The strong magnetic coupling between each coil is described, and the control degrees of freedom are presented. This study provides a model based on an analytical method and on an equivalent electrical circuit calibrated through experimental results. A trade-off is found between the losses related either to the distribution of the fundamental component of currents or to the switching power converter supply.

Introduction

Electric drive generally involves a power electronics converter associating a DC-link voltage to an electrical machine. In standard EV drive solutions, the voltage rating is between 400 V and 800 V . These voltage values meet the requirements of fast charging and are enabled by the new generation of WBG semi-conductors power components. The implementation of these solutions relies on one high power con-verter to feed the whole drive. However, the resiliency of this conventional structure is limited, and its power rate implies a separate thermal management for the converter and the electrical machine. Integrated solutions are proposed to bring the power electronics closer to the electrical machine. In these Integrated Modular Motor Drives (IMMD), the power electronics components are spatially distributed around the motor [START_REF] Hussein | Design of an Integrated DC-Link Structure for Reconfigurable Integrated Modular Motor Drives[END_REF]. This is facilitated by the use of multiple elementary converters to supply each portion of the motor winding.

Furthermore, such segmented structure enables the use of low voltage semi-conductor devices, that generally present reduced losses compared to their high voltage counterparts. The various control degrees of freedom available in those segmented structure offer promising opportunities in terms of energy management and resiliency.

The patent proposed in [START_REF] Hoang | Electric machine supplied at low voltage and associated multicellular power train[END_REF] extends the segmentation limit of the motor winding to the turn limit and enables the use of several stacked coils inside the same slots. [START_REF] Ben | Electric Traction Chain with Segmented Power Supply[END_REF] details the scientific and technical obstacles related to this innovative structure. Recent works like [START_REF] Zhang | Design and Performance Analysis of Segmented Three-Phase IPMSM for EVs Integrated Battery Charger[END_REF] expose the use of a similar structure for an integrated battery charger. In [START_REF] Zhang | Vibration reduction for dual-branch three-phase permanent magnet synchronous motor with carrier phase-shift technique[END_REF], authors combine this kind of segmentation with external coupled inductors to reduce mechanical vibrations. These references only consider balanced current distribution between the different winding sets during the motoring mode. The present paper addresses this issue in a representative geometry. It provides the modelling tools enabling to minimize the losses thanks to uneven current distribution between the different coils.

Each elementary coil being independently supplied, the whole structure presents new control degrees of freedom, these latter were investigated in the particular case of two coils with strong magnetic coupling [6]. A PWM-induced current ripple modelling method is detailed in [START_REF] Chang | Prediction and Evaluation of PWM-Induced Current Ripple in IPM Machines Incorporating Slotting, Saturation, and Cross-Coupling Effects[END_REF] for conventional synchronous motor. It provides accurate estimation of the current waveforms but is based on an inductance model. Other models are used to understand the PWM-induced losses [START_REF] Chang | Comparative analysis of pwm power losses in ipm machines with different modulation schemes using wide-bandgap-based inverters[END_REF][START_REF] Cittanti | Pwm-induced losses in electrical machines: An impedance-based estimation method[END_REF] but are not suited for taking into account the new degrees of freedom. Moreover, the influence of the air-gap on the current ripple induced by the Pulse Width Modulation (PWM) in the aforementioned stacked coils has not been yet studied. This paper emphasizes the link between the control parameters, namely the differential voltages, and the current ripple by modelling two stacked coils within the same magnetic circuit. Measurements are performed to validate the fundamental electrical model and the high frequency model, which estimates the PWM-induced losses over a wide range of operating points. The rest of this paper is structured as follows. Section 2 introduces the concept of the segmentation of ordinary winding into several stacked coils and the case of study. Then, section 3 describes the analytical model applied to the fundamental frequency. Finally, section 4 evaluates the aforementioned trade-off between PWM-induced losses and the losses related to the fundamental distribution of currents in the stacked coils. The last section provides conclusions and perspectives.

2 Stacked coils sharing the same magnetic circuit

Generic Problem

Fig. 1 shows a generic geometry symbolizing one phase within a concentrated winding arrangement. The conventional winding presents a number of electro-conductive turns written as N t . These turns are wound around one stator tooth in a concentrated way as shown in Fig. 1a. This winding is then split into n electrically independent coils that are stacked around the same stator tooth.

In the considered structure, each elementary coil is supplied by a low voltage full-bridge converter. This topology can be seen as an extension of the single tooth inverter integration proposed in [START_REF] Reber | Lowinductive sic h-bridge for direct-inverter-single-tooth-integration (german bmbf public funded research project 'verse')[END_REF]. Whereas the conventional winding Fig. 1a requires a H-bridge supply with a DC-link voltage V DC , each elementary coil Fig. 1b only requires a H-bridge supply with a reduced DC-link voltage of V DC /n in order to meet the same power rating of the whole structure (Fig. 2). The magnetic circuit was initially designed for the conventional winding. Therefore, the ampere-turns are kept constant to guaranty the proper use of the ferromagnetic material. Thus, the current rating remains the same in both structures since the number of turns per coils is equivalent.

Compared to the conventional supply, the proposed structure presents additional control degrees of freedom since each H-Bridge inverter can be controlled independently. The PWM is used to generate the low frequency voltages applied to each coil. Thus, two frequency bands can be distinguished in this problem -The Low Frequency band (LF : 0 to 1 kHz) is related to the fundamental frequency of the electrical machine. It corresponds to the back electromotive force (emf) and the fundamental of current in the winding. -The High Frequency band (HF : > 10 kHz) is related to the PWM signals applied to the winding.

The behaviour of the magnetic circuit and the conductors is different in those two bands. Even if the current magnitude is usually higher in the LF band than in the HF band, [6] has shown that the HF current ripple increase in stacked coils should also be taken into account. Thus two separate models are described, the LF model, based on an analytical method assesses the losses induced by a given current distribution between the stacked coils. Then, it defines the required voltages and so the PWM signals that are injected in the HF model in order to estimate the PWM-induced current ripple.

Study case of two stacked coils

A practical model was built. It is composed of two stacked coils made of 20 turns each wound around a laminated steel magnetic circuit including an air-gap. Its geometry is detailed in Fig. 3 and used for the analytical model in the LF band. The two coils are supplied either by SI MOSFET H-Bridge converters or linear amplifiers, thus enabling to investigate separately the LF and HF phenomena.

LF electrical model and losses assessment

LF electrical equivalent circuit

Generally, each coil k could be represented as a serial resistance r k and leakage inductance l k . One mutual inductance L m is used for all the stacked coils as they are wound around the same stator tooth. The air-gap ag influences the leakage inductance l k of each coils. The coil 1 located close to the air-gap barely has no leakage inductance value whereas coil n embraces all the flux lines that are crossing the slot and thus presents the highest leakage inductance value. The serial DC resistance is the same in all coils as they present the same number of turns and thus the same conductor length. This model rely on the fact that it exists two kinds of flux lines. The first contour C ag corresponds to the flux lines crossing the air-gap vertically (following the y-axis), the other contour C slot is assumed to cross the slot horizontally (following x-axis). These two paths are represented in Fig. 4a. The induction value in the air-gap and in the slot (between the conductor and the air-gap) can be approximated by:

B ag = µ 0 • i 2 • t ag and B slot = µ 0 • i w slot (1)
where t ag is the air-gap thickness and w slot is the width of the slot. Thus, for a given position of the conductive turn following y-axis, the inductance can then be computed. The analytical results are compared with the results of a Finite Element Analysis (FEA) method in Fig. 4b, both models are based on the parameters described in Tab. 1. The analytical model represents accurately (≈ 5% error) the influence of the conductor position on its inductance value according to the FEA results.

Application to the case of study

The analytical model is used to identify the terms of the LF electrical modelling in the case of two stacked coils (Fig. 3). Considering the computed inductance values of one conductor within the slot (Fig. 4b), the worst case in terms of imbalance corresponds to take the extreme values for the self inductances of stacked coils 1 (Y = 30 mm) and 2 (Y = ). The two flux paths taken into account leads to consider that all the flux lines embraced by a conductor close to the air-gap are also embraced by any other conductor. In the case of two stacked coils, this leads to set l 1 to zero (L 1 ≈ L m ) and to admit that the leakage inductance of the coil 2 is the difference between the self inductances :

l 2 = L 2 -L 1 = 40 µH while L m = 180 µH.
By supplying successively each stacked coil with a linear amplifier ensuring sinusoidal current, it is possible to identify the inductance and resistance values of the electrical equivalent circuit (Fig. 5). The measured voltages and currents are represented in Fig. 6 and the electrical parameters are reported in Tab. 2. The resistance analytical value is computed taking into account the conductor length and section of the practical model. The results are in close agree- ment, experimental values of inductance are slightly higher than analytical values since the model neglects the inductive part of conductors outside the 2D slot geometry. The resistance measurement used here is not accurate since it is based on the temporal delay observed between u 1 and u 2 in the experimental results (Fig. 6), however, the results are in the same order of magnitude. The analytical model provide the parameters of the LF equivalent electrical circuit with an appropriate accuracy and enables to estimate the voltages that are ought to be applied to the stacked coils for any current distribution (i 1 , i 2 ). These voltages are then applied to the model that computes the PWM-induced current ripple and losses.

LF losses assessment

The expression of the induction field values in the slot and air-gap can be used to assess the losses in the two materials, namely the copper conductors and the laminated iron core. Even when multiple conductors are considered, the same flux line paths (Fig. 4a) can be taken into account to compute the induction values between each conductors as in [START_REF] Hussein | Design of an Integrated DC-Link Structure for Reconfigurable Integrated Modular Motor Drives[END_REF]. Then, within the coils, the magnetic field is ruled by the following equation:

∆ ⃗ H = σ µ 0 • ∂ ⃗ H ∂t ( 2
)
where σ is the copper conductivity. The problem is supposed to have its only dependence following the y-axis. The influence of the position of the conductor within the slot width (following x-axis) is neglected. Equation (2) can then be projected into a 1D problem:

∂ 2 H x ∂ y 2 = σ µ 0 ( jω)H x . (3)
By computing the magnetic field values within each conductors, the current density induced by the flux linkage can be obtained. The losses in the copper are then computed by integrating the square of the current density on the copper volume.

Having the induction values in any position within the slot, the induction in the iron core can be obtained through the conservation of the induction flux. A Steinmetz model is then used to assess the core losses P Core in the geometry for any current distribution. The LF losses are then decomposed into three terms:

P LF = P Copper,DC + P Copper,supp + P Core ( 4
)
where P Copper,supp corresponds to the additional losses causes by the leakage flux in the slot through the conductors, and

P Copper,DC = R • (i 1,RMS ) 2 + (i 2,RMS
) 2 corresponds to the losses related to the DC resistance terms of the stacked coils. This latter is not related to the analytical model and strongly depends on the experimental setup serial resistance.

In the LF band, both copper conductors and iron core have frequency dependent behaviour, respectively through the terms P Copper,supp and P Core . These two terms also depend on the current distribution among the several stacked coils. The 20-turns coils are supposed to be arranged in two rows within the slot as Fig. 3a shows. Then, for each stacked coils, 10 conductors having twice the section of a practical conductor are used to represent the 20 turns. Among the LF losses, the terms P Copper,supp and P Core are computed with the analytical model represented in Fig. 4a. The term P Copper,DC is adjusted, taking the estimated value of resistance (30 mΩ ) into account.

The two coils are each powered by a dedicated linear amplifier. In this configuration, it is possible to impose the current magnitude in each coil and measure the overall losses related to the fundamental electrical frequency f e . For several values of fundamental frequency f e , the total peak current in the slot is fixed at 6 A and several operating points are tested in terms of current distribution between the two stacked coils. The current imbalance varies from -100 % (i 1 = 0 A) to +100 % (i 2 = 0 A). The losses computed from (4) are compared to measured results in Fig. 7.

The measured and estimated results are consistent, showing that a positive current imbalance (i 1 > i 2 ) could lead to minimize the LF losses in the magnetic circuit and its stacked coils. At low fundamental frequency, Fig. 7b shows that the experimental optimal current imbalance is mostly due to an eventual resistance disparity and could set to a non zero value. However, when the fundamental frequency increases, the optimal current imbalance clearly shifts toward positive values (i 1 > i 2 ) since the current in coil 1, because of its proximity to the air-gap, induces no induction field in the slot around the conductors of coil 2. However, this result is obtained under a linear power supply and has to be confirmed under PWM excitation. In order to reproduce the PWM-induced current ripple in stacked coils, the LF model is not sufficient. In the HF band, the copper conductors and the iron core present non negligible skin effect. Thus the inductive term is degraded and an additional impedance may be added in parallel to each inductive term to reflect the HF losses in those materials. This is visible on Fig. 8. The impedance z k represents the additional losses in the copper, which is directly related to the HF leakage inductance l ′ k . Impedance Z m symbolizes the core losses related to the HF mutual inductive term L ′ m . Considering the geometry and the frequency of the PWM excitation, an FEA method would require a refined mesh leading to high computation time. Admittance analysis can be used to observe the frequency dependent behaviour of the coils but is not accurate to estimate the current ripple and PWMinduced losses [6]. Therefore, experimental results are used here for the identification of the HF model parameters.

Practically, each of these terms depend on the LF operating point where the switching period occurs. However, in order to simplify this aspect and make rapid conclusions, this dependency is neglected. It is then assumed that this HF model is valid for any switching period based on the hypothesis of magnetic material linearity.

Voltages and HF currents are measured on the practical model with its two H-bridge supply (Fig. 9). The model described in Fig. 10 is calibrated on experimental waveforms of PWM-induced current ripples due to duty-cycle differ- ences (∆ α = α 1α 2 ̸ = 0). The obtained model precisely reproduces the current ripples and provides an estimation of losses (Fig. 9). Using this modelling, it is possible to compute the HF losses due to any duty-cycle configuration (α 1 , α 2 ).

In this way, the voltages provided by the LF model are applied to the stacked coils in order to meet a given current distribution. The duty-cycles used in each H-bridge are determined for any switching period. At each switching period, the voltages signals are generated from a centred PWM method. The harmonic component of these signals are then treated by the HF model in order to reproduce the PWMinduced current ripple. For each LF current distribution, the required duty-cycle values are used to assess the average HF losses over a fundamental period 1/ f e . The combination of these two models thus enables to estimate the PWM-induced losses required to reach a given current distribution. The total losses are computed with P tot = P Copper,DC + P Copper,supp + P Core + P HF [START_REF] Zhang | Vibration reduction for dual-branch three-phase permanent magnet synchronous motor with carrier phase-shift technique[END_REF] where P HF is the PWM-induced term averaged over one LF fundamental period. The terms P Copper,supp and P Core are the same than in the LF model (Fig. 7). The losses in inverters are included in P Copper,DC . This latter is computed with an equivalent resistance of 220 mΩ . It enables the comparison with the losses measured on the DC-link of the two inverters in the experimental setup.

Synthesis

For a given Ampere-turn level in the slots, the degrees of freedom offered by the stacked coils enable to consider several current distribution. The measured current waveforms are reported in Fig. 11 showing that the control of the current distribution requires small duty-cycle differences inducing negligible current ripple modifications. The corresponding total losses seen form the DC source are presented in Fig. 12. Similarly to the LF losses (Fig. 7), the minimum of total losses including the PWM-induced losses in the coils and the inverter losses shifts toward positive values when the fundamental frequency f e increases. It exists a non zero differential currents that minimizes the total losses in the two stacked coils. These minima are reported in Fig. 12b. Thanks to the high coupling between the two stacked coils, small duty-cycle differences (< 5 %) are required to ensure any LF current imbalance. Therefore, the term P HF presents a negligible variation over the different values of current imbalance. The shift of the optimal current imbalance toward positive values is mainly due to LF phenomena represented by P Copper,supp and P Core .

Conclusion

The stator winding segmentation studied in this article enables to increase the modularity of the segmented drive and to reduce the voltage rating of the elementary converter. However, the winding configuration obtained presents several coils stacked in the same slots, having then different relative position to the air-gap and strong magnetic couplings. The present papers describes the phenomena involved in these segmented supply and proposes an estimation of the losses depending on the current distribution enforced in the stacked coils. The example of two stacked coils wound around the same magnetic circuit is used to describe the modelling process and an experimental setup provides the necessary information to build the PWM-induced losses model. Finally, non zero differential currents references are investigated in order to mitigate the total losses in the magnetic circuit and conductors. This study is reinforced by taking into account the losses in the several H-bridge converters of the experimental setup. It fully justifies the use of differential currents to properly supply the stacked coils segmented structure at high fundamental frequencies.

Fig. 1 :

 1 Fig. 1: Cross section of the proposed winding segmentation within a generic magnetic circuit.

Fig. 2 :

 2 Fig. 2: Segmented supply of the stacked coils.

Fig. 3 :

 3 Fig. 3: a) Geometry of the practical model with two stacked coils and b) the corresponding practical model.

Fig. 4 :

 4 Fig. 4: a) Analytical model of a generic conductor among a slot and b) comparison of analytical and FEA based inductance value.

Fig. 5 :

 5 Fig. 5: LF electrical equivalent circuit.

Fig. 6 :

 6 Fig. 6: Voltages and currents measured when only a) coil 1 or b) coil 2 is supplied. 6 A peak current at f e = 200 Hz.

Fig. 7 :Fig. 8 :4

 78 Fig. 7: a) Measured and estimated LF losses versus current imbalance for several fundamental frequency values f e and b) Current imbalance minimizing LF Losses.

TimeFig.

  Fig. Measured (colored) and estimated (black) PWMinduced current ripple for ∆ α = 4 %.

Fig. 10 :

 10 Fig. 10: PWM-induced current modelling of the two stacked coils. The values are calibrated through experimental measurements at fixed duty-cycle.

  Fig. 11: Measured currents in the stacked coils with LF peak total current : 6 A while maintaining a) equal duty-cycles or b) equal fundamental currents (∆ α < 2.5 %).

  Fig. 12: a) Measured and estimated total losses P tot versus current imbalance for several fundamental frequency values f e and b) Current imbalance minimizing the overall losses.

Table 1 :

 1 Geometric parameters of the magnetic core and conductors

	Parameter	Symbol	Value
	Air-gap thickness	t ag (mm)	0.5
	Laminated steel thickness	e t (mm)	0.5
	Number of laminated steel sheet	n t	35
	Magnetic section	S m (mm²)	350
	Average magnetic lenght	l m (mm)	130
	Slot depth	d slot (mm)	30
	Slot width	w slot (mm)	10
	Vaccum permeability	µ 0 (µH/m)	4π/10
	Copper conductivity	σ (S/m)	59.6 • 10 6
	Conductor section	S c (mm²)	1.3
	Conductor length	l c (m)	2.5

Table 2 :

 2 Parameters of the electrical equivalent circuit

	Parameter	Analytical value Experimental value
	r 1 = r 2	30 mΩ	40 mΩ
	L 1 ≈ L m	180 µH	200 µH
	l 2 = L 2 -L m	40 µH	50 µH


François Legeais
email: francois.legeais@univ-rennes1.fr

Roger Lewandowski
email: roger.lewandowski@univ-rennes1.fr

Continuous boundary condition at the interface for two coupled fluids

Keywords: MCS Classification: 76D07, 35J20, 65N30 Stokes equations, coupled problems, variational formulation, numerical simulations

We consider two laminar incompressible flows coupled by the continuous law at a fixed interface Γ I . We approach the system by one that satisfies a friction Navier law at Γ I , and we show that when the friction coefficient goes to ∞, the solutions converges to a solution of the initial system. We then write a numerical Schwarz-like coupling algorithm and run 2D-simulations, that yields same convergence result.

Introduction

We consider the two coupled fluids problem with a rigid lid assumption, given by two 3D stokes equations, (1.3) for i = 1, 2, where the velocities (u 1 , u 2 ) = (u 1 (x h , z 1 ), u 2 (x h , z 2 )) are decomposed as u i = (u i,h , w), u i,h = (u i,x , u i,y ). Moreover, x h ∈ T 2 , where T 2 is the two dimensional torus, which means that we consider horizontal periodic boundary conditions. The interface Γ I is given by Γ I = {(x h , 0), x h ∈ T 2 }, the boundaries Γ i are given by Γ 1 = {(x h , z + 1 ),

-ν i ∆u i + ∇p i = f i , ∇ • u i = 0, (1.1) u 1,h|Γ I = u 2,h|Γ I , w i|Γ I = u i • n i|Γ I = 0, (1.2) u 1|Γ 1 = u 2|Γ 2 = 0,
x h ∈ T 2 }, Γ 2 = {(x h , z - 2 ), x h ∈ T 2 }, z 1 ∈ J 1 = [0, z + 1 ], z 2 ∈ J 2 = [z - 2 , 0]
, where z + 1 > 0 and z - 2 < 0. The coefficient ν i > 0 is the viscosity of the fluid i, p i its pressure. The main characteristic of this problem is the continuity boundary condition (1.2), which is natural and physical [START_REF] Batchelor | An introduction to fluid dynamics[END_REF], and usually considered for free interfaces [START_REF] Lannes | A stability criterion for two-fluid interfaces and applications[END_REF]. Notice that the rigid lid assumption we consider is reasonable for laminar coupled flows, as well as for large scales. In this paper we adress the question of the existence and uniqueness of a weak solution to Problem (1.1)-(1.2)-(1.3), given as the limit of "frictional solutions", for which we can write a numerical Schwarz-like algorithm. More specifically, we approach this problem by the following problem

-ν i ∆u i + ∇p i = f i , ∇ • u i = 0, (1.4) ν i ∂u i,h ∂n i | Γ I = -α(u i,h -u j,h ), w i|Γ I = u i • n i|Γ I = 0, (1.5) u i|Γ i = 0, (1.6)
i, j = 1, 2, in which the continuity condition (1.2) is replaced by the Navier law (1.5) where i = j. Similar problems have been already studied before, see [START_REF] Bernardi | A model for two coupled turbulent fluids. I. Analysis of the system[END_REF][START_REF] Jeffrey | Partitioned time discretization for atmosphere-ocean interaction[END_REF][START_REF] Connors | Partitioned time stepping for a parabolic two domain problem[END_REF][START_REF] Lions | Mathematical theory for the coupled atmosphere-ocean models[END_REF], and the existence and uniqueness of a weak solution is guaranteed. We aim to investigate how Problem (1.4)-(1.5)-(1.6) approaches Problem (1.1)-(1.2)-(1.3) when the friction coefficient α goes to infinity. Such question has already been adressed in [START_REF] Tapia | Stokes and Navier-Stokes equations with Navier boundary conditions[END_REF] for a single fluid, where it is proved that the corresponding solution strongly converges to a solution to the corresponding Stokes (Navier-Stokes) equations with a no slip boundary condition when α → ∞. We show in this paper the convergence in H 1 space type of the solution of (1.4)-(1.5)-(1.6) to a solution of (1.1)-(1.2)-(1.3) (see Theorem 2.1). As we shall see, numerical simulations are easily carried out by (1.4)-(1.5)-(1.6) thanks to a Schwarz-like coupling algorithm, that does not work for (1.1)-(1.2)- (1.3). This method has already been successfully implemented for coupled problems, see for example [START_REF] Tayachi | Design and analysis of a Schwarz coupling method for a dimensionally heterogeneous problem[END_REF]. The note is organized as follows. In the first part we set the functional framework and then we prove the convergence result, namely Theorem 2.1. In the second part, we describe our algorithm and show some numerical results in the 2D case. In particular we check the numerical convergence of the algorithm.

Convergence analysis 2.1 Energy balance

This section is devoted to the derivation of the main a priori estimate, which is standard. Let (u 1 , u 2 ) be any enough smooth solution to Problem (1.4)-(1.5)- (1.6). Taking the scalar product of equation (1.4) i by u i in integrating over T 2 × J i over yields by (1.6) i , because of the periodic boundary conditions in the x -y axes, the incompressibility condition and

u i • n i = 0 at Γ I , ν i T 2 ×J i |∇u i | 2 -ν i Γ I ∂u i,h ∂n i = T 2 ×J i f i • u i , giving by (1.5), ν i T 2 ×J i |∇u i | 2 + α Γ I u i,h • (u i,h -u j,h ) = T 2 ×J i f i • u i .
Summing up the two equalities yields the following energy balance, (2.1)

ν 1 T 2 ×J 1 |∇u 1 | 2 + ν 2 T 2 ×J 2 |∇u 2 | 2 + α Γ I |u 1,h -u 2,h | 2 = T 2 ×J 1 f 1 • u 1 + T 2 ×J 2 f 2 • u 2 .

Functions spaces, variational formulation

Let W i = {u ∈ C ∞ (T 2 × J i ), u |Γ i = 0, u • n i | Γ I = 0, ∇ • u i = 0}, equipped with ||u|| i,1 = ||∇u|| L 2 (T 2 ×J i )
which is indeed a norm due to the condition u |Γ i = 0. Let W i denotes the completion of W i with respect to this norm,

(2.2) W = W 1 × W 2 , W 0 = {(u 1 , u 2 ) ∈ W, u 1,h|Γ I = u 2,h|Γ I a.e. in Γ I }.
We equip W with the scalar product, for any

U = (u 1 , u 2 ), V = (v 1 , v 2 ) ∈ W , (2.3) Λ(U, V) = ν 1 T 2 ×J 1 ∇u 1 • ∇v 1 + ν 2 T 2 ×J 2 ∇u 2 • ∇v 2 .
The space W 0 is the kernel of the form L : (u 1 , u 2 ) → u 1,h|Γ I -u 2,h|Γ I , which is continuous by the trace theorem. Therefore W 0 is a closed hyperplane of W . Let P denotes the orthogonal projection over W 0 , and Φ = (φ 1 , φ 2 ) a unit orthogonal vector to W 0 , so that

W ⊥ 0 = vect Φ. Definition 2.1. (weak solution) A couple U = (u 1 , u 2 ) ∈ W is a weak solution to Problem (1.4)-(1.5)-(1.6) when ∀ V = (v 1 , v 2 ) ∈ W , (2.4) Λ(U, V) + α Γ I (u 1,h -u 2,h ) • (v 1,h -v 2,h ) = T 2 ×J 1 f 1 • v 1 + T 2 ×J 2 f 2 • v 2 = (F, V).
Throughout the rest of the paper, we assume that

f i ∈ L 2 (T 2 × J i ), i = 1,

Convergence

Let U α = (u α 1 , u α 2 ) ∈ W be the solution of (1.4)-(1.5)-(1.6). We study in this section the convergence of the familly (U α ) α>0 when α → ∞, proving the following result.

(2.5) ∀ V = (v 1 , v 2 ) ∈ W 0 , Λ(U, V) = (F, V).
Moreover, the solution of (2.5) is unique.

Proof. Let U α = (u α 1 , u α 2 ) ∈ W = W 1 × W 2 be
the solution of (S 1 , S 2 ). We first show that the familly (U α ) α>0 is bounded in W . We have, by (2.1), (2.6)

U α 2 W + α Γ I |u α 1,h -u α 2,h | 2 = T 2 ×J 1 f 1 • u α 1 + T 2 ×J 2 f 2 • u α 2 ,
which yields

(2.7) U α 2 W ≤ T 2 ×J 1 f 1 • u α 1 + T 2 ×J 2 f 2 • u α 2 .
We deduce from Poincaré and Cauchy-Schwarz inequalities that (U α ) α>0 is indeed bounded in W . Therefore, we can extract a subsequence (U αn ) n∈I N (α n → ∞ as n → ∞) which converges weakly in W to some U ∈ W . Moreover, by the trace theorem and usual Sobolev compactness results, the corresponding traces are strongly convergent in L 2 (Γ I ). As by (2.6) lim n→∞ tr(u αn 1,h -u αn 2,h ) = 0 in L 2 (Γ I ), then U ∈ W 0 . Finally, take V ∈ W 0 in (2.4) as test, so that the boundary term vanishes. By passing to the limit in this case when α → ∞, we obtain that U is a weak solution to (1.1)-(1.2)-(1.3). Uniqueness is straightforward, which in addition garanties that the entire familly does converge to U.

It remains to show the strong convergence. Let λ α ∈ IR, be such that U α = P U α + λ α Φ (Φ being given in section 2.2). We first show the strong convergence of (P U α ) α>0 to U by taking P U α as a test in (2.4) which gives, by using the orthogonal decomposition of

U α , (2.8) Λ(U α , P U α ) = P U α 2 W = (F, P U α ),
since the boundary term on Γ I equals to zero by orthogonality. Therefore (P U α ) α>0 is bounded in W 0 , then converges weakly -up to a subsequence (keeping the same notation)to a limit W, strongly in L 2 ((T 2 × J 1 ) × (T 2 × J 2 )). Taking V ∈ W 0 in (2.4) as test, noting that in this case Λ(U α , V) = Λ(P U α , V), and passing to the limit when α → ∞, we see that W is solution of the problem (1.1)-(1.2)-(1.3), hence W = U by uniqueness, and the entire sequence converges. Therefore, passing to the limit in (2.8) yields

lim α→∞ ||P U α || 2 W = (F, U) = Λ(U, U) = ||U|| 2 W ,
which, together with the weak convergence, ensures the strong convergence as claimed. To conclude, it remains to prove that lim α→∞ λ α = 0. By the energy balance (2.6), we have

(2.9) lim sup α→∞ ||U α || 2 W + lim sup α→∞ α Γ I |u α 1,h -u α 2,h | 2 = (F, U).
However, we have

Γ I |u α 1,h -u α 2,h | 2 = λ 2 α Γ I |φ 1,h -φ 1,h | 2 . Therefore, since (||U α || 2 W
) α>0 is bounded, by (2.9) we have (αλ 2 α ) α>0 is bounded, which can happen only if λ α → 0 as α → ∞, concluding the proof.

3 Numerical simulations

Algorithm

We solve the problems (1.4)-(1.5)- (1.6) for large values of α, with a coupling Schwarz like algorithm, using the software Freefem++, for solving 2D Stokes problems by the finite element method. Our algorithm is set as follows.

Step 1: We solve the problem on the upper part which gives a first value u α,0

1 .

-ν 1 ∆u α,0 1 + ∇P

(0) 1 = f 1 , ∇ • u α,0 1 = 0, (3.1) ν 1 ∂u α,0 1,h ∂n 1 | Γ I = -αu α,0 1,h , (3.2) u α,0 1,h | Γ 1 = 0, u α,0 1,h • n 1 = 0. (3.3)
This velocity allows us to solve the problem on the lower part, and to calculate the velocities step by step, up and down.

Step 2: We calculate u α,n 2 and u α,n+1 1 by solving

-ν 2 ∆u α,n 2 + ∇P (n) 2 = f 2 , ∇ • u α,n 2 = 0, (3.4) ν 2 ∂u α,n 2,h ∂n 2 | Γ I = -α(u α,n 2,h -u α,n 1,h ), (3.5) u α,n 2,h | Γ 2 = 0, u α,n 2,h • n 2 = 0. (3.6) and -ν 1 ∆u α,n+1 1 + ∇P (n+1) 1 = f 1 , ∇ • u α,n+1 1 = 0, (3.7) ν 1 ∂u α,n+1 1,h ∂n 1 | Γ I = -α(u α,n+1 1,h -u α,n 2,h ), (3.8) u α,n+1 1,h | Γ 1 = 0, u α,n+1 1,h • n 1 = 0. (3.9)
Note that we are able to prove the stability of this algorithm, and numerical simulation confirm the convergence (see table below). Problem (1.1)-(1.2)-(1.3) cannot be solved in a similar way. Indeed, the interface conditions

(3.10) u (n) 2,h | Γ I = u (n) 1,h | Γ I , u (n+1

Simulation results

We take z To check the numerical convergence of the method, we study the error term U α,n+1 -U α,n L 2 . On the left, for a given α, we have the first value of n for which U α,n+1 -U α,n L 2 < 10 -3 . The method is always converging, and the convergence is almost instantaneous for large α (> 10 6 ).

Theorem 2 . 1 .

 21 The familly (U α ) α>0 strongly converges in W to a weak solution U = (u 1 , u 2 ) ∈ W 0 of Problem (1.1)-(1.2)-(1.3) when α → ∞, in the sense:

) 1 ,

 1 h | Γ I = u (n) 2,h | Γ I ,imply that the sequences (u(n) 1,h ) n and (u (n) 2,h) n are constant on the interface Γ I which doesn't allow any iterations on the coupling algorithm.

+ 1 =

 1 50 , z - 2 = -5, L = 100, ν 1 = ν 2 = 1 and the source F

Figure 1 :

 1 Figure 1: Domain and mesh

Figure 2 :

 2 Figure 2: L2 norm of the errors and rate of convergence for n=100

Figure 3 :

 3 Figure 3: Velocities, f 1 = (1, -1) = f 2 , n = 9, α = 10 9

  2. The existence and the uniqueness of a weak solution to Problem (1.4)-(1.5)-(1.6) that satisfies the energy balance (2.1) is straightforward by the Lax-Milgram Theorem for any given α > 0. Notice that work remains to be done about the pressures, by a suitable adaptation of a De Rham like theorem in this framework, which is an open problem.


Fabien Siron

Dumitru Potop-Butucaru

Robert De Simone

Damien Chabrol

Amira Methni

The synchronous Logical Execution Time paradigm

Real-Time industrial systems are not so much of those that have to perform tasks incredibly fast, but in a timepredictable manner; they rather focus on meeting previously specified timing requirements in a provable way. Consequently, time must be taken into account from the very start of the design. However, exact timing constants may not be available yet in early design stages as they may depend on the target. In answer, formalisms based on the Multiform Logical Time have been introduced to abstract real-time durations. The Synchronous-Reactive (SR) approach introduced a discretized abstraction of time on which computations happen logically instantaneously. Contrary to SR, Logical Execution Time (LET) mandates to specify the actual logical duration a task has to fulfill. This allows a more efficient compilation, at the price of a lower expressiveness. Classical LET (i.e. as introduced in Giotto/TDL) sticks to uniform pseudo-physical time, i.e. based on one logical clock mapped to the real-time. In this paper, we introduce a new paradigm called synchronous Logical Execution Time (sLET) that builds upon both SR and LET paradigms. It keeps the idea of logical durations coming from the LET paradigm, while having logical instants based on logical clocks. This extends the expressivity of LET, as time is totally abstracted as sequences of events. The various schedulings provide physically timed versions that, while having distinct non-functional properties (in terms of performance mostly), remain mutually functionally equivalent (in the logical time realm). A particular instance, where computations are executed "in a single instant", and then time is advanced (as in classical event-driven simulation), can lead to a direct translation into synchronous formalisms (in our case Esterel). We started inquiring how this could open new ways of verification and analysis on PsyC programs.

I. INTRODUCTION

The design of embedded control software calls for stringent real-time constraints due to their permanent interaction with the physical environment. Such real-time systems are usually safety-critical because of the context in which they are used (avionic, automotive . . . ). Therefore, their designs should be verified with the highest possible level of confidence. Various formalisms based on the concept of Multiform Logical Time have been introduced to abstract real-time durations which are usually not known at the design phase, as they might be targetdependent. Then, specific analysis, called schedulability analysis (or time safety analysis) ensure that physical computations satisfy their logical time constraints. Among those formalisms, the Logical Execution Time (LET) paradigm [START_REF] Kirsch | The logical execution time paradigm[END_REF] brings a compromise between the strong ex-pressiveness of the synchronous-reactive approach (SR) [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] and the efficiency of traditional task scheduling. Contrary to the SR model, computation takes time and is bounded by a fixed logical duration usually known before the computation happens. However, in SR model, because computations are considered instantaneous, time can be refined using multiple logical clocks. LET is usually based on a unique clock representing the progress of time. Consider the following temporal interval where the bounds are described by sync: sync 1 ms; f(); sync 1 s;

Depending on the semantics of the language, the interval might last one second or up to the next second tick. The former semantics, usually used by LET based languages, would mean that both ms and s refer to the same clock (i.e. 1 s being just a short for 1000 ms) while the latter semantics, corresponding to synchronous languages, mean that both ms and s are logical clocks related with an affine relation.

In industrial projects, while classical synchronous languages such as Lustre or Scade fit perfectly the need to describe the functional behavior of the system, Synchronous LET, which is implemented by the industrial language PsyC [START_REF] David | Safety properties ensured by the oasis model for safety critical real-time systems[END_REF], actually focus on a different level of the system life-cycle, namely, the software integration level. This approach is quite classical (see [START_REF] Caspi | From simulink to scade/lustre to tta: a layered approach for distributed embedded applications[END_REF], [START_REF] Bliudze | From model-based to real-time execution of safety-critical applications: Coupling scade with oasis[END_REF] and [START_REF] Forget | A Multi-Periodic Synchronous Data-Flow Language[END_REF]) and divides the functional design of the system from its non-functional design, that is, how different functional components (either designed with synchronous languages or manually) are integrated. This multilayer approach is depicted in figure 1. This is, however, still a challenge today to verify that properties specified at system-level (i.e. in the specification) still holds at integration level.

System-Level Specification (e.g. SysML)

Functional Design (e.g. Scade)

Integration level (e.g. PsyC) Fig. 1: Typical Software Life-Cycle of industrial real-time systems

In this paper, we introduce a new paradigm called synchronous Logical Execution Time (sLET) that builds upon both the SR model and LET. It keeps the idea of logical durations coming from the LET paradigm while having logical instants based on logical clocks. This extends the expressivity of LET as time is totally abstracted as sequences of events. Thus, it inherits 1) the simple compilation and fine-grained schedulability analysis coming from LET and 2) the synchronous semantics that is well-suited for formal verification.

After an overview of existing Multiform Logical Time models, we introduce sLET as an extension of LET in section II. To illustrate sLET, we then give an overview of the formal semantics of a subset of the industrial language PsyC [START_REF] David | Safety properties ensured by the oasis model for safety critical real-time systems[END_REF] in section III, developed by the company Krono-Safe, and then, we give translation rules to the synchronous language Esterel in section IV. While PsyC is compiled using techniques similar to the ones of LET, formal verification techniques coming from synchronous languages such as Esterel could then be re-used. Such approach is illustrated in section V where PsyC and Scade interact in a classical use-case showing a flight control system.

II. MULTIFORM LOGICAL TIME MODELS

Synchronous languages introduce the notion of Multiform Logical Time, through logical clocks, so they could better be called Polychronous time. A logical clock counts the successive ticks of any relevant event. Still, their operational semantics expands the behaviors in terms of a grandmother clock, the discrete reaction step measuring the instant. On the other hand, initial LET formalisms rely on such a totally ordered discrete time, and introduce quantitative durations measured in it; although they may use pseudo-physical unit names (milliseconds, nanoseconds,...) for it, the compliance with actual physical implementation is only a wish (or a requirement), that will have to be checked later when the latter becomes available. We now extend on these notions, and the abstraction they allow for fast and human-legible logical design, cleanly split from the later physical implementation (and without back and further adjustements hopefully, a main purpose of these formalisms). The purpose is to position the sLET formalism as combining multiform logical time and the ability of user-defined logical clocks, together with quantitative durations in which a certain amount of computation can be spread in any instants of a given interval, providing it does not overflow it.

A. The Synchronous-Reactive Model (SR)

The Synchronous-Reactive model, implemented by synchronous languages [START_REF] Benveniste | The synchronous languages 12 years later[END_REF], totally abstracts execution time to focus on logical instants, allowing both determinism and concurrency. For that, computations react simultaneously and instantaneously with the tick of a global common logical clock as shown in figure 2a, which can be further refined to give multiple logical clocks. The synchronous hypothesis, then, ensures that if every computation is bounded by the next reaction, then physical time can be safely ignored as shown in figure 2b. Combining synchrony and concurrency allows a very expressive formal model while keeping it very simple [START_REF] Benveniste | The synchronous languages 12 years later[END_REF]. Hence, SR model is well adapted for formal verification. SR model has been introduced in languages such as Esterel [START_REF] Berry | The Constructive Semantics of Pure Esterel[END_REF] which is mainly imperative and Lustre [START_REF] Pilaud | Lustre: A declarative language for programming synchronous systems[END_REF] which is declarative. The compilation of synchronous languages is, however, quite complex. First, instantaneous communication makes the compilation for parallel platforms (i.e. multicore, distributed) quite hard and can produce causality issues. Secondly, execution time is usually limited to the reaction time of the system. Thus, applications with non negligeable execution time (i.e. also called the long period problem) can be rejected during the schedulability analysis. Nonetheless, languages based on the synchronous model such as Prelude [START_REF] Forget | A Multi-Periodic Synchronous Data-Flow Language[END_REF], can adress this problem, but are usually limited to specific patterns such as a multi-periodic synchronous model.

B. The Logical Execution Time Model (LET)

The Logical Execution Time paradigm abstracts physical time through a sequence of instants, similarly to the synchronous approach. However, contrary to the latter, computation takes time and is not considered instantaneous. For that, each computation must fit in a logical interval, called LET interval [START_REF] Kirsch | The logical execution time paradigm[END_REF]. Furthermore, communications are only made on the boundaries of LET intervals. Inputs are read at their beginning and outputs are made visible to other tasks at their termination (see Figure 2c). As communication is only done on predefined instants, computations behave like they always take the same time, which is the LET interval duration (see Figure 2d). As such, the determinism property is ensured.

The LET paradigm has been initially introduced in the Giotto [START_REF] Kirsch | The logical execution time paradigm[END_REF] language in which each task correspond to a periodic LET interval. The language can also express global modes that allows to change from a given set of task frequency to an other. Timing Definition Language (TDL) [START_REF] Pree | Modeling with the timing definition language (tdl)[END_REF] extends Giotto via a decomposition in concurrent modules. Each of them defines a set of task and their frequency, similarly to Giotto, but a module can also change its local mode independently of other modules. Giotto has also been extended with events in the xGiotto language [START_REF] Ghosal | Event-driven programming with logical execution times[END_REF] through the mechanism of event scoping. In all these approaches, the bounds of LET intervals are either defined using chronometric pseudo-physical durations or events (i.e. in xGiotto). This paper proposes a more abstract and homogeneous model of time, that is, the use of multiple logical clocks.

In a way, LET extends logical time with the concept of logical durations. This allows a precise schedulability analysis due to an increased execution time variability. In turn, this also makes compilation to parallel platforms easy and causality issues of the SR model are avoided. However, as LET forbids

C. sLET as a synchronous extension of LET

Synchronous Logical Execution Time (sLET) extends the classical LET paradigm with the concept of logical clocks. As such, interval bounds can be triggered by clock ticks while in classical LET, interval duration is usually defined by a constant fixed duration. However, as in LET, communication can only happen on interval bounds, which ensure temporal determinism. Being closer to the SR model, sLET is actually fully compatible with the synchronous hypothesis, that is, computations can happen instantaneously on the activation date of the interval as long as outputs are delayed to the end of the interval. This model is actually a generalization of the Psy model introduced in the OASIS framework since the 90's [START_REF] David | Safety properties ensured by the oasis model for safety critical real-time systems[END_REF].

A logical clock, in the sense of Lamport [START_REF] Lamport | Time, clocks, and the ordering of events in a distributed system[END_REF], abstracts time through a series of events called clock ticks. The sLET paradigm adds to LET a finite set of logical clocks C on which LET interval boundaries can be based. We will call clock event an event of the form n × c where n ∈ N * and c ∈ C. This event basically means that we have to wait n ticks of the clock c. An interval termination is then defined with a clock event which may also define the activation of a next interval. As an example, if an interval termination is defined with n t ×c t , then after its activation, it should wait for n t ticks of the clock c t . Note that if all intervals use the same clock, then clock events actually express fixed constant durations with respect to that clock which is similar to the classical LET model. sLET keeps with the LET idea of imposing at an early phase of logico-functional design some fixed interval durations at which boundaries the I/O external behaviors will be supposed to occur. Then there is flexibility on when exactly computations are scheduled, as long as they remain inside these bounds. Next, we will introduce the PsyC languages foundations, bringing syntax to the sLET view. In particular, PsyC allows conditional durations, in alternative "if-thenelse" behaviors. Note however that classical LET, as defined in Giotto, can actually be expressed easily with sLET, that is, using only one logical clock mapped to real-time. The converse is not true. Indeed, even considering that all logical clocks are derived using affine relation from a unique global clock, sLET interval might have different duration depending on the current global state.

III. THE PSYC LANGUAGE: ABSTRACT SYNTAX AND FORMAL SEMANTICS

A. Language description

In this section, we give a brief description of PsyC, a language developed by the company Krono-Safe, as an illustration of the sLET paradigm. This language is implemented in the ASTERIOS product which provides a set of tools to design safety-critical real-time software. Such applications can then be certified at the highest level of criticity for the avionic domain (DAL-A, DO-178C). ASTERIOS and PsyC are inherited from the OASIS and PharOS projects coming from the CEA. Initially, PsyC (for Parallel SYnchronous), in OASIS, was presented as a timed-triggered approach for safety-critical real-time software under a model called Psy [START_REF] David | Safety properties ensured by the oasis model for safety critical real-time systems[END_REF]. This model is actually very close to LET, and synchronous LET generalize both of them. The modern version of PsyC is actually better characterized with synchronous LET due to the existence of multiple logical clocks.

A PsyC application is composed of a fixed set of tasks called agents that are formed of sLET intervals. The content of agents is composed of C code that is extended with special instructions like advance n with c which specifies the boundaries of sLET interval. Informally, its semantics is to advance the logical time of n ticks of a clock c. Such clocks have to be defined in the application as an affine relation with respect to another clock. They are actually two kinds of logical clocks in PsyC, sources, which are defined externally, and clocks, which basically refine sources. In practice, most of the applications define only one source called realtime which is mapped on real time. The Figure 3 shows an example of an agent using multiple clocks and a conditional computation in the PsyC language. A possible timeline is shown in figure 4. Such pattern with multiple clocks allowed by sLET couldn't be expressed, as is, in the classical LET model.

Fig. 3: Example of a PsyC agent

Inter-agent communication is performed through a dedicated channel called Temporal Variable. Multiple agents can read a temporal variable but only one can write into it. According to sLET, inputs and outputs in a sLET interval can only be made on sLET interval bounds. Hence, data emitted by an agent is made visible (i.e. the temporal variable is updated) on the instant described by the next advance instruction. Futhermore, Temporal Variables also have another layer of sampling (defined with a clock) and values are read with the expression $[n]var which denotes the n th last sampled value.

The Figure 5 describes the abstract syntax of a subset of PsyC which is necessary to introduce the formal semantics described in the next section. Only specific PsyC constructions are described and the syntax of C expressions is omitted for simplicity. Also, we add a specific statement skip that does nothing (i.e. in the C syntax, this could represented as a semi- Where n, n 1 ∈ N * , n 2 ∈ N, v is the initial value of the temporal variable and c i are clock identifiers.

Fig. 5: Abstract syntax of our PsyC subset colon). This will be helpful in the semantics to express that a statement doesn't have any work left to do. Futhermore, both the starttime and the source parameter of the agent are omitted because the former could be defined as an advance and the latter could be deduced from the clocks. However, we assume that an agent could only use clocks that are derived from a unique source. Nonetheless, different agents might use different sources.

B. Formal Semantics

In this section, we introduce the formal semantics of a subset of the PsyC agents thanks to term rewriting. Basically, a term can be rewritten sucessively, with respect to a given set of formal rules, to form a symbolic execution that serves as an oracle. The rules are given using the Structural Operational Semantics (S.O.S.) approach introduced by Plotkin [START_REF] Plotkin | A structural approach to operational semantics[END_REF] and later adapted to synchronous languages [START_REF] Berry | The Constructive Semantics of Pure Esterel[END_REF]. The idea behind S.O.S. is to give the rewriting relation of a term with respect to the rewriting of its parts. This is usually defined through a set of inference rules. The behavioral semantics of Esterel extends S.O.S. through the use of an integer encoding the type of the transition, i.e. whether the transition takes time or is instantaneous [START_REF] Berry | The Constructive Semantics of Pure Esterel[END_REF] (i.e. fits entirely inside an instant reaction).

As described above, the semantics of PsyC is described in this paper using two relations, one that makes time progress, and one that is not explicitely timed, and should be considered as conceptually guaranteed to fit in the time interval closed by the next time-progress behavior. This vision can be considered as borrowed from discrete-event simulation models. While the semantics allows many later schedules for the internal behaviors before next time advance (a property exploited by the PsyC compiler to optimize their time allocations according to other criteria), one may also consider the (valid) interpretation where all computations are made in the initial instant, than will allow the translation to synchronous languages of the next section.

Consider that the configuration of an agent (i.e. its state) can be described as a pair E, T where E is the assignements of its local variables and T the assignments of its displayed temporal variables (i.e. its outputs). Based on this, we define an instantaneous transition as following:

E, T s ---→< E , T > s
This transition denotes that a statement s can be rewritten in s instantaneously while the configuration < E, T > is updated to < E , T > with respect to s. Based on this definition, we can define the semantics of the first basic statement, the assignation:

E, T x := exp ---→ E , T skip (assign)
where E = U pdate(E, v, exp).

The assign rule evaluates its expression and updates its local variables accordingly using a function that we call U pdate(). This statement is logically instantaneous because it is executed at the beginning of a sLET interval to read the correct input values. Also, it is reduced to skip as there is no work left to do. The advance rule is a little bit more complicated as this statement makes time progress. Recall however that clocks used in an agent by advance statements are all derivated from a unique source. Thus, we first instantaneously rewrite the advance statement with respect to its parameters (i.e. number of ticks and clock) to a counter of source ticks that we define with the pseudo statement Remains. This pseudo statement takes an absolute duration in parameter (with respect to the agent source) that is computed with the function Duration with respect to the advance parameters, and the current source date.

E, T

advance n with c ---→ E, T Remains(N ) (advance) where N = Duration(n, c, date)

The Remains pseudo statement defined above can then be used to make time progress. For that, we define a new transition that takes one time unit of the agent source (symbolized by the double arrow shape):

E, T s = === ⇒< E , T > s
Remains is then defined with two different rules. The first one, Remains-1, decreases the counter when its greater than 1. When the counter is equal to 1, Remains is rewritten to skip as there is no work left to do and temporal variables are updated. As other agents can only read the latter (and not the agent local environment), this actually means that outputs can only be made visible at the end of a sLET interval, which is actually the semantics of the sLET paradigm.

E, T

Remains

(N ) = === ⇒ E, T Remains(N -1) (Remains-1) if N > 1 E, T Remains(1) = === ⇒ E, T skip (Remains-2)
where T = U pdateOutputs(T, E).

The rules defined above can be composed quite easily. First, when there are executed in sequence, multiple instantaneous transitions can be represented with one big instantaneous transition and when they are followed by a non instantaneous transition, they can be represented with a big non instantaneous transition. In other words, the semantics of an agent can be represented with only big non instantaneous temporal transitions, which is interesting to show the expected temporal behavior of the agent. Second, the rules defined above only show simple statements, but actually, this can be extended very easily to control flow statement (e.g. condition statement) given that the type of transition is propagated correctly.

To illustrate a little bit the semantics above, we show a very basic example of a sLET interval with an output variable called tv. Let's consider the following PsyC statements: tv := 2 ; advance 2 with c2 Let's assume that c2 is a strictly periodic clock with period 2 and offset 0, and the current date of the agent is odd. Then, with respect to the semantics rules, we have the following execution: tv =?, tv =? tv := 2 ; advance 2 with c2

---→ tv = 2, tv =? advance 2 with c2 ---→ tv = 2, tv =? Remains(3) = === ⇒ tv = 2, tv =? Remains(2) = === ⇒ tv = 2, tv =? Remains(1) = === ⇒ tv = 2, tv = 2 skip
This execution shows that the communication is correctly updated at the end of the interval and its duration is actually 3 as it depends on the current date (e.g. if the date was even, then the interval should have lasted 4 ticks). Also, one can notice that such pattern is not possible in classical LET languages (e.g. Giotto, TDL) as they usually can only represent duration on a single clock. IV. CONSEQUENCES

A. Synchronous semantics of sLET via Esterel translation

To show the direct relation between sLET and the SR model, we can encode the semantics of PsyC into a synchronous language like Esterel. We choose Esterel because its syntax is quite close to PsyC, that is, both languages are controlflow and imperative. Basically, the approach is to focus on the logical behavior of synchronous LET and to totally abstract actual executions, i.e. to consider them (logically) instantaneous. Thus, of course, this kind of simulation is not representative of the actual operational behavior where each computation takes time. Nonetheless, sLET (similarly to classical LET) allows multiple valid schedules that are mutually equivalent, that is, that form an equivalence class with respect to logical clock ticks. The synchronous semantics is then actually a particular instance where in a sLET interval: 1) all computations are executed during the first instant, thus inputs are read synchronously to the start of the interval; 2) time is advanced to the end of the interval; 3) outputs are produced synchronously to the end of the interval.

This abstraction should be sufficiant to express properties at the model level as actual computations should not impact temporal properties of a sLET design. While this approach focus on synchronous LET, as the latter is a generalization of the classical LET model, this approach could be adapted easily to other LET languages.

Esterel is an imperative synchronous language control-flow based. In particular, it allows imperative concurrency, that is, handling multiple instruction pointers at a same instant and preemption where a computation can be aborted when a given signal is present. The only communication means available in Esterel is through instantaneous signals which can have a status, present or absent, and a value. A signal which only has a status is said to be pure, otherwise, it is said to be valued. The instructions are basically divided into two categories similarly to PsyC:

• the instantaneous instructions, like emit s which emits a given signal s; • and temporal instructions, like pause or await s which respectively wait for the next instant and for the next instant where signal s is present.

For a full definition of the Esterel language, we invite the reader to consult [START_REF] Berry | The Constructive Semantics of Pure Esterel[END_REF].

In a PsyC application, different elements are actually composed in a synchronous parallel fashion. Assuming that PsyC clock ticks are represented as Esterel pure signals (i.e. either the clock ticks or not), we have to represent clock generation, temporal variable sampling and agent behavior. As explained at the beginning of section III, a PsyC source clock is actually an external logical clock, and is thus represented as an Esterel pure signal input. Clock generation basically refines a clock tick, i.e. a pure signal, based on an affine relation. A temporal variable is represented as a persistent valued signal (i.e. which is always present) and is sampled on the signal corresponding to its clock. And finally, an agent is almost translated as is. The control-flow is translated in its equivalent form in Esterel and basic body patterns can be translated either to sequence or to infinite loop depending on the next body to be executed. We do not cover advanced body patterns in this paper for simplicity but this could be adapted with some more advanced Esterel control-flow patterns. All variables of an agent, either local to the agent or temporal variables, are translated to Esterel Fig. 7: The Flight Control System use-case variables. The advance statement is translated to the await statement followed immediately by an update of the temporal variables, that is, the valued signal of each temporal variable in output of the agent is updated with its local variable value. This is actually the main difference from the synchronous model, communication cannot be updated instantaneously. All the translation patterns described above are sketched in the table 6.

B. Synchronous Observers and Formal Verification

The main interesting consequence of the Esterel translation is to apply techniques and tools coming from the synchronous community. Among them, formal verification has been successfully applied on the synchronous approach due to two main things: 1) the synchronous composition semantics, the state space is usually quite reduced comparing to asynchronous approaches and 2) symbolic model-checking allows to scale verification techniques even more, using BDD or SMT techniques. Futhermore, in synchronous languages, properties that are usually modeled using temporal logics, can also be modeled directly in the language using the socalled synchronous observer approach [START_REF] Halbwachs | Validation of synchronous reactive systems: From formal verification to automatic testing[END_REF]. Observers only observe the state of the application and defines which state is acceptable or not. Classicaly, a signal Error is emitted if an error is detected and verification is reduced to a reachability analysis on this signal.

In PsyC, the typical kind of property that we might want to check applies on the different temporal variables rhythms. Typically, due to the sampled buffered semantics of temporal variable, if the producer and the consumer are not at the same rhythm, produced values can be lost or consulted values might be duplicated. While the second scenario is usually not a problem, the first one might be one. Of course, this kind of property is quite simple, but this opens the possibility to more evolved properties like freshness constraints of temporal variables or end-to-end latencies.

V. USE-CASE

A. Description

As highlighted in the introduction, languages based on LET, such as PsyC, focus on the architecture description at integration level of real-time systems. In other words, they allow to express how multiple functional components connect and what is their real-time behavior. In particular, it is possible to express multi-rate behavior. Typically, in industrial applications, a synchronous language such as Scade, is used to design the functional components of the system and a language like PsyC can be used to specify the real-time software integration of these components.

To reflect this approach, the following use-case is based on an example found in the Scade Suite software. Consider a basic Flight Control Software in which you expect the following behavior:

• given an altitude and a speed command, a flight controller implements control laws for altitude and speed with respect to altitude and speed sensors; • depending on the flight mode, i.e. AutoPilot or Manual, the commands given to the flight controller is either the input commands or the commands of the previous instant; • regulation alarm is raised if the altitude sensor (respectively the speed sensor) is too far from the altitude command (respectively the speed command).

The system architecture is shown in the figure 7. Basically, each component correspond to a function detailed above. Each functional component can either be implemented manually or with a Scade sheet. To illustrate a little bit more this approach we show in figure 8 the Scade sheet corresponding to the mode handling. The mode handling is modeled with a finite state machine with only two states, the AutoPilot state and the Manual state corresponding to the AutoPilot mode and Manual mode described at the beginning of the section. The control laws are based on PID controllers but we will not dive into the details as it is beyond the scope of this paper. For simulation purpose, the example also contains a plane simulation that takes throttle and elevation commands (from the controller), simulates the actuator response as well as the aerodynamics behavior and yield updated speed and altitude sensors values (i.e. to the controller). Running such a simulation in Scade typically ignores real-time constraints and every components run at the same rate. While this is totally justified for simulation, it's is not necessarily the case for compilation. Indeed, in practice, different components might have different real-time behaviors, depending on available resource, as well as different mapping constraints (typically for multicore target). Thus, in the next section we will show how PsyC can be used to integrate seamlessly multi-rate functional components.

B. Architecture Description in PsyC

We choose a very simple mapping where each components described above are directly mapped to an individual agent. Futhermore, each communication signal is mapped to an individual temporal variable. Of course, while only one agent can write a temporal variable, multiple agents can read it. Thus, a temporal variable is not duplicated for each reader component. The temporal behavior of a component in an agent is the following:

In our example, let's consider that we want to ensure that all commands values are taken into account once by the Alarm agent. This functional chain is thus composed of the FlightMode agent, the SpeedTarget and the AltitudeTarget temporal variables, and the Alarm agent. Considering that the corresponding tick events have the same name, we have the following constraints to check to ensure that no values are lost:

FlightMode Write ∼ SpeedTarget ∧ FlightMode Write ∼ AltitudeTarget ∧ SpeedTarget ∼ Alarm Read ∧ AltitudeTarget ∼ Alarm Read
Based on [START_REF] André | Verification of clock constraints: CCSL Observers in Esterel[END_REF], an Esterel observer for Alternate can be derivated. To fit our purpose, alternation strictness (whether two synchronous tick events could be synchronous or not) should be adapted to our communication mode. That is, the transition from event A to B in A ∼ B can be instantaneous but not the transition from B A. The verification of the constraints described above can be written in instanciating the Alternate observer for each of the constraint with the help of parallel composition: In our scenario, this property is verified only when the AutoPilot mode is set as the Alarm agent might read multiple times the same values of the FlightMode agent when AutoPilot is disabled. The observer could be adapted easily to enforce alternation only when the mode is set. Also, note that this property does not enforce that all clocks should be synchronous (i.e. at the same rate). In fact, the agent Alarm might have some offset with FlightMode still without losing any value.

VI. CONCLUSION

This article presents a new paradigm called sLET that builds upon both the synchronous and the LET paradigm. Such unification allows to re-use formal analysis techniques coming from synchronous languages while keeping the simple compilation and efficient schedulability analysis of LET. sLET (as the classical LET paradigm) allows multiple valid schedulings that are mutually equivalent (in the logical time realm), as long as all of the computations fits in their sLET intervals (the PsyC compiler actually using one of these solutions also satisfying actual physical computation durations). The synchronous semantics is then actually a particular instance where computations are executed "in a single instant" and outputs are delayed. Consequently, this shows that one can use analysis based on the synchronous approach to reason on languages based on LET (e.g. Giotto, xGiotto, TDL) at the logical level; in particular, it unveils the possibility to re-use synchronous formal verification techniques for these languages. The sLET paradigm is implemented by the PsyC language which is part of the ASTERIOS framework, produced by KRONO-SAFE. This work opens twofold perspectives: firstly, express more advanced properties, such as end-to-end latencies or data freshness, in a simpler way; secondly, synchronous verification techniques should be adapted more thoroughly to be efficient. In particular, (s)LET intervals are composed of instants that only make time progress and thus, can be optimized for verification. Finally, while our main focus is formal verification, other well known techniques can be adapted to languages based on (s)LET, such as model-in-the loop simulation, automatic test generation or test coverage.

Fig. 2 :

 2 Fig. 2: Multiform Logical Time models: physical vs logical views

Fig. 4 :

 4 Fig. 4: Possible timeline of example in fig 3. A Red frame represent computation A while a blue frame represent computation B.

Fig. 8 :

 8 Fig. 8: Mode Controller of FCS

• The reset function of the component generated from Scade is executed once with a given phase constraint (not specified here) and; • The cycle function of the component generated from Scade is executed indefinitely after the reset function, on a given rate.

As an illustration, the code in 9 describes the temporal behavior of the controller component. The inputs and outputs of the cycle function are transfered with the help of global variables which is the usual calling convention with code generated from Scade. Chosing the right clocks is primary in the integration flow. First, the clocks should be sufficiently fast to ensure the stability of the different functional behaviors (typically control laws). Second, the clocks should be not too fast to ensure that the platform has enough resource to execute the system.

The expected clock periods are described in the figure 7. We consider that FlightController is a heavier task and hence, has a period three times bigger (e.g. 150ms) than the other tasks that run with a period of 50ms.

As explained in section III, PsyC also allows more advanced control flow patterns. As an example, the agent FlightMode may need to relax the deadline when the AutoPilot mode is disabled as specified in 7. This is illustrated in the listing 10. This kind of pattern doesn't have any impact on the functional behavior (i.e. commands are unchanged in manual mode), except during mode transition, and this allows to relax the cpu load.

C. Esterel Translation

The PsyC architecture given above allows efficiant compilation. However, as illustrated, the right clocks has to be used to respect the temporal behavior. For that, one can use the synchronous semantics of PsyC, that is, its translation to Esterel to analyze the application's behaviors. Of course, theorically, such analyze could be made in the Scade application simulating the multi-rate behavior. However, we think that such approach should be made at the integration level. Futhermore, the whole application might not be available in Scade, i.e. either coming from different teams or written manually.

The listing 11 shows the translation of the Mode agent. We remove the function calls for simplicity but not the communication mechanism. In the main loop, altitude is read instantaneously then, some computation is made and the AltitudeTarget is set to a private variable. Finally, after the advance, that is, the await here, this private variable is made visible through the emit statement. Special communication events are also added to show when communication is read and when communication is produced in this module. This will be needed later for verification. The typical properties we might want to check with synchronous observers on such system are clock tick event ordering and clock tick event synchronizations. Typically, in PsyC, if the ticks of two tasks are alternating, and assuming temporal variable is at the speed of the producer, then, no value is lost. In a clock algebra such as CCSL, the predicate c 1 ∼ c 2 states that the ticks of the logical clocks c 1 and c 2 alternate. Such constraints can be easily implemented in Esterel observers as shown in [START_REF] André | Verification of clock constraints: CCSL Observers in Esterel[END_REF]. Thus, verification of such properties is then reduced to a reachability analysis on a signal representing an error state that we will call Error in this section.