Cascaded Centered Moving Average Filters for Energy Management in Multisource Power Systems with a Large Number of Devices

: Hybrid systems constitute one of the solutions for supplying isolated applications. Such systems are classically based on clean energy sources. When the renewable energy sources have intermittent productions, they are associated with storage systems. This makes the system economically more interesting. Economically speaking, hybrid energy systems using multiple energy sources are often expensive and their cost must be optimized. This optimization can be done for the system sizing or for its energy management. However, optimizing one does not guarantee the optimization of the other. Indeed, previous studies optimize either the design and apply it with a simple energy management strategy, or the energy management with predetermined sizing supposed optimized, while minimizing the number of sources that contain the hybrid system. In this paper, an energy management and sizing algorithm, applicable to multisource systems, composed of a large number of sources, is proposed. The method is based on a modiﬁed centered moving average ﬁlters architecture for energy management, which permits one to consider and to automatically balance the forecasting errors in solar and load proﬁles. The energy management is then limited to a small number of parameters, which are the averaging horizon and weight coefﬁcients. It is then possible to optimize, at the same time, the sizing and the energy management of such power systems. The proposed optimization criterion is based on a techno-economic approach, by considering acquisition and operation costs, as well as the ageing of the different devices. The main novelty of this approach is the use of energy management formulation that is able to manage an architecture with a high number of controlled devices. An original formulation of centered moving average ﬁlters also permits one to automatically balance the power bias due to forecasting errors on the renewable resources and the load proﬁle. The method is applied to ﬁve devices, including photovoltaic panels, a fuel cell, two batteries with different technologies (Li-ion and lead-acid) and supercapacitors. M.M.,


Introduction
From the 1920s to the 1970s, a new electrical load-feeding procedure was introduced consisting of separating the loads and supplying them with various generators [1]. This technique brought a lot of advantages, especially in terms of construction cost per kW. Since the 1970s, large generating resources such as tidal stream generators have been added to these systems [1]. On the other hand, the benefits of these applications have been progressively reduced because of environmental laws that have been issued by governments and international organizations. All of this has forced companies to motivate their scientists and engineers to look for new solutions.
One of these solutions was the construction of independent power generating units of small sizes. Then, at the end of this period, these units were developed in such a way as to be able to integrate renewable energy sources, which are mainly based on wind and requirements to be respected by the power profile supplied/stored by the source [37]. Additionally, there is a third indicator that combines the first two. It consists of determining the minimum size of the system that respects the technical constraints and then setting the technology that minimizes the financial indicator [40]. Energy management and sizing are two very connected phases and directly affect the optimization of the total cost of the system. Generally, the authors study the optimization, either of energy management with the preliminary setting of the sources sizing, or sizing devices with a simple management strategy [41]. Compared to previous works, in this article, a new energy management strategy is proposed. It is inspired from frequency separation algorithms [19,42], that allow one both to manage energy flows between sources and storage systems and to determine the sizing of each device, to optimize the total system cost.
The proposed strategy has a generic structure and is based on a frequency separation with centered moving average (CMA) filters, which gives an adapted power profile for each device. An issue of such filter structure is the need for an estimation of the power profile. The CMA filter is a special case of FIR filters, where the weight coefficients are all equal to each other. As it is an acausal structure, it permits one to consider future states and past states to generate the power profile of a source or storage device. Obviously, future states are not known but only estimated. It is the reason why some forecasting errors can appear. In case of overly large forecasting errors, deviations on the power profiles at the output of CMA filters can lead to high power and state of charge deviations, that have to be balanced. In this work, we propose an original version of CMA filters, which automatically compensate the forecasting errors to ensure that they keep the expected dynamics of the power profile supply by each device.
Moreover, this structure can easily be applied to hybrid systems composed of a large number of energy sources or storages, with different technologies, which is an important step forward for such multi-source application. With the same idea, it is also proposed here to extend the concept of cascaded filters to reach a more general architecture of control, with series-parallel filters. It permits one to include a large choice of device technologies, to try to reach an optimal control, while allowing redundancy, which finally improves the reliability of the power system. In addition, the energy management is set by a limited number of parameters, which then also allows the optimization of the sizing of the different devices. In this work, we consider the minimization of the total cost of the multisource power system. This criterion permits one to evaluate the effect of several phenomena such as the lifetime of the different devices, as well as acquisition, maintenance and operation costs.
The paper is organized as follows: Section 1 explains the principle of the proposed algorithm. In Section 2, the algorithm will be applied to a multisource system consisting of solar panels as a renewable source, a PEM fuel cell, two batteries of different technologies and a bench of supercapacitors as energy storage systems. The results obtained, in terms of energy management, sources technologies and sizing will be presented in Section 3. Finally, the main conclusions will be given in Section 4.

