Rémy Rigo-Mariani
email: remy.rigo-mariani@g2elab.inp-grenoble.fr

Sean Ooi
email: seanooi@ntu.edu.sg

Chea Wae

Stefano Mazzoni
email: mstefano@ntu.edu.sg

Impact of the Economic Environment Modelling for the Optimal Design of a Multi-Energy Microgrid

Keywords: Multi-Energy, Optimal Planning, Stochastic Optimization, Linear Programming Introduction I. NOMENCLATURE

The paper discusses the impact of the problem formulation for the optimal design of a multi-energy microgrid.

The study shows that the energy rates affect the output results much more than any other parameter such as the cost of equipment or the assets efficiencies (ratio 1 to 100). The main focus of the paper is then the representation of the economic environment for such a planning problems. Five different modeling approaches are investigated with both deterministic and stochastic methods (scenario-based formulation), as well as constant or increasing rates along the planning horizon. The obtained results show great impacts on the installed capacities (variations from simple to triple for the gas engine rated power), while the objective function (i.e. system cost of ownership) displays small variations less than 5 % with the different sets of hypothesis. The papers then concludes on a recommendation to present the results in terms of optimal areas instead of a single global optimum for such problems.

Multi energy systems have long been proved to enhance the integration of renewable energy sources by providing efficiency and flexibility with the interaction of different energy vectors [START_REF] Jana | Polygeneration as a future sustainable energy solution -A comprehensive review[END_REF]. The optimal design of so-called "multienergy microgrids" has been extensively addressed in the literature. Such problems consist in finding the best capacities of the assets (solar generator, co-generation plant, electrical/thermal storages, etc) with a tradeoff between the capital expenditures and the operating cost estimated along the system lifetime [START_REF] Gamarra | Computational optimization techniques applied to microgrids planning: A review[END_REF]. Such decision processes are subject to a wide range of uncertainties that has to be taken into consideration in the design phaseuncertainties regarding the load level, the renewable generation, the energy prices, the investment costs as well as the models accuracy [START_REF] Mavromatidisa | Uncertainty and global sensitivity analysis for the optimal design of distributed energy systems[END_REF]. Considering those uncertainties may lead to prohibitive computational times is usually tackled with different approaches mixed with simplified linear models of the systems. Typical methods rely on robust optimization with the definition of best/worst case scenarios that represent the volatility of uncertain parameters [START_REF] Gazijahani | Robust design of microgrids with reconfigurable topology under severe uncertainty[END_REF]. Other methods introduce chance constraints while setting a tolerance on the probability to meet specific operating conditions (typically power balance) [START_REF] Zachar | Microgrid/macrogrid energy exchange: a novel market structure and stochastic scheduling[END_REF]. Finally, scenario-based stochastic optimization is the most widely used approach. It define sets of potential scenario (with uniform or normal distribution to represent the uncertainties). The objective function is then formulated as the sum of the optimal results for each scenario weighted by the corresponding probability of occurrence [6] [7]. It is important to note that the majority of those studies consider the uncertainties for the load profiles and renewablebased generation. This paper investigates the relevance of the uncertainties consideration on a generic benchmark and with a design problem that is solved using simple linear programming. Is the system design significantly impacted in terms of installed capacities or objective function (typically cost)? Do those impacts justify the degree of complexity and associated computational time of advanced stochastic optimization techniques? After the problem is presented in Section III, a sensitivity analysis identifies the energy costs (i.e. electricity and gas here) as the most influent parameters. Section IV then investigates different representations of the economic environment before the planning problem is run for different hypothesis and the conclusions are drawn in Section V.

III. CASE STUDY AND MODELING A. Considered System

The work presented in this paper has been led in the framework of the Smart Multi Energy System (SMES) project, whose objective is the implementation of an energy management strategy at a building scale or for a cluster of buildings for the co-operation of electrical and cooling systems. Indeed, this approach facilitates two key transitions necessary for achieving reduced carbon emissions in the long term: 1) it enhances the integration of renewable energy sources and 2) it increases the overall system efficiency by exploiting the interaction between different energy vectors. Figure 1 displays the single line diagram of the considered system and highlights the interface of both thermal and electrical networks at the chiller plant level and with the combination gas engine / absorption chiller. Both thermal and electrical load profiles corresponds to an office/research building located in Singapore with more than two years recorded data in addition to solar radiation measurements [START_REF] Rigo-Mariani | Comparison of optimization frameworks for the design of a multienergy microgrid[END_REF]. For the sake of simplicity and in the absence of seasonality effects on the loads and solar generation, the system operation is simulated along a single representative week. It should be noted that the representation of the simulation period (i.e. system lifetime) is an area of research by itself with the appropriate techniques ( [START_REF] Pfenninger | Dealing with multiple decades of hourly wind and PV time series in energy models: A comparison of methods to reduce time resolution and the planning implications of inter-annual variability[END_REF]) (e.g. downsampling, clustering) which are not in the scope of the proposed work.

B. Linear Modeling for Planning Study 1) Gas Engine Model:

Specific attention is paid to the model of the gas engine depending on the operating point and the rated power. As different engine sizes and operating conditions will be investigated in the course of the optimal planning, there is a need for simple/fast models in order to avoid prohibitive computational times. At first, a quadratic formulation is considered and highlights the strong impact of the operating point and rated capacity on the engine efficiency (Figure 2a) [START_REF] Mazzoni | A Simulation tool for concentrated solar power based on micro gas turbine engine[END_REF]. With the considered model, the gas engine consumption M t ge (in kg/kWh) can then be computed with the fuel lower heating value (LHV = 48,000 kJ/kg) for different sizes P geR and output power P t ge . This model is then furtherly used to implement a linear formulation of the fuel consumption with the operating point and rated capacity of the engine. At first, a least square minimization is applied in order to fit a set of consumption curves for different classes of engines (Figure 2b) with an overall error of 3.5 %. Finally, that curve fitting allows to express the consumption linearly with the engine power and size by utilizing the coefficient denoted α ge and β ge ((1)).

2) 'All in One' Problem Formulation -OPT0

As previously mentioned, the scope of the work presented here is to investigate the impact of the input parameters and problem formulation on the obtained results. The optimal planning of the system is formulated using linear programming in order to avoid long computational times and allow to test several formulations. Typical approaches encountered in the literature tackle the planning problem with bi-level architectures. An inner loop optimizes the system operation while an outer loop investigates different configurations by using evolutionary methods such as genetic algorithm in [START_REF] Li | Microgrid sizing with combined evolutionary algorithm and MILP unit commitment[END_REF] or particle swarm optimization in [START_REF] Moradi | A hybrid method for simultaneous optimization of DG capacity and operational strategy in microgrids utilizing renewable energy resources[END_REF]. Alternative methods take advantage of the convexity of the inner problem (i.e. the system operation) and integrate it in the outer loop using its equivalent KKT conditions [START_REF] Quashie | Optimal Planning of Urban Microgrids with an Energy Management System[END_REF]. The approach considered here lies on a 'All in One' method where both operating and sizing parameters are variables of a single optimization problem. The objective is to minimize the annualized Total Cost of Ownership (TCO). That cost is calculated with the contribution of the capital expenditures (CAPEX) computed with the cost of the assets and considering the system lifetime (Ny = 20 years). C geR ($/kW) is cost of co-generation unit (i.e. gas engine and absorption chiller) that depends on the engine rated power, C pvR is the cost for the installed solar generation, E batR and E thR are the cost of electrical and thermal storages in S/kWh and $/RTh respectively. The operating cost (OPEX) is estimated as the system energy bill, with the charge for purchasing electricity and gas as well as the penalty for the peak power imported from the main grid. The OPEX is computed over a representative week with electrical/thermal load profiles at a resolution of Δt = 30 min. That OPEX is then multiplied by the number of weeks within a year (Nw = 52.18) to finally estimate the TCO over a year ((2)). A full list of used symbols is given in Section I.

 

:

                         
Equations ( 3) and ( 4) define the bounds/constraints for the operating and the sizing variables respectively. For the battery and gas engine, the operations are limited by the rated capacities. Note that, contrary to the battery, the design of the thermal storage dissociates power and energy with two different variables E tsR and P tsR . The solar power generated is computed by multiplying the value of the installed capacity P geR by a normalized representative weekly profile P t pvN .

Also note that an overhead is specified to limit the power imported from the main grid at 50 % more than the base electrical load (for electrical devices except the chiller plant). That value allows the problem feasibility with the chiller plant supplied by the grid in case there is no battery or gas engine.

                                 0 kW
