Rémy Rigo-Mariani

Voon Keck

Jan Ling

Maciejowski

A Clusterized Energy Management with Linearized Losses in the Presence of Multiple Types of Distributed Generation

Keywords: Distributed generation, MILP, CO2 emissions, energy storage, controllable load, network partitioning

This paper presents an optimal management (OM) strategy for distributed generation (DG) planning whose objective is the CO2 reduction for the power generation on the Jurong Island in Singapore. Different DG resources are investigated with solar generation, energy storage, small gas turbine as well as controllable loads in addition to the centralized generation already in site. Each of those resource is modeled in an optimal scheduling method that allow to test a great number of different DG configurations (i.e. type/size/site). The paper mainly focuses on the OM implementation that would allow avoiding prohibitive computational times. At first it is achieved thanks to a linearization of the line losses and a modified DC power flow. When simulating the system over a representative day. Then a generic clustering method is developed as well as a sequential optimal management (S-OM) with the zonal and nodal models is implemented. Different validation tests are performed as well as sets of simulation with different DG configurations. The optimal DG planning procedure will be part of further developments.

INTRODUCTION

The work presented in this paper focuses on reducing carbon emissions due to the power generation on Jurong Island in Singapore. That island hosts around 40 % percent of the total generation capacity in Singapore with thirteen identified Combined Cycle Gas Turbine (CCGT) units aggregated in four power plants [START_REF] Rigo-Mariani | A generic method to model CO2 emission Performances of combined cycle power plants for environmental unit commitment[END_REF][2]. A great part of the generated energy is consumed on the island by refineries, oil storage facilities and other petrochemical activities while the surplus is exported to the Singaporean mainland. The integration of distributed generation (DG) to the existing system is investigated as a way to reduce the emissions compared to the current base scenario in which only the centralized generation (i.e. the four power plants) is considered. The advantage of DG assets is to bring more diverse resources as well as more flexibility with various additional benefits such as loss reduction, improved reliability or the deferment of grid investments [START_REF] Prakash | Optimal sizing and siting techniques for distributed generation in distribution systems: A review[END_REF] [START_REF] Jordehi | Allocation of distributed generation units in electric power systems: A review[END_REF]. In such a framework typical DG planning problems aim at finding the best types, sizes and allocations of the assets that minimize the considered objective function (i.e. CO2 emissions for a representative day here) [START_REF] Paliwal | Planning of grid integrated distributed generators: A review of technology, objectives and techniques[END_REF]. The main concern when facing such problems is to avoid the computational burden while estimating the performances of many different DG configurations. The complexity of these DG planning studies can be disaggregated in three parts. Firstly, the numbers and types of the considered distributed resources directly impacts the amount of possible configurations to test. Simplest case studies investigate the placement of single generic PQ units with exhaustive search or analytical approaches [START_REF] Acharya | An analytical approach for DG allocation in primary distribution network[END_REF] [START_REF] Ghosh | Optimal sizing and placement of distributed generation in a network system[END_REF]. More complex problems consider multiple generic PQ units [START_REF] Nayeripour | Multi-objective placement and sizing of DGs in distribution networks ensuring transient stability using hybrid evolutionary algorithm[END_REF] or multiple types of resources with both renewable and storage units [START_REF] Prenc | Distributed generation allocation based on average daily load and power production curves[END_REF] [START_REF] Vigano | Methodology for the optimal siting and sizing of storage systems in distribution networks[END_REF]. This paper investigates different DG technologies, including photovoltaic generation (PV), energy storage (ST) and distributed CCGT units (DU). The presence of controllable loads (CL) is also studied which is not commonly done in the literature. The second degree of complexity for DG planning studies comes from the size of the modeled power system which obviously increases the number of possible DG configurations (i.e. types/sizes/sties). Typical systems encountered literature display a wide range for their number of buses from less than 10 [START_REF] Abbasi | Optimal DG allocation and sizing in presence of storage systems considering network configuration effects in distribution systems[END_REF] to over 100 [START_REF] Popovic | Placement of distributed generators and reclosers for distribution network security and reliability[END_REF] while the Jurong Island system is modeled here with a little more than 200 nodes. Thirdly, the difficulty to solve those planning problems is enhanced when longer time horizons are considered to compute the objective function. Indeed, many studies only refers to a single snapshot power flow with one load level, to compute the cost, losses or reliability index corresponding to every DG configuration estimated [START_REF] Abbasi | Optimal DG allocation and sizing in presence of storage systems considering network configuration effects in distribution systems[END_REF] [13] [START_REF] Vallem | Siting and sizing of distributed generation for optimal microgrid architecture[END_REF]. Some works run simulations over a whole year with successive independent AC optimal power flow (AC OPF), as in [START_REF] Ochoa | Distribution network capacity assessment: variable DG and Active networks[END_REF] [START_REF] Porkar | Optimal allocation of distributed generation using a two-stage multi-objective mixedinteger-nonlinear programming[END_REF]. However, when storage units are considered, as it the case in this paper, the introduction of time related constraints (e.g. the state of charge computation) implies that the whole time horizon needs to be represented in the optimization problems and successive AC OPF cannot be applied. Thus, authors often refer to simplified DC OPF with the extensive use of linearized formulations in order to avoid prohibitive computation times [START_REF] Pandžić | Near-optimal method for siting and sizing of distributed storage in a transmission network[END_REF] [START_REF] Stich | Enhancing the integration of renewables by trans-border electricity trade in ASEAN[END_REF] (See table in Appendix A for a summary of the reviewed studies). With time series profiles considered, DG planning studies can be seen as systematic optimization problems The goal is to find the best configuration for the DG assets (i.e. design loop) whose operation over a given time horizon (i.e. management loop) corresponds to the best values of the objective. Several approaches are identified in the literature, with different couplings between the management and design loops [START_REF] Fathy | On the coupling between the plant and controller optimization problems[END_REF]. Bi-level optimization is often used where the DG sizes/sites are part of an outer loop problem while the management is solved in an inner loop (Figure 1a). Typically, the design variables generated in the outer level are used as input parameters for the inner level with two integrated optimization problems. Such an approach is presented optimize DG management scope of this paper in [START_REF] Prenc | Distributed generation allocation based on average daily load and power production curves[END_REF] [START_REF] Ugrandi | Optimal wind turbine sizing to minimize energy loss[END_REF] [START_REF] Celli | Optimal integration of energy storage in distribution networks[END_REF] with the DG management over a day performed by dynamic programming or successive AC OPF while genetic algorithms try to find the best design in terms of DG sites/sizes. A second class of works refer to "all in one" approaches where both design and management variables are part of a single optimization problem (Figure 1b), usually using LP or MILP formulation [START_REF] Pandžić | Near-optimal method for siting and sizing of distributed storage in a transmission network[END_REF] [START_REF] Stich | Enhancing the integration of renewables by trans-border electricity trade in ASEAN[END_REF]. However, the convergence of the global problem might not be ensured in a reasonable time. The architecture considered in this paper corresponds to a classical iterative method which consists in several estimations of the management loop for different configurations and attempts to improve the design (Figure 1c). [10] [22] use such a method with successive AC OPF for the management strategy while different designed are generated with Monte Carlo simulations or an analytical approach.

This paper mainly focuses on the implementation of the DG management strategy within an iterative architecture for DG planning. Also note that the optimal operation of the assets is jointly considered with the schedule of the centralized resources on sites. DG integration will change the generation mix of the Jurong Island along with more flexibility provided by storage units and the ability to perform load control. Thus the management loop will be implemented in the form of a typical unit commitment problem (UC) with the existing power plants and the addition of DG resources whose connection buses and sizes are input parameters. Note that today the generations companies in Jurong Island bid in the national Singaporean market before the Power System Operator runs the UC [2]. In this work, the assumption is made that all the power plants and DG assets are centrally owned and managed by a vertical entity that could be the already existing "Jurong Island Community" for instance. A specific attention is paid to the computational time that has to be small enough in order to test a great number of DG configurations (i.e. types/sizes/sites). It is achieved here in two steps that are the main contributions of the proposed approach.

Firstly, losses within the lines are linearized using a modified DC OPF. Unlike conventional iterative approaches [START_REF] Chen | Jump and shift method for multi-objective optimization[END_REF] [24], the method developed here is faster and reference independent. It lies on the linear coefficients extracted from preliminary sets of simulations with several run of AC OPF and it is denoted as OM (Optimal Managemen). Secondly, a clustering approach allows fast computational times while dealing with a reduced number of variables. In DG planning studies, clustering methods are often used to reduce the input information and the long time horizons are then represented by a limited number of representative time steps [START_REF] Ochoa | Distribution network capacity assessment: variable DG and Active networks[END_REF]. Other efficient clustering approaches consist in grouping the generations assets by type of resources [START_REF] Zhang | Unified unit commitment formulation and fast multi-service LP model for flexibility evaluation in sustainable power systems[END_REF] or by units displaying similar characteristics [START_REF] Palmintier | Heterogeneous unit clustering for efficient operational flexibility modeling[END_REF]. In this paper the clustering procedure is considered to perform a network partitioning based its topology while computing the "electrical distances". In particular, that zonal decomposition lies in a partitioning method that automatically ensures the cluster connectedness to improve the convergence of the clustering method. Finally, a sequential management procedure with zonal and nodal representations (denoted as S-OM) is implemented in order to ensure the convergence of the management strategy in a faster computational time. It lies on successive runs of the previously OM within each cluster and for the overall system zonal representation.