Proposed Energy Management Approach
Algorithms based on frequency separation and wavelet transform inspire the proposed algorithm. Figure 1 illustrates its general architecture based on a series-parallel architecture of blocks called "device controllers". Each of them is responsible for calculating the reference power (the set point) of one power source. The algorithm is based on centered moving average (CMA) filters [13]. The idea is to consider a predicted profile of the requested power over a defined future time horizon. Based on this predicted profile, the algorithm determines the profile, the sizing and the technology of a power device. By then cascading these CMA filters, it becomes easy to distribute the power profile on all the sources and storage devices of the multi-source system. A classic cascading architecture consists of organizing the components from the slowest (first stages) to the fastest (last stages). From the output of one stage, the residual power is computed and used as input of the next stage (see Figure 1). It is common to find research works using these cascaded Energies 2021, 14, 3627 4 of 21 CMA filters, but considering only one component per stage. A first evolution that is proposed in this paper is to consider the possible connection of several components on one stage. Each component is then associated to a CMA filter with a given time horizon and given weighting coefficients. The only constraint on weighting coefficient is that the sum for one stage is equal to unity: for one stage is equal to unity: With the stage number and the total number of components for this stage. The serialization of the blocks has several advantages. Firstly, it allows the easy management of energy in hybrid systems composed of several and different natures of sources. Secondly, thanks to the notion of residual power, the power demanded by the customer is completely provided. Additionally, it guarantees the non-divergence of the states of charge of the energy storage systems over the operation time. Finally, it allows an automatic compensation behavior of prediction errors that can affect the predicted input of the algorithm when the hybrid system is turned on. These last two findings are mathematically demonstrated in the following sections. In addition, the paralleling of the controllers allows the possibility of managing energy between sources of identical or similar technologies.  With k the stage number and N the total number of components for this stage. The serialization of the blocks has several advantages. Firstly, it allows the easy management of energy in hybrid systems composed of several and different natures of sources. Secondly, thanks to the notion of residual power, the power demanded by the customer is completely provided. Additionally, it guarantees the non-divergence of the states of charge of the energy storage systems over the operation time. Finally, it allows an automatic compensation behavior of prediction errors that can affect the predicted input of the algorithm when the hybrid system is turned on. These last two findings are mathematically demonstrated in the following sections. In addition, the paralleling of the controllers allows the possibility of managing energy between sources of identical or similar technologies. Now that the global architecture of the energy management is presented, it is possible to focus on the structure of a CMA controller. Figure 2 represents the internal structure of such a controller. It consists of four main functions: a CMA filter [13], a power limit supervision function, an energy limit supervision function and a state of charge controller for storage devices [42]. Now that the global architecture of the energy management is presented, it is possible to focus on the structure of a CMA controller. Figure 2 represents the internal structure of such a controller. It consists of four main functions: a CMA filter [13], a power limit supervision function, an energy limit supervision function and a state of charge controller for storage devices [42]. is a weighting coefficient, the input power of the stage , the time horizon of the CMA filter, , are the maximum and minimum power limits of the device , the balance power to compensate the measurement errors, * and are respectively the setting and actual state of charge of the storage device .
The previous architectures correspond to the optimization phase for which the algorithm will need only the future estimate of the requested power to estimate the optimal distribution, as well as the optimal dimensions and technologies of devices. On the other hand, after the installation of the system, a second entry will be added, which represents the data measured in the past of the requested power. This allows for the inclusion of forecast errors that may be present.

Centered Moving Average Filter
The moving average filter is the fundamental function of each device controller [33,42,43]. The output of a filter represents the reference power of the energy source associated with it. When the system is turning on, the reference power is calculated in real time from the data measured in the past of the requested power (measures), noted _ , and predicted future data of the same requested power, noted _ (to highlight the predicted nature of this power profile). The next equation gives a mathematical expression of such a CMA filter [33,34]: with ( ) the output power of the filter at time , the calculation horizon of the filter , a weighting coefficient. This filter can be seen as a convolution between the input power Pink and a centered rectangular window with a total length of Hki. The main interest of such a filter is the small number of degrees of freedom (weighting coefficient and filtering horizon) and the possibility of taking into account, in a balanced way, the power profile in the future and in the past, and then without adding any delay in the generated output power profile. However, it is then obvious that such a filter is acausal and it must be modified to be implemented on real applications. Then, a first modified version consists of separating the history of the measurements and the predictions, such as: The previous architectures correspond to the optimization phase for which the algorithm will need only the future estimate of the requested power to estimate the optimal distribution, as well as the optimal dimensions and technologies of devices. On the other hand, after the installation of the system, a second entry will be added, which represents the data measured in the past of the requested power. This allows for the inclusion of forecast errors that may be present.