E P E                
Equation ( 5) reflects the typical constraints for energy storage operation with the state of charge (SOC) limit and energy conservation over the time horizon, whereby the final value at the end of the simulated period is equal to the value at the start. However, compared to typical formulation, each side of the constraint is multiplied by the battery capacity in order to avoid non-linearity (the capacity is also a variable of the considered 'all in one' problem). Note that a similar set of constraints has to be considered for the thermal storage. Equation ( 6) computes the cooling generation for both the chiller and absorption chiller based on their coefficient of performances and the conversion ratio between RT and kW. The cooling power generation from the absorption chiller is also associated to the electrical consumption of the pumps computed with a fixed coefficient (β pcr abch ).

                                                               tR   / 3.
Q P Q            
Finally, similar to every energy dispatch problem, equation [START_REF] Narayan | Risk-averse stochastic programming approach for microgrid planning under uncertainties[END_REF] ensures the balance between supply and demand for both thermal and electrical energy vectors.

P Q Q Q Q Q                    C.

Preliminary Sentsitivity Analysis

The planning problem previously defined is formulated in MATLAB using YALMIP and solved with CPLEX 12.8.1 (16 threads in parallel, 32 GB RAM, 2.7Hz processor). The implemented models are then used to perform a simple sensitivity analysis in order to estimate the TCO variations with different input parameters. The sizes of the equipment remain constant with values at the middle of the range specified in (4). The impacts of energy prices, the cost of the different assets and their efficiencies are then studied using a method similar to a Taylor development around the average point for every parameter (in the range specified in Table I). The electricity prices are based on weekly data for real time prices in Singapore and vary between 85 and 145 % of the original values in the course of the sensitivity analysis. The price for the gas is considered constant under the entire simulated period. Figure 3 displays the results obtained when varying the specified parameters from their middle point (i.e. 0.5 on the x-axis). Gas and electricity prices are the criteria that have the most significant impact on the TCO computed (Figure 3a) as also observed in [START_REF] Mavromatidisa | Uncertainty and global sensitivity analysis for the optimal design of distributed energy systems[END_REF]. This is to be expected since the OPEX represents around 90 % of the estimated TCO in the investigated scenarios with 10 % corresponding to the annualized CAPEX. That ratio is reflected when comparing Figure 3a&3b with the effect of the equipment cost ten times smaller than the impact of energy prices (magnitude of 1 % compared to magnitude of 10 %). The contribution of the equipment efficiencies in the TCO variations is even smaller with values around 0.1 % in Figure 3c.

IV. SCNENARIO-BASED STOCHASTIC OPTIMIZATION

As observed in the previous section, the economic environment (i.e. rates, policies/regulation) is an important input parameter for the cost-driven design of multi-energy systems. Thus, the main motivation of the work presented in the paper is to investigate different approaches to incorporate the various economic parameters, rather than implementing high accuracy models for the different equipments. Note that conventional robust optimization for microgrid operation and planning focuses on capturing the uncertainties for the load and/or renewable generation profiles. In this paper, that class of uncertainties is implicitly embedded in the whole representative week considered for the planning problem (which may not be long enough if other geographical areas or types or loads are considered).

A. Stochastic Energy Prices-OPT1

A classical approach is firstly considered to tackle the prices uncertainties and corresponds to scenario-based stochastic optimization. Different values (or levels) are considered for the energy rates (π t e , π t g and π p ) with different probabilities (Figure 4a). Tests with two, three and five levels are performed (normalized range as in section III.C), and the set of potential scenarios S contains all the possible combinations of prices in every cases. Obviously the total number of scenarios increases significantly when more levels are explored for each parameter (Figure 4b), which may lead to prohibitive computational times. Each scenario has a probability of occurrence λs computed as the product of the probabilities for each price parameter. The planning problem objective function is then reformulated following (8) with the OPEX computed along all the potential scenarios and weighted by their probabilities. Finally, all the operating variables are defined along the set of potential scenarios S (e.g. P t,s ge ) and the CAPEX calculation remains unchanged as it is not affected by the energy rates.

                      

B. Multi-Period Optimization -OPT2

When the system design is expected to have an operational lifetime in the range of decades (tpypically 20-30 years), it is not realistic to adopt constant values for the different energy rates. Thus, another way to model the economic environment consist in discretizing the planning period into sub-periods p representing a set P. Similar to the previous problem formulation, the OPEX computation is affected (( 9)) and the operating variables are defined along the set P. For each period p, the OPEX is then divided by the number of years in the sub-period so that the final TCO remains on a yearly basis. Note that the price for peak power is projected to be constant over the planning period in the optimization and that the gas and electricity rates are projected to increase by 28 % and 15 % respectively (www.eia.gov).

 

               C. Multi-Period Stochastic Optimization -OPT3
Another possible method to model the economic environment could consist in combining OPT1 and OPT2. Thus, different set of price profiles are considered along the periods p and with different probabilities of occurrence (Figure 5), similar to the scenario-based approach implemented in [START_REF] Alharbi | Stochastic optimal planning of battery energy storage system for isolated microgrids[END_REF]. In the considered test case, four periods of five years each are considered in the planning over 20 years. The OPEX computation is re-formulated (( 10)) and the operating variables are defined over the sets T, P and S, which tends to significantly increase the problem complexity if too many potential scenarios are investigated (i.e. high number of different levels for the parameters).

 

, , , ,

                   

D. Stochastic Energy Prices and Equipment Costs -OPT4

Another investigated problem formulation that captures the economic uncertainties consists in considering stochastic values for the equipment cost as well as for the energy prices. Similar to OPT1 and OPT3, a set of different likely scenarios s is defined depending on the values of electricity/gas/peak rates (3 parameters) and the price level for the different types of assets (4 parameters). For a given size of the assets, the CAPEX is calculated according to [START_REF] Li | Microgrid sizing with combined evolutionary algorithm and MILP unit commitment[END_REF] with the equipment cost taking different values and with the corresponding scenario probability λs. The OPEX calculation is the same as in [START_REF] Rigo-Mariani | Comparison of optimization frameworks for the design of a multienergy microgrid[END_REF]. Note that, compared to OPT1, the problem formulated here includes a larger number of scenarios for the same number of price levels as there are more parameters involved in the uncertainties consideration.

CAPEX C P C P C E C E              E.

Obtained Results

Several optimizations are run following the different formulations described in the previous subsections. The results are displayed in Table II with OPT0 denoting the base case problem presented in Section III. As every formulation is different, the obtained values in terms of TCO are obviously different in every cases and vary from 978 k$ to 1,110 k$. As highlighted by the sensitivity analysis, those variations are mainly explained by different representations of the economic context. It is important to note that all the cases are mathematically different and do not correspond to a same problem. The different formulations then do not have to be compared in terms of TCO but rather with regard to the obtained value for the size of the different assets. Those different simulations have to be interpreted a different ways to tackle the planning issues. In terms of optimal sizes of the equipment, the PV is systematically set at the maximum allowed capacity (upper bound at 500 kW) as it represents immediate energy/cost savings with low investment costs compared to the other assets. The size of the battery displays small variations around 130 kWh on average depending on the cases while the results for the thermal storage varies more significantly from 100 to 283 RTh. The most noticeable outcome is the gas engine optimal capacities covering 60 % (301 -929 kW) of the specified range 0 -1 MW depending on the problem formulations. This specifically shows the necessity to pay the most attention to the modeling of the economic environment, which is likely to be more uncertain than load or solar radiation profiles, especially when planning decades ahead. However, a stochastic representation that allows to fully capture the uncertainties over the planning period may result in higher complexity with intractable problems due to the increased number of variables and constraintsespecially if too many scenarios are investigated -see OPT3 and 3 levels for energy prices and equipment cost. The obtained results show a great variety of solutions with small variations of the objective function (i.e. the yearly TCO here). That outcome is explained by the nature of the planning problem with the objective function in terms of yearly TCO and a small contribution of the equipment cost compared to the operating expenditures. Figure 6 displays small variations of the TCO around the optimal solution. The gas engine and energy storages capacities can vary significantly (> 50 % of the specified range) without varying the TCO by more than 5 %. This implies that under the specified equipment cost and energy prices, the capital expenditures are somewhat compensated by the OPEX savings. Ultimately, that 'flatness' of such long term planning problem with linear formulation implicitly guarantees the robustness of the results.

V. CONCLUSIONS

