H Kanchev
email: hkanchev@tu-sofia.bg

V Lazarov
email: vl_lazarov@tu-sofia.bg

B Francois
email: bruno.francois@ec-lille.fr

F Colas

Keywords: Distributed generators, energy efficiency, energy management system, microgrid, optimization, renewable energy, power generation planning, smart grid, storage

In order to take full advantage of distributed generators, an evolution of the classical power system organization and management is also necessary. An aggregator of a residential urban electrical network can be considered by the Distribution System Operator as a stakeholder, which is able to control a cluster of local generators and loads with technical constraints for the connection with the remaining distribution grid and commercial contracts with outer electrical producers. This paper is focused on the design of the Microgrid Central Energy Management System which relies on a day ahead operational planning and an online adjustment procedure during the operation. A dynamic programming based algorithm is derived to solve the Unit Commitment Problem with a multiobjective function in order to reduce the economic cost and CO 2 equivalent emissions. The proposed energy management system is implemented into a SCADA and tested by using a hardware-in-the-loop simulation of the urban network. Economic and environmental gains are evaluated.

I. INTRODUCTION

esearch activities are more and more headed towards solutions for satisfying the ever growing energy demand. In order to ensure a continuous development in a sustainable way, a considerable portion of the electrical energy has to be generated by Renewable Energy Based Generators (REBG). One of their main drawbacks is the nonconstant nature of the primary energy source (solar irradiation, wind, …). Hence, the increasing of renewable generators penetration into the energy mix could cause difficulties for system operators in matching the power production and demand thus degrading the quality of power supplied to the customers and further causing disruptions in power supply [START_REF] Bhadane | Investigation for causes of poor power quality in grid connected wind energy -a review[END_REF], [START_REF] Robyns | Electricity Production from Renewable Energies[END_REF]. Moreover further investments have to be made in conventional generators to create additional power reserve for compensating the cyclic and stochastic nature of renewable energy.

Alternative to grid reinforcement could be a restructuring of the power system architecture and an increase of the share of Distributed Energy Resources (DER) that generate electricity at a local scale. Hence Micro Gas Turbines (MGT) based Combined Heat Production (CHP) and PV generators play an essential role for domestic small scale electricity generation. In this paper, only electrical production is considered.

A first restructuring strategy is to abandon feed-in tariffs and favor self-consumption of home produced electricity through incentives as experimented in Germany [START_REF] Ernst | Grid integration of distributed generation[END_REF]. But this energy policy restricts the energy sharing for neighborhood customers and limits the energy security and the prevention of over investment in production plants by effect of expansion.

A further strategy is to transform actual PV generators into controllable Active Generators (AG) in order to offer new flexibilities for energy management of electrical networks. Thanks to embedded storage technologies and a dedicated local control algorithm, this generator is able to deliver prescribed power references, power system services and can be dispatched to the distribution system operator.

But the current electrical transmission and distribution networks are rather passive and centralized from the supervision point of view. This situation makes difficult the coordination of DER in the grid [START_REF] Ernst | Grid integration of distributed generation[END_REF]. In order to coordinate power generation in an optimal way and to improve efficiency, reliability, security, the Smart Grid (SG) organization has to incorporate distributed intelligence and interactive communication at all levels of the electric network [START_REF] Gaviano | Challenges and integration of PV and wind energy facilities from a Smart Grid point of view[END_REF], [START_REF] Li | Smart Transmission Grid: Vision and Framework[END_REF].

A step towards the SG is to integrate locally REBG, conventional generators and loads in clusters called microgrids. These microgrids may be operated in islanded or connected mode (with the distribution grid) and can also provide ancillary services to the grid [START_REF] Fang | Smart Grid -The new and improved power grid: A survey[END_REF][START_REF] Farhangi | The path of the smart grid[END_REF][START_REF] Guerrero | Advanced control architectures for intelligent microgrids -part I: decentralized and hierarchical control[END_REF]. They must be locally aggregated and controlled by a Microgrid Central Energy Management System (MCEMS) (fig. 1). This aggregation can be considered by the Distribution System Operator (DSO) as a stakeholder, which is able to locally control a cluster of generators and flexible loads.

The main problem to overcome is to match locally the power demand and the production in an optimal way while minimizing the use of non-renewable energy sources for electricity generation and decreasing economic costs during operation (OPEX) [START_REF] Zhang | Robust Energy Management for Microgrids With High-penetration Renewables[END_REF]. Contributions of this paper are practical methods and solutions to design two stages of the MCEMS taking into account environmental and economic tasks. The first proposed stage is a day ahead operational planning for the minimization of CO 2 emissions and fuel consumption of an urban microgrid. This algorithm solves the Unit Commitment Problem (UCP) by the means of a dynamic programming with predictions of the available energy from PV generators, the power demand from the loads and the State Of Charge (SOC) of batteries inside active PV generators. Based on these data, this stage gives power references for micro gas turbines while maximizing the use of PV generators and reducing equivalent CO 2 emissions and the economic operating cost.

The second proposed stage is implemented during the day and consists in reducing variations due to the power uncertainty (from the PV production and load demand). This stage enables to retrieve day ahead calculated optimal economic costs and/or CO 2 emissions. An adjustment algorithm corrects, each 30 minutes, references coming from the day ahead operational planning, if changes in forecasted values occur. This adjustment algorithm is based on a sequential quadratic programming method. Both energy management functions can be implemented thanks to a communication network and the obtained complexity reduction of algorithms coming from the consideration of a cluster of generators and consumers.

This paper is organized as follows. In section II, the concept of AG is recalled. Then general functions for energy management of an electrical system are presented and ordered into the MCEMS for implementation in section III. As emissions and costs of primary energy come from micro gas turbines, manufacturer characteristics are derived to obtain a model useable for optimizations in section IV. Hence, the proposed day ahead operational planning is detailed in section V by adapting the formulation of the UCP to the studied power system, by mathematically expressing constraints and detailing the application of dynamic programming for the solving. Section VI is focused on the integration of an adjustment procedure during the day in order to compensate deviations from the day ahead forecasted data of loads and PV power. Finally implementation of this MCEMS into a supervisory control and data acquisition (SCADA) of an urban microgrid is exposed in section VII. Results from various tested optimization tasks are given to compare obtained CO 2 emissions and costs according to the use of distributed MGTs alone, integration of dispatched PV based active generators with the proposed operational planning and adjustment algorithm.

II. CONCEPT OF AG

Considered AG are based on a 3,6 kW PV generator, 106 Ah batteries and a 160 Farad supercapacitor bank (fig. 2). All components are connected to a common inner DC bus, which is interfaced to the electrical network though a three-phase inverter. Batteries are used to create a long term energy reserve and ultra-capacitors to supply power with very high dynamics. The inner instantaneous power balancing and power dispatching among the PVsource and storage units according to the storage level capacity and to the specific requirements/limitations of each source are performed by a Local Controller (LC), whose functions are detailed in [START_REF] Marnay | Optimal Technology Selection and Operation of Commercial-Building Microgrids[END_REF]. Real and reactive power references are received from the MCEMS. SOC of batteries and sensed powers at the connection point are sent to the MCEMS.

III. MICROGRID ENERGY MANAGEMENT SYSTEM

The MCEMS must assign power references and also other appropriate control signals to the DER units, conventional production units and controllable loads [START_REF] Kanchev | Energy management and operational planning of a microgrid with a PV-based active generator for Smart Grid Applications[END_REF]. The goals of the MCEMS algorithm are:

-to ensure uninterruptable power supply to the loads, -to maximize renewable generators in the energy mix, -to minimize the economic costs and the CO 2 equivalent emissions of gas turbines by setting their power references such, that they produce the minimum pollution and have minimum startups and shutdowns.

Many difficulties arise. First, the REBG power is variable and meets rarely the power request from the MCEMS. Moreover, there are numerous strategies for managing multiple gas turbines as auxiliary power sources for supplying loads (as example, a strategy may be to use the MGTs having the highest maximum power output as a priority source). Finally, other turbine constraints exist as the minimum power set point, the response time …

The microgrid management is analyzed through various functions that can be classified in a time scale. The long-term energy management elaborates a 24 hour-ahead operational planning including:

-the REBG production forecast taking into account the time dependency of the prime source, environmental impacts and cost of generation [START_REF] Ahn | Power scheduling of distributed generators for economic and stable operation of a microgrid[END_REF],

-the management of non-sensitive loads that may be disconnected/shed according to the supervision requirement, -the provision of an appropriate level of power reserve capacity according to the electricity market and the load demand forecast, -the maintenance intervals.

The medium-term energy management operates during the day and includes:

-the adjustment of forecasts for the power available from REBG and power demanded by the loads, -the adjustment of the long-term operational planning, based on deviations in the above mentioned forecasts from those predicted 24 hour-ahead.

The short-term power balancing is performed in the local controllers and includes the primary RMS voltage regulation and the primary frequency control [START_REF] Li | Participation in the frequency regulation control of a resilient microgrid for a distribution network[END_REF].

The long term operation schedule and the energy management can be mathematically expressed as an UCP. Due to the complexity of the problem, the required computation time may vary according to used optimization tools but the dynamic programming approach remains a good compromised choice [START_REF] Longenthiran | Formulation of Unit commitment (UC) problems and analysis of available methodologies used for solving the problems[END_REF].

IV. CHARACTERIZATION OF MICRO GAS TURBINES

A. Assessment of MGT Fuel Consumption

Fuel consumption of MGT can be assessed by using their partial load efficiency characteristics [START_REF] Boicea | Optimal operation of a microturbine cluster with partial-load efficiency and emission caracterisation[END_REF]. The energetic efficiency between thermal energy (F MGT_i in kWh thermal ) supplied to the gas turbine combustion chamber and electric energy output ( E MGT_i in kWh electric ) is defined as:

i MGT i MGT i F E _ _   (1)
i is the unit number (we will consider 3 micro gas turbines). The efficiency characteristic is a nonlinear function depending on the partial load ratio (fig. 3):

