Chaoyun Wang
email: c.wang]@efficacity.com

Damien Faille
email: damien.faille@edf.fr

Lilia Galai-Dol
email: [l.galai-dol]@efficacity.com

Gilney Damm
email: gilney.damm@centralesupelec.fr

Optimal Power Flow for a Multi-Energy Vector MicroGrid

Keywords: Optimal operation and control of power systems, Control of renewable energy resources, Smart-Cities, Renewable Energy System Modeling and Integration, Distributed Energy Resource Integration and Coordination

The paper introduces a combined modeling that can allow the Optimal Power Flow of a Multi-Energy Vector merging electric and thermal networks. This topic is motivated by the goal of introducing larger shares of renewable energy in the overall energy consumption, which is hindered by the intermittent nature of these renewable energies. To tackle this problem, it is acknowledged that almost half energy consumption in urban environments is used for thermic objectives. For this reason, a possible solution is to use in parallel electric and thermic networks, called the Multi-Energy Vector. It is then first proposed a combined model, and then aiming at cost optimization, it is performed an optimal power flow, such as to supply references for lower-level controllers, attaining the desired objectives, and keeping all variables inside their operational margins. The proposed scheme is applied to both electric and thermal networks of the under-construction Moulon quarter of Paris-Saclay University. The simulation results illustrate the economic gains that can be made with this joint operation.

INTRODUCTION

Current environmental concerns have motivated a large number of works in electric grids, and in particular the concept of MicroGrids. MicroGrids are portions of the electric grid that may be isolated or not from the main electric network. The main objective of such MicroGrids is to attain balance on produced and consumed energy, while keeping all variables inside operational bounds. Produced energy within environmental objectives is mostly produced by renewable energy sources (renewables). Such balance may be reached in all time or in average over a time span. The difference being that in the first the grid could work in island mode (disconnected to the main grid) while in the second case, this main grid is supplying or absorbing the lack/excess production.

Either case is difficult to fulfil while relying on renewables. Their intermittent nature makes balance difficult while respecting constraints. It is then necessary to use large amounts of storage, or to oversize the electric grid to relax constraints. At the same time, it is important to remark that roughly 40% of consumption, for residential and commercial buildings, is used for heating/cooling. In some countries this amount may be even larger. In addition, thermal dynamics can provide a natural cheaper storage, so if one works with both electric and thermal energies, it would be possible to attain the same objectives of consumption and well-being, in a more optimized way, and it would be easier to fulfil grids' constraints, the socalled multi-energy vector.

There are few articles discussing combined thermal and electrical models. A series of works, in particular [START_REF] Geidl | Optimal power flow of multiple energy carriers[END_REF] and [START_REF] Geidl | Integrated Modeling and Optimization of Multi-Carrier Energy Systems[END_REF] proposed the notion energy hub. Each energy hub contains three basic elements: direct connections, converters, and storage. All kinds of real facilities can be modelled as energy hubs. Then the model of an energy hub is designed, as well as the transmission networks, for both electricity transmission lines and heating transmission pipelines, with a set of nonlinear equations. [START_REF] Liu | Combined analysis of electricity and heat networks[END_REF] and [START_REF] Liu | Combined analysis of electricity and heat networks[END_REF] proposes a combined analysis of integrated networks. The author made at first separated then combined analysis of electrical power flow and district heating network, then both systems are combined in large state vectors, and the overall system is linearized and the analysis is carried out solving in time its Jacobian matrix with the Newton-Raphson algorithm. Others have focused on the mix electricity and gas, [START_REF] Fang | Dynamic optimal energy flow in the integrated natural gas and electrical power systems[END_REF] and Martínez-Mares and Fuerte-Esquivel (2012) in a similar manner.

In the present work, similarly to the works above, it is presented a combined model for an electrical and a thermal grid. The main contributions of the present paper can be stated as the development of this model, followed by a combined analysis and, as a step forward from the literature, these results are used for creating an Optimal Power Flow calculation [START_REF] Chatzivasileiadis | Cartographie des réseaux aériens d'Enedis[END_REF] such as to minimize costs and power losses, and respect limits in power flows. In the best knowledge of the authors, this is the first work proposing an optimal power flow for combined Electro-Thermal grids. These first results can also be used to maximize renewable energy use, respecting the grid constraints, and also to plan peak reduction in overload times of the day. In the future these models will be enhanced aiming at using predictive controllers such as to improve the optimal use of the combined grid.