While optimal design of multi-energy systems under uncertainties has been widely addressed in literature, this paper aims at a step back in order to estimate the relevance of complex problem formulation. After a sensitivity analysis is performed, it is showed that the energy prices are , by far, the most impacting parameter, which is consistent with the conclusions drawn in [START_REF] Mavromatidisa | Uncertainty and global sensitivity analysis for the optimal design of distributed energy systems[END_REF]. This paper investigates different formulations in order to model the economic environment with different energy rates as well as uncertain equipment costs. Ten different optimization runs are performed with various assumptions and choices of models. The problem complexity tends to increase drastically when high numbers of stochastic scenarios are explored, which leads non tractability, even with the simple linear models considered here. The obtained results display significant differences in terms of installed capacities while the objective remain somewhat constant due to the chosen hypothesis to represent the price uncertainties (close to a normal law, centered on the average value compared to uniform law for instance). That outcome is also explained by the nature of the planning problem itself with wide 'flat' areas in the search space regarding the objective function. Thus, in addition to the robustness with respect to the input parameters (e.g. the prices here), such planning problem should also be tackled in terms of robustness of the results. Those results should then be expressed as intervals/sub-search space and not only vectors, and the search for a unique global optimum with the all the sources of uncertainties may not be the most suitable approach. The idea would be to include new constraints in the optimization in order to return the widest 'clusters' with small variations around the cluster optimal value of the objective. Those should be part of further developments that would require customized problem formulations and/or solving methods.

Figure 1 .Figure 2 .

 12 Figure 1. Considered Multi Energy System

Figure 3 .

 3 Figure 3. Sensitivity Analysis-a) Energy pricesb) Cost of equipmentc) Efficiencies

Figure 4 .Figure 5 .

 45 Figure 4. Stochastic prices-a) Price levelsb) Set of scenatios

Figure 6 .

 6 Figure 6. Small variations of TCO around the best solution

TABLE I

 I

				.		PARAMTER RANGE FOR SENSITIVITY ANALYSIS
	Param.		Range					Param.		Range
	π t g			0.5 -0.9 $/kg				C batR		800 -2000 $/kWh
	π t e π t p			0.85×π t e -1.145 π t e 40 -60 $/kW		C ts η bat		100 -200 $/RTh 0.88 -0.98
	C geR		800 -1600 $/kW		η ts		0.85 -0.95
	C pvR		750 -1750 $/kWp		COP abch 1.05 -1.15
	SOC	bat	E	batR	0 bat SOC	E	batR	t	t P	bat bat	bat	t P	bat	100	t
									i	1
	0 bat SOC	E	batR	t	t P	bat bat		bat	t P	bat	100		t SOC	bat	E	batR
					i	1
	0 bat SOC	E	batR	T	t P	bat bat		bat	t P	bat	100		0 bat t SOC	E	ba
					t	1

TABLE II .

 II OBTAINED RESULTS FOR DIFFERENT PROBLEM FORMULATIONS


Judy Najnudel
email: judy.najnudel@ircam.frremy.muller@ircam.fr

Rémy Müller

Thomas Hélie
email: thomas.helie@ircam.fr

David Roze
email: david.roze@ircam.fr

A POWER-BALANCED DYNAMIC MODEL OF FERROMAGNETIC COILS

This paper proposes a new macroscopic physical model of ferromagnetic coils used in audio circuits. To account for realistic saturation and hysteretic phenomena, this model combines statistical physics results, measurement-driven refinements and port-Hamiltonian formulations that guarantee passivity, thermodynamic consistency and composability according to both electric and thermal ports. As an illustration, the model is used to simulate a passive high-pass filter. Different types of audio inputs are considered and simulations are compared to measurements.

INTRODUCTION

Ferromagnetism is frequent in analog audio: it is involved in transducers (dynamic microphones, loudspeakers), tape recorders, coils and transformers. As major non-linearities arise from ferromagnetic components (saturation, hysteresis), the need of refined models is critical to accurately simulate behaviors in circuits.

Since the 1980s, a large body of empirical models have been proposed, among them the Jiles-Atherton model [START_REF] David | Theory of ferromagnetic hysteresis[END_REF], the Gyrator-Capacitor model [START_REF] David | Gyrator-capacitor modeling: a better way of understanding magnetic components[END_REF][START_REF] Chen | Gyrator-capacitor simulation model of nonlinear magnetic core[END_REF], or the Preisach model [START_REF] Bertotti | Dynamic generalization of the scalar preisach model of hysteresis[END_REF]. But very few have a strong physical meaning [START_REF] Sergey E Zirka | On physical aspects of the Jiles-Atherton hysteresis models[END_REF] and those retaining some energetic interpretation [START_REF] François-Lavet | An energy-based variational model of ferromagnetic hysteresis for finite element computations[END_REF] either lose major phenomenological properties or are heavy to implement [START_REF] Lobue | Entropy and entropy production in magnetic systems with hysteresis[END_REF]. As a consequence, preserving the model passivity (no artificial hidden sources of energy) comes with a price -computation time.