MAX i MGT i MGT i P t P t _ _ _ ) ( ) (   (2)
P MGT_i (kW) is the generated MGT electric power, P MGT_i_MAX (kW) is the rated MGT electric power.

Based on (1) and ( 2) the consumed fuel thermal energy for a 30 minutes operation (=1800s) at constant electrical output power is obtained (fig. 4): By considering a 0,04 €/kWh thermal cost for the consumed gas, the operational cost of each generator (i) is obtained from (3) as a function of the generated electric power P MGT_i (t):

   i MAX i MGT i i MGT P t t F _ _ _ ) ( ) (  (3)
)) ( ( _ t P f C i MGT i  (4)

B. Assessment of MGT Emissions

The CO 2 equivalent emissions are calculated by applying to the pollutant gases different weights corresponding to their global warming potential. For the assessment of emissions, the masses of the three exhaust gases, NOx, CO and CO 2 are evaluated in g/kWh as a mathematical function of the generated useful power [START_REF] Canova | Assessment of the Emissions due to Cogeneration Microturbines under Different Operation Modes[END_REF]:

    ) ( ) ( ) ( _ _ _ MAX i MGT i x i MGT x x P t t E t m   (5)
x  (mg/kWh electric ) is the emission factor (also called specific emissions) for the pollutant gas x to produce the generic useful electrical energy output they are absorbed in the earth's atmosphere slower than CO 2 . This means that these gases also contribute to the greenhouse gas effect. Global warming potential has been estimated, according to [START_REF]Intergovernmental Panel on climate change (IPCC)[END_REF]: 1 g of NO x has been considered equivalent to 298 g of CO 2 and 1 gram of CO equivalent to 3 g of CO 2 . The sum of the three characteristics (CO 2 , CO and NO x ) according to the MGT partial load ratio represents the CO 2 equivalent emissions of each MGT [START_REF] Canova | Assessment of the Emissions due to Cogeneration Microturbines under Different Operation Modes[END_REF].

In order to apply a multiobjective optimization procedure for a tradeoff between pollutant emissions and consumed fuel price, a price per ton of CO 2 equivalent emissions will be considered. According to economical researches in [START_REF] Abadie | European CO 2 prices and carbon capture investments[END_REF], the historical peak prices for trading a ton of CO 2 emission quota on the european market for industrials is equal to 30 euros. Practical and numerical applications in this paper will be developed with this price. According to a survey of the International Emissions Trading Association in 2013, an increase of emission prices in the European Union Emissions Trading Scheme to 10÷20 EUR/ton is highly probable by 2020 [START_REF]GHG market sentiment survey -2013[END_REF].

V. DAY AHEAD OPTIMAL OPERATIONAL PLANNING

A. Formulation of the Unit Commitment Problem

The 24 hour ahead operational planning is discretized in 48 time steps (t) of 30 minutes and power references are considered constant during each step. For all time steps the operational planning consists in selecting MGT to be used, determining the instant they should be committed and calculating the optimal power references of each generator. Power references of the three studied MGT are gathered in a vector:

x(t)= [P MGT 1 (t), P MGT 2 (t), P MGT 3 (t)] ( 6
)
i  is the state of each MGT during each time step (1 if the plant is running or 0 if the plant is shut down). Boolean MGT states are gathered in a vector:

u(t)= [δ 1 (t), δ 2 (t), δ 3 (t)] (7)
The general objective of unit commitment is to minimize the total operating cost of an electrical system, while satisfying all system constraints.

Power forecastings

Operational planing of MGT

u(t) x(t) (t) J  (t)

Weather forecast

Fig. 6. Scheme of the day-ahead optimal operational planning As the power industry goes to new restructured forms, the UCP must be adapted and applied to small DG clusters including REBG. So in this paper the UCP is used to formulate and solve our objective functions for cost minimization and CO 2 equivalent emissions reduction (fig. 6).

B. Objective functions

The CO 2 equivalent emissions of each generator and the cost of the consumed gas are expressed as a non-linear function of its power output and respectively called

)) ( (

_ _ 2 t P CO i MGT i and )) ( ( _ t P C i MGT i .
Penalties for startup and shutdown of the units are also considered. If one unit is shut down in the current time step and will run in the next time step, a startup penalty is applied. If one plant will not be committed at t+1 and it is running at t, a shutdown penalty is applied. The penalties avoid switching on and off the units because it increases the emissions and shortens the lifetime of the units. In this study, start-up penalty is considered equal to the consumed fuel cost during 5 minutes operation at full load. Shutdown penalty is considered equal to 2.5 minutes operation at full load. The startup and shutdown penalties are expressed by cost and CO 2 functions:

)

) ( ), 1 ( ( _ _ t t C i i i c pe    and )) ( ), 1 ( ( 2 _ t t C i i i co pe    .
The two objective functions to be minimized are defined as:

      48 1 3 1 _ _ _ )) ( ), 1 ( ( )) ( ( ). ( ) ( t i i i i c pe i MGT i i C t t C t P C t t J    (8)       48 1 3 1 _ 2 _ _ _ 2 2 )) ( ), 1 ( ( )) ( ( ). ( ) ( t i i i i co pe i MGT i i CO t t C t P CO t t J    (9) ) ( _ t P i MGT
is the generated power. Maintenance and management costs are outside of the scope of this paper.

C. Non-linear constraints

Day ahead planning strategies use renewable energy based production and load forecasts. With N active generators and M micro gas turbines, the power balancing between the loads ( L P ) and the generators ( n AG P _ and i MGT P _ ) in each time step must be performed with a maximum use of the "clean" PV energy. This is expressed as an equality constraint:

0 ) ( ). ( ) ( ) ( ) ( 3 1 _ 1 _         t P t t P t P t i i MGT i N n n AG L   (10)
The micro gas turbine loading level has to be higher than 50% of the MGT's rated power for improving efficiency (fig. 3) and reducing CO 2 equivalent emissions (fig. 5). Moreover the power in reserve must be equal to or larger than 10% of the generator rated power. The corresponding inequality constraint ( ) (t  in fig. 6) is expressed as:

i MGT i MGT i MGT P P P max_ _ _ max_ _ % 90 % 50   (11)
A last group of constraints refers to the microgrid operation mode [START_REF] Kanchev | Smart monitoring of a microgrid including gas turbines and a dispatched PV-based active generator for energy management and emissions reduction[END_REF], [START_REF] Kanchev | Unit commitment by dynamic programming for microgrid operational planning optimization and emission reduction[END_REF]. The constraints differ from one mode of operation to another (day/night, PV power available or not, AG's battery state of charge) and are taken into account in part E.

The MCEMS does not control the SOC of embedded storage units in various AG since it is managed by the corresponding local controller.

D. Application of the Dynamic Programming

Several approaches can be applied to implement an optimization procedure to solve the UCP by minimizing equ. [START_REF] Marnay | Optimal Technology Selection and Operation of Commercial-Building Microgrids[END_REF] and [START_REF] Kanchev | Energy management and operational planning of a microgrid with a PV-based active generator for Smart Grid Applications[END_REF]. The optimal solution of the Bellman's recursive equation for all time steps is used to construct the optimal solution of the overall problem [START_REF] Kanchev | Unit commitment by dynamic programming for microgrid operational planning optimization and emission reduction[END_REF], [START_REF] Hu | A new tree structure for weighted dynamic programming based stereo algorithm[END_REF].

For one state of MGT's (u(t)), the operational cost is expressed as: -the cost of electricity production during the time step [(t-1), t]

  3 1 _ )) ( ( ). ( i i MGT i i t P C t  (12)
-and the cost during the previous time step taking into account the transition cost due to the start or stop of generators

        3 1 _ _ )) ( ), 1 ( ( )) 1 ( , 1 ( )) ( ), 1 ( ( i i i c i pe t t C t u t F t u t u Tr   (13)
At step t , the unit commitment problem formulation for the studied system can be expressed in the form of the following recursive dynamic programming equation:

)) ( ), 1 ( ( )) ( ( ). ( )) ( , ( 3 1 _ t u t u Tr t P C t t u t F i i MGT i i       ( 14
) t Pref_MGT_24h _1 Pref_MGT_24h _2 Pref_MGT_24h _3 Pref_MGT_24h _1 Pref_MGT_24h _2 Pref_MGT_24h _3 Pref_MGT_24h _1 Pref_MGT_24h _2 Pref_MGT_24h _3 Pref_MGT_24h _1 Pref_MGT_24h _2 Pref_MGT_24h _3 Pref_MGT_24h_1 Pref_MGT_24h _2 Pref_MGT_24h _3 Pref_MGT_24h _1 Pref_MGT_24h _2 Pref_MGT_24h _3 Pref_MGT_24h _1 Pref_MGT_24h _2 Pref_MGT_24h _3 1 48 2 3 … … … … … … … … … Tr 1 (48)
Tr n (48)

Tr 1 [START_REF] Ernst | Grid integration of distributed generation[END_REF] Tr 2 [START_REF] Ernst | Grid integration of distributed generation[END_REF] Tr m [START_REF] Ernst | Grid integration of distributed generation[END_REF] Tr p ( 2)

… … Tr 1 (2) F 48 F 3 F 2 F 1 … … Fig. 7
. Principle of optimal path by dynamic programming

The optimal solution of the overall problem is obtained by selecting the optimal values of control variables in the previous equation for all time steps recursively from t=48 to t=1 as depicted in fig. 7 ( [START_REF] Singhal | Dynamic programming approach for solving power generating unit commitment problem[END_REF], [START_REF] Kanchev | Algorithm for operational planning and emissions reduction of a microgrid comprising gas turbines and PV-based active generators[END_REF] and [START_REF] Kanchev | Environmental and economic optimization of microgrid long term operational planning including PVbased active generators[END_REF]). In a same way, a similar formulation can be found to evaluate CO 2 emissions:

)) ( ), 1 ( ( )) ( ( ). ( )) ( , ( 3 1 _ 2 t u t u Tr t P CO t t u t F i i MGT i i        (15) With         3 1 _ 2 _ )) ( ), 1 ( ( )) 1 ( , 1 ( )) ( ), 1 ( ( i i i i co pe t t C t u t F t u t u Tr   (16)

E. Strategy for maximisation of renewable penetration

In order to give a priority to AG, their power references are first calculated. Then, power references of MGT's will satisfy the remaining power balancing if necessary. According to the available PV power and stored energy in batteries, three cases have been considered (fig. 8). During the day, if the predicted PV power ( ~), then power references of N active generators are set to limit the total generated power (all AG's are closed and have the same PV power sizing):

N t P t P h Load n h AG f ) ( ) ( 24 _ _ 24 _ _ Re  (17)
The energy surplus is automatically saved in batteries by the local control system of active generators and can be estimated: And power references of MGT must balance the remaining power:

1800 . ) ( ) ( ) ( ) 1 ( 1 _ 24 _ _ Re 24 _ 1 _ 1 _                 N n n h AG f h Load N n n Bat N n n
      N n n h AG ref h Load i i h MGT ref t P t P t P 1 _ 24 _ _ 24 _ 3 1 _ 24 _ _ ) ( ) ( ) ( (20)
Then the dynamic programming algorithm considers all possible states of MGT and, for each one, power references of MGT are calculated by a quadratic optimization algorithm (satisfying previous equation).

During the night, if the stored energy in batteries is enough to feed the predicted load demand, MGT are switched off and power references of AG's are set to limit the total generated power:

N t P t P h Load n h AG f ) ( ) ( 24 _ _ 24 _ _ Re  (21)
Calculation of battery SOC can be refreshed:

) ( ) ( ) 1 ( 24 _ 1 _ 1 _ t P t E t E h Load N n n Bat N n n Bat         (22)
Otherwise, power references of MGT's must balance the remaining power:

) ( ) ( 24 _ 3 1 _ 24 _ _ t P t P h Load i i h MGT ref    (23)
And, again, the dynamic programming algorithm considers all possible states of MGT's and their power references are calculated by a quadratic optimization algorithm.

VI. ON LINE ADJUSTMENT

The AG is capable of maintaining a prescribed power reference (

n h AG ref P _ 24 _ _
), received from the MCEMS, in the limits imposed by the actual SOC of the batteries. According to the daily predictions of the available power from the PV (

~

), a day ahead power production planning for the AG and for micro turbines is determined by the MCEMS. Sometimes the real situation (weather conditions, power demand by loads) could differ from the forecasted values. As a 10% power reserve is scheduled in power references of MGT ( equ (3) ), instantaneous deviations due to the power unbalancing can be handled by the primary controllers of local controllers. Current forecasting techniques allow us to have an updated forecast every 30 minutes. In practice, the current forecasted loads (

P

) is expressed as:

1 _ 1 _ 24 _ 1 _ ~      t PV t h PV t PV P P P (24) 1 _ 1 _ 24 _ 1 _ ~      t Load t h Load t Load P P P (25)
During the day the on line adjustment of the MCEMS consists in modifying AG power references each 30 minutes so that primary power reserves of MGT are not used and so that their produced power are as close as possible to the optimal operational point as calculated one day ahead. Deviations are taken into account by modifying power references of the generators according to the new situation:

1 _ 1 _ 24 _ _ 1 _ _ ~      t PV

P

is the AG power reference, calculated one day ahead by the MCEMS. A correction in power references for the micro gas turbines in the system may be also necessary. If the power deviation from forecasted values is larger than 10%, the algorithm has to switch on a turbine that was not planned to be run at t+1 in order to compensate the lack of PV power.

VII. TEST OF THE MCEMS AND APPLICATION

A. Implementation in a SCADA and Used Test-Bed

In industry, SCADA systems encompass the transfer of data between a central host computer and Transmission Control Protocol / internet Protocol (TCP/IP) remote terminal units (RTUs) or programmable logic controllers (PLCs) embedded in generators, gather and display information in a logical and organized fashion, carry out necessary analysis and control, execute sequential commands of remote devices, and transfer the information back to a central site.

The implementation of our proposed MCEMS in a SCADA needs to add centralized remote control and TCP/IP network. So, a test in real scenarios is required to evaluate the performance of the developed solution and to detect necessary adaptations before it will be implemented in the distribution system. For security reason, this test cannot be done directly in real situations onto a power system under operation. Moreover precise energy scenarios have to be played in order to valid the correct interactivity between the MCEMS and RTUs. So a test-bed has been developed and is based on a real time simulator (RTS), which is connected to the laboratory grid through a 30kW power amplifier (Fig. 9). The RTS is also connected to the SCADA, which executes the MCEMS (coded in the PcVue software), via a MODBUS communication network. The three micro gas turbines and the loads are simulated in RTS and the other equipment are real: the communication network, one active generator and the SCADA system. A real prototype of the AG is also connected to the grid and monitored by the SCADA (Fig. 10). The feedback loop of local control systems is closed by implemented current/voltage measurements at the clamp. Models of residential loads, but also RTU, local controllers and models of gas micro turbines are embedded in the RTS and executed in real time. This test-bed enables to reproduce identical demand scenarios under various sun conditions in order to validate the day ahead operational planning and the adjustment algorithm and to test various optimization goals.

MODBUS communication network

B. Presentation of the Application

The studied urban microgrid includes residential loads, two 30kW MGT, one 60kW MGT and twelve homes with 3,6 kW PV based AG's with embedded storage (fig. 11). All power generators and electrical loads are locally connected. So line losses and voltage drops can be ignored. As twelve generators must be considered and are very closed in the urban network (so they receive the same solar irradiation), the measured output power of the real single active generator is multiplied by 12. The microgrid central controller measures the microgrid state variables and dispatches power references (refreshed by the adjustment algorithm every half of an hour) to micro sources through a communication network. LC receive these power set points. In the same time they send various data, as example the sensed power production at the coupling point.

C. Experimental Results

First, a basic dispatching strategy without PV generators and optimization has been experimented with the load demand forecast (fig. 12). It consists in setting MGT power references proportional to the generator rated power while the added power references corresponds to the remaining power required by loads. Then PV based AG's are considered and their calculated global power references from the strategy for the maximization of the renewable penetration are shown on fig. 13.

The dynamic programing based optimization procedure has been performed with three different objective functions: minimization of the CO 2 equivalent emissions from the microgas turbines, the consumed fuel and a tradeoff between these two functions. Results in table I show that any of the three objective functions causes about 10% reduction in the total operational cost (the cost of CO 2 equivalent emissions and the cost of consumed fuel), compared to the same system without optimization of the operational planning.

CO2 cost Fuel cost Total cost Without optimization

50€

153€ 203€

Mono-objective (CO 2 )

35€ (-30%) 147€ (-4.2%) 181€ (-10.7%)

Mono-objective (Fuel)

36€ (-28.7%) 146€ (-4.6%) 183€ (-10.1%)

Multi-objective

35€ (-29.3%) 147€ (-4.1%) 182€ (-10.3%)

There is not a great difference in the overall cost between the different objective functions because, in fact, the CO 2 equivalent emissions and the consumed fuel are not really independent functions: the fuel consumption is used in the calculation of the CO 2 emissions, which is then used in the calculation of CO 2 equivalent emissions. Thus minimizing either the CO 2 equivalent emissions or the consumed fuel is enough for approximating very well the optimal system operation.

A statistical analysis of sensed powers from MGT has been done without optimization and with an optimized operational planning under the same scenario (figS. 14 and fig. 15). Results prove that the optimization algorithm se5ects power references such that MGT operate more frequently in the domain between 0.7 and 0.8 p.u. Operation in this region causes less CO 2 equivalent emissions, as depicted in fig. 5.

VIII. CONCLUSION

A day ahead operational planning has been proposed to calculate power references of DER in order to optimize CO 2 emissions and the operational economic cost. A dynamic programming approach is used for solving the unit commitment problem in the studied microgrid. Based on the load and PV production forecasts, the MGT characteristics, startup and shutdown penalties and the SOC of AG batteries, the algorithm searches by recursion the optimal path through all possible system states in the 24 hour-ahead operational planning. Deviations in forecasted values are compensated by an adjustment algorithm that recalculates power references for the generators if the difference between day-ahead forecasted values and one hour ahead forecasts is greater than the 10% power reserve. This algorithm enables to retrieve day ahead optimized values while satisfying the power balancing.

This MCEMS has been implemented into a SCADA and tested to manage an urban microgrid simulated in real time. Four different scenarios have been considered: without optimization of the operational planning and optimization with three different objective functions: minimization of the CO 2 equivalent emissions from the micro-gas turbines, minimization of the consumed fuel or a tradeoff between the two. Results demonstrate that using any of the three objective functions, a reduction of not less than 10% in the system total operational cost is achieved, compared to the same system without optimization of the operational planning.

Fig. 1 .

 1 Fig.1. SG based on the concept of microgrid integration

Fig. 2 .

 2 Fig. 2. Scheme of an AG including various storage technologies

Fig. 3 .Fig. 4 .

 34 Fig.3. Efficiency characteristic of MGT 2 (30kW)

Fig. 5 .

 5 Fig. 5. CO 2 equivalent emission characteristic of MGT No2 (30kW)

Fig. 8 .

 8 Fig. 8. Power reference calculation and dispatching.

  the required power of the loads (

Fig. 9 .

 9 Fig. 9. Schematic diagram of the Power Hardware In the Loop.

Fig. 10 .

 10 Fig. 10. Overview of our SCADA system.

Fig. 11 .

 11 Fig. 11. Schematic diagram of the urban microgrid

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Day ahead load forecast (kW) and PV power forecast in MPPT (kW)

Fig. 14 .

 14 Fig. 14. Occurrence of MGT 2 power set points without optimization

TABLE I DAY

 I AHEAD OPERATIONAL PLANNING RESULTS

Currently he is working with the Faculty of Electrical Engineering, Technical University and is head of the laboratory on Renewable Energy. His research interests includes electrical machines with electronic switching, small electrical machines and renewable energy sources (electrical aspects).


Dr Jean-Yves Yves Dieulot
email: jean-yves.dieulot@polytech-lille.fr

Frederic Colas
email: frederic.colas@ensam.eu

PhD Lamine Chalal
email: lamine.chalal@univ-lille1.fr

Genevieve Dauphin Tanguy
email: genevieve.dauphin-tanguy@ec-lille.fr

PhD Frédéric Colas

Professor, PhD Geneviève Dauphin- Tanguy

Professor A Kugi

Frédéric Dieulot

Lamine Colas

Geneviève Chalal

Dauphin-Tanguy

Fréderic Colas

Geneviève Dauphin-Tanguy

Event-triggered variable horizon

Keywords: Economic model predictive control, power system control, supervisory predictive control, renewable hybrid system, hardware in the loop, variable receding horizon Economic model predictive control, power system control, supervisory predictive control, renewable hybrid system, hardware in the loop, variable receding horizon

published or not. The documents may come

Introduction

Model Predictive Control (MPC) is a control algorithm that minimizes an objective function over a receding horizon. It is well suited for multivariable or large-scale systems which exhibit large time-constants, long time-delays or middle-range predicted disturbances, and is known to be robust with respect to forecast uncertainties and unmodelled dynamics [START_REF] Scattolini | Architectures for distributed and hierarchical model predictive control -a review[END_REF]). MPC has been applied in the field of hybrid power cell supervision. It allows to deliver an overall smooth power when intermittent and renewable power sources (wind turbines, solar panels) are combined with conventional sources and storage units (batteries, supercapacitors, etc.). Hence, a predictive supervisor helps to promote the integration of the renewable sources into the grid [START_REF] Valenciaga | Supervisory control for a standalone hybrid generation system using wind and photovoltaic energy[END_REF] , [START_REF] Guerin | Supervisory control design for systems of multiple sources of energy[END_REF], Chalal et al. (2012) , Liu et al. (2011a), [START_REF] Qi | Distributed supervisory predictive control of distributed wind and solar energy systems[END_REF]). Generating a smooth active power in a short time window, in a time scale around 10-30 seconds, is mandatory for a power cell (see a list of demonstrators in [START_REF] Daguzan | Managing intermittency and storage: Update of the French demonstrators[END_REF]). To set the context of this supervision in the power generation frawework, the sampling period of local level power converters controllers is below the millisecond [START_REF] Steinke | Switching frequency optimal pwm control of a three-level inverter[END_REF]). On the contrary, tertiary regulation has a time scale around 15 minutes. The supervisor provides the reference trajectories for the controlled units (in this case a microturbine and a storage unit). Moreover, unlike other algorithms, predictive controllers can consider the planned power reference and a short-time forecast of the intermittent sources powers which generally improves the results. When the predictive criterion is based on economic considerations (e.g. tariff policies, fuel or CO 2 emission cost, storage unit depreciation), it is very easy to tune the supervisor. Moreover, extra components can be integrated into the cell. It is sufficient to consider their dynamical models and to add their economic contribution to the objective function, which can depend on the operating regime as in the case of wind turbines (Chalal et al. (2012), [START_REF] Qi | Distributed supervisory predictive control of distributed wind and solar energy systems[END_REF] ). Hence, a predictive approach allows for an (economically) optimal and modular way to tune a hybrid power cell supervisor. On the contrary, rule-based algorithms need to be fully designed anew in case of component modification or cell reconfiguration. The choice of the receding horizon results of a compromise between the intermittent power prediction accuracy, which decreases with the horizon length, and the benefits of using a long-range control algorithm. Typically, in (Chalal et al. (2012)), this horizon was fixed around 10 seconds. However, whenever an abrupt change of the renewable power occurs, a possible false prediction can hamper the results given by the MPC algorithm. As an example, this could happen for solar panels when it quickly clouds over, or for wind turbines with squally conditions. In this case, it is possible to improve the predictions by using a shorter horizon. The so-called Variable-Horizon model predictive control (VH-MPC) which varies the horizon length according to a constrained optimization problem is well-suited to handle time-varying disturbances. The VH-MPC algorithm drives the system state to a closed set in finite time irrespective of bounded disturbances [START_REF] Michalska | Robust receding horizon control of constrained nonlinear systems[END_REF] , [START_REF] Scokaert | Min-max feedback model predictive control for constrained linear systems[END_REF], [START_REF] Richards | Robust variable horizon model predictive control for vehicle maneuvering[END_REF] , [START_REF] Shekhar | Robust variable horizon mpc with move blocking[END_REF]) . It is possible to reduce the complexity of the algorithm by setting the variable horizon to predefinite values according to the onset of a disturbance. The horizon is kept to a nominal value when the prediction is accurate enough. Otherwise, it is triggered to a shorter value. This avoids important prediction biases, when the occurrence of a large disturbance is detected [START_REF] Eqtami | Eventtriggered strategies for decentralized model predictive controllers[END_REF]). Note that this policy is quite different from multi-rate predictive control which triggers during a discrete step or uses adaptive step control algorithms (see e.g. [START_REF] Kowalska | An approach to variable time receding horizon control[END_REF] and references therein). To sum up, such an event-triggered predictive controller has not been used so far for the supervision of a hybrid power cell. This paper presents the design of an Economic Supervisory Predictive (ESP) controller, which original criterion integrates tariff policies, storage units and conventional sources costs, whereas other papers generally propose a standard quadratic criterion. Next, it is shown how predictions can be integrated into the supervisor. Compared with existing papers, the economic tradeoff between the use of a longer horizon which may improve the overall cost and the accuracy of the renewable sources power prediction is evaluated, which shows the relevance of using a variable receding horizon strategy. The design of such an event-triggered variable horizon predictive supervisor, coupled with a detection algorithm is addressed next. Eventually, the algorithm will be applied on a hybrid cell including a storage unit, solar panels and an emulated gas turbine. Real-time experiments in the framework of a Hardware In the Loop platform support the relevance of the proposed approach.