The rest of the paper is organized as follow. Section 2 presents the management strategy run over a representative day.

The mathematical formulation of the OMLL method is described with a particular attention paid to the power plant operation constraints for the centralized generation. The models for the different types of DG assets are detailed as well as the modified DC OPF constraints that integrates the proposed approximations for the line losses. Then Section 3 introduces the partitioning algorithm that performs the network clustering and the sequential zonal/nodal algorithm is presented. Finally, Section 4 gives a set of different results that allows to validated the proposed approaches (OM and S-OM) as well as the developed models.

Different DG configurations added to the existing centralized assets are investigated. However, it is important to remind that the paper focuses on the management loop only (Figure 1c). The design strategy that would aim at estimating the best configurations will be part of further works.

OPTIMAL MANAGEMENT WITH LOSS LINEARIZATION

A. Jurong Island Power System

As already mentioned the case study is the transmission grid in the Jurong Island of Singapore. Figure 2a displays the corresponding model with the topology of the high voltage lines/buses as well as the position of the four existing power plants.

Note that the 208 buses/219 branches model was built based on assumptions and does not necessarily represent the real topology.

In particular, line impedances are computed from cables manufacturers datasheets and the geographic coordinates of the identified load points. Lines ratings are estimated based on the load at each node with four types of consumers are identified: refineries, oil storage, chemical plants and offices. The voltage connection of each of those consumers depends on their rated power following the rules of connection in Singapore. The centralized generation on site consists in thirteen CCGT units displaying five different technologies (denoted from A to E). Those thirteen units are among the four power plants as follow: Keppel (KPL): 2 units type A and two units type B -Pacific Light Power (PLP): 2 units type C -Sembcorp (SMB): 1 unit type A and 2 units type D -Seraya (SRY): 2 units type C and 2 units type E. Previous works focused on the modeling of the CO2 emissions of the CCGT technologies based on their operating points and efficiency (ranging from 53 % to 62 %) [START_REF] Rigo-Mariani | A generic method to model CO2 emission Performances of combined cycle power plants for environmental unit commitment[END_REF]. Their nominal operating points and part load operations were simulated to finally compute their specific CO2 emissions (SCE in kg/kWh) (Figure 2a). An environmental unit was developed to consider those "carbon costs" in a CO2 minimization while supplying a representative load profile. The optimal management developed in this paper lies on the same strategy with more complexity coming from the addition of different types of DG assets to the existing generation portfolio. Figure 2b shows the a the daily deterministic load that have to be supplied with a 30 minute time step while aggregating the load for all the four types of consumers.

B. Environmental Unit Commitment Objective

As already mentioned the management strategy is implemented in the form of a typical unit commitment (UC) with both centralized power plants and distributed resources. A typical MILP formulation is consider here with the objective to minimize those emissions while fulfilling a set of traditional constraints such as minimum up/down times and considering startup and shut-down carbon costs [START_REF] Van Den Bergh | A Mixed-Integer Linear Formulation of the Unit Commitment Problem[END_REF]. It is important to note that the units contributing to the emissions are the centralized generators in existing power plants (index CG) as well as the smaller distributed units (index DU) that are considered in the DG planning framework. Mathematically the formulation is the same for those two types of generators. The distinction is made as the centralized assets will remain the same (size/site) while the smaller units portfolio will change depending on the different DG configurations that might be investigated. Thus equation (1) allows to compute the total CO2 emissions for both centralized and distributed thermal units located at the system buses (bus set B) and along the simulated day (time set T) discretized at an half hourly time step Δt. The computation considers the on/off (u b,t CG , u b,t DU ) and startup (v b,t CG , v b,t DU ) conditions of the units at bus ensures that each units operates within its acceptable range if it is on. Conventional ramping limits and up/down-times constraints do not appear here for the sake of simplicity. A full list of the symbols used in the optimal management strategy is given in Table 1.

                                               T t B b DU t b DU b C c DU t b c DU b c DU t b DU b CG t b CG b C c CG t b c CG b c CG t b CG b v SU t p a u C0 v SU t p a u C0 obj  
     C c CG t b c CG t b CGmin b CG t b p u P P , , , , and      C c DU t b c DU t b DUmin b DU t b p u P P , , , ,   CGmax b c CG t b c p p , , , 0   and  DUmax b c DU t b c p p , , , 0     CG t b CGmax b CG t b CG t b CGmin b u P P u P , , ,     and DU t b DUmax b DU t b DU t b DUmin b u P P u P , , ,      

C. Centralized Power Plants OperationConstraints

Unlike most UC studies in the literature, the different CCGT units are here aggregated among distinct power plants that display several generators with the same technology (i.e. the same parameters). All the generators k within a same power plant being connected to the same bus, the optimization may lead to unrealistic cases: two units k and k' of the same technology in a same plant opposite deviations (up/down ramps) between two successive time steps while the sum of their output power remains constant. With the considered mathematical formulation, such "distorted profiles" would correspond to the same emissions as a case where units k and k' would generate a constant power. Thus additional variables and constraints are introduced for the centralized generation in order to obtain more realistic operation of the power plants [START_REF] Boyce | Handbook for cogeneration and combined cycle power plants[END_REF]. Binary variables ∆u b,t CG{k}+ and ∆u b,t CG{k}-allow the identification of the up/down ramp direction of each concerned units k at bus b at time t. Those ramp direction are computed with constraint (5) with ∆u max CG being the maximum expected power deviation (set to a high value)

and P b,t CG{k} the output power of the power plant unit k at bus b at time t. That power has minimum and maximum value and is subject to conventional operating constraints related to the piecewise linearization of the carbon costs. Logical constraint [START_REF] Acharya | An analytical approach for DG allocation in primary distribution network[END_REF] ensures that the deviations can be either positive or negative. Finally ( 7)-( 8) allow to give priority (arbitrarily) to unit k' over another unit k'' if k' and k'' display the same technology within the same power plant. Note that only one of the two constraints in ( 8) is necessary, the other being implicitly fulfilled with [START_REF] Acharya | An analytical approach for DG allocation in primary distribution network[END_REF].

             } { , max } { 1 , } { , } { , max k CG t b CG k CG t b k CG t b k CG t b CG u P P P u P   1 } { , } { ,       k CG t b k CG t b u u    } ' ' { , } ' { , k CG t b k CG t b P P    1 or 1 } ' ' { , } ' { , } ' ' { , } ' { ,             k CG t b k CG t b k CG t b k CG t b u u u u  

D. Modified DC OPF for Line Losses Linearization

This subsection describes the introduction of the line losses estimation in the developed optimal management. That estimation is necessary if the impact of different locations of the resources is investigated. AC power flow equations provide an exact computation of those losses when DG planning studies only considers single load level. However, a DC power flow should be considered here to avoid prohibitive computational times while optimizing over a daily time horizon. Conventional loss estimation in DC power flow lies in the computation of the marginal loss factor and generation shift factor matrixes [START_REF] Li | Fully reference-independent LMP decomposition using reference-independent loss factors[END_REF]. The approaches require the integration of an offset that represent an initial guess of the total system loss which is not necessary available. Iterative approaches such as [START_REF] Li | DCOPF-Based LMP Simulation: algorithm, comparison with ACOPF, and sensitivity[END_REF] allow to overcome that problem but might unnecessary increase the computational time with several runs needed. The Jump and Shift iterative method [START_REF] Chen | Jump and shift method for multi-objective optimization[END_REF] will be further considered as a base case to validate the method developed in this paper. Thar method lies in a modified DC power flow that discriminates the flow directions F l,t + and F l,t -in line l at time t before performing a linearization of the corresponding losses δF l,t + and δF l,t -. A preprocessing step is needed in which several AC OPF are run with random load levels and using a dedicated power flow package (Matpower here). A piecewise linearization of W blocks (W = 1 and W = 2) is then performed for each lines using a regular least square method (Figure 3).

  t l t l F F , , , , , ,   t l t l F F   MW MW max , 2 l w f  max , 1 l w f  F l w   , 1  F l w   , 2 
transformer definition (i.e. start/end buses). Finally, line losses are computed following [START_REF] Abbasi | Optimal DG allocation and sizing in presence of storage systems considering network configuration effects in distribution systems[END_REF] with the slope coefficients α w,l δF estimated in the preprocessing phase for each block w of every line l.

, , ,

     W w t l w t l f F and , , ,      W w t l w t l f F   max , , , 0 l w t l w f f    and max , , , 0 l w t l w f f      , , , ,       W w t l w F l w t l f F    and , , , ,       W w t l w F l w t l f F     