In this paper, we propose a new nonlinear model of ferromagnetic coil that is physically-based, passive, modular (allowing electric and thermal connections) and with a reduced complexity (few state variables and parameters). As it is built on statistical physics results on magnets, it is thermodynamically consistent. It also inherits macroscopic characteristics (hysteresis and its conditioned activation) from underlying microscopic phenomena (metastability and phase transition). This lumped-element model is used to simulate a passive high-pass filter. The circuit modeling relies on Port-Hamiltonian Systems [START_REF] Maschke | An intrinsic Hamiltonian formulation of network dynamics: Non-standard Poisson structures and gyrators[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF] (PHS) that structurally fulfill the power balance. Simulations are based on numerical methods [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach[END_REF] that preserve this guarantee in the discrete-time domain.

The paper is structured as follows: Section 2 first presents some short recalls on PHS. Section 3 develops a primary model derived from statistical physics. This model exhibits saturation and hysteresis but does not take into account some other phenomena, such as non-homogeneities, thermal fluctuations and eddy currents. Section 4 refines the primary model with a polynomial interpolation based on measurements of a Fasel inductor. In section 5, the final nonlinear model is implemented to simulate a passive high-pass filter.

REVIEW OF PORT-HAMILTONIAN SYSTEMS

The following modeling relies on Port-Hamiltonian systems [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF], under a differential-algebraic formulation [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach[END_REF]. A dynamical system is represented as a network of: (i) storage components of state x and energy E (x), (ii) dissipative components described by an efforts law w → z(w) that dissipates the power P diss = z (w) ⊺ w ≥ 0 for all flows w, and (iii) connection ports conveying the outgoing power Pext = u ⊺ y for inputs u and outputs y. The flows f and efforts e of all the components are coupled through a skew-symmetric interconnection matrix J = -J ⊺ :

  ẋ w y   f = J   ∇E(x) z(w) u   e . (1)
Such systems satisfy the power balance P stored +P diss +Pext = 0 where P stored = ∇E(x) ⊺ ẋ denotes the stored power. Indeed, P stored + P diss + Pext = e ⊺ f = e ⊺ Je is zero since e ⊺ J e = (e ⊺ Je) ⊺ = -(e ⊺ Je) due to the skew-symmetry of J .

All models herein will be formulated as (1).

PRIOR THEORETICAL MODEL

Macroscopic model of a ferromagnetic core

This section presents a bi-stable core model rooted in the mean field Ising model [START_REF] Ising | Beitrag zur theorie des ferromagnetismus[END_REF][START_REF] Gordon | On the theory of the Ising model of ferromagnetism[END_REF][START_REF] Strecka | A brief account of the Ising and Ising-like models: Mean-field, effective-field and exact results[END_REF][START_REF] Utermohlen | Mean field theory solution of the Ising model[END_REF][START_REF] Néel | Magnetism and local molecular field[END_REF]. Using statistical physics, Ising derives a macroscopic scalar state (the core magnetization) from a microscopic representation of the core (a set of normalized atomic magnetic moments, which can be either up or down), at a given temperature T . For simplicity, additional assumptions are: a homogeneous, isochoric (constant volume V ) and closed (constant number of atoms N ) ferromagnetic crystal with one (local) magnetization axis and periodic boundaries (typically, a toroidal geometry often found in audio circuits [START_REF] Wm | Transformer and inductor design handbook[END_REF]). In this section, there is no external magnetic field (issue addressed in section 3.2).

DAFx.1

Proceedings of the 23 rd International Conference on Digital Audio Effects (DAFx-20), Vienna, Austria, September 8-12, 2020

Macroscopic quantities and laws

Following [START_REF] Barry | The two-dimensional Ising model[END_REF] with our assumptions, the core internal energy is

E = N α 1 2 m 2 -m tanh m θ , (2)
where parameter α (≈ 5 × 10 -21 J for transition metals) denotes the exchange energy between one moment and its nearest neighbours [START_REF] Stuart | Direct exchange in ferromagnets[END_REF][START_REF] Michael | Lattice statistics in a magnetic field, I. A two-dimensional super-exchange antiferromagnet[END_REF][START_REF] Hasegawa | Considerations on double exchange[END_REF], and where variables m and θ are average intensive quantities (homogeneous over the body) that statistically characterize the core configuration at a macroscopic scale:

• m ∈ [-1, 1]
is the mean normalized magnetic moment: m = ±1 if all moments are aligned in the same direction, and m = 0 if no particular direction is favored;

• θ = T /Tc ∈ R+ is the reduced temperature relative to the core Curie temperature Tc [START_REF] Gs Rushbrooke | On the Curie points and high temperature susceptibilities of Heisenberg model ferromagnetics[END_REF]: if θ < 1, there are multiple equilibria m (ferromagnetism), and only m = 0 (disordered moments) otherwise (paramagnetism).

Note that the core parameters α and Tc are related through the Boltzmann constant k b = 1.38 × 10 -23 J.K -1 as α = k b Tc.

A measure of the number of possible microscopic states (atomic moments) consistent with the core macroscopic configuration is given by the entropy [START_REF] Edwin | On the rationale of maximum-entropy methods[END_REF], which is found to be

S = N k b f m θ with f (χ) = ln(2 cosh χ)-χ tanh χ, χ ∈ R.
(3) This statistical entropy coincides with the thermodynamic entropy for a core in internal thermodynamic equilibrium (possibly timevarying at macroscopic scale). This variable is extensive (proportional to N ) and quantifies the macroscopic "order degree" of the core, on which phase transitions and hysteresis depend.

In addition to E and S, a third extensive variable is introduced, namely, the total magnetic flux of the core (of volume V )

BV = B V, ( 4
)
where B is the magnetic flux density. For the core, B is related to the core magnetization M = m Ms through B = µ0 M where µ0 is the vacuum magnetic permeability and Ms is the saturation magnetization (see Table 1).

Choice of state and energy function

We choose to express the core energy E as a function of the state

xcore = [BV , S] ⊺ , (5)
so that, in (1), the flow ẋcore accounts for the time variation of extensive quantities (to balance with quantities external to the core, or equilibrium-establishing) and, concomitantly, the effort accounts for intensive quantities (shared with the exterior at the core interface, or equilibrium-determining). Choosing extensive energy variables over intensive co-energy variables stems from two arguments. The first one is physical: except for linear, mono-variate components, constitutive laws derived from the co-energy are not equal to those derived from the energy, and should be handled with care. The second is numerical: solving an ODE by integration instead of differentiation is generally preferable [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF]. This function is derived from (2), in which m/θ and m are expressed with respect to S and BV using (3-4) and noting that f is smooth, even on R and strictly monotonic 1 (so invertible) on R + . Its formula expressed w.r.t. ( 5) is given by (see Fig. 1),

Ecore [BV , S] ⊺ E0 = 1 2 BV BV s 2 - BV BV s tanh f -1 S S0 , (6) with core constants E0 = N α, S0 = N k b and BV s = V µ0Ms. The energy gradient (effort) is ∇Ecore = [Hcore, Tcore] ⊺ , (7)
where, omitting variables in functions, the core internal magnetic field is

∂ Ecore ∂B V = E0 BV s   BV BV s -sign (BV ) tanh f -1 S S0   := Hcore (8)
and the core temperature is

∂ Ecore ∂S = E0 S0 BV BV s / f -1 S S0 := Tcore. (9)
Fig. 1 shows that as S increases, the core goes from two ordered (aligned moments) meta-stable equilibrium states to one non-ordered (no alignment) stable equilibrium state: it exhibits a phase transition (from ferromagnetic to paramagnetic). Table 1 recaps the physical quantities involved and their units.

Connection to coil and external electrical ports

Ideal coil model

The coil is considered to be linear. Choosing BV as its state variable, the coil energy is

E coil (BV ) = B 2 V 2µ0V , (10)
and its derivative with respect to BV is the coil magnetic field H coil (BV ) = BV /(µ0V ).

Symbol Quantity S.I. units

N atoms nb. dimensionless α nearest neighbours exchange energy kg.m 2 .s. -2 m norm. magnetic moment dimensionless M magnetic moment A.m 2 M magnetization A.m -1 H magnetic field A.m -1 B magnetic flux density kg.s -2 .A -1 µ0 vacuum magnetic permeability kg.m.s -2 .A -2 Φ magnetic flux linkage kg.m 2 .s -2 .A -1 n coil turns nb. dimensionless k b Boltzmann constant kg.m 2 .s -2 .K -1 T temperature K S entropy kg.m 2 .s -2 .K -1 V = Al volume = section x length m 3
Label L coupled core and coil

Coupled system

To express the coupled system L as a PHS, one needs to determine the relations between the core and coil flows and efforts. Using the extensivity of the total magnetic moment ML [START_REF] Daniel | Statistical physics: An introductory course[END_REF], one gets

ML = Mcore + M coil ⇔ BV L = BV core + BV coil , (11)
since BV = µ0M. Differentiating Eq. ( 11) with respect to time, one obtains the relation between the flows:

Ḃ V L = ḂVcore + ḂV coil . (12)
Conversely, the coil and core share their efforts, namely,

HL = H coil = Hcore. (13)
Fig. 2 represents the coupling as a series connection.

For any fixed entropy S and for all x = [BV coil BV core ] ⊺ , we introduce the function ES : x → Ecore(BV core , S) + E coil (BV coil ) (total energy of the system). With these notations, the core and coil coupling can be expressed as the constrained Dirac structure [START_REF] Van Der Schaft | Port-Hamiltonian differential-algebraic systems[END_REF] 

 ẋ 0 y   =   . A B -A . . -B . .     ∂ E S ∂x λ u   , (14)
with

A = [1, -1] ⊺ , B = [0, 1] ⊺ , λ = ḂV coil , u = ḂV L
and y = -HL (dots indicate zeros). This constrained Dirac structure can be reduced to (see also [START_REF] Van Der Schaft | Port-Hamiltonian differential-algebraic systems[END_REF] for more details):

ż y = . Br -Br . ∂ E L ∂z u , (15)
with à such as Ã⊺ A = 0 to eliminate the constraint, Br = Ã⊺ B, z = Ã⊺ x, EL the total energy with respect to z.

Taking à = [1, 1] ⊺ , this yields Br = 1 and z = BV coil + BV core . Therefore, for any given entropy S, the dynamics of the coupled system is that of an equivalent component of state xL = [BV L , S], energy EL(xL) and magnetic field

HL = ∂ E L ∂B V L
. This equivalent component energy can be computed (see [START_REF] Najnudel | Simulation of the ondes Martenot ribbon-controlled oscillator using energy-balanced modeling of nonlinear time-varying electronic components[END_REF] for a detailed derivation) through the expression

EL(BV L , S) = E coil • H -1 coil + Ecore • H -1 core • H -1 coil + H -1 core -1 (BV L , S), (16)
where the symbol • stands for function composition. In practice, all mathematical functions in this expression can be defined as piecewise affine functions (computation of inverse efforts in particular becomes straightforward when analytical expressions are not available, as is the case here).

Connection to external electrical ports

Denoting n the number of turns, l the length of the coil, A its section, Φ the magnetic flux linkage, the magnetic field HL is related to coil current iL through

HL = n l iL, (17)
and the state BV L is related to the coil voltage vL through

ḂV L = Φ nA V = l n vL. (18)
In section 4.4, variables iL and vL will be related to external ports u and y of Eq. (1).

Thermodynamics

In this section (only), we assume that the ferromagnetic coil is put in a isothermal bath (i.e. the exterior is much larger than the coil size), so that the temperature of the system TL is considered constant and equal to the exterior temperature Text during a change of state, supposedly below the Curie temperature. A convenient and classical way to study the energetic behavior of the ferromagnetic coil is to examine how, for all BV L , the energy EL(BV L , S) of the component subject to a constant magnetic field H0, deviates from the energy H0 BV L . The energy deviation of this conditioned component, called the Gibbs free energy [START_REF] Bertotti | Hysteresis in magnetism: for physicists[END_REF], is defined by, for all BV L , S, and all constant-in-time H0 as

GH 0 (BV L , S) = EL(BV L , S) -TLS -H0 BV L . (19)
For any given S, at H0 = 0, two symmetric meta-stable equilibrium states corresponding to GH 0 minima with respect to BV L are available (Fig. 3a, red curve). For H0 ̸ = 0, the symmetry is broken and a previously stable equilibrium state can be made unstable. We suppose now H0 slowly controlled (so it is still considered constant during a change of state). When decreasing H0 from Hmax ≥ 0 (Fig. 3a, green curve) to -Hmax, BV L starts from its initial equilibrium state and follows a trajectory solution of

∂ G H 0 ∂B V L
= 0 (Fig. 3a, black curve), until the minimum degenerates into an inflection point at H0 = Hcoerc (Fig. 3a, solid blue curve). Then, a Barkhausen jump occurs [START_REF] Bertotti | Hysteresis in magnetism: for physicists[END_REF] so that BV L occupies the remaining stable equilibrium state (Fig. 3a, intersection of solid blue curve and left yellow area). Since EL is even with respect to BV L for all S, GH 0 (BV L , S) = G-H 0 (-BV L , S). Therefore, when increasing H0 from -Hmax to Hmax, the Barkhausen jump occurs at -Hcoerc. Consequently, BV L follows a different path depending on whether H0 decreases or increases (Fig. 3a, black curve and arrows), hence the hysteresis (Fig. 3b) between Hcoerc and -Hcoerc.

Thermodynamics laws show that the difference of energy before and after the jump is irreversibly dissipated as heat. Indeed, the first principle of thermodynamics states that the internal energy variation dEL is the work performed on the ferromagnetic coil δW = H0dBV L , plus the received heat δQ = TLδeS where δeS is the variation of incoming entropy and δ denotes an inexact differential [START_REF] Peter T Landsberg | Thermodynamics and statistical mechanics[END_REF]: The second principle of thermodynamics states that the internal heat TLdS is the received heat plus the heat internally produced by irreversible phenomena TLδiS:

dEL = H0 dBV L + TLδeS. (20)
TLdS = TLδeS + TLδiS. ( 21
)
Replacing TLδeS from Eq. ( 21) in Eq. ( 20) yields

dGH 0 = dEL -TLdS -H0dBV L = -TLδiS, (22)
which is consistent with the assertion that the difference of energy is entirely and irreversibly dissipated as heat. Now, let us denote HL the observed effort law such as H0 = HL(BV L ) (definition given in appendix A). Replacing H0 with HL, the entropy production rate δiS/dt is obtained differentiating Eq. ( 22) with respect to BV L and multiplying with ḂV L :

δiS dt = 1 TL HL(BV L ) - ∂ EL ∂B V L (BV L , S) ḂV L . ( 23
)
To model the conversion between excess electro-magnetic power and thermal power, the ideal thermal exchanger r th is introduced (Fig. 4) so that

i th v th = TL δiS dt ( 24
)
where v th is the exchanger voltage and i th its current. Introducing the function

P th : xL → HL(BV L ) -∂ E L ∂B V L (BV L , S) ḂV L , (25)
equations [START_REF] Edwin | On the rationale of maximum-entropy methods[END_REF][START_REF] Daniel | Statistical physics: An introductory course[END_REF] allow to model the dissipation in the PHS formalism:

w = [v th , TL] ⊺ z P th (x L ) (w) = [ P th (xL) v th , - P th (xL) TL ] ⊺ = [i th , - δiS dt ] ⊺ . ( 26
)
The passivity condition P diss ≥ 0 is fulfilled as z(w) ⊺ w = 0.

The complete PHS structure is given in section 4.4.

REFINED MODEL BASED ON MEASUREMENTS

Measurements and observations

As thermodynamically meaningful as the bi-stable model is, it does not capture the variety of phenomena contributing to hysteresis, as measurements on real ferromagnetic coils reveal. To conduct such measurements, a Fasel Red inductor (which can be found in Cry Baby wah-wah pedals [START_REF] Falaize | Simulation of an analog circuit of a wah pedal: a port-Hamiltonian approach[END_REF] for instance) in series with The magnetic field HL and state BV L are then obtained using Eq. ( 17)-( 18) and the relation Φ(t) = t 0 vL(τ ) dτ . As the coil and the core share the same volume V , Eq. ( 11)-( 13) yield the relation BV L = µ0V (HL + M ) from which M is obtained. These measurements (Fig. 5b) lead to two observations.

• First, M has an order of magnitude of 6, whereas HL has an order of magnitude of 1, as expected for soft materials [START_REF] Wm | Transformer and inductor design handbook[END_REF].

• Second, instead of the large jumps predicted by the bi-stable model, one observes a continuous progression, which calls for a model refinement to determine the entropy production rate law responsible for hysteresis.

Model reduction

According to measurements, for this inductor BV core = µ0M V ≫ BV coil = µ0HLV and Ecore ≫ E coil . The influence of the coil on the overall energy of the coupled component is negligible (Fig. 6) and we thus may use

EL(BV L , S) ≈ Ecore(BV L , S). (27)
The dynamics of the coupled system is therefore that of a driven core alone.

Entropy production rate law

In real ferromagnetic cores, domain structure and non-homogeneities [START_REF] Kittel | Physical theory of ferromagnetic domains[END_REF] yield an energy function with not two but multiple local minima. Consequently, multiple Barkhausen jumps give the effort law the shape of a staircase. The Preisach model generates this effort law by computing a statistical mean on a collection of bi-stable systems such as the one presented in section 3, each one representing a domain. This averaging "damps" the large bi-stable jumps.

Here, to obtain a similar result while remaining at a macroscopic level, the hysteresis loop is modeled using a cubic polynomial

P (χ) = p0 + p1χ + p2χ 2 + p3χ 3 interpolating the effort ∂ E L ∂B V L
, and an additional friction term of the form r f ḂV L , r f ≥ 0, to account for thermal fluctuations [START_REF] Brown | Thermal fluctuation of fine ferromagnetic particles[END_REF] and eddy currents [START_REF] Cd | Physical origin of losses in conducting ferromagnetic materials[END_REF]. The coefficients of P are computed through

[p0 p1 p2 p3] ⊺ = X -1 Y
where, given two interpolation data points χ1 and χ2, X and Y are defined as

X =      1 • • • χ 3 1 1 • • • χ 3 2 0 • • • 3χ 2 1 0 • • • 3χ 2 2      , Y = ∂ E L ∂B V L (χ1) ... ∂ 2 E L ∂B 2 V L (χ1) ... ⊺ .
The final hysteresis loop P (BV L ) is thus defined by

P (BV L ) = δBP (δBBV L ) + r f ḂV L , (28)
where δB = sign(dBV L ), and the entropy production rate δiS/dt is:

δiS dt = 1 TL P (BV L ) - ∂ EL ∂B V L (BV L , S) ḂV L , (29)
which is the expression given in Eq. [START_REF] Edwin | On the rationale of maximum-entropy methods[END_REF] where HL has been replaced with P . For a given ferromagnetic coil, such a loop is accurate in a range from saturation approach to saturation and higher, provided that the data points are taken in that range. At lower fields though, a Rayleigh law would be more adequate [START_REF] Bertotti | Hysteresis in magnetism: for physicists[END_REF].

Final model

Finally, Kirchhoff laws on the equivalent circuit shown on Fig. 4, together with Eq. ( 21), yield the PHS in Fig. 7 structured as in Eq.

(1), in which EL is given by Eq. ( 27)-( 6), w and z(w) are given by Eq. ( 26)-( 29), u = [U, δeS/dt] ⊺ , y = [i, -Text].

Model identification with the Fasel inductor

The measurements are taken during an isothermal transformation, so that, replacing S from Eq. (3) in the magnetic field, one can use the expression

∂ EL ∂B V L = E0 BV Ls   BV L BV Ls -tanh BV L BV Ls θ  
for identification. A least squares optimization between the Eq. ( 28) spline model and the measurements, i.e. solving p = argmin p HL -Pp BV L 2 2 with p = [E0, BV Ls , θ, r f ] yields the parameters in Table 2. Figure 8 shows a good match between measurements and the estimated model.

∇E(x) z(w) u TL ∂ E L ∂B V L i th -δ i S dt U δeS dt               ẋ Ṡ . . . -

APPLICATION TO A PASSIVE HIGH-PASS FILTER

Circuit modeling

The ferromagnetic coil model is used to simulate a high-pass filter (Fig. 9). The resistance R is linear of constitutive law vR(iR) = RiR. Kirchhoff laws yield the PHS shown in Fig. 10.

Simulation

Discretization

The state vector x(t) is discretized to x[k] = x(hk) where h = 1/Fs is the sampling step, and we denote δx

[k] = x[k+1]-x[k].
To preserve the passivity of the PHS in discrete time and reduce the energy gradient sensitivity to the state indexing, we rely on the symmetric discrete energy gradient [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach[END_REF].

∇E(x) z(w) u TL ∂ E L ∂B V L i th -δ i S dt vR U δeS dt                   ẋ Ṡ . . . -1 . . 1 ḂV L . . . . -l n l n . v th . . . . - 1
∇E(x, δx)i =        1 nx! δxi π∈P(nx) △i(xπ, δxπ) δxi ̸ = 0 ∂ E ∂x i otherwise ( 30
) where △i(x, δx) = E(x + δxi) -E(x + δxi-1) and δxi = [δx1, ..., δxi, 0, ..., 0] ⊺ . The discrete energy variation is obtained with the chain rule:

δE[k] h = ∇E(x[k], δx[k]) ⊺ δx[k] h . (31)
The PHS of Fig. 10

Newton-Raphson iteration

The interconnection matrix is decomposed as

J = [J stored J diss Jext] ⊺ . We denote ē(x[k], δx[k]) = [ ∇E(x[k], δx[k]) z(w[k]) u[k]] ⊺ , ν = δx[k] and
F : ν → J stored ē(x[k], ν) -ν/h (32)
At each sample k, δx[k] is solution of F (ν) = 0. If F ′ (νi) is invertible and given an initial value ν0 and a relative error ϵr, this solution can be computed iteratively with the update

νi+1 = νi + ∆ν i (33)
where

∆ν i = -F ′ (νi) -1 F (νi), until ∥∆ν i ∥/∥∆ν 0 ∥ ≤ ϵr. The state x[k+1] is then computed using x[k+1] = x[k]+δx[k].

Simulation parameters

The circuit is driven with a sinusoidal voltage whose parameters are given in Table 3, as well as with an instrumental bass sample. The ferromagnetic coil model parameters are those indicated in Table 2. The incoming entropy flow δeS/dt is set so that the ferromagnetic coil temperature stays constant.

Results and comparison to measurements

The circuit is simulated with the non-linear coil model and a simple linear coil model (iL = ΦL/L with L = 840 mH) for comparison. Simulation results on Fig. 11a-11b show a good correspondence between the non-linear model and measurements. Fig. 11c-11d show that the produced entropy is always positive and that the coil temperature stays constant. Spectrograms of the bass sample is shown on Fig. 11e-11f. Sound results on the bass sample can be heard at https://github.com/JNaj/dafx20-ferromag.

CONCLUSION

In this paper, a physical and passive model of ferromagnetic coil has been developed. It is explicit and maintains a reduced number of variables and parameters. First the core and the coil were treated separately, then their coupling, which determines both their electrical and thermal dynamics, was addressed. This lead to the building of an equivalent component, characterized by a well-established state, energy function, and entropy production rate law. A refined entropy production law based on measurements was then proposed.

As an application, this model was used to simulate a passive high-pass circuit. The simulations are in close agreement with measurements, though extensive measurements (a set of different frequencies, amplitudes, waveforms) would be required to validate the model on a broader scale.

Besides these complementary measurements, further work aims to assess real-time performances, and build a transformer model on the same principle by coupling two ferromagnetic coils.

Figure 1 :

 1 Figure 1: Core energy function with respect to BV and S.

  (a) Gibbs free energy G H 0 for H 0 decreasing from Hmax (green curve) to Hcoerc (solid blue curve), at constant temperature, and trajectory of B V L (black curve) for a complete cycle. In the yellow area, two local potential minima coexist but only one direction is possible for B V L to follow (blue arrows). (b) Observed state B V L during a complete cycle, resulting in Barkhausen jumps (blue curve), and theoretical effort ∂ E L ∂B V for B V ranging from min(B V L )to max(B V L ) (red curve). The area between the blue and red curves is the energy dissipated during a cycle.

Figure 3 :

 3 Figure 3: Gibbs free energy GH 0 for decreasing values of magnetic field H0 (3a), and observed state BV L during a complete cycle of magnetic field variations (3b).

Figure 4 :

 4 Figure 4: Voltage-controlled ferromagnetic coil with thermal dissipation.

Proceedings of the 23

 23 rd International Conference on Digital Audio Effects (DAFx-20), Vienna, Austria, September 8Measurements setup (thermal dissipation not represented). (b) Measured H L and M .

Figure 5 :

 5 Figure 5: Measurements setup and results.

Figure 6 :

 6 Figure 6: Core energy and coupled system equivalent energy computed with piecewise affine functions.

Figure 7 :

 7 Figure 7: PHS of the voltage-controlled ferromagnetic coil with thermal dissipation. Dots in the interconnection matrix indicate zeros.

Figure 8 :

 8 Figure 8: Measurements (red curve) and estimated spline model (blue curve).

Figure 9 :

 9 Figure 9: Passive high-pass filter.

( a )

 a Measured and simulated currents. (b) Measured and simulated voltages. (c) Simulated entropy flows. (d) Simulated temperature. (e) Bass sample spectrogram. (f) Filtered bass sample spectrogram.

Figure 11 :Table 3 :

 113 Figure 11: Simulation results.

Table 1 :

 1 Physical quantities and labels.

	HL	H	coil
			ḂV coil
			ḂV L

coil core ḂVcore Hcore Figure 2: Coil and core connection.

Table 2 :

 2 Denoting nx the number of states, P(nx) the set of all possible permutations on the Physical parameters of the model where zi = zi/BV Ls and pi = piB i V

			Estimated
	E0 2.43.10 -5 3.09.10 -7 1.10 6.07.10 4 θ r f BV Ls	p0 0	p1 8.69 0 p2 p3 8.78
				Given
	µ0	k b	n	V	z1	z2
	4π.10 -7 1.38.10 -23 150 2.22.10 -7 -1	1

Ls . The units are S.I. units given in Table

1

.

Indeed, f ′ (χ) = -χ/ cosh

χ ≤ 0 ∀χ ∈ R + .

Proceedings of the 23 rd International Conference on Digital Audio Effects (DAFx-20), Vienna, Austria,September 8-12, 2020

The authors acknowledge the support of the ANR-DFG (French-German) project INFIDHEM ANR-16-CE92-0028.

A. DEFINITION OF THE BI-STABLE MODEL

OBSERVED EFFORT LAW HL

(BV L 0 ) = 0 (Fig. 3b, green cross), and

green plus), one can define HL as:

where [BV Linf , BV L sup ] = [ BV L 0 , BV L 0 ] if dBV L ≤ 0 and [-BV L 0 , -BV L 0 ] otherwise.


T T Q Nguyen

P T Tran

Vincent Debusschere

Christophe Bobineau

Rémy Rigo-Mariani

T P T Tran

Comparing high accurate regression models for short-term

Keywords: Regression model, Seq2Seq model, attention mechanism, Deep multilayer perceptron, Long short term memory, Support vector regression, random forest, short-term forecasting, load forecasting

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Smart grid technologies help collecting and analyzing a huge volume of data, with the objective of optimizing the management of the power system at both the producer and the demand sides. This volume of data can contribute to the optimization of the power usage, considering economic benefits as well as environmental constraints. The forecast of the loads consumption (for both tertiary, commercial and residential buildings) is at the core of most of the optimization tools design to take the good decision in power conservation and management planning [START_REF] Yildiz | A review and analysis of regression and machine learning models on commercial building electricity load forecasting[END_REF].

Load forecasting can be divided into three categories according to the time scale [START_REF] Feinberg | Load forecasting[END_REF].

1) Short-term: From a few minutes to a few days ahead; 2) Medium-term: From a few days to a few months ahead; 3) Long-term: From months to years.