Centered Moving Average Filter
The moving average filter is the fundamental function of each device controller [33,42,43]. The output of a filter represents the reference power of the energy source associated with it. When the system is turning on, the reference power is calculated in real time from the data measured in the past of the requested power (measures), noted P in_ki , and predicted future data of the same requested power, notedP in_ki (to highlight the predicted nature of this power profile). The next equation gives a mathematical expression of such a CMA filter [33,34]: with P ki (t n ) the output power of the filter at time t n , H ki the calculation horizon of the filter ki, a ki a weighting coefficient. This filter can be seen as a convolution between the input power P ink and a centered rectangular window with a total length of H ki . The main interest of such a filter is the small number of degrees of freedom (weighting coefficient and filtering horizon) and the possibility of taking into account, in a balanced way, the power profile in the future and in the past, and then without adding any delay in the generated output power profile. However, it is then obvious that such a filter is acausal and it must be modified to be implemented on real applications. Then, a first modified version consists of separating the history of the measurements and the predictions, such as: whereP ink is the predicted input power of the Stage k. This filter version is now causal, but because of possible prediction errors, an additional compensation stage must be added to ensure that long-term errors are balanced. Practically, without such compensation, the power profile supplied by the storage stages could lead to a state of charge deviations, due to a non-zero mean value of the power.
Here, an original formulation of a causal CMA filter is proposed to automatically balance the prediction errors. This formulation is the next one: The form of the first part of the Equation (2) is not obvious, but is justified by the presence of the estimation errors, which have to be balanced. The rest of the section presents and demonstrates the results obtained with this new CMA filter formulation.
Each filter is only set by the width of its filtering horizon H ki . This horizon is linked to the dynamic possibilities of the source. A larger time horizon leads to lower dynamics at the filter output. Then, slow sources are controlled by CMA filters with large time horizons. It is the reason why it is expected that the first stages have time horizons larger than the next one. Moreover, for the controllers of the same level, the calculation horizons of the filters can be different, but should remain in the same range. As a synthesis and in order to explain the operating principle of the algorithm, we present below a simple example of application through which the mathematical findings will be demonstrated.
If we consider P Load (t) as a function which represents the profile of the power that a residential load will consume over a week andP Load (t) as a forecast of the evolution of this profile. Added to that, we also consider P PV (t) as a function that represents the profile of the power that some solar panels will generate over a week, andP PV (t) a forecast of the evolution of this profile. Two other profiles ∆P and∆ P are also created. They respectively represent the power difference to be satisfied by the multi-source system (supply it if ∆P(t) > 0 or store it if ∆P(t) < 0) and a future forecast of this difference, expressed as follows: These two profiles are supposed to be the two inputs, measured and predicted, respectively, of the algorithm. The expected energy source profile is not really defined. It can come from historical data or simulations. In the worst case, if no historical data exists, it can be replaced by the power expectancy. It is also considered that a multisource system is composed of a slow dynamic primary energy source, a first storage system with a moderate dynamic and a second storage system with a fast dynamic. Figure 3 shows the architecture of the corresponding algorithm. At this stage, two situations may exist. The first case considers no prediction errors, which leads to ∆P(t) =∆ P(t). At first, we suppose that our forecasts of the consumed and generated powers (and therefore ∆P) will coincide perfectly with what will be really consumed and generated. Then, if the actual power ∆P is replaced by its estimated value∆ P, the output of the first CMA filter is given by: As the stage 1 is here composed of a single controller, it leads that a 1 = 1. For the second stage, if the actual power ∆P(t) is also replaced by estimated future power∆P(t) and if we consider a weight factor a 2 = 1, the output of the second CMA filter is given by: Developing Equation (4) permits one to express the output of the filter with two different average calculations: To help to understand this result, we can consider the mean value of signals. Indeed, it is obvious that the average value of CMA filter over an arbitrary horizon H s is equal to the average value of the input signal, over the same horizon. Noting x H s the average value of a signal x(t) over a time horizon H s , it can be deduced that the average value of the output y(t) of CMA filter with horizon H i is equal to the average value of the input u(t): where H s >> H i . If we applied this relation to Equations (3) and (4), it is obtained for the first stage: and then, for the second stage: In other words, only the first energy source provides the average value of the profile of the requested power. Thus, P in 2 , the input power of the second stage, will have a zero-average value on H s . These null values over an operating period guarantee the nondivergence of their states of charge. However, for this first demonstration, it is considered that the future power inbalance between renewable sources and power consumption is perfectly predicted. It is obvious that in reality, this power can have deviations from those that will actually be measured. These deviations are called forecasting errors, classically defined as follows: In order to show the effectiveness of the algorithm in dealing with these forecasting errors, it is assumed that our future predictions of P PV and P Load are affected by errors compared to the powers that will actually be measured. Considering Equation (2), the calculation of the new reference powers gives the following results: With the formalism given by Equation (9), it is obtained: If it is assumed that the error ε has stationary properties at the CMA horizon scale and over, it results that: The average value for any time horizon H s , longer than H 1 , is then: Finally, the average value of the output power calculated by the first stage is equal to the average value of the measured requested power: It results that all the other stages have then a null average value, as for the case without forecasting error. What may be different are the actual power variations for each stage, which can lead to an over-or under-sizing of the device. However, this dynamic forecasting error is distributed among the different stages of the hybrid system while respecting the dynamics of each of them. In summary, even in the presence of forecasting errors, since it has a certain degree of stationarity, the average value of the reference powers of the different energy stages remained zero. This is essential, since it allows potential storage devices to have no divergence of their state of charge. In accordance, their presence modifies the magnitude of the powers supplied or stored. This must be taken into account in the sizing phase. It should be noted that all the results and previous findings remain valid regardless of the number of parallel controllers in a stage. This has been verified.
The average value for any time horizon , longer than , is then: Finally, the average value of the output power calculated by the first stage is equal to the average value of the measured requested power: It results that all the other stages have then a null average value, as for the case without forecasting error. What may be different are the actual power variations for each stage, which can lead to an over-or under-sizing of the device. However, this dynamic forecasting error is distributed among the different stages of the hybrid system while respecting the dynamics of each of them. In summary, even in the presence of forecasting errors, since it has a certain degree of stationarity, the average value of the reference powers of the different energy stages remained zero. This is essential, since it allows potential storage devices to have no divergence of their state of charge. In accordance, their presence modifies the magnitude of the powers supplied or stored. This must be taken into account in the sizing phase. It should be noted that all the results and previous findings remain valid regardless of the number of parallel controllers in a stage. This has been verified.