The proposed results are then applied to the now under construction Paris-Saclay University's Moulon quarter. This choice was motivated by two reasons. First, because information is available for these networks under construction which is not always the case for existing installations. The second is because if results are concluding, they should be implemented in this grid in the future.

ELECTRICAL NETWORK

The electrical grid consists of buses and branches. A bus can be a power source or a load. The buses are connected by cables or overhead lines, which are called branches.

In an N-bus system, a state vector consists of the phase angle θ i , the voltage amplitude V i , active power P i and reactive power Q i at any bus of the system. Basically, the buses are classified as:

• PV bus: when it's connected to a generator, then active power P i and voltage amplitude V i are known and the unknown variables are reactive power Q i and phase angle θ i ; • PQ bus: when it's connected to a load, then active and reactive powers (P i and Q i ) are known, and phase angle θ i and magnitude of the voltage V i at the bus are to be determined; • r (reference) bus: one PV bus is taken as a reference bus which is also called V θ bus. The phase angle θ i is fixed at 0, and the magnitude of the voltage V i is supposed to be known, while the active and reactive power (

P i and Q i ) become unknown We define θ = [θ 1 , θ 2 , ..., θ N ] T , V = [V 1 , V 2 , ..., V N ] T , P = [P 1 , P 2 , ..., P N ] T and Q = [Q 1 , Q 2 , ..., Q N ]
T , four column vectors with N elements for each, to describe the state of buses.

Vectors of bus currents and bus voltages are linked through the algebraic relationship (see [START_REF] Chatzivasileiadis | Cartographie des réseaux aériens d'Enedis[END_REF] and [START_REF] Hennebel | Lecture notes Fonctionnement et réglage des systèmes de transport et de distribution d'électricité en régime permanent[END_REF]):

I bus = Y bus V bus
For a N bus system, V bus is a vector of the complex voltage value at the buses, and the bus admittance matrix Y bus is:

Y bus =     y 11 y 12 • • • y 1N y 21 y 22 • • • y 2N . . . . . . . . . . . . y N 1 y N 2 • • • y N N    
In the same way, the line admittance matrix Y line links the bus voltages to current flows:

I line = Y line V bus
The key problem of power flow calculation is solving a set of non-linear equations because the power is related to the product of voltage and current. There are many algorithms to solve a set of non-linear equations like Gauss-Seidel and Newton-Raphson (N-R).

The injections of active and reactive powers in the network depend on the phase angles and voltage amplitudes of all N nodes. The power injections can be written as:

P i (θ, V ) = V i j V j (G ij cosθ ij + B ij sinθ ij ) Q i (θ, V ) = V i j V j (G ij sinθ ij -B ij cosθ ij ) (1)
where

θ ij = θ i -θ j , G ij and B ij come from the admitance matrix Y bus,ij = G ij + jB ij .
The optimal power flow (OPF) calculation aims to find an optimal dispatch of an electrical system. The optimization process that determines the operation of the least-cost available generators, given (a) the total electric demand, (b) the minimum and maximum operation limits of each generator, and (c) power flows in the network. It can also minimize other parameters like the losses of transmission lines.

In OPF problem, the state vector is x = [θ V ] T , and the control vector is ũ = [P g Q g ] T , so the total sizing vector is

x = [θ V P g Q g ] T .
The OPF problem consists of two parts: cost function and constraints. In the case of minimization of costs, the objective function minimizes the total costs for generating active power:

min x c T P g (2)
where c T is the column vector of the generation price of each bus. In the case of minimization of active and reactive power losses, the objective function would be:

min i,j S i→j + S j→i (3)
subject to:

S g -S d -diag(V )Y * bus,row-i V * = 0 (4)
In the same time, the sizing parameters should be limited in an operation range:

x min ≤ x ≤ x max

THERMAL NETWORK

A thermal network, also called district heating (DH) network, is a system that distributes the heat generated at the sources to all the loads within an area. The heat is transmitted by water, so basically, a DH network is analyzed at first by a hydraulic model, and then by a thermal model. The analysis with these two models can be done in a decomposed method, which means that the system is analyzed step by step, or in an integrated method, which means that the individual hydraulic and thermal (H-T) analyses are combined in one step.