Most of the researches has focused on short-term load forecasting (STLF), because it is critical for real-time power generation, operation and economic transactions. Besides, it is difficult to achieve an acceptable accuracy for short-term load forecasting. In fact, smart grids technologies intend to collect sensors data at a higher sample rate, targeting real-time monitoring to obtain an improvement in the prediction result as Institute of Engineering Univ. Grenoble Alpes well. STLF is also complex because of the variety of the dataset. Thus, it is challenging, for the predicted model, to capture relevant variations. For instance, buildings energy consumption do not follow clear patterns, rules or cyclic elements. The output of the prediction algorithm (the power consumption) is a real number that is influenced by many input factors such as temperature, humidity, wind speed , most of which are exogenous and quite hard to forecast.

Two main approaches are considered currently in forecasting techniques: physical based and data-driven based ones. The physical based techniques, such as EnergyPlus or Ecotest [START_REF] Mocanu | Optimized parameter selection for assessing building energy efficiency[END_REF], are based on thermodynamics relationships to build the energy model. The accuracy of the predicted values depends on the accuracy of the simulated model and its inputs. However, some details of data are not available in simulation time. This leads to significant uncertainty in the forecasting. On the contrary, the data-driven based technique do not require physical detailed models and the forecasting process is performed based on historical data only.