That piecewise linearization is integrated in a conventional DCPF formulation [START_REF] Seifi | Electric Power System Planning -Issues, Algorithms and Solutions[END_REF] where line flows are computed with the bus admittance matrix B bus and the bus instantaneous voltage angles θb,t [START_REF] Popovic | Placement of distributed generators and reclosers for distribution network security and reliability[END_REF]. Finally, the power balance at each bus is computed according to [START_REF] Gozel | Optimal Placement and Sizing of Distributed Generation on Radial Feeder with Different Static Load Models[END_REF] as the contribution of the in the branches arriving at the bus or departing from it. Noted that similarly to the method developed in [START_REF] Li | Fully reference-independent LMP decomposition using reference-independent loss factors[END_REF] the losses within a line are equally distributed at the start and end buses as additional loads.

                 1) bus (slack 0 to from defined line with , 1 , ' ' , ' , , b T t b'' b' l F F t t b t b t l t l    bus ' b' b' B                                                        b l t l t l t l t l b l t l t l t l t l F F F F F F F F P b,t to lines , , , , from lines , , , , 2 2 2 2
     

E. Integration of DG Resources 1) Conventional DG assets

This section introduces the modeling for the different DG assets in the developed optimal management with linearized losses (OMLL). Firstly, PV distributed generation at bus b and time t P b,t PV is considered with different sizes of installations simulated using a normalized profile for a 1MWc plant in Singapore that was generated with the PVwatts tool from the National Renewable Energy Laboratory [32]. The management strategy considers the possibility to curtailed the PV generation with the variable 𝛿P b,t PV-( 14) in order to avoid potential line congestion is case of oversized installation [START_REF] Alnaser | Optimal sizing and control of energy storage in wind power rich distribution networks[END_REF]. The second type of considered DG resources corresponds to distributed generic energy storage (ST) devices. A typical linear formulation is used while discriminating the charge (P b,t ST-) and discharge (P b,t ST+ ) powers for the units at bus b at time t. A power efficiency η ST = 97.5 % is introduced to compute the state of charge SOCb,t following (15) with the rated energy E b STmax and the initial charge level SOCb,t=0 (typically 50 %) [START_REF] Pandžić | Near-optimal method for siting and sizing of distributed storage in a transmission network[END_REF]. The set of constraints ensures to fulfill the power/energy limitations and enforces the storage devices to come back at their initial SOC at the end of the simulated day. The last type of conventional DG resources that is considered correspond to small CCGT. The operating parameters of those distributed units (DU) derive from the technology A with value scale down for smaller rated power. The operating constraints of those assets are similar to the ones implemented for the centralized thermal units and have already been discussed in section 2.A.

PV t b PV t b P P , , 0                           t i ST i b ST ST ST i b STmax b t b t b P P E t SOC SOC 1 , , 0 , , 100     0 , , , , , , 100 0 , , 0          t b T t b t b STmax b ST t b ST t b SOC SOC SOC P P P  

2) Controlable loads

This paper considers load controllability (CL) as a potential DG resource. Many strategies are investigated in the literature when dealing with load control ranging from highly centralized demand response methods to the fully distributed transactive energy framework [ Kok 16]. In this study a centralized control is implemented as the assumption is made that a vertical utility (Jurong Island Community here) owns and operates all the assets on site. In particular, it implies that the customers on site agrees to change their consumption pattern accordingly to central controls which might lead to some privacy issues and be challenging when heavy industrial customers are concerned. Incentives policies should then be considered but it is not in the scope of this paper. Distributed energy storage already provides an example of full continuous flexibility with controllable charge/discharge powers. Here the load control strategy relies on discrete load increments or decrements that would correspond to the time shifting of chemical processing by the considered consumers. The idea is to provide a different type of flexibility.

Also, in order to better manage the complexity, loads are aggregated within N distinct clusters using a network partitioning strategy that will be discussed in the next section. The actual load control will then be applied to clusters and not individual customers. Each of those cluster n agrees in advance upon a global load level ∆P n ΣCL and duration T n ΣCL for the control phase. If a power decrease (u n,t ΣCL-) or increase (u n,t ΣCL+ ) for the cluster n is triggered at time t (i.e. variables set to 1 by the management), all the loads in cluster n (i.e. b ∈ n) will adapt their profiles according to a predefined individual level ∆p b CL as presented in Figure 4.

Equation [START_REF] Pandžić | Near-optimal method for siting and sizing of distributed storage in a transmission network[END_REF] ensures that power uprate and derate within each cluster do not occur simultaneously while (18) allows controlling the shift phase duration with only one load shifting occurring each day ((

)). 1 , ,       CL t n CL t n u u                   1 , 1 , , 1 CL n T t t i CL i n CL n CL t n CL t n u T u u and                 1 , 1 , , 1 CL n T t t i CL i n CL n CL t n CL t n u T u u   T u u CL n T t CL t n T t CL t n            , , 19

 

3) Total bus power balance

The power balance at each bus is computed according to [START_REF] Ugrandi | Optimal wind turbine sizing to minimize energy loss[END_REF] as the difference between both centralized and distributed generation at bus ad its instantaneous load (P b,t L ). Note that incidence matrices are needed to map the DG assets in the grid and test different configuration. Also the equality between constraints ( 13) and ( 20) allow to compute a Lagrangian that directly aggregates energy, congestion and loss costs (carbon costs here) where other approximation of the losses would require additional post processing. The MILP formulation of that optimal management (denoted OM) is finally written using YALMIP toolkit [START_REF] Lofberg | YALMIP: a toolbox for modelling and optimization in MATLAB[END_REF] and solved with CPLEX 12.7.

                   

GRID CLUSTERING AND REDUCED COMPUTATIONAL TIME WITH ZONAL APPROACH

A. Computational Time Reduction

Preliminary runs of the OM converged for different configurations of the DG assets in the considered system. However, in some cases the computational time exceeded the arbitrary limit set to 5 minutes. If a couple of minutes/hours computation would not be an issue for a day ahead scheduling problem they can become prohibitive if the OM has to be run successively a great number of times as it is the case in the DG planning iterative strategy chosen here (i.e. hundreds/thousands runs). A grid clustering approach is then considered in this section in order to furtherly reduce the computational time of the implemented OM. As already mentioned in the introduction, such clustering methods are often used to reduce the amount of input information while reducing a whole year to a sample of representative time step. Other approaches consist in aggregating the generation unit of same type or displaying close characteristics. In this paper, the clustering lies in a network partitioning (bus clustering) and the considered power system is disaggregated in smaller microgrids. The previously developed OM can then be run faster in

B. Electrical Network Partitioning

The clustering method applied to the Jurong Island power system consists in a regular K-means algorithm [START_REF] Jain | Data clustering: 50 years beyond K-means[END_REF] hybridized with a standard genetic algorithm similarly to [START_REF] Blumsack | Defining power network zones from measures of electrical distances[END_REF]. The first step of the network partitioning is to define the way to compute the "electrical" distance between the buses. A typical approach consists in considering the line impedance with the bus impedance matrix Z bus . A modification is introduced in order to encourage the buses to be aggregated accordingly to their voltage level if the pre-defined number of clusters is high enough. A penalty coefficient λ (set to 10 3 ) is considered and the distance db'b'' between two buses b' and b'' is finally computed following [START_REF] Celli | Optimal integration of energy storage in distribution networks[END_REF] if there is an existing branch between the buses. With highly sparse matrices for radial systems the cluster connectedness might not ensure when running a typical K-means [START_REF] Sanchez | Multi-attribute partitioning of power networks based on electrical distance[END_REF]. Thus the distance between two buses b' and b'' with no direct link (i.e. connection line/transformer) needs to be modified. The shortest path in terms of impedance is then considered using Dijkstra's algorithm. That formulation leads to a distance matrix with no null coefficients. The cluster connectedness is automatically guaranteed as the distance computation with the Dijkstra's algorithm lies on the summation of serial impedances.

          b'' b' b b b'' b' b b V V d V V d if if ' ' ' ' ' ' bus ' b' b' bus ' b' b' Z Z   
The idea of the hybridized K-means is to use the genetic algorithm to generate different sets of clusters centers (i.e. n buses defined as centers clusters Kn). Regular K-mean clustering is then performed for each set of centers before estimating the clustering performances which is computed as the aggregation of two indexes (( 22)). On one hand the clustering tightness index (CTI in [0,1]) measures the extent to which each bus is proximate to the other members in the same cluster Then the cluster size

b B b K b b b B b n n                          , , 1 2 )
( # ln ) ( # ln ' ' ' ' ' ' ' ' ' ' ' '

 

C. Modified Optimal Management

Figure 6a shows the results for a grid partitioning with ten clusters n each of them containing buses b and lines l that are subsets of the original B and Ldenoted b ∈ n and l ∈ n. Once the network partitioning is done, the overall system can be modeled using a zonal representation where each cluster is identified with an equivalent bus b n ∈ B n and the different buses are connected thgouth lines l n ∈ L n that are a subset of L (Figure 6b). The zonal representation obviously corresponds to a much simpler system with a 208 buses / 210 branch system being reduced to a 10 buses / 10 branch model. Then that subsection describes a 3 steps method in order to take advantages of the zonal representation to lower the computation time of the OM when different configurations of DG assets are investigated.

1) Step 1 : Optimal management over the zonal modelConventional DG assets