Additionnal Control Loops and Supervision
The power distribution is not affected by the state of the different devices. Due to unbalance and inaccuracy in the state measurement and control, it is necessary to maintain the state of charge (SoC) of storage devices toward a reference value. Then, taking into account the dynamics of the energy storage systems, a proportional control of the SOC using CMA filters is used, as illustrated in Figure 2. The SoC can be easily obtained through charge counting or open circuit voltage methods. The reference power sent to the source is obtained by summing the output power of the filter and a part of the deviation of SOC from its setting point. The idea is classical in energy management with frequency separation and considers that the average SOC value of a stage is ensure by the previous (and slower) stage − 1. The equations that define this control are as follows: and the power correction is obtained with: Figure 3. Architecture of the algorithm corresponding to the hybrid system taken as an example.∆P and ∆P are the expected and actual unbalance power respectively,P out k and P out k are the expected and actual output power profile of the stage k respectively,P in k and P in k are the expected and actual input power profile of the stage k respectively, H k are the time horizon of the CMA filter k.

Additionnal Control Loops and Supervision
The power distribution is not affected by the state of the different devices. Due to unbalance and inaccuracy in the state measurement and control, it is necessary to maintain the state of charge (SoC) of storage devices toward a reference value. Then, taking into account the dynamics of the energy storage systems, a proportional control of the SOC using CMA filters is used, as illustrated in Figure 2. The SoC can be easily obtained through charge counting or open circuit voltage methods. The reference power sent to the source is obtained by summing the output power of the filter and a part of the deviation of SOC from its setting point. The idea is classical in energy management with frequency separation and considers that the average SOC value of a stage k is ensure by the previous (and slower) stage k − 1. The equations that define this control are as follows: (15) and the power correction is obtained with: where H k the calculation horizon of the filter, SOC k+1 a history of the SOC of the source whose SOC is to be regulated,ŜOC k+1 a prediction of the future state of charge, K corr is the proportional gain and SOC * k+1 the reference of the state of charge. The proportional Energies 2021, 14, 3627 9 of 21 gain K corr is calculated according to the technical characteristics of the source whose SOC is to be regulated. It should be noted that in the actual system, the SOC is a measured parameter, and it is not necessary to to calculate it. However, by simulation, these states are just estimated from a charge counting method. Thanks to the Formulation (10), the presence of the forecasting errors generates an automatic rise of the reference powers in order to compensate them. This increase may result in a reference power that is greater than the power limits of the source and therefore may lead to equipment damage or even danger. In order to remedy this problem, the output of each filter is bound by the maximum and minimum limits of the source by saturation bloc, as shown in Figure 2. On the other hand, the saturation of the filter output signal causes a change in its average value over the operating horizon. However, as explained in the previous section, the hypothesis given by Equation (6) is fundamental in achieving the results and findings that came after. Indeed, it can be assumed that the average power over the operating horizon of the filter output is the actual one. Any inequality between these two averages leads to a divergence of the SOCs of the energy storage systems and consequently the instability of the system.
In order to remedy this problem, a block is added just at the saturation output, which detects the saturations if they exist, accumulates the suppressed powers and adds them, if the accumulated power is positive, or subtracts them, if the accumulated power is negative, of the reference power of the source as soon as possible. In this way, the total energy supplied on the operating horizon remains unchanged, even in the presence of saturation, and consequently, the average value on this same horizon always remains equal to that of the measured input of the filter.
Finally, the moving average filter can cause a power profile on a source involving the exceeding of energy limits. To solve this problem, an energy supervision loop has been added. At each time step, and before sending the command to the source, a saturation is applied on the power, to avoid energy overruns.