Plenty of regressive models exist in short-term building load consumption forecasting. They are developed in four main categories [START_REF] Yildiz | A review and analysis of regression and machine learning models on commercial building electricity load forecasting[END_REF]:

1) Support vector regression (SVR); 2) Artificial neural networks (ANN): Bidirectional fuzzy neural network (BFNN), feed forward neural network (FFNN), multilayer perceptron (MLP), long short-term memory (LSTM); 3) Decision tree (DT): Random forest (RF), classification and regression trees (CART); 4) Statistics methods: Multiple linear regression (MLR), auto-regressive integrated moving average (ARIMA), knearest neighbor (kNN). Statistic and SVR techniques are widely used to find out the most accurate models for short-term load forecasting. The computations are realized with real data-sets with hourly resolution, but the results are not usually optimistic. The performances of SVR techniques are generally better than statistic ones, as illustrated in [START_REF] Edwards | Predicting future hourly residential electrical consumption: A machine learning case study[END_REF], [START_REF] Massana | Shortterm load forecasting in a non-residential building contrasting models and attributes[END_REF] for instance. ANN, such as FFNN and BFNN, are also widely implemented, though machine learning levels with one hidden layer only present an accuracy around 70 % [START_REF] Yun | Building hourly thermal load prediction using an indexed ARX model[END_REF].