The first step consists in running the OM method previoulsy describe for the zonal representation of the model. The mathematical formulation does not need to be modified, excpet that the original buses and branches sets (B and L) have to be subsittued with Bn and Ln in the equations. However, a modification of some of the input paramters is required in order to aggregate the load, the PV generation as well as the rated power/capacity for the storage systems and ditributed CCGT units following [START_REF] Chen | Jump and shift method for multi-objective optimization[END_REF]. Thus each of those resources will be represented by a single equivalent asset at every bus bn. Nota that the centralized units on site are not aggregated and, as previously, control variables will be considered for each of them. Also the controllable loads were defined to be aggregated in clusters. Thus the power magnitude of the control phase at each equivalent bus bn corresponds to the cluster predefined control level ∆P n ΣCL . Finally, the power balance constraint at each bus (20) can be rewritten as in [START_REF] Li | DCOPF-Based LMP Simulation: algorithm, comparison with ACOPF, and sensitivity[END_REF] for better understanding.

                               
    DU t b CL t n CL t n CL n ST t b ST t b PV t b PV t b L t b n b CG t b
                     
The OM run for the defined zonal model is much faster that for the nodal representation of the whole system with many less variables. That number of variables in the OM is never fixed and depensd on the considered DG configuraton the number of distributed assets for each ressrouces. The ouputs of that first step are the control variables for the centralized and aggregated DG assets as well as the power flows in lines Ln. Those flows time profiles denoted F l n ,t +* and F l n ,t - * define the energy exhange (i.e. import/export) between the cluster (i.e. theequivalent buses bn). However the obtained results for the zonal model would diverge from the ones computed for the complete system as losses and congestion issues within each cluster are not considered.

2) Step 2 : Optimal management within each cluster

A second step in implemented in order to refine the results inside every cluster with the OM run successively for each of them. Once again the mathematical formulation of the optimization problem does not require major modifications. The same equations are considered with subsets of B and L conrresponding to every cluster (i.e. considering b ∈ n and l ∈ n). The DG assets are not aggregated and their nodal position and individual control variables are considered when computing the power balance constraint at each bus. However, the bus injection constrainst [START_REF] Gozel | Optimal Placement and Sizing of Distributed Generation on Radial Feeder with Different Static Load Models[END_REF] with the line import/export is modified when considering the buses b ∈ n that have a direct connection line with another cluster (Ln). The equation is modified following [START_REF] Zhang | Unified unit commitment formulation and fast multi-service LP model for flexibility evaluation in sustainable power systems[END_REF] in order to ensure that the cluster n export the amount of energy predicted by the previous steps with the power flows profiles F l n ,t +* (for ln leaving cluster n) and F l n ,t - * (for ln arriving at cluster n). On the contrary, imported flows F l n ,t + (for ln arriving cluster n) and F l n ,t -(for ln leaving cluster n) remain free variables subject to constraint [START_REF] Palmintier | Heterogeneous unit clustering for efficient operational flexibility modeling[END_REF] in order to ensure the convergence of the OM while compensating clusters losses or any deficit of local generation. A term is also added to the objtective function of the OML with a penalty coefficient λ (typically λ = 10 6 ) to keep the import of energy from the other cluster as low as possible and encourage the local generation (for both centralized and decentralized assets) (( 27)).

                                                                        n b l n b l t l t l n b l t l t l t l t l t l t l n b l t l t l t l t l n n n n n n n,t b F F F F F P F P F F F F P to n to * , , from , * , , , , , from n , , , , 2 2 2 2
    (25)            n F F n F F t t
                         t n b l t l n b l t l n n P P obj obj to , from ,   
Once the second step is performed, the net load for each cluster n (P b n ,t Net ) can be computed considering the load and optimal generation profiles at each buses (b∈ n) as well as the losses computed over the cluster lines ((28)).

                                 n l W w t l w t l w F l w DU t b CL t n CL t n CL b ST t b ST t b PV t b PV t b n b CG t b n b L t b Net t b f f P u u p P P P P P P P n n n n n , , , , , , , , , , , , , , ,   


3) Step 3 : Optimal management within each cluster

A final OM is run for the zonal representation while considering the net load profiles for the clusters (buses bn). the power balance constraint at each bus ( 20) is modified folllwong [START_REF] Boyce | Handbook for cogeneration and combined cycle power plants[END_REF]. No DG assets are considered in this step as their operation is already imlied in the computation of the cluster net load. That last runs aims at balancing the system losses and the cluster import/export with the centralized CCGT units. The sequential optimal managemet strategy is denoted as S-OM and is summarize in Figure 7.

B. Line Losses

Another set of simulations aims at validating the estimation of the line losses compared to an iterative reference method (Jump and Shift method here [START_REF] Chen | Jump and shift method for multi-objective optimization[END_REF]). Serval runs of the OM wit 1 block (W =1) and 2 block (W = 2) are performed for the Jurong Island systems as well as for the IEEE 30-bus and 57-bus test cases. As previously, those test run do not consider any DG connected to the system. Also the possibility to "correct" the ouputs of the OM is considered with a sequential AC OFP. Like for the Jump & Shift method, those OPF are run for every time step while updating the bounds for the units generation around the power profiles returned by the OM (a tolerance of 5 % is considered here). Table 2 displays the obtained results computational times are reduced using the OM with linearized losses. The error in the objective function (i.e. CO2 emissions) is very low, under 0.5 %, while the active losses difference is around 5 %. Running a sequential AC OPF after the OM (with W = 2) allows adjusting the results and the obtained solution is even slightly better than the one returned by the Jump and Shift method. The simulations also allow the analysis of the power flowing through the lines along the day for the different systems. Results

for the IEEE 30-bus test case in Figure 9 show consistent power flows in the branches computed with the OM (W = 2) before and after the application of a sequential AC OPF (Figure 9a) (MW for UCLL2 and MWA after AC OPF). Greater differences are observed when considering active losses with deviations either positive or negative (Figure 9b). However, it is noticeable that the rank of the least efficient branches is fully respected which is important when allocation problems are investigated.

Indeed, the UCLL will prove able to properly discriminate two configurations (in terms of losses) where the same DG unit is connected to two different buses.

C. Computatinal Times

The next validation tests investigate the computational times for different DG configuration while comparing the performance of the firstly developed ON and the sequential version (S-OM). At first, twenty random distributions of PV installations and distributed CCGT units (DU) are tested with different capacities and allocations. The results displayed in Figure 10a show good consistency while computing the emissions with the two optimization procedures. The recorded computational times are similar.

The S-OM is slightly faster while the regular OML encounters one case where the arbitrary 5 minutes limit of time is reached.

The OM also exceeds the time limit for one case when distributed storage is added to the investigated DG configurations. With the introduction of storage assets, the S-OM tends to be slower than the OM (when the convergence is reached) while the results in terms of emissions are still consistent (Figure 10b). When controllable loads are considered with PV and DU the emissions computed with the two methods display greater deviations. Especially the rank of the "cleanest" DG distribution is not fully respected any more (Figure 10c). However, the regular OM reaches the maximum computational time in most of the cases and the associated results might be sub optimal. The S-OM overcomes that difficulty, allowing a faster convergence with lower CO2 emissions in most cases.

D. DG configurations 1) Base Case Scenario

The last set of simulation focuses on the performances of different DG configuration regarding the objective of CO2 reduction compared to a base case scenario in which only the centralized generation on site is considered while running the OM (with W = 2 here). As already displayed in Table 2, the simulated day corresponds to 26116 tons of CO2 emitted, with an equivalent grid emission factor of 413 kg/MWh. That value is close to the official 431 kg/MWh estimated by the Singaporean energy authority for the year 2015 [2]. The difference can be explained by the 5 % of energy generated from old, dirtier oil power plants in Singapore, which is included in the official statistics. Regarding the plant generation the bulk load is provide by the cleanest units (type A in Keppel and type C in Seraya) while the dirtier generator are started up to supply the peak demand (Figure 11b). The system displays very low losses equals to 0.35 % when computed over the entire day (Figure 11a). Note that only higher voltage lines are considered here. The total grid losses (transmission/distribution) in Singapore is estimated around 2 % according to the International Energy Agency. Those low values are explained by the low impedance underground cables on site as well as by the short distances between the nodes.

2) Conection of DG assets on a test feeder

A firs set of simulations is performed considering the 20 MVA rated feeder connecting the buses 165-168 (Figure 2a).

Once again, the tests do not imply any optimal DG size/site. Here the idea is to validate the problem formulation for an arbitrary chosen feeder. Different DG resources are successively connected to the nodes and results are estimated regarding the overall CO2 emissions and losses reduction within the feeder. For the base case (S0) the feeder is loaded according to the buses consumption profiles along the day for a total of 4.3 MWh losses. Then a scenario (S1) is simulated with the connection of a 30 MWc PV plant to node 168. That peak power is at an arbitrary high value (greater that the feeder capacity which is not realistic) in order to validate the generation curtailment in case of important surplus of production. In practice, most of the generation is self-consumed by the local loads and the surplus flows upstream to the main grid (Figure 12a). Losses within the feeder are significantly reduced (more than 25 %) as a result of the shortened distance between the generation and the load point (Table 3).

