Yann Rotrou
email: rotrou@enstimac.fr

Thierry Sentenac
email: sentenac@enstimac.fr

Yannick Le Maoult
email: lemaoult@enstimac.fr

Pierre Magnan
email: magnan@supaero.fr

Jean Farre
email: jfarre@supaero.fr

Jean Farré

Near infrared thermography with silicon FPA -Comparison to MWIR and LWIR thermography

Keywords: silicon FPA, thermography, radiometric model

est destinée au dépôt

Introduction

For various technological reasons, the first infrared thermographic systems to be developed were scanning cameras. They came equipped with a cooled high cutoff wavelength materials (InSb or HgCdTe) mono detector, typically operated at liquid nitrogen temperature. Such systems present good thermal sensitivity characteristics, and are often still considered as the reference of thermography in military and laboratory applications. Latter infrared measurement systems were developed with InSb, PtSi and HgCdTe detectors arrays. Such infrared systems do offer better spatial resolution, but suffer from spatial response non uniformities. Due to cooling requirement, these systems remain relatively costly. Applications call for development of larger arrays and of detectors operating at ambient temperatures or so. Recently, new uncooled microbolometer cameras, designed especially for industrial process control (surveillance, predictive maintenance, automotive…), presenting good uniformity characteristics at lower costs, have emerged. In this context, video FPA with silicon detectors seem to offer interesting prospects, and can be used to obtain good temperature measurement. They operate at ambient temperature, are relatively cheap, light weight, compact, and do offer high spatial resolution. In some industrial applications, temperature measurements have already been demonstrated with such cameras, above 800°C. (Meriaudeau et al., 2003).

The purpose of this paper is to show that accurate temperature measurements from 400 to 700°C (end even 1000°C) can be achieved with silicon FPA cameras, thanks to a specific radiometric model, and a very short calibration process. Then a comparison with traditional thermographic cameras is provided.

The first part deals with the radiometric modeling using only one set of parameters for all integration times. Starting from the classical model used for IR scanning cameras, the influence of operating in the NIR spectral range is presented. Some necessary corrections to be efficient at low temperatures are included, and finally a specific radiometric model is established.

In the second part, a very short radiometric calibration procedures is proposed and compared to the ones used classically. Accurate results over broad temperature range have been obtained. Few points of calibration are needed and the parameters are identified without initialization.

In the third part, a comparison of the performances between silicon FPA1 operating in the NIR and two well-known commercial infrared scanning cameras 2 , operating in the MWIR and LWIR, is realized. Intrinsic parameters (accuracy, fidelity, NETD, sensitivity) are presented. Finally true temperature measurement is introduced and sensitivity to the emissivity of the observed material is also discussed.

Radiometric modeling

The radiometric model of a detector describes the relation between the incident flux Ф and the output signal d I . The apparent temperature T -the temperature of a blackbody emitting the same incident flux Ф -is used.

First, one of the "traditional" radiometric models designed for infrared cameras is reminded. It is valid for a single detector with a linear photo response. Then we explain why this model is not totally suited to silicon FPA in the 300 -1000°C temperature range, and we propose a specific radiometric model. Last, we detail the corrections, required necessary to obtain performances comparable with IR cameras.

Traditional infrared radiometric model