HYBRID CELL ECONOMIC SUPERVISORY CONTROL

Power cell, HIL platform and component models

The hybrid renewable power plant embeds a photovoltaic system (108 modules BP solar 3160 with a power of 160 W each, connected to a 3-phase grid via a Fronius IG30 inverter), a storage unit connected to the grid via electronic converters, and a gas microturbine. This microturbine is emulated by a real time digital simulator RT-Lab which is able to run real-time models on a multi-CPU computer (Fig. 1). A power amplifier generates a real power signal from a model (e.g. the microturbine). The simulated sources can be blended with real power measurements, achieving the Power Hardware In the Loop principle. Another interest of the real-time simulator is the ability to emulate and reuse a previously measured power (typically here a solar power profile). Hence, a fair comparison of different algorithms is possible (see Chalal et al. (2012) for additional details).

The storage unit allows to absorb the fluctuations of the renewable source power. A simplified model was used to represent the battery storage [START_REF] Chungpaibulpatana | A simulation model for predicting the performance of a solar photovoltaic system with alternating current loads[END_REF]):

P bat = 1 τ p s + 1 P batref (1)
where P batref , P bat are the storage unit reference and real powers. It is assumed that the storage unit State Of Charge (SOC) varies (in the working range) as the integral of the power. In case of State Of Charge saturation, if t 0 P bat dτ > E max or t 0 P bat dτ < E min , then P bat = 0. E max and E min are respectively the maximum and minimum of stored energy. τ p = 5 s is the time constant of the battery.