Adding a 1 MW/5 MWh storage to the same bus (S2) does not improve the performance much with a net load slightly reduced during short periods (Figure 12a). A cluster of controllable loads is then considered (S3) with 1MW shifted at nodes 165, 166 and 167 for 3 hours (Figure 12b). As expected the load is increased when the PV peak occurs. It leads to less energy exported to the grid, which corresponds to reduced losses in the feeder. Once again the benefits are not significant and many more loads should be controlled to get a real impact on CO2 reduction. A last simulation is performed with the addition of a 5 MW distributed CCGT unit at bus 165 (S4). The optimal profiles returned by the OM method logically correspond to full-time operation of the distributed clean CCGT unit (Figure 12b) to lower the overall emissions.

3) Impact of DG size and site

Additional simulations are run with increasing DG capacities. For each type of DG resources, the installed capacity is distributed among five arbitrary buses (Figure 2a). Storage units are given with a power/energy ratio equal to one and a single cluster aggregates the load controlled at the test buses. As previously observed, photovoltaic and distributed CCGT units have the most significant impact on the carbon cost, with reduced system losses (Figure 13a). Both emissions and losses tend to decrease linearly as more capacity is installed. There is no major improvement when only controllable loads or energy storage are considered, even with greater capacities involved. Indeed, those two resources only allow shifting consumption through time and they cannot directly contribute to carbon reduction by themselves when a given amount of energy has to be consumed along the day. However, a slight improvement is observed with capacities above 15 MW. In those cases, the resources provide enough flexibility to allow a small benefit from the arbitrage between dirtier and cleaner CCGT centralized units. Typically, generation from the cleanest units is increased at night time and the reduced consumption during the day allows minimizing the time period where dirtier technologies are committed. The gain in terms of tons of CO2 is then directly linked to the carbon cost difference between the different type of units. With only one type of fuel considered here (i.e. natural gas) the costs are very close to each other, which explains the almost non observable improvement. Note that the efficiency attached to storage units lead to increased losses even if the overall carbon cost is lower than in the base case scenario (Figure 13b). The DG resource siting is then investigated by distributing a same amount of DG capacity (30 MW) among different nodes. Twenty random sets of buses considered -sets of one, three, five and ten buses. As previously observed the carbon emissions decrease proportionally to losses when PV units are investigated (Figure 14a). However, the PV siting does not have any impact on the system performance, with a variation less than a ton (i.e. < 0.01 %) over the twenty simulated scenarios. It becomes a little more significant when allocation of decentralized CCGT units is studied. In particular, distributing the capacity among several buses seems to correspond to slightly improved performance (around 1 %) (Figure 14b). system in order to reach the total distributed capacities (denoted STtot and PVtot). Note that the considered capacities are unrealistic regarding the available space on Jurong Island and that those simulations should be considered as a partial sensitivity analysis. As already noticed in the previous simulations the improvements in terms of CO2 reduction are not significant. Indeed the results displayed in Table 4 show around 1 % reduction for every 100 MWp of PV installedwhich is in the range of the total installed capacity currently installed in the whole Singapore. The impact of the storage is even more limited less than 0.5 % reduction for 500 MWh! The results are the same when storage units are connected at the same buses as the PV plant and with the same capacities (PVtot = STtot in Table 4). Those low improvements are explained by the nature of the considered Jurong Island power system. Especially, there is no energy mix regarding the centralized generation with all the CCGT units fueled with natural gas. As displayed on Figure 2b the performances of the units in terms of CO2 emissions are close. Thus if DG allow reducing the generation of the "dirtiest" unit, improvement will not be significant. Also the rated power of the centralized unit (around 400MW each) is high regarding the load level (3GW). It means that important levels of DG should be installed in order to cover the contribution of a single unit and avoid starting it which might be unrealistic as observed in Table 4. The uselessness of storage device can be more specifically explained with the optimal scheduling problem that does not consider any time varying parameters (e.g. time of use or real time prices). The objective function is only the amount of CO2 for the power generation and the interest of using storages to perform arbitrage thought time is limited here. The main contribution of storage unit would be to shave the overall peak demand in order to reduce the generation of the "dirtier" unit. However, once again the absence of energy mix with similar generator limit the expected improvements in terms of CO2 savings.

CONCLUSIONS

The method developed in this paper presents an optimal energy management of a power system in the presence of different types of DG resources and with the centralized generation already installed on site. The case study refers to the Jurong Island in Singapore. A particular attention is attached to the computational time of the procedure as that management loop aims at being integrated in a DG planning method that would require to test many configurations of the DG assets. At first, the use of a modified DC OPF implemented in the form of MILP problem avoids computational complexity. Especially, the method lies in a piecewise linearization of the lines losses that has been validated over several systems (Jurong Island and IEEE test cases) and compared to another method with AC OPF constraints taken as the reference. In a second time a clusterised version of the optimal management is developed once a generic method is developed for network partitioning. Validation us showed that the clusterized method allow to reduce the computational time and the convergence before an arbitrary limit of five minutes. The implementation of the DG optimal planning is not in the scope of the paper. However, some simulations have been performed in order to estimate the improvement in terms of CO2 reduction depending on the DG types, sizes and sites. If the obtained results validated the implemented models and constraints, they did not how significant improvements compared to a bases case scenario in which only the centralized generation on site is considered. That is directly link to the nature of the Jurong Island power system with very low losses and only gas fueled units (i.e. no generation mix). Considering objective function with other parameters such as operating costs or carbon taxes might create more space for improvements. That should be part of further works should as well as the integration of the management strategy in a design procedure in order to find the best DG configuration.

Figure 1 :

 1 Figure 1: Coupling between design and management loopsa) integratedb) all in one -c) iterative

Figure 2 :

 2 Figure 2: Jurong Island Power Systema) Grid topologyb) CO2 emissions of centralized CCGT technologiesc) aggregated load profiles

b

  at time t as well as the corresponding base "carbon costs" (C0 b CG , C0 b DU ) and startup emissions (SU b CG , SU b DU ). The piecewise linearization is run along C, the set of linear blocks, in which the generated power (p c,b,t CG , p c,b,t DU ) is attach to a cost (a c,b CG , a c,b DU ).

  of CCGT unit at bus b tons SU b CG ,SU b DU startup CO2 emissions of CCGT unit at bus b tons a c,b CG ,a c,b DU ooperating CO2 emissions of units bus b operating in block c tons/MWh p c,b CGmax , p c,b DUmax upper bound of CCGT unit power at bus b operating in block c MW P b CGmin , P b CGmax min/max output power of centralized CCGT unit at bus b MW P b DUmin ,P b DUmax min/max output power of distributed CCGT unit at bus b MW ∆P max CG maximum ramp value for the centralized CCGT units MW f w,l max maximum power flow in block w of line l at time t MW F l max maximum flow in line l MW α w,l δF loss factor in block w of line l -B bus B×B bus admittance matrix p.u. P b,t PV available PV generation at bus b at time t MW P b ST max , E b ST max rated power/energy for storage at bus b MW/MWh η ST generic storage power efficiency (typically 0.975) -∆P n ΣCL , T n ΣCL power magnitude and duration of control phase of cluster n MW/h ∆p b∈n CL power magnitude for CL at bus b in cluster n MW P b,t L load at bus b at time t MW Carbon costs linearization requires sets of typical variable and constraints [28] in order to compute the output power of the CCGT units (P b,t CG , P b,t DU ) based on their on/off status, minimum output power (P b CGmin , P b DUmin ) and their operation within each linear blocks c (2). The contribution in each block is subject to upper bounds (p c,b CGmax , p c,b DUmax ) ((3)) and constraint (4) finally

Figure 3 :

 3 Figure 3: Line Loss LinearizationSimilarly to the linearization of the CO2 emissions new sets of variables and constraints allow the loss estimation with f w,l,t + ,f w,l,t -the positive/negative instantaneous line flow in block w and their corresponding upper bounds f w,l max . Similarly to the piecewise linearization of the carbon cost function, W intervals are defined to compute the power flowing thought the branch ((9)). Equation[START_REF] Vigano | Methodology for the optimal siting and sizing of storage systems in distribution networks[END_REF] implicitly allows the avoidance of line congestion ( 0 ≤ F l,t + , F l,t -≤ F l max ). Note that for transformers, only unidirectional flows are possible meaning that either the positive or the negative flow has to remain null depending on the

Figure 4 :

 4 Figure 4:Control of loads at buses b' and b'' in cluster n

Controllable load 1

 1 of cluster k at bus b' Controllable load 2 of cluster k at bus b'' microgrids with a reduced numbers of variables while additional operations should ensure to maintain boundary conditions at the points of connection between the clusters.

Figure 5

 5 Figure 5 the use of the genetic algorithm allow obtaining good performances with all the indexes close to unity for ten clusters or more.

Figure 5 :

 5 Figure 5: Clustering resultsa) one run of K-meansb) hybridized K-meanc) partitioning for 10 clusters

Figure 6 :

 6 Figure 6: Network partitioning for N = 10 clustersa) nodal representationb) zonal representation

  with centralized generation and export to mainland 66 kV hub transformer (mono directional flow)set L n , subset of L line (bi directional flow)set L n , subset of L clustersequivalent buses, set B n a) b) subsets of buses B and lines L in each cluster n

Figure 8 :

 8 Figure 8: Sembcorp plant operationa) inactive constraintsb) active constraints

Figure 9 :

 9 Figure 9: IEEE 30-bus results-a) line powerb) line losses

Figure 10 :

 10 Figure 10: Comparison between OM and S-OMa) PV and DUb) PV, DU and BATc) PV, DU and CL.