Let us consider an IR camera observing a black body at temperature T. The system is composed of a filter and a classical objective -their spectral transmittances are noted ) ( λ τ f and respectively ) ( λ τo -, and a detector with spectral quantum efficiency ) ( λ ηe . We denote ) ( λ τa the spectral transmittance of the atmosphere. Then the output signal of a single detector with a linear photo-response is:

∫ = λ λ λ η λ τ λ τ λ τ d T L ti K T I e o f a d ). , ( ). ( ). ( ). ( ). ( . . ) ( 0 0 [1]
with ti the integration time, and where:

1 ) . exp( . ) , ( 2 4 3 0 - = - T C C T L λ λ λ [2]
is the Planck's law expressed in photons/(s.m 2 .sr.µm), where C2 and C3 are the second and the third radiation constants (C2 = h.c/k and C3 = 2.c). h is the Planck's constant (h = 6.62 10 -34 J.s), c is the speed of light (c = 2.99 10 8 m/s) and k is the Boltzmann's constant (k = 1.38 10 -23 J/K). The spectral support of ) ( λ ηe , characteristic of the detector (0.4 -1.1 µm, 3 -5µm or 8 -12 µm) is finite, and as a result the range of integration is also finite and. Lastly, to simplify, we note:

) ( ). ( ). ( ). ( ) ( λ η λ τ λ τ λ τ λ e o f a W ≡ [3]
Next, the introduction of an effective wavelength e λ allows to establish a model without preliminary knowledge of the function W(λ), and avoids the difficulty of extracting T from equation It is known [START_REF] Kostkowsky | Theory and Method of Optical Pyrometry[END_REF] that strictly speaking, the effective wavelength depends on the black body temperature. A constant C is introduced to take into account this variation:

1 ) exp( . . ) ( 0 - = T B C ti A T I d [6]
This equation describes the "classical" model with three parameters (A, B, C).

Our approach

For MWIR and LWIR scanning cameras, this radiometric model [6] is applied with a single integration time ti , and a constant effective wavelength e λ . The set of parameters (A, B, C) deals with one configuration of the camera. Typically, the temperature range from 300 to 700°C can be covered with two or three different aperture values, and so the camera is calibrated with 2 or 3 sets of parameters (A, B, C). The same approach can be used with silicon FPA, where measurement over very small temperature ranges above 300°C [START_REF] Sentenac | Evaluation of a charge-coupleddevice-based video sensor for aircraft cargo surveillance[END_REF] or over small temperature range above 800°C (Meriaudeau et al., 2003) are realized with few integration times, associated with the same number of parameters sets.

Here we want to use a great number of integration times to optimize the sensitivity of the detector, as will appear in section 4.2.1. See [START_REF] Rotrou | Demonstration of nearinfrared thermography with silicon image sensor camera[END_REF] for the on-line control of the integration time. In the spectral range3 (λ < 1.1 µm) of the silicon FPA, the luminance reaches quicker with the temperature in the near infrared wavelength range than in infrared ones. The ratio of the spectral luminances at 300 and 1000°C is about one million in the NIR (see figure 1). To measure such different fluxes, it is impossible to use only three diaphragms or integration times. Many more integration times must be calibrated, typically between 100 µs and 10 s. But this leads two main questions:

-How can we assure that each configuration gives the same temperature? This is a well-known question in infrared thermography (see for example [START_REF] Pron | Focal plane array infrared cameras as research tools[END_REF].

-How can we easily and quickly establish a relation between three variables ( Id , T , and now ti ) ? In other words, how can we calibrate the system for each integration time without one calibration procedure per integration time ? To solve the first issue, we include the integration time in our radiometric model considering the photo response of the system, and we use a single set of parameters (A, B, …) for the whole temperature range. This method allows to establish the radiometric model of the section §2.3, with an effective wavelength properly dependent on the temperature, and not on the configuration (integration time).

Finally, all FPA exhibit array non uniformities that affect the temperature performances. Moreover, the non-linearity of the photo response and/or the smear effect also degrade performances. So in the section 2.4, we quantify the influence of each of these phenomena, and propose some corrections. Note that the temperature of the detector is not controlled. This involves a variation of W(λ) (see [START_REF] Sentenac | Temperature correction of radiometric and geometric models for an uncooled CCD camera in the near infrared[END_REF]. However we consider that the temperature of the detector (typically in the 20 -30°C range) does not vary significantly, and so that it does not affect W(λ).

NIR radiometric model description

NIR spectral range and the Planck's law

Firstly we can see on figure 1 that less energy is emitted in the near infrared spectral range than in others. As a result, the minimum temperature of detection (the lower temperature that we can measure) of our system is higher than these of other systems. It is typically around 300°C [START_REF] Sentenac | Temperature correction of radiometric and geometric models for an uncooled CCD camera in the near infrared[END_REF].

We can also note that in our spectral and temperature ranges, we have:

K m T . 3000 . µ λ <<
. So the Planck's law can be replaced by the Wien's one, without introducing any significant difference (< 2 -15 ). In this context, the "1" can be neglected in the equation [6], and the coefficient C can be combined with A in a single constant; C can not describe the variation of e λ with the temperature. That fact has been confirmed by practical identifications of (A, B, C) parameters of equation [6]. It shows an important problem of conditioning, and does not provide a model fitted to the whole temperature range. So to take into account that, in the NIR, the density of energy moves to lower wavelengths when the temperature increases (i.e. the effective wavelength decreases with the temperature), we need to introduce another method; this is done in the next part.

Extended effective wavelength

As shown in figure 2, we can not use a constant effective wavelength in our temperature and spectral ranges. We propose to use the extended effective wavelength x λ , first introduced by (Saunders, 1997) in the field of pyrometry. Then ... ) (

∫ ∞ = 0 4 3 ) ( . λ λ λ d W C k w [8] ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ = ∫ ∫ - λ λ λ λ λ λ λ d T L W d C W C T T x ). , ( ). ( . . ). ( ln . ) ( 1 0 4 3 2 [9]
1 2 2 1 0 + + + = T a T a a T x λ [10]
The number of terms that we should use depends on the amplitude of the temperature range, for a given -Next the error between the reference values and those obtained using [7] and [10] with a given number of terms is evaluated.

We note that for a temperature range of 300 -1000°C, the use of a linear (first order) expression for ) 1 ( 1 T

x λ leads to a maximal error of 1.6% on x λ , or 21°C of error at 1000°C. With a second order expression, the error is less than 0.2%, or 2.5°C at 1000°C, which seems to be less important than other sources of imprecision. So we consider that when the system is equipped with a filter and used on a very large temperature range, the second order expression is sufficient. For a 400 -700°C range, the first order is sufficient as we will see further.

To conclude :

-To realize accurate measurements between 300 and 700°C, we use the relation

[7], with ) (T x λ
given by the first ordered relation [10].

-This relation is correct for a single detector with a linear photo response. These two assumptions will be discussed in the two next sections, regarding respectively CCD and CMOS silicon FPA.

Additionally, we have realized an efficient identification procedure of ( 1 0 , , a a k w ), from data ) ( T Id and without knowing of ) ( λ W , which is presented in the section 3. Keep in mind that least-square method suffer from problem of conditioning and need a lot of calibration points.

CCD silicon FPA corrections

In this part we examine phenomena that should be corrected to obtain accurate measurements.

Correction of the non linearity of the photo-response

It has been assumed that the output signal is given by [11], i.e. that it is proportional to the energy received during the integration time. In other words, the output signal is proportional to the integration time and to the incident flux.

∫ = λ λ λ λ d T L W ti K I d
). , ( ). ( . . 0 0 [11] CCD cameras that we studied do not feature a strictly linear photo-response, as we can see on the figure 3. This is due to the fact that commercially "video" cameras are tools for visualization, and not radiometric measurement. To obtain better visualization and particularly to minimize effects of saturation and blooming, detectors are equipped of an anti-blooming system. When the detector is highly illuminated, these drains positioned near the photosensitive area drives the excess signal charges into the substrate. This reduces the effects of blooming, but introduces some non-linearity. (see [START_REF] Stevens | Photoresponse Nonlinearity of Solid-State image Sensors with Antiblooming Protection[END_REF] for a typical modeling of the non-linearities). Other causes of non-linearity are the "Floating diffusion nonlinearity", because the value of the floating diffusion capacitance depends on its amount of charge, or "MOSFET non-linearity", related to the on-chip amplifier nonlinearity. Both have similar effects.

The proposed approach to take the non linearity into account consists in determining the photo-response of the detector, and generating a function of correction:

)

( 1 0 Id f I d - = [12]
where f is the photo-response of the system, Id is the read output signal and 0 d I is the corrected output signal, proportional to the integration time and the incident irradiance. The photo-response f is determined by preliminary measuring the digital output signal camera d I versus the integration time ti . Note that during this measurement, the intensity of the hot source in front of the detector is constant. It is chosen to allow to record output signal on the whole range, from about three times the noise standard deviation (low integration time) to detector saturation (high integration time). We observe a field of linearity FoL (typically from three times the noise standard deviation to the half saturation output signal) where d I is proportional to .

ti By extension we define for all ti :

) ( ) ( 1 0 Id C Id Id f I d = = - [15]
This correction is coherent with [START_REF] Stevens | Photoresponse Nonlinearity of Solid-State image Sensors with Antiblooming Protection[END_REF]. We observe that the function

Correction of the inter dependence of the detectors and smearing

Four phenomena make the output signal of two or more pixels of the detector array being cross-correlated: the blooming, the smearing, the optical reflections and the diffusion of the generated charges in the substrate of the detector.

The blooming results from of the overexposure of one or many pixels, leading to vertical or horizontal streaks. We adapt during the measurements the integration time to avoid the overexposure of any detector, and so blooming does not affect our images.

The smearing (see figure 5) is due to the reception of photons during the transfer and/or the read out of the generated charges. It is characterized by the appearance of streaks which depends on the incident irradiance and the number of detectors highly illuminated in the column or the row, but does not depend on the integration time. Its physical causes depend of the architecture of the detector; they are described in [START_REF] Theuwissen | Solid-state imaging with charge-coupled devices[END_REF]. The better way to avoid smearing is to use a mechanical shutter, which protect the detector against incident radiation during the transfer and the reading of the charges. But the use of very short integration times (down to 100 µs) prevent from using mechanical shutters. Then we should post-process the images, and two main corrections are possible. The first one consists to block the light for some rows (if the smear creates vertical streaks). Then the variation of offset of these "black" detectors is totally due to smear, and this variation is subtracted to all detectors of the column. But this system is not implemented in commercial cameras and so we use another algorithm, based on a technique proposed by (Powell et all, 1999). These algorithms have been studied for frametransfer detector, and we work to adapt them for the used interline architecture.

It is sometimes possible to suppress the smearing by reducing the aperture of the objective, but this increases the minimal temperature of detection. Finally, the smear decreases the amplitude of the temperature range that we can observe, even if we can propose some algorithm of correction when the offset due to smear is not the preponderant signal.

The reflections in the objective, also called "ghost images" (see [START_REF] Tarel | Calibration radiométrique de cameras[END_REF] and figure 6) cause a lower error that we can neglect. With traditional objectives, the coefficienet of reflection reaches typically 2%, for all wavelengths. Lastly, the charges which are generated by NIR long wavelength photons are created deep in the substrate and diffused in neighbors detectors until they are captured. This leads to a loss of spatial resolution which is presented in [START_REF] Sentenac | Temperature correction of radiometric and geometric models for an uncooled CCD camera in the near infrared[END_REF].

Non uniformity correction

As we use multi-detectors arrays, we have to take into account that all the detectors are not strictly identical. We characterize the matrix with Fixed Pattern Noise, and Photo-Response Non Uniformity. The point is to eventually correct particular detectors, in order to apply a single radiometric model all over the matrix. We use Calibration based Non Uniformity Correction and proceed like (Healey et al., 1994).

For CCD cameras, we generally do not observe significant variations of properties. For example the figures 7 represents the FPN and the dark current matrix of the EC1380 sensor. Even if we observe a strange vertical streak on the left part of the image, the difference of intensity is less than 3 Gray Levels and may eventually be neglected. The PRNU is measured for different wavelengths and does not present particular pattern.

Figure 7. FPN [Gray Level] and dark current [GL/s] map of the EC1380.

In fact there are only a little bit of detectors which show a particular behavior. For example, 99.8% of the detectors of the matrix admit a dark current less than 4 Gray Levels/s. (see the repartition in figure 8). Then the 0.2% of others detectors can be defined as "wrong" detectors or often "hot" detectors. Their values are not used, and can eventually be replaced by an average value based on their neighbors.

Figure 8. Spatial repartition of the dark current (EC1380).

Radiometric model, for a silicon CCD FPA: summary

The different steps for setting up the radiometric model, for a silicon CCD FPA, can be summed up as follows:

-Photo-response non linearity correction.

-Correction of the smear effects. The value of the 'wrong' detectors are suppressed and eventually replaced by an average of the values of their neighbors.

-Application of the radiometric model [7] and [10].

CMOS Silicon FPA corrections

Smearing effect can also be avoided by using CMOS based silicon FPA, intrinsically resistant to this phenomenon. We evaluate this architecture as an alternative to CCD silicon FPA.

Using CMOS FPA in measuring high temperatures presents an important advantage dealt with the smearing. As there is no transfer of charge with this technology, the offset due to smear affects only the detector highly illuminated, and not all the detectors of its column. So obviously no streak does appear. The correction is individual, and the importance of smearing does not cumulate the effects of all the highly illuminated detectors of a column. We can observe a typical smear effect (the offset increases with high flux) for higher temperatures than for CCD. The main disadvantages of the CMOS FPA concerns their non uniformities. This need a precise characterization. Key point are that the dark currents are higher than these of CCD FPA (that is a well-known fact presented for example in [START_REF] Magnan | Detection of visible photons in CCD and CMOS: A comparative view[END_REF], that their distribution is larger and varies quicker with the temperature detector (see figure 9). So the dark currents should be modeled and taken into account, and notably during the measurements of of low temperature (longer integration times). And last the photo response non uniformity (PRNU) must be taken into account (see figure 10), because the differences of sensibility can reach up to 20%, which causes 25°C of error at 1000°C. Unhappily we have noted that the PRNU is highly dependant of the wavelength, and even if the received energy increases quickly with the wavelength, it is not sure that the behavior of the system is accurately described by the PRNU at a given wavelength.

CMOS silicon FPA are promising for measuring high temperature, but the temperature of the detector need to be controlled, and the characterization of the FPA is more specific.

Calibration procedure

In this part, the calibration procedures are detailed. The silicon FPA one is a very short process, with accurate performances and an identification of parameters without initialization.

Calibration set up

The experiment takes place in a dark room at a stable temperature. It is composed of a reference source viewed directly by the camera. The used reference source is a P-1200 blackbody manufactured by Land Infrared. Its emissivity is 0.998 and its aperture is 54 mm. The studied temperature range is 400 -700°C.

The infrared cameras are manufactured by Agema. 0ne is the 780 camera, equipped with an InSb mono-detector (3 -5 µm) cooled at 77 K and the other is the 880 LWB camera with an HgCdTe mono-detector (8 -12 µm), cooled at 77 K with a Stirling cycle. Their usual frame rate is 6 images per second. The CCD silicon FPA is the Prosilica EC1380.

The alignment blackbody-camera is provided by an optical bench. The distance between blackbody-camera is fixed at 1.2 meter.

Traditional calibration process

The calibration of IR cameras, or silicon FPA in other works (for example Meriaudeau et al., 2003) are characterized by many measurements at different temperatures. A non linear identification needing a precise initialization of the set of parameters is applied.

Here, for the infrared cameras, we do acquisitions from 400 to 750°C, with an increment of 25°C, for each f-number. Each of these acquisitions (Id, T) is an average on one hundred images, to remove temporal noises. The results are shown in figure 11 andsection

4.1.
Here is a set of parameters for the MWIR camera: N = 10; A = 3.00 GL; B = 3.36e3 1/(µm.K); C = 1.78.

Specific calibration process

For silicon FPA, the first step is the determination of the photo response of the detector, as described in §2.3.1.

Then, only three measurements at high temperatures (600, 650 and 700°C) are needed to find w k and the parameters which define ) ( T x λ . So the calibration procedure is very short. For each temperature, we acquire a temporal average of the intensity on two regions of interest, one corresponding to a uniform region of the black body, and a second which is on the same columns, in a non illuminated region; that allows to estimate the offset due to smear. More precisely, we acquired for each temperature different values of Id for different integration times; this is a very short process, which allows to estimate more precisely the output signal and provides a response in Gray Levels per second (see figure 12).

To determine the value of the set of parameters, we have shown the efficiency of the following process: over a small temperature range, we know that T a a

x 1 0 1 + = - λ .
For a given amplitude of the temperature range, this approximation is more true when the temperature is higher. So we consider the three higher points (Id, T), here at 600, 650 and 700°C. These three points provide three relations [17] and so values of the three parameters w k , 0 a and 1 a : For the 400-700°C temperature range, this step is sufficient. The 400°C measurement only provides the confirmation of the fitness of the model (see next section).

Its parameters are: k w = 2.12 e11 GL/s; a 0 = 1.10 e6 m -1 and a 1 = -3.02 e7 m -1 K -1 .

Comparison of performances

This section describes the results of the comparison between two IR cameras and the CCD silicon FPA. These three cameras provide a 12 bit signal per pixel. First, we introduce some criteria to characterize the parameterized radiometric models. Next, intrinsic performances such as sensitivity and Noise Equivalent Temperature Difference (NETD) are compared. Finally, we introduce some information dealing with the true temperature measurements.

Evaluation of the radiometric models

We define three criteria intended to evaluate the performances of the radiometric models : the trueness, the precision, and the maximal error. We evaluate and compare these criteria for each camera. Note that such criteria are not only used to compare different systems, but also serve in study and identification of potential "problems". It is indeed possible to estimate this criteria on only a part of the set of calibration points and that may show some interesting phenomena; for example the non linearity of the photo response was pointed out by noting that the trueness was smaller for higher output signals.

Definitions

We consider the set of calibration data (Id, T) and we note E[x] the mean value of x. Otherwise we note T ˆ the estimation of T provided by the radiometric model, knowing Id.

Then the trueness J of the parameterized model is defined [START_REF] Perruchet | Estimer l'incertitude[END_REF]:

] [ 1 T T E J - = [19]
The precision F and the maximal error max

E are defined as: ] ) [( 1 2 T T E F - = and ) max( max T T E - = [20]

Results

The next table summarizes results for different configuration of the three cameras. For the three cameras, the considered set of data is the mean output signal on the region of interest. The on-line adaptation of the integration time (see [START_REF] Rotrou | Demonstration of nearinfrared thermography with silicon image sensor camera[END_REF] for a description of the algorithm) allows to optimize the sensibility of the system, but leads to use a great number of integration times.

Model

The sensitivity of the system is traditionally defined as (

0 d I Id = for IR cameras): T T Id T S ∂ ∂ = ) ( [21]
For infrared systems using the radiometric model [6], it can be analytically computed: This is a maximal sensitivity which is not really obtained. Indeed if we consider that 0 d I is not the read output signal but a corrected signal due to the non linearity of the photo response, we should take into account this photo response f. The relation :

) ( . . ) ( ) ( 0 0 0 0 0 0 0 T S I f T I I f T I f T Id T S d d I d d I d T d T ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂ = ∂ ∂ = [24]
shows that the loss of sensitivity is easily linked to the loss of sensitivity of the photo-response. Here the higher loss is 7%.

If there are smear effects, the sensitivity indirectly decreases; indeed, as the smear creates an artificial offset, we should decrease the integration time to not saturate, thus 0 Anyway, it appears that silicon FPA are promising and can present a better sensitivity than infrared systems. We should lastly note that the sensitivity varies quickly with the temperature; that means that it is not constant in an image with thermal gradients, and that we can not assure that very good sensitivity for all the image if it is very contrasted.

Generalized NETD

The Noise Equivalent Temperature Differences NETD is defined as the ratio of the standard deviation noise to the sensitivity. It is computed for Signal Noise Ratio equal to 1 and a temperature of 30°C. As 30°C is not a detectable temperature with silicon detectors, we can not compare the systems.

We propose to generalize the notion of NETD by introducing the temperature of the black body that we observe. So the GNETD -G for generalized -, at the black body temperature T, is the ratio of the standard deviation noise to the sensitivity, both considered for the temperature T. In figure 14, we present result considering that the standard deviation is constant and unitary. Found values for IR detectors are comparable with values given by the manufacturer.

Regarding the emissivity of a real material, we know that the measurement is less sensitive to the uncertainty of the emissivity -or its spatial variation -, at lower wavelengths. Indeed considering a monochromatic measurement, we see that the For a given state of surface ( ε d ), the uncertainty dT is about 10 times less important in the NIR than in the LWIR (see figure 15). A study published in [START_REF] Rotrou | Demonstration of nearinfrared thermography with silicon image sensor camera[END_REF] shows that the error on the NIR image of a mould is less than 3°C (between 300 and 520°C). During the cooling of the mould, gradients of temperatures appear due to its shape. These gradients (few dozen of degrees) are shown on the NIR image. On the LWIR image, the differences of apparent temperature are higher (more than one hundred degrees) and so they are not due to the cooling of the mould. The observed gradients come from gradients of emissivity.

True temperature measurement

In this last part we introduce some notes about the measurement of true temperature. Indeed the radiometric model allow to obtain black body temperature, or the incident flux. Some external factors must be taken into account to access to the true temperature of the observed object: the emissivity of the material, the temperature of the air, the shape of the object which determines the relation between radiance and emitted fluxes.

Conclusion

In this paper, we present a general radiometric model of a silicon FPA. We compare its performances with those of two commercial MWIR and LWIR cameras.

One of the originalities of the work is that only one set of parameters is determined for all integration time. The study leads to an easy-to-use and effective calibration procedure, without tedious measurement of W(λ).

The comparison with two commercial infrared cameras shows that silicon FPA properties are promising. Thanks to the real time control of the integration time, the sensitivity of our system reaches more than 60 Gray Levels by degree and offers values higher than the sensitivities of infrared cameras. The Generalized NETD of these three systems, introduced as a new metric, are comparable, around 0.1°C. Additionally, as it works at shorter wavelengths, our system offers the advantage of being less sensible to the emissivity of the material. So its use is particularly recommended, for example, in the study of non uniform materials (e.g. forming moulds, metallic steels with degraded state of surface…), and other applications requiring good spatial resolution.

Finally, we demonstrate with this work that silicon FPA complement the large range of thermographic systems commercially proposed. Further investigations are pursued, with CMOS silicon FPA.

Figure 1 .

 1 Figure 1. Planck's law (semi logarithmical scale).

  temperature, and that w k is by definition a constant.By comparison with the relation [1], we obtain the following definitions:

  This model offers an important advantage: only one parameter x λ depends on the temperature. In this context we can numerically determine k w and which is not possible by introducing the temperature influence in relation[5].

Figure 2 .W

 2 Figure 2. Left: Measurement of W(λ). Right: ) ( T x λ

  evaluate the precision that we can obtain for a given number of terms, we measure the relative ) is calculated; thus providing the reference output signals.

Figure 3 .

 3 Figure 3. Error of non linearity, here denoted transfer rate vs. output signal (EC1380).

  on the incident spectral distribution.We see on the figure3that the non-linearity reaches 7%, what causes an error of 10 K at 1273 K. We should correct such a deviation.

Figure 4 .

 4 Figure 4. Photo response of the CCD EC1380 camera: read output signal Id vs. corrected output signal 0 d I , in Gray Levels.

Figure 5 .

 5 Figure 5. Observation of a circular black body and smearing. On the left: first illuminated image in which a half spread is observed; on the right: a entire spread, observed without mechanical shutter.

Figure 6 .

 6 Figure 6. Observation of a black body and of its "ghost" on the bottom right part of the image.

Figure 9 .

 9 Figure 9. Evolution of the spatial repartition of the dark current for a CMOS camera, with the temperature of the detector (20 and 29°C).

Figure 10 .

 10 Figure 10. PRNU of the CMOS CV1280 camera, at λ=880 and 930 nm.

Figure 11 .

 11 Figure 11. Experimental results for the infrared cameras.

Figure 12 .

 12 Figure 12. Normalized calibration curve of the EC1380 camera, for an integration time of 1 second.

Figure 13 .

 13 Figure 13. Sensitivity [Gray Levels / °C].

Figure 14 .

 14 Figure 14. GNETD supposing the noise is constant and unitary.

Figure 15 .

 15 Figure 15. On the left, observation of a real surface with LWIR camera and, on the right, NIR camera. On the left, we see more spatial variations of emissivity than temperature gradients. On the right, isotherm curves are really observed.

The CCD silicon FPA is the Prosilica EC1380. The CMOS silicon FPA studied in section

2.4 is the Prosilica CV1280.2 The IR cameras are presented in section 3.

To avoid the influence of solar radiations, we can work in complete darkness or use a filter.Here we consider that we work in complete darkness, in the 0.7 -1.1 µm range.


may come

Actinic fluxes in Titan's atmosphere, from one to three dimensions Application to high latitude composition Sébastien Lebonnois, Dominique Toublanc


R Saidi

J C Olivier

M Machmoum

Eric Chauveau
email: eric.chauveau@eseo.fr

Energy management and sizing algorithm applied on a hybrid power system suppliying an isolated residential application

Keywords: energy management, centered moving average filters, forecasting errors, power forecasting

This paper presents an energy management algorithm to be embedded for a residential application disconnected from the power grid. The multisource system studied includes solar panels as a renewable energy source, a fuel cell as a secondary source, three batteries and a bank of supercapacitors as storage systems. The proposed algorithm uses as input data, future estimates obtained from power consumption and meteorological forecasting data, and historical of the load power and the renewable power obtained from measurements. It is assumed that these two input powers are imposed and uncontrollable. In the case where the future estimates have errors compared to what has actually been measured, a mathematical approach shows that the algorithm is able to compensate these forecasting errors by sharing them between the different sources of the system while respecting their different characteristics. In addition to the optimal energy distribution, the algorithm gives the optimal size of each source and storage devices. In this work, the total cost of the system is chosen as the criteria to be optimized. The simulations are carried out over a year, with a time step of 1 second and in the presence of significant forecasting errors. The results obtained are particularly convincing and make it possible to validate this energy management strategy.

I. INTRODUCTION

Nowadays, multisource production systems represent a major challenge for the energy production and transport industries. Thanks to these systems, the supply of isolated applications has become possible. In addition, these systems are often based on one or more sources of clean energy such as renewable energies or hydrogen, which makes their studies more and more interesting as the world tends to decrease the exploitation of conventional energy production systems such as power plants or combustion engines, because of the harmful effect they generate on the environment. These multisource systems can contain several energy sources of different characteristics. Indeed, the presence of several sources in a multisource system increases its lifetime and its reliability and on another hand avoids the complete shutdown of the system in case of malfunction of one of its sources [START_REF] Bellaouar | Fiabilité, maintenabilité, Disponibilité[END_REF].

Moreover, during the lifetime of the power system, one or more devices can be subject to replacement, or even substituted by new more efficient technologies. Then, it may be interesting that the architecture of these multisource systems was flexible and adaptable. The distribution of the energy between the different devices and the fed load follows a set of laws that are called energy management strategies (EMS). But the coupling of several energy devices in the same production system makes the arbitrage and the management of energy flows more complex. This is partly one of the reasons why several authors limit their system to one or two sources associated with a single storage organ, thus making the system less complex and with limited degrees of freedom. In fact, many authors use dynamic programming (DP) as EMS [START_REF] Eng | Rule based energy management strategy for a series-parallel plug-in hybrid electric bus optimized by dynamic programming[END_REF]. It is a very efficient technique and well adapted to energy management. It consists of sampling the power margin of the devices in n states and sampling the time axis in k stages and it is up to the algorithm to determine the optimal system state at each stage. As a result, the more the number of sources increases, the more the number of states multiplies and a multidimensional resolution grid is created, thus making the search time of the optimal paths very long, especially as soon as the number of paths to be optimized exceeds two or three. Other authors use fuzzy logic (FL) to share energy between the sources of the system [START_REF] Arcos-Aviles | Energy management strategy for a gridtied residential microgrid based on fuzzy logic and power forecasting[END_REF], [START_REF] Erdinc | A wavelet-fuzzy logic based energy management strategy for a fuel cell/battery/ultracapacitor hybrid vehicular power system[END_REF]. For this technique, a large number of sources cause a large number of membership functions, which also makes difficult the definition of the different rules. Other works such as [START_REF] Hredzak | Model predictive control of a hybrid battery-ultracapacitor power source[END_REF] are based on predictive control. This method is particularly effective for taking into account resource forecasting and usage profiles. This technique requires linear modeling of each device to be effective, which means that some important non-linear phenomena have to be neglected. Finally, the frequency separation is one of the most used EMS thanks to its ability to control sources by power profiles that go with their respective dynamics [START_REF] Erdinc | A wavelet-fuzzy logic based energy management strategy for a fuel cell/battery/ultracapacitor hybrid vehicular power system[END_REF], [START_REF] Bourdim | Efficient Energy Management Strategy for Fuel Cell Ultracapacitor Hybrid System[END_REF], [START_REF] Uzunoglu | Modeling and analysis of an FC/UC hybrid vehicular power system using a novel wavelet based load sharing algorithm[END_REF]. The work proposed in this article is in the framework of this last technique. Indeed, the proposed strategy is based on the cascading and paralleling of Centered Moving Average (CMA) Filters that generate all the reference profiles to the system [START_REF] Arcos-Aviles | Energy management strategy for a gridtied residential microgrid based on fuzzy logic and power forecasting[END_REF], [START_REF] Saidi | Energy Management Strategy for Hybrid Power Systems based on Moving Average Filters and Power Forecasting[END_REF]. The principle of the strategy appears applicable on a system composed of many energy sources. But instead of considering only cascaded filters, with a classical decomposition from low to high frequency ranges, we propose here to extend this idea to more general control architecture, with serie-parallel filters. This new control architecture permits to easily include different technologies of devices, on the one hand to ensure an optimal technical solution, and on the other hand to allow redundancy which improve the global reliability of the system. In addition to the optimal management of energy, the proposed algorithm allows to size these devices in such a way to minimize the total cost of the multisource system. The choice of the total cost of the system as the variable to be optimized, is explained by the fact that with this cost, we can evaluate the effect of several phenomena that the multisource system undergoes as for example aging, lifetime, maintenance, replacement and acquisition. In other words, the final cost of the system will reflect the total cost of one kWh produced. This paper is organized as follows: Section II describes the architecture of the multisource production system. Section III presents the proposed algorithm applied to the multisource production system. Section IV presents the objective of optimization chosen in this study. The simulation results using predicted power profiles are presented in Section V. Finally, Section VI presents main conclusions of this paper.

II. MULTISOURCE SYSTEM DESCRIPTION

Figure 1 shows the architecture of the studied hybrid system. It includes photovoltaic (PV) panels as a renewable energy source, a fuel cell (FC) as a secondary energy source, two batteries Li-ion, a lead acid battery and a bank of supercapacitors (SC) as energy storage systems (ESS). With the exception of PV panels, dimensions of all sources are not fixed. For the renewable source, 10 of PV panels are used. The DC / DC converters associated with PV panels and the FC are unidirectional Boost converters. On the other hand, for ESS, Buck-Boost bidirectional DC / DC converters are used to enable their charging and discharging phases. A negative power corresponds to a charging phase whereas a positive power corresponds to a discharge phase. The representation of sources in this work has remained restricted on classical empirical models. The studied habitat is supposed located at Saint-Nazaire, France. For the estimation of the PV power profile , a solar profile generator has been developed and is presented in [START_REF] Bouabdallah | Robust sizing of a stand alone multi-sources power system[END_REF]. This generator based on the Marcov matrix, generates an estimation of the radiation profile according to the geographical coordinates of the studied location (latitude / longitude), the inclination angle of PV panels and a sampling step. The PV power profile used in this study corresponds to the inputs shown in Table 1: In addition, an estimated domestic load demand profile was generated by a consumption profile generator algorithm that operates according to probabilities of load appearance along the day. This profile corresponds to the power consumption of 5 persons. It is assumed that the habitat does not include an electric heater but rather a refrigerator, lighting, an oven, a washing machine, a dishwasher and multimedia.

Figure 2 shows the used load profile on three different time scales. These representations show that the more the time horizon considered decreases the more the dynamics of the load profile stabilizes. In other words, we can consider that the more the energy source dynamic is slow, the more its reference power on the reduced time horizons is stable and also vice versa. From these last considerations, an algorithm based on a multi-time horizons operation is developed.

III. ENERGY MANAGEMENT ALGORITHM

The proposed energy management algorithm is based on cascading/paralleling of several blocks called devicecontrollers, as shown in Figure 3. The total number of controllers is equal to the number of sources in the system. Each controller calculates, in real time, the reference power that a corresponding energy source must supply/store. In fact, at each instant t, it has at its input a set of estimated future power values and a set of past power values taken from a data history [START_REF] Bouabdallah | Robust sizing of a stand alone multi-sources power system[END_REF]. Based on these two future and past inputs, it gives the good reference power that its corresponding source has to supply/store at the present moment t. The proposed algorithm, illustrated in Figure 3, has for main input, two vectors ∆ P and ΔP. Both represent the difference between the power demanded by the load and the power produced by PV panels over a simulation horizon :

= _ - _ 1 = _ - _ 2
Fig. 3.Architecture of the proposed energy management algorithm.

The internal structure of a device controller is shown in Figure 4. It mainly comprises a Centered Moving Average (CMA) filter [START_REF] Arcos-Aviles | Energy management strategy for a gridtied residential microgrid based on fuzzy logic and power forecasting[END_REF], [START_REF] Pascual | Energy management strategy for a renewable-based residential microgrid with generation and demand forecasting[END_REF] and safety devices (maximum powers (block in green), exceeding the SoC limits (block in red)). At each instant t, each filter calculates a centered average of its input vectors (future estimated vector and past measured vector) at a coefficient * near, over a well-defined time horizon * . The mathematical relationship between the inputs and the output of each CMA filter at instant t is described as follows:

* = ! * ! * " 2 $_ -$_ % &' * ( ) + ! * ! * " $_ &+ * ( % ) 3
The output power expressions of filters are integral functions that connect the future to the past. In other words, these expressions relate past measured data that are always true, to predicted future data that may be erroneous. Taking into account the hypothesis (4), it becomes possible to say that the average of the output signal of each controller on , is equal to that of its measured signal input:

< > . = 0< > 1 . 2 < (4)
This observation means that 34 * , the output power coming out of the first controller will have a mean value equal to the input vector , what makes the average of $_ , the input power of the second controller, null and consequently, its outputs * ensure that their respective States of Charge do not diverge along . The output power of a filter represents the first step in the construction of the reference power of the energy source. Each filter differs from the others by its time horizon * . Indeed, the value of the time horizon of a filter depends on the dynamics of the source at which its output power is intended. The larger the time horizon, the smoother the filter output and vice versa. For this, the reference powers of slow dynamic sources must come from device-controllers containing filters operating with wide time horizons, while this time horizon can be smaller for fast sources. The coefficients : * reflect the maximum and minimum power level that the corresponding source can reach. Since the values of * and : * affect directly the output power of the device-controllers, the pairs (: * , * ) play the role of the optimization variables in order to minimize the objective function described in the next section. The main objective of the algorithm is to manage the measured unbalanced power ∆P to ensure the equality between the power demand and the production. In other words, if the unbalanced power is positive, the controllable devices will supply power to the system, and on the other hand, if it is negative, the excess of power will be stored in the storage systems. A difficulty is to propose control architecture able to handle a large number of devices (sources and storages) with very different characteristics (in power and energy densities).

To illustrate this problem, we propose to use an algorithm, which is inspired by frequency separation methods. Indeed, the cascade structure of the algorithm presents 4 stages according to the dynamics of the controlled energy sources. As for classical frequency separation, the first stage considers very high specific energy and low specific power (very slow source). On the other hand, the last stage (4 th ) consider very high specific power and low specific energy (very fast sources). Each other stage will be used to intermediary devices. Indeed, among the five sources, the FC is the only source that is characterized by a high specific energy and a low specific power. As shown in Figure 3, this energy source will be controlled by the device-controller 11 placed on stage 1. Unlike FC, SC are storage devices that are characterized by high specific power and low specific energy. Then, their reference power must be a very high frequency profile. For this, this power is ensured by stage 4, which represents the residual power. Batteries are energy storage systems that are positioned between these last two devices. Without loss of generality, it is considered here one Lead-acid battery pack and two Li-ion batteries. According to the Ragone diagram [START_REF] Lee | Analogical understanding of the Ragone plot and a new categorization of energy devices[END_REF], lead-acid batteries have a lower specific energy than Liion batteries. For this, the device-controllers of the stage 2 ensure power of the two Li-ion batteries. The stage 3 is then devoted for the Lead-acid battery. Actually, as the proposed algorithm uses predicted values of load consumption and also solar radiation, these values may have deviations from the measured values. This difference is called forecasting error (FE) , defined as follows:

= ∆P + ; 5 Indeed, with keeping the same assumption (4), by substituting the estimated power unbalance with equation ( 5), it appears an automatic error compensation term in the expression of each reference power. In fact, the profile of the forecasting error ( -∆P) is automatically shared between the different sources while respecting the dynamics of each one. Equation ( 6

E6

Where 77 * is the optimal horizon of the FC corresponding filter.

Indeed, the automatic correction provided by the FC consists of increasing or decreasing its reference power according to the sign of the average value of the error over [START_REF] Saidi | Energy Management Strategy for Hybrid Power Systems based on Moving Average Filters and Power Forecasting[END_REF], [START_REF] Saidi | Méthode générique de gestion de l'énergie des systèmes multisources par filtres à moyenne glissante[END_REF]. As for the power profiles of the batteries and SC, the same phenomenon occurs. In fact, averaged terms of forecasting errors appear in the expressions of their reference powers. What is interesting in these results is that the sharing of this automatic compensation between the 5 sources of the system is done while respecting the dynamics of each source and also while keeping the average annual power of each ESS null, which guarantees the non-divergence of its SoC even in presence of forecasting errors. This previous finding is mathematically demonstrated in [START_REF] Saidi | Energy Management Strategy for Hybrid Power Systems based on Moving Average Filters and Power Forecasting[END_REF] and [START_REF] Saidi | Méthode générique de gestion de l'énergie des systèmes multisources par filtres à moyenne glissante[END_REF].

IV. OBJECTIVE OF OPTIMIZATION

The couples (: 77 , 77 ), (:  respect the dynamics of the five sources of the system.

 respect the power limits of the five sources.

 respect the states of charge(SoC) of the storage systems.

 correspond to a minimum total cost of the system. This cost is explained below. In this study, optimization is applied to the total cost of the system. This cost includes four components: the total initial cost, the total operation cost, the total replacement cost and the total maintenance cost. Given that these different costs are not evenly distributed over time and since the different sources do not have the same lifetimes, it will then be necessary to unify the period of time over which these costs will be calculated. In this case, the notion of annualized cost (cost per year) is used.

A. Devices lifetime

A degradation function ∆ GH is considered for the FC. This function depends on its provided energy over the simulation horizon 9 [START_REF] Depature | IEEE VTS Motor Vehicles Challenge 2017 -Energy Management of a Fuel Cell/Battery Vehicle[END_REF]:

∆ GH = " ∂ ) J K 7 Where ∂ = L M 8NKK O1 + Q RS_TUV A 34 * -34_$ W ²Y
where ∂ K , α are coefficients and 34_$ W the FC nominal power (80% of 34_W Z ). The degradation function ∆ GH is expressed between 0 (beginning of life) and 1 (end of life).

Once the degradation of the FC has been calculated, its lifetime is estimated as follows:

[\ 34 = 1]^:_/∆ GH =
The operation of a FC requires hydrogen. The consumption of this one is proportional to the energy provided by the FC. The volume 2 _ $$ ( 8 ) of the consumed hydrogen over a year is calculated as follows [START_REF] Zhang | Optimization for a fuel cell/battery/capacity tram with equivalent consumption minimization strategy[END_REF]: 2 _ $$ = < 34 * > a b * 9 * 2.7778. 10 'g / 3 * η 34

With < 34 * > a b the average power supplied by the FC along the horizon 9 and η 34 the FC efficiency.

For batteries, ∆ jkl , a function describing the evolution over time of their life. This evolution is proportional to the energy that the battery exchanges compared to the maximum energy that it can exchange throughout its lifetime. Generally, a battery is able to perform a few thousand full charge/discharge cycles throughout its life. ∆ jkl is expressed as follows [START_REF] Bouabdallah | Robust sizing of a stand alone multi-sources power system[END_REF]:

∆ jkl t = 1/ n o 6 " pn qZrs $tu p . K ) 8
With n qZrs $tu the energy exchanged by the battery during and n o 6 the total energy that it is able to exchange during its life. This energy is usually obtained by multiplying the storable energy by the maximum number of cycles. Once the equality |n qZrs $tu |=n o 6 is verified, the battery reaches its end of life. The degradation function ∆ jkl is expressed between 0 and 1. The estimated lifetime of the battery can then be calculated:

[\ 5 6 = 1]^:_/∆ jkl =
∆ vH , the degradation function of the SC bank [START_REF] Hmam | A cycle-based and multirate approach for power system simulation application to the ageing estimation of a supercapacitor-based ferry[END_REF], [START_REF] Trieste | Techno-economic optimization of a supercapacitor-based energy storage unit chain: Application on the first quick charge plug-in ferry[END_REF]. ∆ vH depends on the state of charge of the bank (SoCsc) :

∆ vH = " ∂ ) . K 9 Where ∂ = K,K7w
8NKK * x * 8Nw * 2 7 ,w y9 4 r 6 '7

Thus, the estimated lifetime of the bank is deduced as follows:

[\ 94 = 1]^:_/∆ vH =
The renewable source generally has the longest lifetime among the sources of a multisource system. For this, the lifetime of a multisource system is generally equal to that of the renewable source present in this system. The lifetime of PV panels is estimated at 25 years.

B. Devices total costs

Taking into account the previous lifetimes, the annualized costs of each device are calculated [START_REF] Bouabdallah | Robust sizing of a stand alone multi-sources power system[END_REF], as shown in table 2.

Table 2. Annualized costs of each device of the studied multisource system.

PV Panels

Annualized Initial cost

= $ 6_ _ $$ = = $z_ + 0.4 * = $z_ /25 Annualized maintenance cost = W $_ _ $$ = 0.01 * = $ 6_94_ $$ Annualized total cost = o 6_ _ $$ = = $ 6_ _ $$ + = W $_ _ $$ Fuel cell Annualized Initial cost = $ 6_34_ $$ = = $z_34 /[\ 34 Annualized operation cost = }~_34_ $$ = 2 _ $$ * _ W8_ Annualized replacement cost = •u~_34_ $$ = 0.04 * = $ 6_34_ $$ Annualized maintenance cost = W $_34_ $$ = 0.04 * = $ 6_34_ $$ Annualized total cost = o 6 RS€TT = = $ 6_34_ $$ + = }~_34_ $$ += •u~_34_ $$ + = W $_34_ $$ Battery Annualized Initial cost = $ 6_5 6_ $$ = = $z_5 6 /[\ 5 6
Annualized replacement cost

= •u~_5 6_ $$ = 0.03 * = $ 6_5 6_ $$ Annualized maintenance cost = W $_5 6_ $$ = 0.03 * = $ 6_5 6_ $$ Annualized total cost = o 6_5 6_ $$ = = $ 6_5 6_ $$ + = •u~_5 6_ $$ + = W $_5 6_ $$ SC bank Annualized Initial cost = $ 6_94_ $$ = = $z_5 6 /[\ 5 6 Annualized replacement cost = •u~_94_ $$ = 0.03 * = $ 6_5 6_ $$ Annualized maintenance cost = W $_94_ $$ = 0.03 * = $ 6_5 6_ $$ Annualized total cost = o 6_94_ $$ = = $ 6_94_ $$ + = •u~_94_ $$ + = W $_94_ $$
Where = $z_ the investment cost of PV panels, = $z_34 the investment cost of the FC, = $z_5 6 the investment cost of a battery, = $z_94 the investment cost of SC bank and _ W8_ the price of one cubic meter of gaseous hydrogen. Unlike the rest of the components, the initial cost of PV panels includes a cost of installation and acquisition in addition to the purchase cost. Installation and acquisition costs are estimated at 40% of the purchase cost.

C. Annualized total cost of the system Finally, the total annualized cost of the system is given by [START_REF] Pascual | Energy management strategy for a renewable-based residential microgrid with generation and demand forecasting[END_REF]:

= o 6_ • _ $$ = = o 6_ _ $$ + = o 6_34_ $$ + = o 6_5 67_ $$ += o 6_5 6 _ $$ + = o 6_5 68_ $$ + = o 6_94_ $$ 10
This cost represents the objective function to be optimized.

D. Constraints of the problem

As explained in Section 3, the values of time horizons depend on the dynamics of the sources. For the generator 11, in order to avoid the negative values in the generated power (the FC is a unidirectional source) and also to obtain a smoothed profile, the horizon 77 * must be wide. It is chosen to be of the order of days. Li-ion batteries are assumed to do a charge or discharge cycle in less than half a day. For this, horizons will not be well adapted to the capacity of the SC bank (too much energy). Concerning the coefficients, we consider that they must be positive to avoid the risk of reversing the sign of the power at the input of a filter (which causes charge instead of discharge and vice versa). The coefficient must also be less than 1 in order not to amplify the profiles at the input of the filters and consequently ask to supply or store an unnecessary power. Taking all this into account, the following constraints are defined: On stages 1 and 3, there is no power sharing in these two stages and therefore : 77 * and : 87 * are set at 1.

Added to that, the dimensions of ESS that will be given by the algorithm must respect the following constraints:

20% ≤ ‡ˆ= 5 6 ≤ 80%, 0 ≤ ≤ 11 10% ≤ ‡ˆ= 94 ≤ 100%, 0 ≤ ≤ 12
Where ‡‰= Š 6 = -1/n 5 6 " 5 6 * )  n 5 6 (J) represents the total energy of the battery.

And

‡ˆ= 94 = -1/n 94_‹ Z " r * ) 

Where n 94_‹ Z = 1/2 = 94 2 94_‹ Z ² n 94_‹ Z (J) the maximum energy that can be exchanged by the SC, = 94 (F) its total capacitance and 2 94_‹ Z (V) the maximum voltage at its terminals. In fact, through the constraint [START_REF] Lee | Analogical understanding of the Ragone plot and a new categorization of energy devices[END_REF], the deep charges/discharges are prevented because they accelerate the degradation of the batteries. On the other hand, SCs are devices that support deep charges/discharges. The lower limit of the constraint ( 12) must be strictly positive in order not to generate an infinite current of the SC. It is fixed at 10%.

V. SIMULATION RESULTS

Using the two estimated profiles cited in section 2, solving the previous optimization problem gives the following results: The power profiles associated to the previous values correspond to a total annualized cost of the system equal to 28044 €, in other words, equal to 7.25 € per kWh.

In a first simulation, it is assumed that the future estimates coincide perfectly with the powers actually measured. In this case, Figure 5a shows the evolution of <= * over a day (in blue). The slow dynamics and the absence of sudden peaks and changes in power are the main characteristics of this profile. Figures 5b and5c show the daily evolutions of the reference powers of Li-ion batteries (in blue). These two profiles are each characterized by the presence of two phases of charging and discharging each lasting about 12 hours. In fact, during the day, the batteries are charged by the excess of power that is created by solar panels. Figure 5d shows the daily power profile of lead-acid battery (in blue), which presents significant variations compared to those of Li-ion batteries. The figure 5e shows the residual power for the SC bank (in blue). It contains power peaks over very short periods of time, which means high powers and low energies. Then, in a second simulation, forecasting errors are considered according to the equation [START_REF] Zhang | Optimization for a fuel cell/battery/capacity tram with equivalent consumption minimization strategy[END_REF]. Indeed, future estimates of ΔP are assumed to be lower than what will actually be measured: = + 10% , 0 ≤ ≤ 13

Figure 5a shows the new power profile of the FC over a day (in red). This profile presents higher powers than those found in the first simulation, while keeping the same rate over time. This result is explained by equation ( 6), which reflects the evolution of <= * over time. This new profile comes from the sum of the profile without FE (blue) and an averaged term of forecasting errors. As for the power profiles of the batteries and SC, the same phenomenon occurs. The reference powers of the three batteries and the SC bank are shown in the figures 5b, 5c, 5d and 5e respectively. Figure 6 shows the evolution of SoC of batteries and SC in both cases of simulation. In the presence of FE, SoCs exceeded safety limits because of the additional compensating powers. For this reason, the dimensions given by the optimization algorithm are the minimum dimensions to use. In the sizing phase, forecasting errors should be taken into account.

VI. CONCLUSION

This work presents an energy management strategy used to optimize energy distribution and sizes of the sources in a multisource production system feeding an isolated habitat. The studied system is composed of solar panels as a renewable source, a fuel cell, three batteries and a bank of supercapacitors as storage devices. It has a cascading architecture of several blocks called device-controllers, which operate mainly based on centered moving average filters coupled with safety and supervision functions. The proposed algorithm uses, as input data, a future estimate and also a measurement history of the load power and the photovoltaic power. In addition, the mathematical approach of the algorithm shows that it is able to compensate the forecasting errors that can appear as a result of bad estimations. The simulation results showed the effectiveness of the proposed algorithm to respect the different power and energy constraints of the sources and also behaving with forecasting errors. As a perspective, the algorithm is to be generalized by adopting an architecture composed of 4 stages with N controllers per stage.

Fig. 1 .

 1 Fig. 1.Architecture of the studied multisource system.

Fig. 2 .

 2 Fig. 2. Representation of the requested power profile on different time scales.

Fig. 4 .

 4 Fig. 4.Internal structure of a generator.

  also will have a zero-average on . The same analysis is done for the two stages that come after. The zero-averages of 5 6 7 * , 5 6 * , 5 687 * and 94

  ) reflects the evolution of <= * over time in presence of FEs. It presents an averaged term of forecasting errors. In other words, the FC automatically compensates the 77 * -average value of this error over the horizon :

Fig. 5 .

 5 Fig. 5.Daily reference power profiles in presence of FE (red) and without FE (blue).

Fig. 6 .

 6 Fig. 6.Evolution of ESS's SoC in presence of FE (red) and without FE (blue).

Table 1 .

 1 Inputs of the solar profile generator..

	Inputs	Value
	Position (decimal degrees) Saint Nazaire, France X: 47,283329 Y: -2,2
	PV panels inclination (deg)	50
	Sampling step (s)	1
	Hs:Simulation horizon (year)	1

Table 3 .

 3 Optimization results.

	( 7,7 * ,: 7,7 * ) ( ,7 * ,: ,7 * ) ( , * ,: , * ) ( 8,7 * ) * ,: 8,7 Maximum power of the FC	(11.57 days, 1) (11.98 h, 0.47) (6.32 h, 0.53) (6.66 min, 1) 343.4 W
	Li-ion battery 1 min size	12008 Wh
	Li-ion battery 2 min size	14677.5 Wh
	Lead acid battery min size SC bank	6375 Wh = 94 =50.65 F, Pmax=3206 W