This model is implemented using a dedicated test bench with supercapacitors. These consist of 6 Maxwell modules in series, the characteristics of each being 48 V, 160 F. The power reference is generated through a power amplifier. The supercapacitors are connected to the grid (Fig. 2) .

The gas microturbine power is generated by a model which features a simplified representation of a turbine in closed loop. The supervisor does not manage the inner control loops and only provides the power references. A gas turbine consists of a gas compressor, turbine, a recuperator, a permanent magnet synchronous generator, and a power electronics interface. In this case, the microturbine works in operating mode. Start-up is discarded and only a turbine speed PID control loop is considered (see details in Chalal et al. ( 2012)). The torque C m is a function of the rotor speed N and of the fuel flowrate W f :

C m (t) = 1.3(W f (t) -0.23) + 0.5(1 -N(t)) (2)
All variables are expressed Per Unit, allowing scale-up, according to the model in [START_REF] Pai | A design and operation of power converter for microturbine powered distributed generator with capacity expansion capability[END_REF]. The term 0.23 accounts for the need for a significant fraction of rated fuel to support self-sustaining, no-load conditions. One has the following transfer function:

N(s) = 1 T l s + 1 C m (s) (3) W f (s) = 1 T f s + 1 a cs + b V P (s) (4)
where the fuel control valve reference position V P is:

V P (t) = 0.77N(t).F c (t) + 0.23. ( 5
)
and the speed controller is:

F c (s) = 25 0.05s + 1 (N ref (s) -N(s)) (6)
where F c is the per unit fuel command per unit speed, the turbine power P mtg is the product of speed and torque. T l , T f , a, b, c are respectively the turbine rotor time constant, fuel system time constant, and fuel system transfer function parameters. This model is used for validation and implemented into the real-time simulator. In case of "stiff" turbines, it is possible to simplify the closed-loop model of the microturbine. The model of the microturbine reduces finally to a fourth-order linear model, which is better fitted for Model Predictive Control design:

P mtg (s) = 1 (T 1 s + 1)(T 2 s + 1)(T 3 s 2 + T 4 s + 1) P mtgref (s) (7)
where T 1 = 0.05 s, T 2 = 0.15 s, T 3 = 0.22 s, T 4 = 0.27 s. This coarse linear model shows a good agreement with the full model in the region where the microturbine operates when embedded in the hybrid cell. Simulations can be found in Chalal et al. (2012), showing results which are consistent with that found in the literature [START_REF] Martinez | Tools for analysis and design of distributed resourcespart iv: Future trends[END_REF]).

Economic Supervisory Predictive control

The design of an Economic Supervisory Predictive controller, presented in Chalal et al. ( 2012) is recalled. The objective consists of generating the power references of the controllable elements belonging to a hybrid power cell (in this case the storage unit and the microturbine powers P bat , P mtg ). This is achieved by minimizing a global cost over a receding horizon H, notwithstanding the fluctuations of the renewable power source (here the solar arrays, see Fig. 3).

The time-dependent cost function can be split up into the gas consumption cost C f uel , carbon dioxide emissions tax C emissions , the battery cycling cost C cycling , and a penalty for not reaching the reference power P d imposed by the grid manager.

J(t) = t+H t α(P d (τ ) -Ppv (τ ) -P bat (τ ) -P mtg (τ )) 2 + ... (8) + C f uel (τ ) + C emissions (τ ) + C cycling (τ ) dτ
or in a discrete-time fashion:

J = k α(P d,k -Ppv,k -P bat,k -P mtg,k ) 2 + C f uel,k + C emissions,k + C cycling,k (9)
α is a penalty term imposed by the grid manager. k corresponds to the k th time period over the time interval [t, t + H].

The delivered power P is the sum of all the components' powers. It incorporates the intermittent solar power predictions Ppv . In our case, these predictions are supposed to be linear functions of time. Their parameters are estimated by a Least-Square method at every instant t. The predictive control strategy embeds the dynamical closed loop models of the controlled sources' powers and the corresponding operational costs. The cost function is derived only from technical and economic considerations or policies. One can consider marginal costs as in [START_REF] Edlund | Hierarchical model-based predictive control of a power plant portfolio[END_REF]); however, over the short prediction horizon, the parameters are considered to be constant (e.g. fuel pricing does not vary). It is also possible to enrich the criterion with further considerations on power quality. For example, one objective of the MPC controller is to smooth the output power and follow a reference power P d in the paper. The output power should thus be smoother than if the renewable source was directly connected to the grid. Moreover, in a full scale operation, such a hybrid power cell could participate to frequency support: the reference power could contain a ∆P d term which comes from a droop in frequency as depicted in (De-Brabandere et al. ( 2007)) Another possibility would be to add an economic penalty in the criterion. The cost function of a gas microturbine is assumed to be a linear function of the power P mtg (Dufo-Lopez and Bernal-Agustin ( 2005)) :

C f uel (t) = C F [BP N mtg (t) + AP mtg (t)] (10)
where C F is the fuel price (e/l), A = 0.246 l/kW and B = 0.08415 l/kW are the fuel curve coefficients, P N mtg is the rated power (kW), P mtg is the generator output power (kW), C f uel is the fuel consumption cost (e).

The emission cost C emisions is assumed to be proportional to the CO 2 emissions weight, with:

C emisions (t) = 0.017 * E(P mtg (t)) (11)
The volume of these emissions is modelled as a quadratic function of the generator power P mtg [START_REF] Johansen | A qualitative approach to economic-environmental dispatch[END_REF], where the A i , i = 1, 2 are constant coefficients:

E(P mtg (t)) = A 0 + A 1 P mtg (t) + A 2 P 2 mtg (t) (12)
The operational cost of the microturbine is a quadratic function of the microturbine power. Battery cycle lifespans vary according to their type and monitoring policy. A mean cycle, with full charge and discharge, is considered, given by the ratio of total cost of the storage unit and the lifespan. cycle cost = cost of the battery number of expected cycles (13)