Figure 11 :

 11 Figure 11: Power profiles for the base case scenario-a) load and generationb) plants operation

Figure 12 :

 12 Figure 12: Feeder loada) S0, S1 and S2b) S3 and S4

Figure 13 :Figure 14 :

 1314 Figure 13: Impact of DG size at 5 arbitrary test buses-a) CO2 emissions for different capacitiesb) losses for different capacities

Table 1 :

 1 List of the used symbols in the OMLL

	Sets :
	b ∈ B	set of buses	-
	t ∈ T		set of time steps	-
	c ∈ C	set of linear blocks for the CCGT carbon emissions	-
	l ∈ L		set of lines	-
	w ∈ W	set of linear blocks for line losses estimation	-
	n ∈ N	set of grid cluster (with buses and lines)	-
	b n ∈ B n l n ∈ L n	set of equivalent buses in zonal model set of lines in in zonal model (subset of L)	--
	Variables :
	u b,t CG ,u b,t DU	on/off status of CCGT unit at bus b at time t {0,1}	-
	v b,t CG ,v b,t DU	startup of CCGT unit at bus b at time t {0,1}	-
	p c,b,t CG ,p c,b,t DU	CCGT unit power at bus b operating in block c at time t	MW
	P b,t CG , P b,t DU	power output of CCGT unit at bus b at time t
	P b,t CG{k}	power output of a unit k in power plant at bus b at time t	MW
	∆u b,t CG{k}+ , ∆u b,t CG{k}-F l,t -+ , F l,t δF l,t + , δF l,t -	up/down ramp of unit k in power plant at bus b at time t {0,1} positive/negative power flows in line l at time t positive/negative power flows in line l at time t	-MW MW
	f w,l,t +	, f w,l,t
	P b,l Σ		power balance of at bus b at time t	MW
	δP b,t PV-

positive/negative power flows in block w of line l at time t MW θb,t voltage angle at bus b at time t p.u. curtailed PV generation at bus b at time t MW P b,t ST-,P b,t ST+ charge/discharge powers of storage at bus b at time t MW SOC b,t state of charge of distributed storage at bus b at time t % u n,t ΣCL-, u n,t ΣCL+ control for decrease/increase of cluster n at time t {0,1} -Parameters :

Table 2 :

 2 Obtained Results for Jurong Island System

		Jump Shift	OM (W = 1)	OM (W = 2)	OM (W = 2) + AC OPF
	Jurong Island Power System -13 units
	Tons of CO2	26119	26116	26114	26118
	Losses (MWh)	219	211	206	219
	CPU Time (sec)	63	6	8	13
	IEEE 30 buses -6 units
	Tons of CO2	1701	1681	1698	1699
	Losses (MWh)	59	36	52	55
	CPU Time (sec)	38	3	7	10
	IEEE 57 buses -7 units
	Tons of CO2	11136	11036	11113	11132
	Losses (MWh)	309	275	286	300
	CPU Time (sec)	51	9	23	28

Table 3 :

 3 DG connected to the tested branch

		S0	S1	S2	S3	S4
	CO2 (tons)	26114	26056	26055	26054	26029
	Losses (MWh)	4.3	3.0	2.9	2.7	1.8

Table 4 :

 4 CO2 emisiosn for fdiffenrt PV and ST configurations

	PVtot	STtot	0 MWh	100

MWh 200 MWh 300 MWh 400 MWh 500 MWh

VALIDATION RUNS ANS SAMPLE RESULTS

A. Power Plant operations

In this section 4, validation tests for the implemented methods (OM and S-OM) are presented along with preliminary results for different DG configurations. Note the DG planning procedure is not in the scope of that paper. Thus the DG configurations that will be considered are not optimal and are only use to estimate the sensitivity of the objective function for different scenarios.

The first set of simulations is run with no DG assets connected to the system. The objective is to verify the validity of the constraints in Section 2.C. that ensure appropriate operations of centralized power plants displaying several units of the same type. Figure 8 shows the results for the generation of the Sembcorp power plant with its two units of type D when the constraints are inactive (Figure 8a) and activate (Figure 8b). As expected, the introduction of the power plant constraints allows to smooth the generating profiles while the value objective function remains unchanged (i.e. tons of CO2). Note that the overall power plant generation is consistent with for the greatest part of the day. The differences observed for the last two time steps are compensated by the same type of units in another power plant. Also note that the avoidance of multiple ripples in the units' power profile could have been done with the introduction of ramping cost that are not considered here.

Run OM for zonal model with aggregated DG capacities in each cluster Sets : B n , L n Objective : (1)minimize CO 2 Constraints : (2)-( 19), ( 23)-( 24)

Run OM for each cluster with disaggregated DG and cluster import constraints Sets : b∈n, l∈n (buses and lines in cluster n) Objective : (27) -minimize CO 2 /energy import Constraints : (2)-( 19), ( 25)-( 26)

Run OM for zonal model with clusters net load, no DG assets Sets : B n , L n Objective : (1)minimize CO 2 Constraints : (2)-( 19), ( 29)


Noureddinecherrad

Djamelselloum

Sophie Tingry

Modeling the control of the desorption rate of hydrogen released from the adsorption storage bed to supply a fuel cell

Keywords: Hydrogen, DHR modeling, Excess amount, Desorption, Fuel cell, Power Nomenclature

The storage of hydrogen by adsorbent materials has been addressed by several researchers. These materials can serve as a reservoir of hydrogen at a very lower temperature and the emptying operation for fuel cells is simply heating the adsorbent bed. However, in order to maximize the ability of the adsorbent materials to meet the instantaneous hydrogen demand, adequate knowledge of desorbed hydrogen flow rate (DHR) must be investigated. The objective of this study is to model the control of the DHR induced by heating. The results show that the excess of hydrogen stored in the adsorbent material can be completely released at room temperature and the DHR increases with temperature. A solution is proposed for stabilizing the DHR, which consists in controlling the fuel cell supply section and, consequently, the power to be produced.

Introduction

Most of the energy used is produced from fossil fuels, which requires a lot of time to be renewable.

Although these energies are less expensive, they must meet the challenge of replacing them with renewable and environmentally friendly energies [START_REF] Abdalla | Hydrogen production, storage, transportation and key challenges with applications: A review[END_REF]. Power to hydrogen is a promising solution for storing variable Renewable Energy (RE) to achieve a 100% renewable and sustainable hydrogen economy [START_REF] Dawood | Hydrogen production for energy: An overview[END_REF]. For this reason, hydrogen appears as a potential alternative to fossil fuels and it became one of the most abundant elements on the planet [START_REF] Sinigaglia | Production, storage, fuel stations of hydrogen and its utilization in automotive applications-a review[END_REF]. In the ever growing demand of ___________________________________________________________________________ Page 3 of 19 future energy resources, production reaction of this fuel has attracted much attention among the scientific community [START_REF] Sharma | Hydrogen adsorption on pristine and platinum decorated graphene quantum dot: A first principle study[END_REF] and it can be obtained from many sources such as water, biomass, natural gas, ethanol, etc. [START_REF] Sinigaglia | Production, storage, fuel stations of hydrogen and its utilization in automotive applications-a review[END_REF].

Although hydrogen is considered to be one of the most promising green fuels, its efficient and safe storage and use still raise several technological challenges [START_REF] Czakkel | Static and dynamic studies of hydrogen adsorption on nanoporous carbon gels[END_REF]. The reason why the production and the storage of hydrogen have recently been the subject of particular attention [START_REF] Czylkowski | Hydrogen-enriched gas production from kerosene using an atmospheric pressure microwave plasma system[END_REF], because the exploitation of this gas as a clean fuel has many advantages [START_REF] Md Arshad | Preparation of activated carbon from empty fruit bunch for hydrogen storage[END_REF]. Its impact on the environment is greatly reduced and comes from a clean and sustainable renewable energy source. Hydrogen will gradually replace non-renewable energies for clean energy purposes [START_REF] Oriňáková | Recent applications of carbon nanotubes in hydrogen production and storage[END_REF], as in fuel cells technologies for vehicles matrices with the expansion of H2 filling infrastructures in Europe [START_REF]Germany adds more hydrogen stations in Wolfsburg, Ingolstadt[END_REF].

However, the treatment of hydrogen as a fuel must have the same degree of risk as other conventional liquid and gaseous fuels and, since this gas operates at high pressure, its existence in a confined space presents a risk of explosion [START_REF] Shirvill | Effects of congestion and confining walls on turbulent deflagrations in a hydrogen storage facility-part 1: Experimental study[END_REF]. Therefore, the storage of hydrogen for commercial purposes with high gravimetric density is a major task [START_REF] Sathe | First-principles study of hydrogen storage in metal functionalized [4,4]paracyclophane[END_REF][START_REF] Barthelemy | Hydrogen storage: Recent improvements and industrial perspectives[END_REF]. Safe storage of large capacities is an essential step in the exploitation of hydrogen as renewable energy [START_REF] Edalati | Design and synthesis of a magnesium alloy for room temperature hydrogen storage[END_REF].