Optimization Process
In this study, the considered application which illustrates this energy management strategy is a habitat, with photovoltaic panels already installed and whose owners want to be disconnected from the grid to switch to self-consumption. Thus, the sizes of the renewable source and load profiles are input parameters that are initially imposed. At this stage, by knowing the size of the renewable source and the geographical position of the studied habitat, a profile of the power that will be produced by this source over a year is estimated. Additionally, considering technical and sociological criteria (size of the habitat, number of people living there, loads, etc.), a profile of the power that will be consumed over a year is estimated.
From the PV production and load consumption estimations, an estimated unbalance power profile∆P is created. The actual profile ∆P is then built as and when. Based on the dynamics and magnitude of the profile∆P, architecture of a multisource system should be defined. This definition consists of specifying the types of sources or energy storage systems, there number and the parameters of the energy management (i.e., the power profile distribution according to the dynamics of each device). This joint sizing energy management optimization is made possible thanks to the very simple form of energy management based on CMA filters. There are only two additional optimization variables per device. Then, based on a predicted load and solar profile unbalance∆P and using an optimization algorithm, all the pairs (a ki ; H ki ) associated with the different CMA controllers of the series-parallel structure (see Figure 1) are optimized at the same time as the sizing. To help to understand the role of these filter parameters, Figure 4 illustrates their effects on a Ragone diagram. It is done for a given power profile and for different values of weight coefficient (a ki ) and time horizon (H ki ). This figure shows that if the horizon H ki becomes larger, the energy of the profile increases and its maximum power decreases. In other words, the reference power of slow dynamics sources (also known as energy sources) must use controllers with CMA filters operating with wide horizons. As for a ki , it allows one to adjust the power and energy demand to the size of the device. Indeed, the higher the coefficient, the larger the source size. It can be seen as a parameter that reflects the size of the source, unlike H ki , which rather reflects its dynamic.
sizing. To help to understand the role of these filter parameters, Figure 4 illustrates their effects on a Ragone diagram. It is done for a given power profile and for different values of weight coefficient ( ) and time horizon ( ). This figure shows that if the horizon becomes larger, the energy of the profile increases and its maximum power decreases. In other words, the reference power of slow dynamics sources (also known as energy sources) must use controllers with CMA filters operating with wide horizons. As for , it allows one to adjust the power and energy demand to the size of the device. Indeed, the higher the coefficient, the larger the source size. It can be seen as a parameter that reflects the size of the source, unlike , which rather reflects its dynamic. Here, the objective function is the total cost of the power system. It combines the total investment, operation, replacement (which depends on the lifetime of the components and therefore their power profiles) and maintenance costs: Each cost of Equation (17) depends on many assumptions and on the techno-economic hypothesis. For instance, the replacement cost of a device depends on its health degradation rate, which needs a lifetime function, linked to the operation modes of the device (number of full-charge and discharge cycles for a battery, total supplied energy for the fuel cell, etc.). A full description of these costs and the degradation functions for each device are presented in Section 3.
As explained, a time horizon value depends on the capabilities of the connected device. For the first stage controllers (1 ), in order to obtain a smooth profile, the horizons * must be large enough. It is chosen to be around at least a few days. The horizons * are in the order of minutes. This last horizon is chosen to ensure a balance between the characteristics of the connected device and the residual power profile which must be adapted to sources with low specific energy. Then, the following constraints could be defined: Here, the objective function is the total cost of the power system. It combines the total investment, operation, replacement (which depends on the lifetime of the components and therefore their power profiles) and maintenance costs: Each cost of Equation (17) depends on many assumptions and on the techno-economic hypothesis. For instance, the replacement cost of a device depends on its health degradation rate, which needs a lifetime function, linked to the operation modes of the device (number of full-charge and discharge cycles for a battery, total supplied energy for the fuel cell, etc.). A full description of these costs and the degradation functions for each device are presented in Section 3.
As explained, a time horizon value depends on the capabilities of the connected device. For the first stage controllers (1i), in order to obtain a smooth profile, the horizons H * 1i must be large enough. It is chosen to be around at least a few days. The horizons H * 3i are in the order of minutes. This last horizon is chosen to ensure a balance between the characteristics of the connected device and the residual power profile which must be adapted to sources with low specific energy. Then, the following constraints could be defined: The optimization algorithm tries to find the right combination of pairs (a ki ; H ki ), which, at the same time, respects the previous constraints and optimizes the predefined objective function. The optimization steps are described below. For each pair (a ki ; H ki ) of a tested combination, the following steps are performed: Calculation of the corresponding power profile according to Equation (2). b.
Calculation of the maximum energy and the maximum power of this profile. c.
Based on the results of the step b and on a database containing specific powers and specific energies of the various technologies whose source ki can be built with a power sizing and an energy sizing of this source, these are calculated for each technology. These dimensions represent the minimum sizes to be respected for the source ki to be able to supply/store the energy and power calculated in b. For each technology, the largest dimensioning (between that in power and that in energy) is considered. d.
Based on a database containing approximate lifetimes (in hours, in complete cycles, etc.) of the different technologies whose source can be built with, and on the power profile calculated in step a, an approximate lifetime of the source is calculated for each possible technology. e.
Based on the results of step c and on a database containing the provided Wh price for the different technologies whose source ki can be built with an initial investment cost of the source is calculated for the different possible technologies. f.
Using the lifetimes calculated in step d, a replacement cost is assigned to each technology. g.
A maintenance cost is associated with each technology. h.
Based on lifetimes calculated in d and costs e, f and g, the total cost (investment + replacement + maintenance) over the lifetime of the multisource system, of each technology that the source ki can be built with, is calculated. The power system lifetime is assumed to be equal to the lifetime of the renewable source. i.
The technology that corresponds to the lower cost is chosen as the optimal technology of the source ki. We go back to step c to know the minimum dimension that must be considered for this technology.
These steps are done for all couples at each iteration. In the end, the algorithm gives an optimal solution for each device, through the optimal couple (a ki ; H ki ) which corresponds to an optimal power profile calculated by Equation (2).