It is assumed that locally, a portion of a cycle cost corresponds to a linear fraction of the total cost of one cycle. The State Of Charge Q varies linearly:

C cycling (t) = λ dQ(t) dt (14)
where λ can be calculated from equation ( 13).

In addition to the cost determination, operational constraints on the power amplitude or power change rate [START_REF] Liu | Supervisory predictive control of standalone wind/solar energy generation systems[END_REF]) are considered. The load level of the microturbine should not reach less than 50% for efficiency reasons:

15 kW ≤ P mtg (t) ≤ 30 kW (15)

Inner-loop constraints on the fuel flowrate F mtg yield in turn:

0 ≤ F mtg (t) ≤ F max (16)
The battery operational constraints consider power and capacity limitations:

-5 kW ≤ P batref (t) ≤ 5 kW.

(17)

0.5 ≤ Q(t) ≤ 0.8 (18)
and

| dQ(t) dt ≤ 0.1| (19)
The overall predictive control problem can be put under the following form: min

P batref (k),P mtgref (k) J(k) s. t.C(k)∀k equations( 1, 7) (20)
where C(k) are the operational constraints. Since the overall criterion is a quadratic function of the states (i.e. the powers), the overall optimization problem is a (QP) quadratic program and can be solved easily. As recommended for most predictive control algorithms, only the first value of the sequence of references provided to the supervisor is fed to the controlled units.

VARIABLE HORIZON PREDICTIVE SUPERVISOR

Event-triggered Supervisory Predictive control

The choice of the prediction horizon H is fundamental for the success of the MPC supervisor. When the horizon is too long, short-term predictions of the solar power become less accurate. If this horizon is too short, the predictive algorithm shows poorer results because of the slow reactivity of the controlled sources (Chalal et al. (2012)). Previous works have shown that the economic criterion can be improved if the predictions are embedded into the supervisor. As was said before, a linear trend was chosen to represent the solar power prediction. For every time t, the relative prediction error is

ǫ s = (1/H) t+H t |P solar (τ ) -P predicted solar (τ )| P solar (τ ) dτ (21)
Table 1 provides an example, for a standard solar power profile, of the mean value of the relative prediction error MRAPE ǫ s as a function of the horizon. One can see that the forecast errors increase with the prediction horizon, which is quite obvious. The following proposition considers the effects of a time-linear drift in the output. It is assumed that this drift is not compensated by the predictive controller. This case may happen, for example, when there exists a mismatch in the trend of the solar power predictions. The slope of the linear drift is supposed to depend on the value of the horizon and to account for the precision of the predictions. It is assumed to be constant over a time interval [0...T ]. When the prediction is perfect, it will be assumed that the mean economic criterion is a decreasing function of the horizon H. Indeed, in this case, previous results have shown that using a longer receding horizon yields better results for the supervision of steady solar power P pv profiles (Chalal et al. ( 2012)).

Lemma 1. (Leibniz integral rule) Let f (x, t) be continuous over a time interval [a(x), b(x)], then

∂ b(x) a(x) f (x, y)dy ∂x = f (x, b(x)) ∂b(x) ∂x -f (x, a(x)) ∂a(x) ∂x + b(x) a(x)
∂f (x, y) ∂x dy

Proposition 1. Consider a dynamical SISO state-space model ẋ = f(x, u), ŷ = h(x) and a smooth control sequence which minimizes, at each time t, the criterion J(t, H, ŷ) where

J(t, H, w) = 1 H t+H t r(τ )dτ, t ∈ [0, T ], r(t) = |y ref (t) -w(t)| + g(x, u).
It is assumed that the real output y obeys, over the time window [0, T ] :

ŷ(t) -y = a(H)t, a(H) > 0.
Then, one has

1 T T 0 ∂J(t, H, y) ∂H dt = 1 T T 0 ∂J(t, H, ŷ) ∂H dt + a(H)/2 + (T /2 + H/2) ∂a(H) ∂H Proof.
Using the lemma, one differentiates the criterion without linear drift

∂J(t, H, ŷ) ∂H = 1 H r(t + H) - 1 H J(t, H, ŷ)
The partial derivative of the criterion with drift is now

∂J(t, H, y) ∂H = 1 H r(t + H) + 1 H a(H)(t + H) - 1 H J(t, H, ŷ)- 1 H 2 t+H t a(H)τ dτ + 1 H t+H t ∂a(H) ∂H τ dτ = ∂J(t, H, ŷ) ∂H + 1 H a(H)(t + H) - 1 2H (2t + H)a(H) + (t + H/2) ∂a(H) ∂H = ∂J(t, H, ŷ) ∂H + (t + H/2) ∂a(H) ∂H + a(H)/2
One can now integrate over [0, T ], which gives the final result

The proposition shows that when the mismatch on the solar power trend is important, the benefits of using a larger horizon are backed by the term which depends on the partial derivative of the slope of the linear trend. It is difficult to have a clear account of the mean value of the criterion without drift. However, as was said before, it is likely to be a decreasing function of the horizon. For unsteady profiles, the precision can be insufficient and a short time horizon will be preferred. Hence, a heuristic controller is proposed, using a decision variable d v :

if d v < κ then set controller with H = 3 s (22) if d v > κ then set controller with H = 11 s. ( 23
)
where κ is an appropriate detection threshold. The values of H are chosen according to preliminary simulation results for steady and unsteady profiles.

Detection methods of solar power change

Different online algorithms allow to detect a change in time series. Some algorithms may detect small variations (e.g. the CUSUM algorithm). However, only high variations are sought in the presented case. This advocates for simple detection policies like Shewhart control charts which are recalled briefly in [START_REF] Basseville | Detection of Abrupt Changes -Theory and Application[END_REF] ). Consider independent variables y i with a probability density p θ , where the value of θ changes from θ 0 to θ 1 at some instant, and let the decision function be S n = n i=1 ln p θ 1 (y i ) p θ 0 (y i ) . The algorithm reads simply θ = θ 0 if S n < κ else θ = θ 1 , where κ is a "wellchosen" threshold. Consider a Gaussian distribution of mean µ, standard deviation σ, and p θ = 1 σ √ 2π e (yµ) 2 2σ 2 . When a change in mean occurs, one has

S n = n i=1 µ 1 -µ 0 σ 2 (y i -µ 1 +µ 0
2 ). After a few manipulations, the detection algorithm corresponds, at time k, to:

1 n nk i=n(k-1)+1 y i -µ 0 ≥ κ σ √ n , (24)
In this case, some adaptations are needed. The algorithm will be used to detect a change of slope in the solar power signal following a five step procedure:

• differentiate the signal • at each sample k, estimate recursively the linear trend for the solar power derivative

• use the detection algorithm with a threshold factor κ = 2

• Set a relaxation time of 200 s to enforce the controller to keep a small horizon for some time after the detection of disturbance to avoid frequent switching

• Reinitialize after detection Note that differentiation of a signal amplifies the noise. Of course, it would be possible to use low-pass filters instead of changing the horizon; however, changing the horizon allows for a more explicit and discarding policy for handling severe disturbances.

EXPERIMENTAL RESULTS

Previous experiments in (Chalal et al. (2012)) have shown that a ESP controller could outperform a rule-based controller used as a standard supervisor for a hybrid cell , and which is not detailed here. The final cost is reduced by more than 30 %. The main goal of the current experiment is to show the relevance of using an event triggered variable strategy for this supervisory predictive control. The solar power profile which is considered is given in Fig. 4. This profile is quite steady, let two negative power disturbances appearing after a time of 450 seconds. After disturbance detection, the event-triggered algorithm will switch the horizon from a high value fit for steady operations, that is H = 11 s, to a shorter horizon H = 3 s. Note that there will be a relaxation time to avoid frequent switching, as that the horizon is set back to H = 11 s only at a time of 750 seconds. The switching instants are shown in Fig. 4. The κ threshold is chosen by a trial-and-error method as twice the standard deviation of the signal and allows for a rapid detection. An interesting aspect of the HIL architecture presented in Fig. 1 is that this real solar power profile can be stored and later emulated by power amplifiers, and then reused to allow a fair comparison between the event-triggered strategy, a short-horizon strategy and a long horizon strategy. This allows clearly to a quantitative comparison of different strategies and to provide a fair account of the benefits of using an event-triggered variable horizon predictive algorithm, which would not be possible if different solar power profiles were used for each algorithms. Providing the overall cost for each strategy would be, however, not very useful to the control practitioner, as the weighting of the components of the economic criterion are subjected to a modification of the economic context (e.g. a rise in fuel pricing or a change in tariffs). Table 2 summarizes the different costs (Emission, fuel, cycling and total costs) for a short, long and a variable horizon for the solar profile given in Fig. 4. Before drawing a finer comparison for each strategy, it is necessary to investigate the influence of the horizon value on the behavior of the microturbine and storage unit powers. Specifically, explaining the behavior of the different physical components is important to understand the benefits of using a time-variable horizon. Fig. 5 shows the profiles of the overall power produced by the cell for a short, a long and a variable horizon strategy. In each case, the optimization of the economic criterion realizes a tradeoff between tracking (corresponding to a tariff policy) and other costs (battery cycling and turbine workout). Hence, the power tracking is never perfect. As could be expected, when using a short horizon, tracking the total power is not very good when the solar power is steady, because predictions are carried out over a short time-window. At the onset of a disturbance, however, the reaction is faster, and prediction errors are not amplified (contrary to the use of a long time-window), even if one can notice power peaks (overshoots)with larger amplitudes. The overall tracking cost is thus not very good for both steady horizon strategies, as can be seen in Table 2. The moving horizon strategy combines the advantages of the two methods, that is, ensure the same tracking performances than the long horizon strategy over steady conditions, and being more reactive during unsteady conditions.

Fig. 6 shows a comparison of the storage unit real-time State of Charge for the three strategies. The microturbine power shown in Fig. 7 exhibits an inverse behavior to that of the battery power profile. For a long horizon, when the solar power is steady, the storage unit and microturbine powers keep an almost constant value. When the solar power is unsteady, prediction errors will cause the microturbine and storage unit to be called for to maintain the total power close to the reference, hence dramatically increasing cycling and fuel costs. When the horizon is kept to a low value, the MPC performance is poorer when the solar power is steady, degraded tracking performance leads also to an increase of the control effort. When solar power fluctuations occur, prediction errors are not so important as in the long horizon case, the control is more reactive, leading to an improvement in the monitoring of the microturbine. The storage unit is more called for to smooth peaks or rapid variations of the solar power, hence enduring larger SOC variations which will increase the cycling cost.