Data-sets in STLF present a significant variation, which increases the uncertainty of the predicted values. In recent work, ANN techniques with deep learning models were showing promising results. However, their performances on benchmark data was still to improve. Moreover, the size of the training data-set and the influencing features have to be adapted to such a deep learning based model for relevant results (not too small and if possible adaptive) [START_REF] Wang | LSTMbased short-term load forecasting for building electricity consumption[END_REF].

In this paper, we focus on comparing and analyzing the most recent techniques with the highest accuracy. These methods are SVN, RF, MLP, LSTM and an attention-based LSTM. The objective is to identify the best model for the STLF of smart buildings and local energy communities (or smart districts) and evaluate factors impacting the accuracy of the prediction, with in mind a future usage of the data in an energy management system, i.e. based on an optimization loop. These models are trained on a real data-set with and without taking into account of weather conditions and time-related features.

The paper is structured as follow. First, the general theoretical models of the five techniques are presented. The analysis of the features impacting the predicted values and the choice of the parameters to train the models are explained in Section II. Then, the data-set used to compare the forecasting techniques is describe in Section III. Next, the discussion on the experimental results is presented in Section IV. Finally, the conclusions is drawn at the end of the paper.

II. FORECASTING TECHNIQUES

A. Definitions and performance metrics

The generic regression model for energy forecasting based on time series analytic is expressed in [START_REF] Yildiz | A review and analysis of regression and machine learning models on commercial building electricity load forecasting[END_REF], where {x t-T +1 , . . . , x t } is the feature of d-dimensional time series in the n th moving window, including T time steps before (t + 1), and ŷt+1 ∈ R is the prediction of the output values at time (t + 1). ŷt+1 = F (x t-T +1 , . . . , x t )

The objective function is to minimize the distance D between the true values y and the predictive values ŷ, as expressed in [START_REF] Feinberg | Load forecasting[END_REF], where t can be any value or time period, and D is any distance.

We evaluate the accuracy of the techniques with two criteria: the mean absolute percentage error (MAPE) and the root mean squared error (RMSE).

MAPE = 1 N N i=1 |y i -ŷi | y i × 100 % (3) RSME = N i=1 (y i -ŷi ) 2 N ( 4
)
The MAPE determines the average percentage of error per prediction, defined in [START_REF] Mocanu | Optimized parameter selection for assessing building energy efficiency[END_REF], and the RMSE focuses on the deviation between each predicted value and its corresponding true value, defined in (4), with N the number of elements of the time series.

B. Implemented algorithms 1) Multilayer perceptron: MLP is a model of feed-forward neural network that is composed of at least three layers: input, hidden and output layers. Information are moved only one direction, "forward", i.e. from the input to the output, as shown in Fig. 1. The input layer is used to enter the input information, {x t-T +1 , . . . , x t }, thus the number of unit is a product of the number of features and the width of the window. The output layer is only a single unit: the output value y t+1 . In the hidden layers, the nodes rely on a nonlinear activation function (e.g. relu, tanh or sigmod) to map the weighted inputs to the output of nodes. A 2-hidden layers MLP is formulated in [START_REF] Massana | Shortterm load forecasting in a non-residential building contrasting models and attributes[END_REF], where W corresponds to a weight matrix and b is a bias parameter vector.

-1 +1 ˆ 1 1 2 2 . . . . . .
                   Z 0 = flatten(x t-T +1 , . . . , x t ) Z 1 = W 0 .Z 0 + b 0 A 1 = f 1 (Z 1 ) Z 2 = W 1 .Z 1 + b 1 A 2 = f 2 (Z 2 ) ŷt+1 = W 2 .A 2 + b 2 (5)
"flatten" is an operator that creates a vector with one row from a matrix. Deep MLP implementations present more than one hidden layer, which is more advantageous to learn complex functions F [START_REF] Kuo | A high precision artificial neural networks model for short-term energy load forecasting[END_REF]. However, the subsequent increase in parameters requires a larger data-set to train the model.

2) Long short-term memory: LSTM is an artificial recurrent neural network (RNN) in which the output from previous step is fed as input to the current step, relying namely on a "feedback" connection. This network is capable to learn longterm dependencies and remember information for prolonged periods of time with the help of a memory cell (i.e. the previous time step). The architecture of LSTM is a chain of cells in which the information can be removed or added using gates, as illustrated in Fig. 2.

For the t th cell, the input includes the current input information x t , the previous information of the state of the cell h t-1 and the memory cell c t-1 . The information flows through three gates: the update gate Γ u , the forget gates Γ f (used to update the cell memory c t ) and the output gate Γ o (used to update the state cell h t-1 ). The LSTM principle is defined in [START_REF] Yun | Building hourly thermal load prediction using an indexed ARX model[END_REF], where σ and tanh are the activation function. LSTM is a special type of RNN that solves the problem of the vanishing gradient in the training process and that captures much more long-term dependencies, justifying its preferred selection in the last few years. Furthermore, deep LSTM, stacking at least two LSTM layers as presented in Fig. 3, addresses correctly the estimation of complex functions, as in our case load consumption in smart buildings.

                       ct = tanh(W c (h t-1 , x t ) + b c ) Γ u = σ(W u (h t-1 , x t ) + b u )) Γ f = σ(W f (h t-1 , x t + b f )) Γ o = σ(W o (h t-1 , x t + b o )) c t = Γ u .c t + Γ f .c t-1 h t = Γ o . tanh(c t ) ŷ = W s .h t + b s (6) where W c , W u , W f , W o , W
ˆ Cell Cell Cell Cell Cell Cell ... ... -+1 ℎ -+1 -+2 -1 -+2 ℎ -+2 ℎ -+1 -+1 -+2
ℎ′ -+1 ℎ′ -+2 ℎ′ Fig. 3. Principle of the long short-term memory implementation with more than two layers.

3) Attention-based long short-term memory: At-LSTM is based on the attention mechanism, a part of the Seq2Seq model [START_REF] Sehovac | Deep learning for load forecasting: Sequence to sequence recurrent neural networks with attention[END_REF], which has demonstrated interesting qualities in multiple domains, particularly in language translation because of its capacity to capture important part of the input sequence. This improves the correlation between the input and the output data so that the model prediction is more accurate. Technically, the attention mechanism is applied to the top layer of an LSTM, as presented in Fig. 4. The implementation of the At-LSTM is presented in [START_REF] Wang | LSTMbased short-term load forecasting for building electricity consumption[END_REF], where H is the features extracted by the prediction model at the current step t, such as H = [h t-T +1 , ..., h t ]; e is the content-based function of H; α denotes the attention weights of the features of H; r is the context vector of the attention layer; h is the attention vector that produces the predicted output and V α is a parameter learned with the rest of the system.

           ht = tanh(W c (r t , h t )) r t = i=1 T α ti .h i α ti = exp(eti) i=1 n exp(eti) e ti = V T a tanh(W a (h t , H)) (7)
4) Random forest: RF represents a type of learning methods that combine prediction from multiple models to increase accuracy and control over-fitting. RF techniques aggregate many DT, using bagging. An individual tree in this technique is built on a random subset of the data-set, with random features that can be split on at each node, as presented in Fig. 5. The average prediction from all the individual regression trees helps improving the predictive accuracy. Hence, RF is efficient on large data-sets.

5) Support vector regression: SVR positions an optimal hyperplane by maximizing, under constraints, the distance between the predicted value ŷt+1 and the desired output y t+1 . The principle of this method is illustrated in Fig. 6 and expressed in [START_REF] Kuo | A high precision artificial neural networks model for short-term energy load forecasting[END_REF], where w is the normal vector of the hyperplane weight vector, ε a margin of error, and b a bias parameter.

min 1 2 w 2 such that | y t+1 -w, flatten(x t-T +1 , . . . , x t )) -b |≤ ε, ∀t (8)
+ Φ( ) -Φ( ) 0 Φ( ) + - 0 Φ -1
Fig. 6. Principle of the Support vector regression technique.

The standard SVR implementation is linear. Using alternative kernel functions (like the polynomial (poly), the Gaussian radial basis function (RBF), or the Sigmoid) to replace the inner product ., . in [START_REF] Kuo | A high precision artificial neural networks model for short-term energy load forecasting[END_REF], creates a nonlinear SVR, much more convenient to analysis complex data. Furthermore, a soft margin can be employed by adding slack variables ε * t , ε t ≥ 0 to relax the optimization as expressed in [START_REF] Sehovac | Deep learning for load forecasting: Sequence to sequence recurrent neural networks with attention[END_REF], where Φ is a kernel function and C is a regularization parameter.