We consider a simplified model of a DH network for which a section is described in Figure 1. The supply line transports the water from the source to the loads. Each load i extracts a fraction of the flow ṁqi at the temperature T s,i to heat the load. The water at the output T o,i is mixed with the water in the return line resulting in a temperature T r,i . The flow is supposed to be imposed by local controllers actuating on the pumps and valves. The model represents the mass and enthalpy balance equations.

The enthalpy depends mainly on the temperature H = c p T , where c p is the specific heat of water kJ/kg/ o C.

The model considers two kind of losses. Pressure or head loss due to the friction in the pipe ant heat loss between the pipe and the environment. We define the difference with the ambient temperature T a by the expression:

T = T -T a
Basically, the model consists of six sets of equations:

• Continuity equations: A ṁ = ṁq • Loop pressure equations: Bh f = 0 • Head loss equations: h f = K. * ṁ. * | ṁ| • Heat power equations: ϕ = c p ṁq . * (T s -T o ) • Supply temperature equations: C s T s = b s • Return temperature equations: C r T r = b r
where ṁ is a vector whose components are the mass flow rate in pipes, ṁq is the mass flow rate at the source and loads. The topology network is defined in the two first equations by the A and B matrices which correspond to the flow balance equations for the nodes and the pressure loss equations for loops. The third equation is the relation between the flow and the head losses along the pipes. The pipe friction coefficients vector K is component wise multiplied (. * Schur product) by the square of the mass flow rate vector. The heat power equation links the temperature with the thermal power temperature vectors defined by

T o = [T o,1 T o,2 • • • T o,N ] T and T s = [T s,1 T s,2 • • • T s,N ] T .
The supply and return temperature are calculated by energy balance that can be reformulated by linear equations where C s and C r depend of the network topology and the flow in the pipes.

In the case of Fig. 1, the steady state equation for the supply and return temperature are given by the following expressions where Ψ j is the heat loss in the pipe j that depends on the flow, the heat exchange coefficient (supposed constant) and the length of the pipe.

T s,i+1 = Ψ j T s,i [ ṁj-1 -ṁj Ψ j ] = T r,i T r,i+1 = ṁqi T o,i (5)
Ψ j = e - λL j cp ṁj
The six sets of equations can be reformulated as a computation of a zero of a nonlinear function, the total mismatch ∆F which consists of the heat power equation the loop pressure equation, supply temperature equations and return temperature equations:

∆F =    ∆ϕ ∆h f ∆T s ∆T r    =    c p A ṁ. * (T s -T o ) -ϕ B(K. * ṁ. * | ṁ| C s T s -b s C r T r -b r    (6)
These equations build an integrated H-T model and are used in the derivation of N-R algorithm with the iterative form:

x k+1 = x k -J -1 k ∆F x k (7)
where k is the iteration; J k the Jacobian matrix. The state vector x contains all the unknowns: ṁ for all pipes, T s and T r for all load and source nodes:

x =   ṁ T s T r   (8)

Optimal hydraulic-thermal model calculation

The objective of an optimal hydraulic-thermal (OH-T) calculation is to find an optimal dispatch of a district heating system. The optimized criterion can be the total cost of heat generation and the total cost of electricity consumption by the pumps among others. Just like optimal power flow (OPF) calculation in the electrical network, this is also an optimization problem with nonlinear constraints.

Similar as integrated H-T calculation, the unknown variables comprise ṁ, T s,load , T r,load . The difference is that T s,source become no longer known variables, but the results of optimization. More generally, we use the heat power of the source nodes ϕ source as the sizing parameters and define their upper and lower bounds. Thus, the total sizing vector should be x = [ ṁ T s T r ϕ source ] T .

Fig. 2. Flow chart of combined electro-thermal model 4. OPTIMAL COMBINED ELECTRO-THERMAL NETWORK

In real application, the electrical grid and DH network are not independent. There are many coupling components (e.g., CHP units, heat pumps, electric boilers and circulation pumps). In this paper, the coupling components are the circulation pump that compensates the friction of the networks pipes and thermo-refrigerating heat pump (TFP) located at the station. In our study the TFP is a device that can convert electric power in heat or cold with a given efficiency.

The calculation of the hybrid thermal and electrical power flow is done mainly in three steps as described in the Fig. 2 where the known and computed variables are defined for each stage. First is the initialization of system: the known values, topology of network, linkage nodes between electrical and thermal network are defined at first. Then, the problem is solved in an iterative way. Finally, the results are put out.

In the paper we also propose an Optimal Power Flow for the combined electrical and thermal network. The objective of an optimal combined electro-thermal calculation is to find an optimal dispatch of a combined model. In this part, key parameters linking the thermal and electrical networks have been connected for the provision of flexibility services. The objective was to study the suitability of such modelling for control and optimization of power, while fulfilling grid constraints.

The optimization was carried out using the Matlab Optimization toolbox. Future works will be focused on the comparison of optimization methods, such as Kuhn Tucker, to verify which would be better fit for this problem.

In the combined optimization, the sizing vector will be:

x = [θ V P g Q g ṁ T s T s,load φ source ] T
In a combined optimization, the objective is to minimize the total cost, including the electrical generation and thermal generation. Therefore, the cost function is: where c e is the vector of the cost of electrical generation for each electrical bus, c t is the vector of the cost of thermal generation for each thermal source node.

min x c T e P g + c t φ source (9)
For a radial DH network, the constraints are:

c p A ṁ. * (T s -T o ) -φ = 0 B(K. * ṁ. * | ṁ|) = 0 C s T s,load -b s = 0 C r T r,load -b r = 0 (10)
For an electrical network, the constraint is:

S g -S d -diag(V )Y * bus,row-i V * = 0 (11)
where S g and S d represents the complex power of the generator and demand. The combined optimization should satisfy all these constraints. The difference is that, in the combined analysis, the electrical power demand S d includes not only the consumption of each load bus, but also of the circulation pumps and of the TFPs:

S d = S d,load + S pump + S T F P
Besides, the upper and lower bounds are kept:

x min ≤ x ≤ x max

PARIS-SACLAY'S MOULON QUARTER

Modeling of Electrical Network

Paris Saclay is a district in construction in the south of Paris for which a high energy efficiency is aimed. In this section we explain how we build a realistic model of the district.

To calculate load powers on this new quarter, we have the area and type of each building and data of consumption per square meter for each type of building.

Then, we correspond the positions of the source stations and the substations obtained with the data from the local distribution system operator ENEDIS (see Enedis ( 2019)).

We consider that each group of living building is fed by a Medium to Low Voltage station, and each large tertiary or education building is fed by an independent sub-station, Figure 3. Here the Central Installation is the heat source and District sub station corresponds to the substation with the TFP. The image of Google Maps is shown in Figure 4, from which we can derive the electric grid structure as shown in Figure 5 The Electrical/Thermal coupling components are the circulation pump in the heat source and TFPs in the substation. The coupling is calculated as introduced above. Since the heat power consumption of each load node is fixed, P T F P , the electrical active power consumed by the TFPs, is known.

Optimization of Thermal Network

For the thermal network in Moulon, the thermal power comes from two parts: geothermal power P ther,geo and gas power P ther,gas , as shown in Figure 6.

The geothermal energy can be considered with no cost, so we will focus on gas consumption. The geothermal heating power (in MWt) can be calculated as: Thus, P ther,gas can be calculated as follows:

P ther,gas = φ source -P ther,geo

where φ source is the total heat power of the source node (MWt).

In the thermal optimization, it is necessary to minimize not only the thermal generation, but also the electrical generation. That is to say, there are 2 variables to minimize: electrical power consumption of circulation pump P elec,pump and thermal power generation with gas P ther,gas . In order to solve this multi-objective optimization problem, a scalar parameter a is used. Thus, the cost function should be:

min x aP elec,pump + (1 -a)P ther,gas (14)
where 0 < a < 1. When a is set to 0 (to 1) we minimize only the electric consumption of the pump (gas consumption). A compromise is obtained with intermediary value of a. For smaller a, the weight of P elec,pump is lighter than on P ther,gas , that will lead to smaller value of P ther,gas and a larger value for P elec,pump . Two values of parameter a are tried: a = 0.5 and a = 0.98.

where both P g and P source,gas are power generation, and tariffs c e and c t are for electricity and gas respectively.

Taking c e = 49.50Euro/M W he (see Selectra ( 2019)), and c t = 63.61Euro/M W ht (see Jechange ( 2019)), the result of the optimization is a total cost of 2184.81Euro/h.

CONCLUSIONS

In this paper it is proposed a Multi-Energy Vector Optimal Power Flow merging electric and thermal networks. Their combined model is derived, and then it is proposed a modelization adapted to a cost optimization for both energies at the same time. The proposed scheme is applied to the grid of the under construction Moulon quarter of the Paris-Saclay University. The results show how the economic gains can be optimized with this joint operation and will be implemented to the network when completed. This approach will be enhanced in future works to also include a gas network. In addition, this configuration is well fitted for Model Predictive Control (MPC) approach in order to provide a predictive optimal power flow. Future works will develop a MPC framework such as to optimize the combined network with different time horizons. In the same way, different optimization methods, such as Kuhn Tucker, will be tested to select the most adapted for the Optimal Power Flow for Combined Electro-Thermal grids.

The main idea is first to cope with future variable prices during the day, and second to use the thermal inertia as energy storage, such as to reduce electric power peaks while fulfilling grid constraints.

Fig. 1 .

 1 Fig. 1. District heating Network

Fig. 3 .

 3 Fig. 3. Electrical network modeling on Google Maps of Moulon

Fig. 4 .

 4 Fig. 4. Thermal network in Moulon

Fig. 6 .

 6 Fig. 6. Principle of thermal source in Moulon P ther,geo = ṁq c p (T s -T r ) (12) where ṁq is the mass flow rate in the node (kg/s), c p is the specific heat capacity M J/kg/ o C, T s is the supply temperature of geothermal center (30 o C in the figure) T r is its return temperature (10 o C in the figure, in the following calculation T r is a result of optimization).

Table 1 .

 1 The results are shown in Table1.H-T optimization of MoulonAnother way to transform the multi-objective problem is to minimize the cost taking into account the prices. This is a partial optimisation in the sense that only the economic cost is taken into account and no other aspects like the CO 2 emission or use of renewable energy. This combined optimization is obtained by an economic cost function in Moulon as: min

	a	0.5	0.98
	P elec,pump	101 kWe	41 kWe
	P ther,gas	17511 kWt 19490 kWt


Fengdi Li

Giao T M Nguyen

Cédric Vancaeyzeele

Frédéric Vidal

Cédric Plesse

Healable ionoelastomer based on poly(ionic liquid)

Keywords: poly(ionic liquid), ionoelastomer, vitrimer

ABSTRACT: Ionoelastomers are elastomer based on poly(ionic liquid) (PIL) and have been introduced recently as liquid-free and stretchable ionic conductors. [1][2][3] These ionically conducting materials find their applications in all-solid stretchable and wearable electrochemical devices. In this presentation, we report a novel PIL based on imidazolium and bis(trifluoromethanesulfonyl)imide synthesized from commercially available copolymer rubber of polyepichlorohydrin and poly(ethyleneoxide). The PIL bearing allyl functional groups was then crosslinked with dynamic boronic ester crosslinkers through thiol-ene photoaddition. The resulting ionoelastomers are healable via vitrimer chemistry of boronic ester transesterification. The effect of crosslinker and additional poly(ethylene glycol) based dangling chains acting as free volume enhancer has been studied on thermomechanical properties and ionic conductivities. At optimal conditions, the resulting ionoelastomer is soft with a Young modulus of 0.2MPa, stretchable up to 300% elongation, and exhibits an ionic conductivity of 1.6  10 -5 S•cm -1 at 30 °C. In addition, the vitrimer properties of the material have been demonstrated with Arrhenius-type stress relaxation, healability, and recyclability. The recycled materials were able to retain their original mechanical and conductive properties. The reported material has allowed obtaining sustainable high-performance polymer materials with high ionic conductivity and by combing the facile processability of thermoplastics with the performance advantages of permanent crosslinked thermosets via vitrimer chemistry.

Scheme 1 :

 1 Scheme 1: Scheme for the formation of ionoelastomer and the dynamic exchange of boronic ester function