To sum up, Table 2 presents the global results obtained for this experiment. Fuel cost is nearly the same for a short and a long horizon; however, it has been shown, through the previous analysis, that the underlying grounds were, in each case, quite different. Making the most of the two strategy by using an adequate switching strategy, the moving horizon algorithm has been able to save up to 15 % fuel and 14.5 % emissions. The battery cycling cost for a moving horizon strategy results in a weighting average between the long horizon periods, where it keeps a lower value, and the short-horizon period, when the SOC of the storage unit fluctuates much more. As a result, the total cost savings can reach 12 %, the shorter horizon being a reference, or 10 % with respect to a long horizon strategy. In any case, the moving horizon strategy proves to be effective.

Conclusion

A hybrid cell, which consists of a storage unit, a gas microturbine and solar panels is controlled with an Economic Predictive Supervisor. This strategy relies on an economic criterion. This criterion is a function of the grid management policy, the fuel, gas emission and battery cycling costs, all of which have been presented. The overall cost, which is non standard, is optimized over a receding horizon, feeding the battery and storage unit with proper references. The implementation and tuning of such a predictive supervisor is quite easy, since it involves, in this case, the optimization of a criterion at each step time, with any constraints and components of the criterion being based on physical or financial models or considerations. There is no arbitrary weight left to the choice of the user. It has been shown that this controller outperforms a rule-based supervisor by 30 % However, the results may depend strongly on the prediction horizon of the predictive supervisor. Indeed, when the solar power is well predicted, results are better for a long-range predictive supervisor because the dynamics can be controlled over a wider time-window. When the solar power is unsteady, the predictions are less accurate and using a shorter horizon avoids to deteriorate the performance. Theoretical results show that the choice of a receding horizon for the predictive supervisor results of a tradeoff between performance and prediction accuracy. A new moving horizon strategy has been designed which takes advantage of the two situations. When solar power fluctuations occur, the horizon is switched to a lower predefined value. The switching instants are detected through the use of Shewhart charts. Experimental results are obtained using a real-time digital simulator, which allows to emulate virtual components, to blend them will real devices and use the same solar power profiles to get a fair comparison of algorithms. These results have shown the relevance of such an approach, as the overall cost savings reach 12 % with respect to a short horizon and 10 % with respect to a long horizon strategy. Fuel savings are by 15%. This approach remains very simple as only very simple decision algorithms are used along with the model predictive supervisor, and, as the time period is 1 second, no very powerful computation devices are required. A generalization of these results, including a study of the sensitivity to the different cost parameters, will be the topic of future research. Indeed, the methodology can easily be scaled up to a real world application, as welcoming an extra component needs only to know its dynamics and its economic contribution.

Figure 1 :

 1 Figure 1: Hybrid Cell and Power Hardware In the Loop architecture

Figure 2 :

 2 Figure 2: Storage unit test bench

Figure 3 :

 3 Figure 3: Simplified control scheme

Figure 4 :

 4 Figure 4: Solar power profile used for real-time experiments

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Total hybrid cell power under Predictive supervision with constant and moving horizons

Table 1 :

 1 Mean relative absolute prediction error for a standard solar power profile as a function of the prediction horizon

	Horizon H [s] 3	5	7	9	11
	MRAPE [%] 0.4 1.7 4.5 8.6 13.9

Table 2 :

 2 Cost function components as a function of the predictive supervisor horizon value

		H = 3 [s] H = 11 [s] Moving Horizon
	Emission cost (e)	0.115	0.111	0.098
	Fuel cost (e)	1.617	1.569	1.376
	Battery cycling cost	0.181	0.079	0.118
	(e)
	Tracking cost (e)	0.047	0.035	0.039
	Total cost (e)	1.86	1.794	1.631
	Savings (%)	100	96	88

Table 3 :

 3 List of Acronyms

	CUSUM	CUmulative SUM control chart
	ESP	Economic Supervisory Predictive
	HIL	hardware In the Loop
	MRAPE	Mean Relative Absolute Prediction Error
	SOC	(Battery) State Of Charge
	VH-MPC Variable Horizon Model Predictive Control

Acknowledgements

The authors would like to thank the anonymous reviewers for their indepths comments and M. Mohit Makkar for fruitful help with the manuscript.

Response to Reviewers, Control Engineering Practice, Manuscript CONENGPRAC-D-14-00018R1 Editor in Chief: Pr. A. Kugi Please find the following answers to the reviewers' comments. Yours Sincerely Jean-Yves Dieulot Editor: The article is suitable for publication, but the author is requested to consider the further comments of the reviewers. Specific changes are not required but some issues have been raised for the author's further consideration in the final manuscript.

R The authors would like to acknowledge the in-depth and valuable work of the reviewers that helped greatly to improve the paper. We replied to some queries by reviewer #2 and added a few related comments in the paper.

Reviewer #2: This article has clearly been improved in this revision. I do have a few additional comments.

Introduction:

-Regarding economic costs, This might more complicated for renewables like wind turbine turbine were operation influences the costs of energy. R. We totally agree. In the PhD work of Dr. Chalal (unpublished results), we have shown that it was possible to compute the economic cost of a wind turbine, but that one should consider this cost in a different way according to the operating regime, i.e. whether the turbine was operated in MPPT mode or in a sub-optimal (controlled) mode.

Variable Horizon Predictive Supervisor

-page 14: detection based on differentiated signals can be problematic in real world applications as it amplifies the noises. R. In general, this is true. In our study, detection was found effective as noises were small enough.

Reviewer #3: Most of my comments have been satisfactory considered and found errors have been removed.However, I do not think that all weak places have been completely improved. Therefore, my recommendation to publish the papers without further amendments is valid only under the assumption that the others reviewers will be satisfied with all that they asked from the authors to improve. My attitude ls also influenced by the fact that the article may not have any extraordinary scientific value, but it represents for practice useful application of theory. All points in the CEP review guidelines are now tolerable fulfilled and this is the most respectful aspect for publishing in Control Engineering Practice journal Reviewer #4: This is a revision of a previously submitted manuscript, with extensive changes by the authors. Attention has been given in particular to the motivating sections which were queried in a previous version as the particular power system application was less well motivated than the control theory itself. The theory and application of the variable horizon control itself appear sound.

*Response to Reviewers

R. We agree with reviewers #3 and #4 and we think that this paper, while not very theoretical, will be useful to control practitioners. Thanks again to reviewer #3 for his very careful reading of the paper.


Nyoman Juniarta
email: nyoman.juniarta@loria.fr

Miguel Couceiro

Amedeo Napoli

Chedy Raïssi

Sequential pattern mining for analyzing visitor trajectories

Hecht Museum is an archaeological museum in Haifa, Israel.

In this work, we present our work on the grouping and mining of trajectory patterns of visitors in this museum.

Within the framework of CrossCult project, we are working on a dataset of the trajectories of 254 visitors in Hecht Museum. Each trajectory contains a list of visited items.

An example is illustrated in the table below. Based on his/her movement, a visitor can be grouped as one of four defined behaviors: ant, grasshopper, butterfly, and fish, as illustrated in the figures below.

Each trajectory is modeled as a sequence of items.

The trajectory in

Conclusion

We apply two algorithms to obtain two lattices:

• MFCS (mining frequent contiguous sequences) [START_REF] Buzmakov | On mining complex sequential data by means of FCA and pattern structures[END_REF] (102, 202, 301, 402) → (102, 202, 301)

• MRGS (mining rare general sequences) [START_REF] Codocedo | A proposition for sequence mining using pattern structures[END_REF] (102, 202, 301, 402) → (102, 202, 402)

Jumping emerging pattern (JEP)

From the two lattices, there are some concepts whose extent contains visitors from the same label.

The intent of such concepts is a JEP.

Our results highlight some interesting patterns that may define visitor behaviors. This can help museum researchers to analyze and evaluate the placement of items and the visiting styles. Furthermore, these patterns can be analyzed to build a recommendation system for future visitors.

Table I

 I

		. An example of one visitor trajectory
	Start	Finish	Item	Room
	14:10:09	14:11:14	102	1
	14:12:45	14:20:19	402	4
	14:22:10	14:25:42	407	4


Ludivine Javourey-Drevet

Stéphane Dufau

Núria Gala

Jacques Ginestié

Johannes C Ziegler

IMPACT OF SIMPLIFICATION -ORIGINAL TEXTS -SIMPLIFIED TEXTS

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Digital as a teaching tool: Can the comprehension and fluency of reading be increased by the simplification and the nature of texts for young students?

Digital as a teaching tool: Can the comprehension and fluency of reading be increased by the simplification and the nature of texts for young students?

Ludivine Javourey-Drevet Written comprehension is fundamental for academic, professional, and personal success. According to the "simple" model of reading proposed by Hoover and Gough (1990), written comprehension is the product of two independent skills: oral comprehension (vocabulary and retrieval of semantic information) and decoding (identification of written words ).

This study aims to investigate whether text simplification improves the fluency and reading comprehension of students aged from 7 to 9 years old. Simplification was carried out at three levels (lexical, syntactic and discursive) for literary and scientific texts.

The first experimental phase allowed to pass cognitive, linguistic and sociocultural tests in individual. For the second phase, a computerized tool was developed on iPad to measure and transmit "online" reading speed and comprehension scores. A cohort of 170 students in 9 classes followed the experiment. The students read independently 20 texts including 10 scientists and 10 literaries, in the original or simplified version. Each text was followed by 5 comprehension questions and 3 suggested answers.

➢ What are the factors (cognitive, linguistic) that determine text comprehension and fluency? ➢ Does reading comprehension of texts vary according to the nature of the text (scientific and literary)?

➢ Can text simplification increase fluency and comprehension, and if so for whom?

Early results suggest that text simplification improves the speed and the comprehension of literary texts but has no effect on the reading of scientific texts. The comprehension of simplified literary texts was superior to the comprehension of scientific texts in original and simplified versions. The same group of students will be followed for 3 years to determine the factors which influence the texts comprehension over time.

The integration of normative data in the computer tool allows the obtention of relatively precise and dynamic measures of the students' level of comprehension. Reading on iBook (autonomy)

Objectives

Introduction Preliminary results

Conclusion

  . A., & Gough, P. B. (1990). The simple view of reading. Readingand writing, 2(2), 127-160.READING FLUENCY ACCORDING TO THE NATURE OF TEXTS AND SIMPLIFICATION NUMBER OF CORRECT RESPONSES ACCORDING TO NATURE OF TEXTS AND SIMPLIFICATION

1,2 , Stéphane Dufau 2 , Núria Gala 3 , Jacques Ginestié 1 , Johannes Ziegler 2

	1 Laboratoire Apprentissage Didactique Evaluation Formation EA 4671
	2 Laboratoire Psychologie Cognitive -UMR 7290
	3 Laboratoire Parole et Langage -UMR 7309


Paavo Salminen
email: phsalmin@abo.fi

Lioudmila Vostrikova

On moments of integral exponential functionals of additive processes