Hybrid Power System Description
The REDD (Reference Energy Disaggregation Data Set) provides a database containing measurements of the total power consumed by six different houses [44]. The data is recorded at a sample time of 1 second. It is the power profile of the House 1, that is used in this study. The load profile includes the consumption of lighting, an oven, a refrigerator, a microwave, a stove, a dishwasher, a washer dryer and kitchen outlets. Figure 5 shows the daily pace of this profile. We consider this profile as the forecast of the total consumption of the house over the coming year,P Load .
We also assume that this house is located in Saint-Nazaire, France, and equipped with 5 m 2 of solar panels. As for the estimated photovoltaic power profile, an estimated solar radiation profileĜ is generated by a solar profile generator [41]. This generator uses a Markov matrix to simulate the cloud cover. Then, from the inclination angle of the PV panels, the geographical coordinates of the studied habitat (longitude and latitude) and a sample step, a solar radiation profile can be generated. Figure 6 shows a daily pace of the profileĜ. Table 1 summarizes the inputs used for the profile used in this study.  We also assume that this house is located in Saint-Nazaire, France, and equipped with 5 m 2 of solar panels. As for the estimated photovoltaic power profile, an estimated solar radiation profile is generated by a solar profile generator [41]. This generator uses a Markov matrix to simulate the cloud cover. Then, from the inclination angle of the PV panels, the geographical coordinates of the studied habitat (longitude and latitude) and a sample step, a solar radiation profile can be generated. Figure 6 shows a daily pace of the profile . Table 1 summarizes the inputs used for the profile used in this study.    We also assume that this house is located in Saint-Nazaire, France, and equipped with 5 m 2 of solar panels. As for the estimated photovoltaic power profile, an estimated solar radiation profile is generated by a solar profile generator [41]. This generator uses a Markov matrix to simulate the cloud cover. Then, from the inclination angle of the PV panels, the geographical coordinates of the studied habitat (longitude and latitude) and a sample step, a solar radiation profile can be generated. Figure 6 shows a daily pace of the profile . Table 1 summarizes the inputs used for the profile used in this study.   The photovoltaic panels available to the consumer are supposedly made of polycrystalline silicon. The efficiency of this material is estimated at 10% and assumed to be constant throughout the operating period. The power produced by the panels is proportional to the solar radiation and is calculated according to Equation (18): with S PV the surface of the photovoltaic panels, η PV their efficiency. From the two previous profiles, an estimated profile∆P is created. It represents the unbalance power that the multisource system is expected to store and/or provide in the coming year. Figure 7 shows one day of this power profile. On the other hand, differences will most likely be found between the actual profile ∆P that will be measured over the coming year and the predicted profile∆P. These differences will represent the forecasting errors. These errors come from a bad prediction of either the load demand or the production of solar panels. As shown in the same figure, a profile ∆P is defined, which is supposed to be the unbalanced power that will be actually measured over the coming year. In this work, the estimates of the unbalance power ∆P are assumed to be 20% lower than the actual profile.
ous profiles, an estimated profile Δ is created. It represents the unbalance power that the multisource system is expected to store and/or provide in the coming year. Figure 7 shows one day of this power profile. On the other hand, differences will most likely be found between the actual profile Δ that will be measured over the coming year and the predicted profile Δ . These differences will represent the forecasting errors. These errors come from a bad prediction of either the load demand or the production of solar panels. As shown in the same figure, a profile Δ is defined, which is supposed to be the unbalanced power that will be actually measured over the coming year. In this work, the estimates of the unbalance power ΔP are assumed to be 20% lower than the actual profile. Figure 7. Estimate of a daily power to be managed by the hybrid system and power that will actually be managed. Now, based on the estimated profile Δ , the architecture of a multisource system must be defined. Indeed, the power profile Δ has power peaks of short durations (at the scale of the second) as between 2:00 a.m. and 3:00 a.m. in Figure 7. To meet these sudden and fast demands, we need a fast dynamic power source. We choose supercapacitors (SC), which is a fast device, with high efficiency in short time energy transfers [45]. According to Section 2.1, the average of the measured input profile Δ will be provided by one (or more) source(s) controlled by the first stage controllers. In other words, it must have a slow dynamic. If the forecasts coincide perfectly with the profile Δ that will be actually measured, this source will supply the annual average power of Δ (equal to 137 W in this case). Given that the value of the power is not very high, a fuel cell (FC) is chosen as an auxiliary source. This choice involves, for the unbalanced power ∆P, a positive average Figure 7. Estimate of a daily power to be managed by the hybrid system and power that will actually be managed. Now, based on the estimated profile∆P, the architecture of a multisource system must be defined. Indeed, the power profile∆P has power peaks of short durations (at the scale of the second) as between 2:00 a.m. and 3:00 a.m. in Figure 7. To meet these sudden and fast demands, we need a fast dynamic power source. We choose supercapacitors (SC), which is a fast device, with high efficiency in short time energy transfers [45]. According to Section 2.1, the average of the measured input profile ∆P will be provided by one (or more) source(s) controlled by the first stage controllers. In other words, it must have a slow dynamic. If the forecasts coincide perfectly with the profile ∆P that will be actually measured, this source will supply the annual average power of∆P (equal to 137 W in this case). Given that the value of the power is not very high, a fuel cell (FC) is chosen as an auxiliary source. This choice involves, for the unbalanced power ∆P, a positive average value, to be compatible with the behavior of a fuel cell (the power cannot be negative). Indeed, if the annual average of∆P is negative, it is necessary to consider another solution, such as the addition of an electrolyzer.
For the second and third stage controllers, while the average power is supplied by the first stage (a fuel cell in this application), the average power of the other stages is null and then the corresponding devices can be storage solutions. Here, it is considered to use two different dynamic batteries, one for stage 2 (the slowest) and one for stage 3 (the fastest). Figure 8 shows the architecture of this energy management. In each stage, there is only one controller, so that: a * 11 = a * 21 = a * 31 = a * 41 = 1 Finally, particular swarm optimization algorithm (PSO) is used to determine the optimum energy distribution through the determination of the optimal values of horizons H 11 , H 21 and H 31 and to obtain the optimal size and technology of each source. The objective function considered in this study is presented in the following section.
the first stage (a fuel cell in this application), the average power of the other stages is null and then the corresponding devices can be storage solutions. Here, it is considered to use two different dynamic batteries, one for stage 2 (the slowest) and one for stage 3 (the fastest). Figure 8 shows the architecture of this energy management. In each stage, there is only one controller, so that: * = * = * = * = 1 Figure 8. Architecture of the algorithm corresponding to the defined hybrid system. Finally, particular swarm optimization algorithm (PSO) is used to determine the optimum energy distribution through the determination of the optimal values of horizons , and and to obtain the optimal size and technology of each source. The objective function considered in this study is presented in the following section.