min 1 2 w 2 +C t ε * t + ε t such that y t+1 -Φ(w, flatten(x t-T +1 , . . . , x t )) -b ≤ ε t + ε * t Φ(w, flatten(x t-T +1 , . . . , x t )) + b -y ≤ ε t + ε * t (9)
III. DATA-SET Many parameters influence the energy consumption, such as the ambient dry bulb temperature, humidity, solar radiance, dew point, cloud cover, wind speed, pressure or scheduling (weekdays/weekend/holidays) and so on [START_REF] Yildiz | A review and analysis of regression and machine learning models on commercial building electricity load forecasting[END_REF]. Considering influencing factors with the energy consumption leads to an increase of the training data-set size, and then elongate the training time. Energy consumption trend is drawn from the same figure with these parameters. Discarding factors presenting lower or identical correlation coefficients (multico-linear) is recommended at this step. Energy consumption data, collected from smart meters of 5 567 households in London, and historical weather data during more than 2 year, from November 2012 to February 2014, are available in Kaggle [START_REF]Smart meters in London[END_REF]. The dew point and UV index are correlated with the temperature and pressure, and the moon phase present insignificant correlation with energy, so they are discarded. Thus, the exogenous parameters are reduced to the temperature, the wind speed, the humidity and time scheduling data.

The correlation of weather parameters is shown in Fig. 7. The temperature and humidity present an inverse correlation to energy. The wind speed presents a low correlation with energy. We also consider two features for the time scheduling: the hour of the day (from 0 to 23) and the day of the week (0 for weekday or 1 for holiday and weekend).

IV. RESULTS AND DISCUSSION

A. Scenarios and modeling

The prediction of hourly energy consumption is performed following four scenarios:

1) All buildings without influencing parameters; 2) All buildings with influencing parameters; 3) One private house without influencing parameters; 4) One private house with influencing parameters. The idea is to compare the performances of the algorithms facing aggregated energy consumption patterns (for instance at the level of a local energy community or a smart district) or directly raw energy consumption curves, with much more variability (at the level of a smart household). Also, considering influencing parameters impacts the needed sensors information, thus indirectly the practical acceptability, replicability, and in the end the installation cost as well as the computational tractability. Those considerations define effectively four scenarios for the comparison, knowing that the usage of the forecasting algorithm lies within an optimization loop of an energy management system (with technical and economic constraints). This optimization layer will probably influence the performance metrics, without being directly considered in this paper.

Training data is aimed to predict the energy for the weekahead. The energy consumption readings from the 1 st to the 7 th of January 2014 are used as test data, the rest is used for training and validation. Note that this week is randomly chosen. Indeed, for comparison between the technologies, we should focus on the manipulations on the same data-set, independently of it it is winter or summer.

Regression neural network (MLP and LSTM) are implemented with two neutron layers, each one built of 32 units. The mean absolute error (MAE) is used as loss function because its gradient varies significantly more than the mean square error (MSE) throughout iterations, as shown in Fig. 8. The Adam algorithm was adopted as the gradient-based optimizer because of its fast convergence rate and lower error ratios. Besides, the mean absolute values are steadily decreasing during training periods, so using dropout layer was ignored. The learning rate is automatically set by Adam.

A single layer LSTM was used for the At-LSTM technique to better compare with the other techniques as well as to insure the best results for this paper.

The RF was implemented to create trees using random selections of features and samples. The random_state parameter allows controlling these random choices. It was fixed throughout the cases study to compare adequately the results after the executions. The SVR used a 4 th order polynomial function as kernel. Data in the past 24 hours (T =24) was used to predict one hour-ahead of energy consumption. As supplementary data, the detailed implementations of the four scenarios are presented in [START_REF]Detail implementation of four scenarios[END_REF].

B. Comparison of the implemented techniques

The result of the prediction of the five implemented techniques in the case of the hourly average load consumption prediction of all buildings and only a single private household are shown in Fig. 9.

With the same data-set, the accuracy of the average energy consumption prediction error MAPE is high for the whole set of houses, reaching up to 96 %. Meanwhile, the prediction for the private house is less accurate, simply due to the increased variability of the load curve. The value of the two performance criteria (MAPE and RMSE) for the four scenario and the five implemented techniques are summarized in Table I. For hourly predictions, the consideration of influencing parameters during the training period do not seems to be adaptive as the accuracy for scenarios 1 is higher (MAPE: ∼4 % vs ∼7 % respectively). This could be linked to the chosen influencing parameters or a de-synchronization in their time stamp. The size of the data-set or the chosen period (hourahead prediction) could be an explanation of this counterintuitive result to be more deeply assessed.

We can say that Deep LSTM is a good technique for load prediction without influencing parameters as its MAPE is the smallest. Meanwhile, At-LSTM shows a better performance in case the influencing parameters are taken into account. RF provides comparable results to LSTM or At-LSTM and is better performing if you consider the RMSE values as the variation of energy consumption is important. The SVR is the worst technique in the four scenarios in this experience. Its MAPE and RMSE values are the highest and its inaccuracy (MAPE) is up to more than 34 % while the At-LSTM ones goes down to ∼23 %. For weekends, its prediction seems to be more accurate (not shown in this paper).

The computational tractability of the five implemented techniques is shown in Fig. 10. The kernel of SVR allows constituting a optimal hyperplane which aims to increase the accuracy of the prediction. In this issue, the training data-set is used with various parameters, the linear kernel can not give a better result than the nonlinear kernel. However, the chosen nonlinear kernel of the SVR makes its training time much longer than the linear one. Besides, MLP, and RF present moderate training times because they perform simple computations (e,g., the MLP model contains many units in which only linear operators are performed). Finally, LSTM presents the longest training time as it requires memory cells and more complex computations in each unit.

V. CONCLUSION

In this paper, we describe and implement five state of art machine learning models and applied them to short-term load forecasting to compare and analyze their performances on two relevant criteria. The idea is to help determine which one could be considered in the optimization loop of an energy management systems, either at the level of a district (local energy communities, or smart districts) or directly at the level of individual households (smart residential buildings). The five implemented techniques are deep multilayer perceptron (deep MLP), deep long short term memory (deep LSTM), attentionbased LSTM, support vector regression (SVR), and random forest (RF), all of them chosen as the most efficient forecasting algorithm of the state of art for energy consumption prediction.

The comparison of the forecasting algorithms is presented with and without the consideration of influencing parameters (weather data for instance), based on real data (relying on a large training data-set). Two level of energy consumption are assessed: the single private household in which the variation of data is significant and the whole data-set, presenting a smart district or a local energy community, in which the variation of the data is much smoother. The prediction without influencing parameters seems to be more accurate in the considered scenarios, which is counter-intuitive at first, but interesting from the point of view of its application to energy management systems. The prediction on the single household data-set is promising, the mean absolute percentage error of the LSTM technique can be less than 4%. The RF, At-LSTM and deepLSTM are more efficient than the deep MLP or the SVR. To conclude, in the context of energy management systems, the At-LSTM is better fitting in cases where influencing parameters are taken into account. Otherwise, deep LSTM is the most interesting, if the long training time is not a problem.

The techniques compared in this paper represent the state of art of STLF in the energy context. To select one for an energy management system, sensitivity studies should be conducted on the training data-set duration as well as its impact on seasonal predictions. Also, a real-time approach could be enforced to improve accuracy in contexts where low-latency online predictions are needed.

  min

  t ; ŷt )

Fig. 1 .

 1 Fig. 1. Principle of the multilayer perceptron with two hidden layers.

Fig. 2 .

 2 Fig. 2. Cells description of the long short-term memory implementation.

  s are the weight matrices parameters, b c , b u , b f , b o , b s are the bias vector parameters, and ct is a new candidate of cell state.

Fig. 4 .

 4 Fig. 4. Attention-based long short-term memory implementation.

Fig. 5 .

 5 Fig. 5. Illustration of a sample random forest tree.

Fig. 7 .

 7 Fig. 7. Correlations of weather parameters and energy in the data-set [10].

Fig. 8 .

 8 Fig. 8. Gradient of the mean square error (MSE) and the mean absolute error (MAE) as a function of the number of the number of iterations.

Scenario 1 :Scenario 2 :Scenario 3 :Scenario 4 :

 1234 All buildings without influencing parameters. All buildings with influencing parameters. One private house without influencing parameters. One private house with influencing parameters.

Fig. 9 .

 9 Fig. 9. Load forecast from the 1 st to the 7 th of January 2014.

Fig. 10 .

 10 Fig. 10. Execution time of the five implemented techniques.

2014-01-01 2014-01-02 2014-01-03 2014-01-04 2014-01-05 2014-01-06 2014-01-07 2014-01-08 Time