Keywords: independent increments, Lévy process, subordinator, Bessel process, geometric Brownian motion AMS Classification: 60J75, 60J60. 60E10 1

For real-valued additive process (X t ) t≥0 a recursive equation is derived for the entire positive moments of functionals

in case the Laplace exponent of X t exists for positive values of the parameter. From the equation emerges an easy-to-apply sufficient condition for the finiteness of the moments. As an application we study first hit processes of diffusions.

Introduction

Let X = (X t ) t≥0 , X 0 = 0, be a real valued additive process, i.e., a strong Markov process with independent increments having càdlàg sample paths which are continuous in probability (cf. Sato [START_REF] Sato | Lévy Processes and Infinitely Divisible Distributions[END_REF], p.3). Important examples of additive processes are:

(a) Deterministic time transformations of Lévy processes, that is, if (L s ) s≥0 is a Lévy process and s → g(s) is an increasing continuous function such that g(0) = 0 then (L g(s) ) s≥0 is an additive process.

(b) Integrals of deterministic functions with respect to a Lévy process, that is, if (L s ) s≥0 is a Lévy process and s → g(s) is a measurable and locally bounded function then

Z t := t 0 g(s) dL s , t ≥ 0,
is an additive process.

(c) First hit processes of one-dimensional diffusions, that is, if (Y s ) s≥0 is a diffusion taking values in [0, ∞), starting from 0, and drifting to +∞ then H a := inf{t : Y t > a}, a ≥ 0, is an additive process.

Of course, Lévy processes themselves constitute a large and important class of additive processes. The aim of this paper is to study integral exponential functionals of X, i.e., functionals of the form

I s,t := t s exp(-X u )du, 0 ≤ s < t ≤ ∞, (1)
in particular, the moments of I s,t . We refer also to a companion paper [START_REF] Salminen | On exponential functionals of processes with independent increments[END_REF],

where stochastic calculus is used to study the Mellin transforms of I s,t when the underlying additive process is a semimartingale with absolutely continuous characteristics.

The main result of the paper is a recursive equation, see [START_REF] Borodin | Handbook of Brownian motion -Facts and Formulae[END_REF] below, which generalizes the formula for Lévy processes presented in Urbanik [START_REF] Urbanik | Functionals of transient stochastic processes with independent increments[END_REF] and Carmona, Petit and Yor [START_REF] Carmona | On the distribution and asymptotic results for exponential functionals of Levy processes[END_REF], see also Bertoin and Yor [3]. This formula for Lévy processes is also displayed below in (18). In Epifani, Lijoi and Prünster [START_REF] Epifani | Exponential functionals and means of neutral-to-right priors[END_REF] an extension of the Lévy process formula to integral functionals up to t = ∞ of increasing additive processes is discussed, and their formula [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF] can be seen as a special case of our formula (15) -as we found out after finishing our work. We refer to these papers and also to [START_REF] Salminen | On exponential functionals of processes with independent increments[END_REF] for further references and applications, e.g., in financial mathematics and statistics.

In spite of the existing closely related results we feel that it is worthwhile to provide a more thorough discussion of the topic. We also give new (to our best knowledge) applications of the formulas for first hit processes of diffusions and present explicit results for Bessel processes and geometric Brownian motions.

Main results

Let (X t ) t≥0 be an additive process and define for 0

≤ s ≤ t ≤ ∞ and α ≥ 0 m (α) s,t := E I α s,t = E t s e -Xu du α , α ≥ 0, (2) and m (α) t
:= m (α) 0,t , m (α) ∞ := m (α) 0,∞ .
In this section we derive a recursive integral equation for m

(α)
s,t under the following assumption:

(A) E(e -λXt ) < ∞ for all t ≥ 0 and λ ≥ 0.

Under this assumption we define Φ(t; λ) := -log E(e -λXt ).

(

) 3
Since X is assumed to be continuous in probability it follows that t → Φ(t, λ) is continuous. Moreover, X 0 = 0 a.s. implies that Φ(0; λ) = 0. If X is a Lévy process satisfying (A) we write (with a slight abuse of the notation)

E(e -λXt ) = e -tΦ(λ) . (4)
See Sato [START_REF] Sato | Lévy Processes and Infinitely Divisible Distributions[END_REF] Theorem 9.8, p.52, for properties and the structure of the Laplace exponent Φ of the infinitely divisible distribution X t .

In particular, (A) is valid for increasing additive processes. Important examples of these are the first hit processes for diffusions (cf. (c) in Introduction). Assumption (A) holds also for additive processes of type (a) in Introduction if the underlying Lévy process fullfills (A).

Remark 2.1. If X is a semimartingale with absolutely continuous characteristics, a sufficient condition for the existence of the Laplace exponent as in (A) in terms of the jump measure is given in Proposition 1 in [START_REF] Salminen | On exponential functionals of processes with independent increments[END_REF].

By continuity we have -from Jensen's inequality -the following result Lemma 2.2. Under assumption (A) m

(α) s,t < ∞ for all 0 ≤ α ≤ 1 and 0 ≤ s ≤ t < ∞.
The main result of the paper is given in the next theorem. In the proof we are using similar ideas as in [START_REF] Carmona | On the distribution and asymptotic results for exponential functionals of Levy processes[END_REF].

Theorem 2.3. Under assumption (A) it holds for α ≥ 1 and 0 ≤ s ≤ t < ∞ that the moments m (α) s,t are finite and satisfy the recursive equation

m (α) s,t = α t s m (α-1) u,t
e -(Φ(u;α)-Φ(u;α-1)) du.

(

) 5
Proof. We start with by noting that the function s : [0, t] → I s,t is for any t > 0 continuous and strictly decreasing. Hence, for α ≥ 1

α s 0 I α-1 u,t dI u,t = α Is,t I 0,t v α-1 dv = I α s,t -I α 0,t ,
where the integral is a pathwise (a.s.) Riemann-Stiltjes integral and the formula for the change of variables (see, e.g., Apostol [START_REF] Apostol | Mathematical Analysis[END_REF], p. 144) is used. Consequently, from the definition of I s,t it follows that

I α s,t -I α 0,t = -α s 0 I α-1 u,t e -X u-du = -α s 0 I α-1 u,t e -Xu du.
Introducing the shifted functional I s,t via

I s,t := t-s 0 e -(X u+s -Xs) du
we have

I s,t = e Xs I s,t = e Xs t s e -Xu du, (6)
and, therefore,

I α s,t -I α 0,t = -α s 0 I α-1 u,t e -αXu du. (7)
4

Notice that the independence of increments implies that I α-1 u,t and e -αXu are independent, and, hence, for all α ≥ 1

E I α u,t = E I α u,t /E e -αXu (8)
Then evoking Lemma 2.2 and ( 8) yield for 0 ≤ α ≤ 1 and 0

≤ s ≤ t < ∞ E I α s,t ≤ E I α 0,t /E e -αXu < ∞. (9)
Assume now that α ∈ [START_REF] Abramowitz | Mathematical Functions, 9th printing[END_REF][START_REF] Apostol | Mathematical Analysis[END_REF]. Taking the expectations in [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF] and applying Fubini's theorem gives

E I α s,t -I α 0,t = -α s 0 E I α-1 u,t E e -αXu du > -∞ (10)
where the finiteness follows from ( 9). Since I s,t → 0 a.s. when s ↑ t we obtain by applying monotone convergence in ( 10)

E I α 0,t = α t 0 E I α-1 u,t E e -αXu du < ∞. (11)
Putting [START_REF] Salminen | On exponential functionals of processes with independent increments[END_REF] and [START_REF] Salminen | Perpetual integral functionals as hitting and occupation times[END_REF] together results to the equation

E I α s,t = α t s E I α-1 u,t E e -αXu du. (12)
Finally, using ( 8) and ( 9) in [START_REF] Sato | Lévy Processes and Infinitely Divisible Distributions[END_REF] and recalling (3) yield (5) for α ∈ [START_REF] Abramowitz | Mathematical Functions, 9th printing[END_REF][START_REF] Apostol | Mathematical Analysis[END_REF]. Since ( 8) is valid for all α and, as just proved, the finiteness holds for α ∈ [1, 2] the proof of (5) for arbitrary α > 2 is easily accomplished by induction.

Corollary 2.4. Let (X t ) t≥0 be a Lévy process with the Laplace exponent as in (4). Then ( 5) with s = 0 and t < ∞ is equivalent to

m (α) t = α e -tΦ(α) t 0 m (α-1) u e uΦ(α) du. ( 13
)
Proof. Put s = 0 in (5) to obtain m (α) t = α t 0 m (α-1) u,t e -u(Φ(α)-Φ(α-1)) du. (14)
Consider

m (α-1) u,t = E t u e -Xv dv α-1 = E e -(α-1)Xu t u e -(Xv -Xu) dv α-1 = e -uΦ(α-1) E t-u 0 e -(X v+u -Xu) dv α-1 = e -uΦ(α-1) E t-u 0 e -Xv dv α-1 = e -uΦ(α-1) m (α-1)
t-u . Subsituting this expression into ( 14) proves the claim.

For positive integer values on α the recursive equation ( 5) can be solved explicitly to obtain the formula (15) in the next proposition. However, we offer another proof highlighting the symmetry properties present in the expressions of the moments of the exponential functional.

s,t = n! t s dt 1 t t 1 dt 2 • • • (15) • • • t t n-1 dt n exp - n k=1 (Φ(t k ; n -k + 1) -Φ(t k ; n -k)) .
In particular, m

s,∞ < ∞ if and only if the multiple integral on the right hand side of (15) is finite.

Proof. Let t < ∞ and consider m (n) s,t = E t s e -Xu du n = E t s • • • t s e -Xt 1 -•••-Xt n dt 1 . . . dt n = n! E t s dt 1 e -Xt 1 t t 1 dt 2 e -Xt 2 • • • t t n-1 dt n e -Xt n = n! t s dt 1 t t 1 dt 2 • • • t t n-1 dt n E e -(Xt 1 +•••+Xt n ) ,
where, in the third step, we use that (t

1 , t 2 , • • • , t n ) → e -(Xt 1 +•••+Xt n ) is symmetric.
By the independence of the increments E e -αXt = E e -α(Xt-Xs)-αXs = E e -α(Xt-Xs) E e -αXs .

Consequently, E e -α(Xt-Xs) = E e -αXt /E e -αXs = e -(Φ(t;α)-Φ(s;α)) .

Since,

X t 1 + • • • + X tn = n k=1 (n -k + 1) X t k -X t k-1 , t 0 := 0, we have m (n) s,t = n! t s dt 1 t t 1 dt 2 . . . • • • t t n-1 dt n exp - n k=1 (Φ(t k ; n -k + 1) -Φ(t k-1 ; n -k + 1)) .
Applying the initial values given in (4) yields the claimed formula (15). The statement concerning the finiteness of m (n) s,∞ follows by applying the monotone convergence theorem as t → ∞ on both sides of (15). Corollary 2.6. Variable I ∞ has all the positive moments if for all n ∞ 0 e -(Φ(s;n)-Φ(s;n-1)) ds < ∞.