Obective of Optimization
The optimal values of the calculation horizons , and of the filters of controllers 11, 21 and 31, respectively, must correspond to:

−
An optimal distribution of energy between the fuel cell, the two batteries and the supercapacitors. The optimal reference power should respect the dynamics, the power and the energy limits of each source. − An optimal size and an optimal technology of each source. − A minimum power system cost.
Here, the optimization is based on the total cost of the power system. It includes the capital, operation, replacement and maintenance costs. As the lifetime of the different devices is not the same, the annualized cost is considered in this work. To estimate the lifetime of the fuel cell, a degradation function Δ is considered. It only depends on the supplied energy over a simulation horizon [12]: where ∂( ) = ∂ 3600 1 + α ( * ( ) − )² with ∂ and α two empirical coefficients and (in W) the nominal power of the fuel cell (80% of ). The lifetime is then calculated as follows:

= ∆ ( ) [in years]
The consumption of hydrogen is linked to the efficiency and to the power supplied by the fuel cell. The volume _ ( ) of the annualized hydrogen consumption is calculated as follows [13]:

Obective of Optimization
The optimal values of the calculation horizons H 11 , H 21 and H 31 of the filters of controllers 11, 21 and 31, respectively, must correspond to: − An optimal distribution of energy between the fuel cell, the two batteries and the supercapacitors. The optimal reference power should respect the dynamics, the power and the energy limits of each source. − An optimal size and an optimal technology of each source. − A minimum power system cost.
Here, the optimization is based on the total cost of the power system. It includes the capital, operation, replacement and maintenance costs. As the lifetime of the different devices is not the same, the annualized cost is considered in this work. To estimate the lifetime of the fuel cell, a degradation function ∆ FC is considered. It only depends on the supplied energy over a simulation horizon H S [12]: where with ∂ 0 and α two empirical coefficients and P FC nom (in W) the nominal power of the fuel cell (80% of P FC max ). The lifetime is then calculated as follows: The consumption of hydrogen is linked to the efficiency and to the power supplied by the fuel cell. The volume V H2_ann (m 3 ) of the annualized hydrogen consumption is calculated as follows [13]: For electrochemical batteries, the health degradation over the time is noted ∆ Bat . It depends on the energy exchanged by the battery. The simplest function is a parabolic relationship between the health degradation and the depth of charge and discharge [23,40,41]. In this case, the health degradation is proportional to the exchanged energy. With this assumption, ∆ Bat is given by [9]: With P Exchanged the power supplied by the battery during the simulation horizon H s and E Tot the maximum energy exchange of the battery. This energy is deduced from the energy capability and the maximum number of cycles. The end of life is finally obtained when the degradation function ∆ Bat reaches unity. The estimated lifetime of the battery can then be calculated: To finish, the degradation function of the supercapacitor bank is noted ∆ SC . This function depends on the state of charge of the supercapacitors (SoC sc ) [14,15]: where Additionally, with µ SC nominal aging of SC per year [%/s] (1.5%/year in this study). Thus, the estimated lifetime is deduced as follows: The renewable sources, such as wind turbine or photovoltaic panels, generally have the longest lifetime among the other devices. For photovoltaic panels, the expected lifetime is 25 years. It is this duration which is considered for the multi-sources power system.
The previous health degradation function permits one to estimate the lifetime of each device [9]. Then, the annualized costs can be calculated from equations summed up in Table 2. The PV panels need additional installation and acquisition costs. These costs are 40% of the acquisition cost [41,46]. Finally, the system cost is given by (22) (22) In addition to the time horizon range for the optimization of the energy management and sizing, the storage devices must respect the following constraints: where E Bat and E SC represent the storage energy of the battery and the supercapacitors, respectively. Supercapacitors can theoretically support large charges and discharge limits.
Here, the discharge threshold is imposed to 5%, to not generate a too high current of the supercapacitor. Table 2. Annualized costs of the multisource power system [40,46].