High strength steel coatings have been used on high pressure gas cylinders containing hydrogen or natural gas. Nevertheless, hydrogen penetration is possible in this coating layer, making the material brittle to withstand the damage caused by high pressures [START_REF] Yoon | Hydrogen damage in 34CrMo4 pressure vessel steel with high tensile strength[END_REF]. To be able to store energy on board, the success of fuel cells requires a compact, lightweight and affordable hydrogen storage system to replace pressurized hydrogen tanks under pressure [START_REF] Sadhasivam | Dimensional effects of nanostructured Mg/MgH2 for hydrogen storage applications: A review[END_REF]. Therefore, the challenge of its storage extremely confined its applications. Current researchers are mostly focused on the development of porous nanomaterials having large specific surface areas for improving the gravimetric densities, which undoubtedly restrain the volumetric densities of hydrogen storage [START_REF] Nabi | Gallium vacancies role in hydrogen storage of single-crystalline GaN hexagonal micro-sheets[END_REF]. The discovery of new hydrogen storage materials has advanced hydrogen storage technology in recent decades [START_REF] Wang | Recent advances in additive-enhanced magnesium hydride for hydrogen storage[END_REF]. In this respect, various developed composite and nanostructured materials are currently considered promising candidates for the storage of hydrogen to supply fuel cells [START_REF] Sadhasivam | Dimensional effects of nanostructured Mg/MgH2 for hydrogen storage applications: A review[END_REF][START_REF] Garroni | Recent progress and new perspectives on metal amide and imide systems for solid-state hydrogen storage[END_REF][START_REF] Li | Nanostructured Materials for Next-Generation Energy Storage and Conversion: Fuel Cells[END_REF][START_REF] Broom | Hydrogen storage materials: the characterisation of their storage properties[END_REF].

In a previous work dealing with on the simultaneous generation and conditioning of hydrogen storage, under minimum pressure by continuous solar adsorption in activate carbon AX-21 to supply a fuel cell [START_REF] Cherrad | Conditioning of hydrogen storage by continuous solar adsorption in activated carbon AX-21 with simultaneous production[END_REF], the adsorption rate of H2 was modeled to be suitable for production rate.

___________________________________________________________________________ Page 4 of 19

The system used adsorbent beds as a storage unit that could play two roles, one as a fuel storage tank and the other as a fuel supply unit, as shown in Fig. 1. The previous work did not deal with the second role of the storage unit, which will be studied in this work.

While the high-density reversible hydrogen storage as a vehicular fuel is a biggest challenge [START_REF] Sathe | Hydrogen storage capacity of low-lying isomer of C24 functionalized with Ti[END_REF] and a reliable, efficient, and affordable medium that can store hydrogen reversibly at ambient temperature and pressure for on-board applications is also current challenge for researchers worldwide [START_REF] Rather | Preparation, characterization and hydrogen storage studies of carbon nanotubes and their composites: A review[END_REF], one of the main technical barriers associated with adsorbent-based hydrogen storage systems are their ability to effectively discharge hydrogen, depending on the real-time demand of the fuel cell [START_REF] Corgnale | Hydrogen desorption using honeycomb finned heat exchangers integrated in adsorbent storage systems[END_REF]. Hydrogen can be evacuated if necessary by increasing the temperature or decreasing the external pressure [START_REF] Jain | Hydrogen storage in Mg: A most promising material[END_REF]. However, the hydrogen release rate is not easy to control because it also varies with the temperature and composition of the materials [START_REF] Graetz | Controlled hydrogen release from metastable hydrides[END_REF]. Graetz and Vajo [START_REF] Graetz | Controlled hydrogen release from metastable hydrides[END_REF] have developed a kinetic equation of rates from a series of isothermal measurements to describe the relationship between temperature, hydrogen release rate and composition for aluminum hydride and metastable hydrides, during operation of a fuel cell at a controlled rate of about 1% wt%/hr. Although different kinetic models of hydrogen desorption for some adsorbent materials have been reported in the literature [START_REF] Kitagawa | Development of a kinetic model of hydrogen absorption and desorption in magnesium and analysis of the rate-determining step[END_REF][START_REF] Ma | Kinetics of hydrogen desorption from titanium hydride under isothermal conditions[END_REF][START_REF] Yang | Improved hydrogen absorption and desorption kinetics of magnesium-based alloy via addition of yttrium[END_REF][START_REF] Yang | Hydrogen absorption and desorption behavior of Ni catalyzed Mg-Y-C-Ni nanocomposites[END_REF][START_REF] Gabis | Kinetics of hydrogen desorption from the powders of metal hydrides[END_REF][START_REF] Terrani | The kinetics of hydrogen desorption from and adsorption on zirconium hydride[END_REF], these models did not consider the control of hydrogen flow rate as a function of real-time demand of the fuel cell.

The objective of this work is therefore to develop a single model based on general equation to control the rate of hydrogen desorption when supplying the fuel cell. In such system, (i) the release of hydrogen from activated carbon AX-21, applied as an adsorbent material, is induced by the rise in temperature, and (ii) the hydrogen flow relies on the pressure difference between the adsorbent storage bed and the hydrogen inlet into the fuel cell. The inlet pressure is considered to be set at atmospheric pressure, and the pressure in the adsorbent bed is modeled using the modified Dubinin-Astakhov equation [START_REF] Richard | Gas adsorption process in activated carbon over a wide temperature range above the critical point. Part 1: modified Dubinin-Astakhov model[END_REF], which proved to be a good analytical tool for describing hydrogen adsorption in a wide range of pressures and temperatures for several adsorbents [START_REF] Sdanghi | Application of the modified Dubinin-Astakhov equation for a better understanding of high-pressure hydrogen adsorption on activated carbons[END_REF] and in supercritical conditions [START_REF] Sdanghi | Modelling of a hydrogen thermally driven compressor based on cyclic adsorption-desorption on activated carbon[END_REF].

The study presents a general mathematical model and a calculation algorithm for estimating the desorbed hydrogen flow rate (DHR) as a function of supply time and temperature. It proposes a condition to stabilize the DHR by controlling the monitoring of the flow section of the hydrogen supply of the fuel cell, thus controlling the generated power and time.

System scheme

To release hydrogen gas to the fuel cell, a desorption process is required by heating the absorbent bed in ambient air at room temperature to increase the pressure of the stored hydrogen [START_REF] Cherrad | Conditioning of hydrogen storage by continuous solar adsorption in activated carbon AX-21 with simultaneous production[END_REF]. In our system, the adsorption is automatically synchronized with an intermittent solar energy source.

Thus, during the day, the solar irradiation is used to heat the beds and to desorb the hydrogen gas, while at night the water of hot fluid tank is used in place of the solar radiation.

The evaluation of heat exchange in the system can be based mainly on the work of Cherrad et al. [START_REF] Cherrad | Effect of heating time of adsorber-collector on the performance of a solar adsorption refrigerator[END_REF][START_REF] Cherrad | Transient numerical model for predicting operating temperatures of solar adsorption refrigeration cycle[END_REF]. The opening of the valve, which connects the outlet of the adsorbent bed and the fuel cell inlet (Fig. 1), induces the flow of gas into the fuel cell by pressure difference. However, this difference becomes stable when the two pressures are equal in the adsorbent bed and at the fuel cell inlet, and does not guarantee the release of all of the excess hydrogen stored hydrogen. The solution to release a maximum of fuel is therefore the increase of the temperature of the adsorbent bed.

Fig. 1. Hydrogen desorption scheme for supplying the fuel cell

Mathematical modeling of DHR

To evaluate and optimize the performance of hydrogen storage systems, it is necessary to take into account the thermodynamic characteristics, in particular temperature and pressure [START_REF] Xiao | Adsorption-desorption cycle thermodynamics for adsorptive hydrogen storage system[END_REF].

The adsorption storage phase occurs at a low temperature (liquid nitrogen temperature of about -195.79 °C) and at a pressure close to that of the atmosphere. Heating the adsorbent bed to room temperature Tamb increases the pressure. The higher the temperature, the more the pressure increases. The density of the bulk gas ρg is estimated using the state equation of an ideal gas at atmospheric pressure:

𝜌 𝑔 = 𝑃 𝑎𝑡𝑚 𝑅.𝑇 (2)
In the case of the activated carbon AX-21, as adsorbent material, the data of the modified Dubinin-Astakhov model parameters are given in Table 1.

Equation ( 3) above means that the pressure in the adsorbent can be controlled for the same constant values in table [START_REF] Abdalla | Hydrogen production, storage, transportation and key challenges with applications: A review[END_REF] and the different parameters which can be modified.

 Hydrogen desorption

The desorption of hydrogen from activated carbon AX-21 occurs at the same time as the fuel cell supply. It can be calculated from the dynamic pressure deduced from the difference of the static pressure between the adsorbent bed and the admission of the fuel cell (pressure Pi). Thus, the desorption rate of hydrogen desorption is expressed by:

𝑉 𝑑𝑒𝑠 = √ 2(𝑃-𝑃 𝑖 ) 𝑀 ℎ .𝜌 𝑔 (4)
The mass flow of hydrogen depends on the flow section S of the valve as shown in Fig. 1 and is described as follows:

𝑞 𝑑𝑒𝑠 = 𝑀 ℎ . 𝜌 𝑔 . 𝑆. 𝑉 𝑑𝑒𝑠 (5)
During desorption, the excess amount of stored hydrogen in adsorbent material decreases with time and can be estimated by the following equation:

𝑚 ℎ 𝑡+∆𝑡 = 𝑚 ℎ 𝑡 - 𝑞 𝑑𝑒𝑠 .∆𝑡 𝑚 𝑑 (6)

Numerical modeling of DHR

The process stops if one of the below conditions is verified:

|𝑃-𝑃 𝑖 | 𝑃 𝑖 ≤ 𝜀 (7)
Or

|𝑚 ℎ 𝑡 -𝑚 ℎ 𝑡=0 | 𝑚 ℎ 𝑡 ≤ 𝜀 (8)
Verification of the condition of equation ( 7) does not mean that all the excess amount of stored hydrogen has been released, because the pressure P (equation ( 3)) depends on the temperature of the adsorbent material and the amount of hydrogen (not yet described). To realize the maximum of the excess hydrogen, the temperature must be increased to increase the pressure until the condition of the equation ( 8) is verified. At this point, no hydrogen can be desorbed from the adsorbent material.

An iterative process with time variation is necessary for equation ( 8) to estimate the change in mass flow of hydrogen during the hydrogen supply of the fuel cell. The computational algorithm for the DHR model is shown in Fig. 2.

Start

Read data: n max , P 0 , α, β, V a , M h , R, Δt, ΔT, P atm , T(t=0), ε, m d , S, m h (t=0), ΔS, q des * t←0, T←(T(t=0)-ΔT) T←(T+ΔT) t←(t+Δt) Equation ( 7)

Equation ( 8)

No

No Equations ( 2), ( 3) and ( 4)

End * Only for a constant mass flow rate q des Equation ( 9)* Equations ( 5)

Recording: q des , t, T and S*

DHR variation

The variation of the desorbed hydrogen flow rate as a function of time was first estimated at the storage temperature T=-195.79 °C and at atmospheric pressure (considered for adsorption storage phase when the liquid nitrogen is used as cooling fluid), for a flow section set a at 1mm 2 and an initial excess of hydrogen stored mh=50 g for 1 kg of adsorbent material having a mass of 100 kg. When the adsorbent bed was heated to ambient air, the pressure increases automatically, while the pressure at the fuel inlet remains atmospheric (Pi=Patm).

Released excess amount

Fig. 4 shows that the excess amount of stored hydrogen is completely released from temperature above -50°C. This observation results from the significant increase of the initial pressure in the adsorbent material with temperature. Higher pressures obtained from ambient temperatures are useful for transporting hydrogen over a long distance. A positive result, seen in Fig. 4, is the ___________________________________________________________________________ Page 11 of 19 complete release of hydrogen at ambient temperature, which minimizes the heating energy supplied to the desorption phase.

Fig.4.

Excess amount of hydrogen released by the adsorbent material and the initial pressure relative to the desorption temperature (S=1mm 2 and initial mh t=0 =50 g/kg).

Desorption duration

Fig. 4 describes the relationship between excess hydrogen released and temperature. However, it may be rather useful to describe the desorption time of hydrogen as a function of temperature. As shown in Fig. 5, increasing the temperature accelerates the release of hydrogen into the fuel cell.

Which requires finding a way to stabilize it according to hydrogen demand of the fuel cell. The minimum desorption time occurs at high temperatures above 100°C, after a total pressure drop.

___________________________________________________________________________

Condition to stabilize the DHR

The above results show that the desorbed hydrogen rate is not stable and may not be useful for the fuel cell, since a constant mass flow rate is required to produce a constant power. Precise mass flow control therefore requires controlling the flow section S of the valve (equation 5) to stabilize the desorption.

For a constant mass flow rate qdes and a variable flow section S, the condition to be verified is given as follows:

𝑆 = 𝑞 𝑑𝑒𝑠 𝑀 ℎ .𝜌 𝑔 .𝑉 𝑑𝑒𝑠 (9)
In this way, the flow section S is controlled to ensure a constant supply of fuel. We have found that the hydrogen mass flow rate depends on the fuel cell supply required to produce the required power. High flow sections in square millimeters have been selected to allow the desorption of hydrogen in its maximum capacity. However, during a long supplying period, it is possible to reduce the desorbed hydrogen rate by decreasing the flow section to the square micromillimeter (μm 2 ). Then, the key to regulating the desorbed hydrogen rate and its duration is the flow section for providing the appropriate hydrogen release rates to a fuel cell.

Control of producing power of fuel cell

The relationship between the amount of hydrogen consumed and the current supplied by the fuel cell is expressed as follows:

𝑞 𝑑𝑒𝑠 = 𝐼 2𝐹
. 𝑀 ℎ [START_REF] Shirvill | Effects of congestion and confining walls on turbulent deflagrations in a hydrogen storage facility-part 1: Experimental study[END_REF] where, I is the current delivered by the cell and F is the Faraday's constant.

For a given voltage V, we can estimate the power Pe produced by the fuel cell as a function of desorbed hydrogen rate using equation [START_REF] Shirvill | Effects of congestion and confining walls on turbulent deflagrations in a hydrogen storage facility-part 1: Experimental study[END_REF] as follows:

𝑃 𝑒 = 2𝑉.𝑞 𝑑𝑒𝑠 .𝐹 𝑀 ℎ (11)
This means that the control of the valve of the control unit (Fig. 1) makes possible to adjust by automatic control the output power. The control unit can operate only by receiving a signal of four parameters as input data, namely: the temperature T and the pressure P in the adsorbent material, the pressure Pi at the admission inlet of the fuel cell and the power Pe selected by the user. A user control panel allows selecting the output power and display the power generation time before turning on the fuel cell. The control system program of the control unit can be based on the modeling of this study.

Conclusion

A hydrogen desorption rate modeling during fuel cell fueling has been proposed. The system process is based on hydrogen storage by adsorption in AX-21 activated carbon material. The system control program is based on the DHR model of this study, in which an iterative process with time variation was presented to estimate the change in hydrogen mass flow supply.

It can be seen that the mass flow rate of hydrogen decreases with time as the pressure difference between the adsorbent bed and the fuel cell inlet decreases. The pressure drop at a lower desorption temperature cannot release all the excess amounts of hydrogen stored. The excess of hydrogen released by the adsorbent material after a total pressure drop is greater at high temperature than at low temperature. However, the excess hydrogen is completely released at ambient desorption temperatures, which is a better result in minimizing the heating energy supplied to the desorption phase.

In addition, the minimum desorption time, after a total pressure drop, corresponds to higher temperatures. The results show that the DHR is not stable, which is why a condition has been proposed to stabilize the desorption rate by controlling the feed section of the fuel cell. However, another study also seems necessary to estimate the thermal energy supplied to the hydrogen desorption process.

F

  for maximum filling of the entire volume of the adsorption space mol/kg

19 

 19 Hydrogen adsorptionFrom the modified Dubinin-Astakhov equation, it can be estimated the excess amount of hydrogen (mh) in kg per kg of adsorbent material as a function of pressure and temperature, as follows.𝑚 ℎ = 𝑀 ℎ . (𝑛 𝑚𝑎𝑥 𝑒𝑥𝑝 [-[ 𝜌 𝑔 . 𝑉 𝑎 )

Fig. 2 .

 2 Fig.2. Computational algorithm of DHR model to estimate the variation of the mass flow of desorbed hydrogen during the supply of the fuel cell.

Fig. 3 (Fig. 3 .

 33 Fig.3(A) clearly shows the decrease in the mass flow rate of hydrogen other the time, resulting from the decrease in the difference between the pressures of the adsorbent bed and the fuel cell inlet , as shown in Fig.3(B).

Page 12 of 19 Fig. 5 .

 195 Fig.5. Duration of desorption as a function of temperature after a total pressure drop (S=1mm 2 and initial mh=50 g/kg).

Fig. 6 ___________________________________________________________________________ Page 13 of 19 Fig. 6 .

 6196 Fig.6. Variation of the controlled flow section versus stable desorbed hydrogen flow rates and supply time at ambient temperature (T=25°C).

Fig. 7 .

 7 Fig.7. Variation of maximum flow section as a function of stabilized mass flow rate for different ambient temperatures.

___________________________________________________________________________ Page 15 of 19 Fig. 8

 198 Fig.8shows the variation decrease of the generated power associated to the diminution of the mass flow rate as a function of the desorption time.

Fig. 8 .

 8 Fig.8. Variation of the produced power and the mass flow rate of hydrogen as a function of the desorption time (T=25°C, S=1mm 2 and V=0.7 V).

Table 1 :

 1 Modified Dubinin-Astakhov model parameters in the case of the AX-21 activated carbon[START_REF] Richard | Gas adsorption process in activated carbon over a wide temperature range above the critical point. Part 1: modified Dubinin-Astakhov model[END_REF].It results that for a given stored excess of hydrogen and a given temperature of the adsorbent bed, the pressure in the adsorbent bed can be deduced from equation (1) as follows:

	Parameter	Value
	nmax	71.6 mol/kg
	P0	1470 Mpa
	α	3080 J/mol
	β	18.9 J/(mol.K)
	Va	0.00143 m 3 /kg
	Mh	2.0159×10 -3 kg/mol
	R	8.314 J/(mol.K)