(16)

Proof. From (15) we have

m (n) t ≤ n! n k=1 ∞ 0 e -(Φ(s;n)-Φ(s;n-1)) ds. ( 17
)
The right hand side of (17) is finite if (16) holds. Let t → ∞ in (17). By monotone convergence, m

(n) ∞ = lim t→∞ m (n)
t , and the claim is proved.

Formula (18) below extends the corresponding formula for subordinators found [START_REF] Urbanik | Functionals of transient stochastic processes with independent increments[END_REF], see also [START_REF] Bertoin | Exponential functionals of Lévy processes[END_REF], p.195, for Lévy processes satisfying (A). It is a straightforward implification of Proposition 2.5.

Corollary 2.7. Let (X t ) t≥0 be a Lévy process with the Laplace exponent as in (4) and define n * := min{n ∈ {1, 2, . . .

} : Φ(n) ≤ 0}. Then m (n) ∞ := E(I n ∞ ) =    n! n k=1 Φ(k) , if n < n * , +∞, if n ≥ n * . (18)
.

Example 2.8. A much studied functional is obtained when taking X = (X t ) t≥0 with X t = σW t + µt, σ > 0, µ > 0, where (W t ) t≥0 is a standard Brownian motion. In Dufresne [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF] and Yor [START_REF] Yor | On some exponential functionals of Brownian motion[END_REF] (see also Salminen and Yor [START_REF] Salminen | Perpetual integral functionals as hitting and occupation times[END_REF]) it is proved that

I ∞ := ∞ 0 exp(-(σW s + µs)) ds ∼ H (δ) 0 , (19)
where H δ 0 is the first hitting time of 0 for a Bessel process of dimension δ = 2(1 -(µ/σ 2 )) starting from σ/2, and ∼ means "is identical in law with". In particular, it holds

∞ 0 exp(-(2W s + µs)) ds ∼ 1 2 Z µ , (20)
where Z µ is a gamma-distributed random variable with rate 1 and shape µ/2. We refer to [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF] for a discussion showing how the functional on the left hand side of (19) arises as the present value of a perpetuity in a discrete model after a limiting procedure. Since the Lévy exponent in this case is

Φ(λ) = λµ - 1 2 λ 2 σ 2 ,
the criterium in Corollary 2.7 yields

E (I n ∞ ) < ∞ ⇔ n < 2µ/σ 2
, which readily can also be checked from (20).

First hit processes of one-dimensional diffusions

We recall first some facts concerning the first hitting times of one-dimensional (or linear) diffusions. Let now Y = (Y s ) s≥0 be a linear diffusion taking values in an interval I. To fix ideas assume that I equals R or (0, ∞) or [0, ∞) and that lim sup

s→∞ Y s = +∞ a.s. ( 21
)
Assume Y 0 = v and consider for a ≥ v the first hitting time

H a := inf{s : Y s > a}.
Defining X t := H t+v , t ≥ 0, it is easily seen -since Y is a strong Markov process -that X = (X t ) t≥0 is an increasing purely discontinuous additive process starting from 0. Moreover, from (21) it follows that X t < ∞ a.s. for all t. The process X satisfies (A) and it holds

E v (e -βXt ) = E v (e -βH t+v ) = ψ β (v) ψ β (t + v) , t ≥ 0, ( 22
)
where β ≥ 0, E v is the expectation associated with Y, Y 0 = v, and ψ β is a unique (up to a multiple) positive and increasing solution (satisfying appropriate boundary conditions) of the ODE

(Gf )(x) = βf (x), (23)
where G denotes the differential operator associated with Y. For details about diffusions (and further references), see Itô and McKean [START_REF] Itô | Diffusion Processes and Their Sample Paths[END_REF], and [START_REF] Borodin | Handbook of Brownian motion -Facts and Formulae[END_REF].

The Laplace transform of X t can also be represented as follows

E v (e -βXt ) = exp - t+v v S(du) ∞ 0 (1 -e -βx )n(u, dx) , ( 24
)
where S is the scale function, and n is a kernel such that for all v ∈ I and t ≥ 0

t+v v ∞ 0 (1 ∧ x)n(u, dx)S(du) < ∞.
Representation (24) clearly reveals the structure of X as a process with independent increments. From ( 22) and (24) we may conclude that

∞ 0 (1 -e -β x )n(u, dx) = lim w→u- 1 -E w (e -βXu ) S(u) -S(w) . (25)
We now pass to present examples of exponential functionals of first hit processes. Firstly, we study Bessel processes satisfying (21) and show, in particular, that the exponential functional of the first hit process has all the moments. In our second example it is seen that the exponential functional of the first hit process of geometric Brownian motion has only finitely many moments depending on the values of the parameters.

Example 3.1. Bessel processes. Let Y be a Bessel process starting from v > 0. The differential operator associated with Y is given by

(Gf )(x) = 1 2 f ′′ (x) + δ -1 2x f ′ (x), x > 0,
where δ ∈ R is called the dimension parameter. From [START_REF] Borodin | Handbook of Brownian motion -Facts and Formulae[END_REF] we extract the following information

• for δ ≥ 2 the boundary point 0 is entrance-not-exit and (21) holds,

• for 0 < δ < 2 the boundary point 0 is non-singular and (21) holds when the boundary condition at 0 is reflection,

• for δ ≤ 0 (21) does not hold.

In case when ( 21) is valid the Laplace exponent for the first hit process X = (X t ) t≥0 is given for v > 0 and t ≥ 0 by

E v (e -βXt ) = ψ β (v) ψ β (t) = v 1-δ/2 I δ/2-1 (v √ 2β) t 1-δ/2 I δ/2-1 ((t + v) √ 2β) , (26)
where E v is the expectation associated with Y when started from v and I denotes the modified Bessel function of the first kind. For simplicity, we wish to study the exponential functional of X when v = 0. To find the Laplace exponent when v = 0 we let v → 0 in (26). For this, recall that for p = -1, -2, . . .

I p (v) ≃ 1 Γ(p + 1) v 2 p as v → 0. (27)
Consequently,

E 0 (e -βXt ) = lim v→0 E v (e -βXt ) = 1 Γ(ν + 1) √ 2β 2 δ/2-1 t δ/2-1 I δ/2-1 (t √ 2β)
=: e -Φ(t;β) .

The validity of (16), that is, the finiteness of the positive moments, can now be checked by exploiting the asymptotic behaviour of I p saying that for all p ∈ R (see Abramowitz and Stegun [START_REF] Abramowitz | Mathematical Functions, 9th printing[END_REF], 9.7.1 p.377)

I p (t) ≃ e t / √ 2πt as t → ∞. (28)
Indeed, for n = 1, 2, . . .

e -(Φ(t;n)-Φ(t;n-1)) = n δ/2-1 I δ/2-1 (t √ 2n) I δ/2-1 (t 2(n -1)) (n -1) δ/2-1 ≃ n n -1 δ/2-1 n n -1 1/4 e -t( √ 2n- √ 2(n-
• lim s→∞ Y s = +∞ a.s if µ > 1 2 σ 2 ,
• lim s→∞ Y s = 0 a.s if µ < 1 2 σ 2 .

• lim sup s→∞ Y s = +∞ and lim inf s→∞ Y s = 0 a.s. if µ = 1 2 σ 2 .

Consequently, condition (21) is valid if and only if µ ≥ 1 2 σ 2 . Since Y 0 = 1 we consider the first hitting times of points a ≥ 1. Consider

H a := inf{s : Y s = a} = inf s : exp σW s + (µ - 1 2 σ 2 )s = a = inf s : W s + µ -1 2 σ 2 σ s = 1 σ log a .
We assume now that σ > 0 and µ ≥ 1 2 σ 2 . Let ν := (µ -1 2 σ 2 )/σ. Then H a is identical in law with the first hitting time of log a/σ for Brownian motion with drift ν ≥ 0 starting from 0. Consequently, letting X t := H 1+t we have

for t ≥ 0 E 1 (e -βXt ) = ψ β (0) ψ β (log(1 + t)/σ) = exp - 2β + ν 2 -ν log(1 + t) σ = (1 + t) - √ 2β+ν 2 -ν /σ , (29)
=: exp (-Φ(t; β)) .

where E 1 is the expectation associated with Y when started from 1 and

ψ β (x) = exp 2β + ν 2 -ν x
is the increasing fundamental solution for Brownian motion with drift (see [START_REF] Borodin | Handbook of Brownian motion -Facts and Formulae[END_REF] p.132). Notice that the additive process X is a deterministic time change of the first hit process of Brownian motion with drift, which is a subordinator. We use now By formula (15) the nth moment is given by

E 1 (I n ∞ ) = n! ∞ 0 dt 1 (1 + t 1 ) -(ρ(n)-ρ(n-1)) ∞ t 1 dt 2 (1 + t 2 ) -(ρ(n-1)-ρ(n-2)) × ∞ t 2 dt 3 • • • ∞ t n-1 dt n (1 + t n ) -ρ(1) =          n! n k=1 (ρ(k) -k) , if n < n * , +∞, if n ≥ n * ,
ρ(n) > n ⇔ n < 2 σ 2 , (32)
i.e., smaller the volatility (i.e. σ) more moments of I ∞ exist, as expected. Moreover, in this case

ρ(n) -ρ(n -1) > 1 ⇔ √ 2n + 2(n -1) < 2 σ ⇔ 2n -1 + 4n(n -1) < 2 σ 2 (33)
showing, in particular, that when σ is "small" there exist "many" n satisfying (32) but not (33).

Proposition 2 . 5 .

 25 For 0 ≤ s ≤ t ≤ ∞ and n = 1, 2, . . . it holds m (n)

0 e

 0 Proposition 2.5 to study the moments of the perpetual integral functional I ∞ = ∞ -Xs ds. To simplify the notation (cf. (29)) introduce ρ(β) := 2β + ν 2 -ν /σ.

  where n * := min{n ∈ {1, 2, . . .} : ρ(n) -n ≤ 0}.Condition (16) in Corollary 2.6 takes in this case the formρ(n) -ρ(n -1) > 1. (30)This being a sufficient condition for the finiteness of m(n) ∞ we have ρ(n) -ρ(n -1) > 1 ⇒ ρ(n) -n > 0. (31)Consider now the case ν = 0. Then

Funding. This research was partially supported by Defimath project of the Research Federation of "Mathématiques des Pays de la Loire", by PANORisk project "Pays de la Loire" region, and by the Magnus Ehrnrooth Foundation, Finland.