Optimisation Results
Based on the estimated profilesP Load andP pv mentioned in the previous section, the optimization results are given in Table 3. Firstly, it is assumed that the real-time measurements made during the coming year, of the photovoltaic power as well as that demanded by the load, will coincide perfectly with what has been predicted (without forecasting errors). The obtained sizing is in accordance with the assumptions done on the expected horizon values for the different devices. For the fuel cell, the filtering horizon is of the order of several days, to obtain a very smooth power profile, which permits the fuel cell to not have to supply high power spikes and to have then a smaller and cheaper sizing. For the batteries, the main one (battery 1, with lead-acid technology) ensures the range of dynamics of the order of the day, while the second (battery 2, with LiFePO4 technology) is used to fill the gap between the battery 1 and the residual power supplied by the supercapacitors. It is not easy to compare these sizing and dynamic results because the use of 4 different sources in the same system. Nevertheless, in [47,48], similar dynamics are obtained for the storage devices (battery and supercapacitors), as for the power-energy ratio of these devices. In [47], an equivalent system is used, but with only the first battery and the supercapacitors. The time horizon to share the power profile between the battery and the supercapacitors is a few minutes, which leads to an equivalent sizing result than in this work. Figure 9 shows daily patterns of the supplied/stored power of the different energy sources. Firstly, it shows the evolution of the power supplied by the fuel cell over one year. Each point of this profile is the result of the sum of about 5 days of ∆P in the past and 5 days points of the∆P, to respect the optimal horizon H * 11 , which is equal to 10 days. Thanks to the large horizon H * 11 , this profile is characterized by a very slow dynamic and a total absence of sudden changes or power peaks. The maximum power of this profile is lower than 200 W, which goes with the sizing mentioned in Table 3. The second and third plots represent, respectively, the evolutions of the reference powers of the second and third stages, which control, respectively, the lead acid and LiFePO4. It is obvious that the dynamic of the power profile for theLiFePO4is higher than those for lead-acid battery. This difference in dynamics is explained by the difference of technologies, result of the optimization. The last plot of Figure 9 illustrates the daily evolution of the residual power, ensured by the supercapacitors. A zoom is made on an hour. Power peaks characterize this profile over short periods.  Figure 9. Daily reference power profiles with (blue) and without (black) forecasting errors.
with predicted profile with actual profile Figure 9. Daily reference power profiles with (blue) and without (black) forecasting errors. Figure 10 shows the annual evolution of the state of charge of the various storage elements. It can be seen that the SoC of lead acid and Li-ion batteries respect the maximum and minimum constraints (22) and (23). On the other hand, the annual evolution of the SOC of the supercapacitor bank, where forecasting errors (FE) are considered, shows over charges. It is due to forecasting errors that lead to a larger power profile for this device. Concerning the influence of forecasting errors, it can be seen that the algorithm automatically modifies the reference powers to compensate it. Indeed, when the predictions are different than the reality, the controllers increased or decreased the power profile until reaching the balance with the actual ∆P power. This process is carried out on all four devices, while always respecting the dynamics of the energy sources.
Nevertheless, the sizing mentioned in Table 3 may no longer be sufficient if the forecasting errors are too large. For this reason, the optimized sizing represents a no margin result. To improve the reliability of the system, oversizing can be considered. In this case, the total costs of the energy sources will increase and therefore the cost of the provided energy. This leads us to say that to have the costs shown in Table 2, it is necessary to optimize our predictions: the better we predict, the closer we are to the minimum cost.    Figure 11 shows a representation on the Ragone diagram of the different power profiles, with and without Forecasting errors. It shows that a relatively high-power with a low energy characterizes the reference profile dedicated to supercapacitors. Conversely, the profile provided by the fuel cell has a relatively low power, while its energy is high. Both batteries are in the middle. The power profile of the LiFePO4 battery is of higher power and energy than the lead acid battery profile.

Conclusions
In this paper, an algorithm for the energy management and sizing of multisource systems was presented. Its particularity is in its generic aspect that allows it to manage energy in multisource systems composed of a large number of sources of different natures To easily analyze these results of this energy management, iso-time lines are added on the figure. It permits one to demonstrate that, even in presence of forecasting errors, the specific energy-power ratio is always respected. It is then clear that the CMA filters approach makes it possible to achieve efficient separation in the dynamics, with good compliance with the specific characteristics of each device.

Conclusions
In this paper, an algorithm for the energy management and sizing of multisource systems was presented. Its particularity is in its generic aspect that allows it to manage energy in multisource systems composed of a large number of sources of different natures and technologies, without resorting to complex resolutions. Its architecture is a series-parallel arrangement of CMA filters. Based on the prediction of the load profile and the renewable power, the output of each filter represents the reference power of a specific energy source. The dynamics as well as the order of magnitude of this power depend, respectively on the calculation horizon of the filter and its weighting coefficient. A mathematical approach has shown an ability of the algorithm to automatically compensate for prediction errors. In addition to energy management, optimal sizing and technology choice are made for each source of the multisource system. To do this, the algorithm uses a database containing technical-economic characteristics of the most widely used technologies of sources and energy storage systems on the market. An application was made on a multisource system feeding an isolated habitat. This system consists of a fuel cell and three storage systems. An optimization was launched to know the optimal distribution of energy as well as the nature and the optimal technology of the devices. The results obtained showed the effectiveness of the proposed algorithm to respect the different power and energy constraints of the sources and to behave in the face of forecast errors.