Xingyu Yan

Dhaker Abbes

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

Uncertainty Analysis for Day Ahead Power Reserve Quantification in an Urban Microgrid Including PV Generators

Keywords: Power reserve scheduling, renewable energy sources, forecast errors, uncertainty analysis, reliability. I

Setting an adequate operating power reserve (PR) to compensate unpredictable imbalances between generation and consumption is essential for power system security. Operating power reserve should be carefully sized but also ideally minimized and dispatched to reduce operation costs with a satisfying security level. Although several energy generation and load forecasting tools have been developed, decision-making methods are required to estimate the operating power reserve amount within its dispatch over generators during small time windows and with adaptive capabilities to markets, as new ancillary service markets. This paper proposes an uncertainty analysis method for power reserve quantification in an urban microgrid with a high penetration ratio of PV (photovoltaic) power. First, forecasting errors of PV production and load demand are estimated one day ahead by using artificial neural networks. Then two methods are proposed to calculate one day ahead the net demand error. The first perform a direct forecast of the error, the second one calculates it from the available PV power and load demand forecast errors. This remaining net error is analyzed with dedicated statistical and stochastic procedures. Hence, according to an accepted risk level, a method is proposed to calculate the required PR for each hour.

to maintain the security and reliability of grids with a high share of renewable generators, primary, secondary and tertiary regulation as well as spinning reserve are now required from renewable generators in more and more grid codes [START_REF] Ipakchi | Grid of the future[END_REF][START_REF] Galiana | Scheduling and pricing of coupled energy and primary, secondary, and tertiary reserves[END_REF]. This operating power reserve should be ideally minimized to reduce system costs with a satisfying security level.

Typically, PV power generation forecasting is needed to optimize the operation and to reduce the cost of power systems, especially for the scheduling and dispatching of required hourly operating [START_REF] Mills | Integrating Solar PV in Utility System Operations[END_REF]. However, the predicted uncertainty associated with forecast cannot be eliminated even with the best model tools. In addition to the load demand uncertainty, the combination of power generation and consumption variability with forecast uncertainty makes the situation more difficult for power system operators to schedule and to set power reserve level. Therefore, the uncertainties from both generation and consumption must be taken into account by an accurate stochastic model for power system management. In addition, forecasting errors from system uncertainty analysis could be used to set power reserve [START_REF] Sobu | Dynamic optimal schedule management method for microgrid system considering forecast errors of renewable power generations[END_REF].

Historically, most conventional utilities have adopted deterministic criteria for the reserve requirement: the operating rules required PR to be greater than the capacity of the largest on-line generator or a fraction of the load, or equal to some function of both of them. Those deterministic criteria are wildly used because of their simplicity and understandable employing. However, I these deterministic calculation methods are gradually replaced by probabilistic methods that respond to the stochastic factors corresponding to the system reliability. Several research works have been focused on calculating the total system uncertainty from all the variable sources. Based on dynamic simulations, the study in [START_REF] Delille | Dynamic frequency control support by energy storage to reduce the impact of wind and solar generation on isolated power system's inertia[END_REF] is focusing on dynamic frequency control of an isolated system and the reduction of the impact due to large shares of wind and PV powers. However this work did not consider other aspects, such as variability and forecast accuracy. A deterministic approach is proposed in [START_REF] Vosa | Revision of reserve requirements following wind power integration in island power systems[END_REF] to analyze the flexibility of thermal generation to balance wind power variations and prediction errors. A stochastic analysis could improve it in order to be able to quantify the power reserve with a risk index. A stochastic model was developed in [START_REF] Sansavini | A stochastic framework for uncertainty analysis in electric power transmission systems with wind generation[END_REF] to simulate the operations and the line disconnection events of the transmission network due to overloads beyond the rated capacity. But the issue is clearly that analysis of the system states, in terms of power request and supply, are critical for network vulnerability and may induce a cascade of line disconnections leading to a massive network blackout. An insurance strategy is proposed in [START_REF] Yang | Insurance strategy for mitigating power system operational risk introduced by wind power forecasting uncertainty[END_REF] to cover the possible imbalance cost that wind power producer may incur in electricity markets. Monte Carlo simulations have been used to estimate insurance premiums for further analysis excesses and so require a significant calculation time.

Our previous works in [START_REF] Buzau | Quantification of operating power reserve through uncertainty analysis of a microgrid operating with wind generation[END_REF][START_REF] Yan | Operating power reserve quantification through PV generation uncertainty analysis of a microgrid[END_REF] showed that forecasting errors from system uncertainty analysis could be used for PR setting.

Following these promising results and experiences, we have carried out further investigations on rigorous methods to quantify the required PR. The task is to calculate it by considering uncertainties from PV prediction and load forecast or with uncertainty estimation. In the second part of this paper, PV power and load uncertainty and variability are analyzed. Then, the artificial neural network based prediction methods are applied to forecast PV power, load demand and errors. In the third part, the Net Demand (ND) Forecasted uncertainty is obtained, for each hour of the next day, as the difference between the forecasted production uncertainty and the forecasted load uncertainty. Two methods are detailed to calculate the ND forecast errors. An hourly probability density function of all predicted ND forecasted errors has been used for the error analysis. In the fourth part, a method is explained to assess the accuracy of these predictions and to quantify the required operating PR to compensate the system power unbalancing due to these errors. Power reserve is obtained by choosing a risk level related to two reliability assessment indicators: loss of load probability (LOLP), and expected energy not served (EENS). Finally, this management tool is proved through an illustrative example.

II.METHODOLOGY

A.PV Power and Load Uncertainty Analysis

The PV power variability is the expected change in generation while PV power uncertainty is the unexpected change from what was anticipated, such as a suddenly cloud cover. The former depends on the latitude and the rotation of Earth, while the latter is mostly caused by uncertainty conditions, such as cloud variations over the PV. The movement of clouds introduces a significant uncertainty that can result in rapid fluctuations in solar irradiance and therefore PV power output. However, the influence of a moving cloud and, hence, the shading of an entire PV site depends on the PV area, cloud speed, cloud height and many other factors. Data from solar installations covering a large spatial extent have an hourly temporal dynamic, while individual zones have instantaneous dynamics as in local distribution networks or micro grids.

The daily operation of a power system should be matched to load variations to maintain system reliability. This reliability refers to two areas:

-system adequacy, which depends on sufficient facilities within the system to satisfy system operational constraints and load demand, and -system security, which is the system ability to respond to dynamic disturbances.

When RES represent a significant part of the power generation, system operating power reserve must be larger to regulate the variations and maintain the security level. This additional power is required to stable the electrical network. Classically this power reserve is provided by controllable generators (gas turbines, diesel plant, ..). Today, the increasing of balancing power reserves leads to a significant increase in the power system operating cost and so system may limit the PV power penetration is due to the variability and uncertainty over short time scales.

There are different ways to manage variability and uncertainty. In general, system operators and planners use mechanisms including forecasting, scheduling, economic dispatch, and power reserves to ensure performances that satisfy reliability standards with the least cost. The earlier system operators and planners know what kind of variability and uncertainty they will have to deal with, the more options they will have to accommodate it and cheaper it will be. The key task of variability and uncertainty management is to maintain a reliable operation of the power system (grid connected or isolated) while keeping down costs.

Energy management of electrical systems is usually implemented over different time scales. One day ahead, system operators have to balance the load demand with electrical generation by planning the starting and set points of controllable generators on an hourly time step. Risks are also considered and thus a power reserve has also to be hourly planned. During the day, unexpected PV power lack is compensated by injecting a primary power reserve. The PV variability can be separated into different time scales associated with different impacts, onto the grid management and costs. Consequently, more capacity to compensate errors in forecasts or unexpected events must be accommodated.

The instantaneous PV power output is affected by many correlated external and physical inputs, such as irradiance, humidity, pressure, cloud cover percentage, air/panels temperature and wind speed. The per unit surface power output is modeled by [START_REF] Fuentes | Application and validation of algebraic methods to predict the behaviour of crystalline silicon PV modules in Mediterranean climates[END_REF]:

) ) 25 ) ( ( . - 1 ( ). ( . . = ) (  t T C t I A t P p r PV  (1)
where is the power conversion efficiency of the module (%), A is the surface area of PV panels (m 2 ), I r is the global solar radiation (kW/m 2 ) and T is the outside air temperature (°C), p C is the cell maximum power temperature coefficient (equal to 0.0035 but it can varies from 0.005 to 0.003 per °C in crystalline silicon).

The PV power, solar irradiance and temperature of our lab PV plant have been recorded during three continuous days (22/06/2010 -24/06/2010) and are presented respectively in Fig. 1. The PV power variability is highly correlated with irradiance, so as to the temperature, while the PV power uncertainty is almost caused by the irradiance change. Sensed PV power data points can be drawn according to sensed irradiance and temperature data points in order to highlight correlations (Fig. 2). The local load consumption demand is also highly unpredictable and quite random. It depends on different factor, such as the economy, the time, the weather, and other random effects. However, for power system planning and operation, load demand variation and uncertainty analysis are crucial for power flow study or contingency analysis. As for PV production, load demand variations exist in all time scales and system actions are needed for power control in order to maintain the balancing.

B.Power Forecasting Methodology 1)PV Power Forecasting

In recent decades, several forecasting models of energy production have been published [START_REF] Espinar | Photovoltaic Forecasting: A state of the art[END_REF][START_REF] Yan | Solar radiation forecasting using artificial neural network for local power reserve[END_REF][START_REF] De Rocha | Photovoltaic forecasting with artificial neural networks[END_REF][START_REF] Mohammed | Probabilistic Forecasting of Solar Power: An Ensemble Learning Approach[END_REF][START_REF] Golestaneh | Generation and evaluation of space-time trajectories of photovoltaic power[END_REF]. For PV power, one method consists in forecasting solar radiation and then forecasting PV power with a mathematical model of the PV generator. A second one proposes to directly predict the PV power output from environmental data (irradiance, temperature, etc.). Statistical analysis tools are generally used, such as linear/multiple-linear/non-linear regression and autoregressive models that are based on time series regression analysis [START_REF] De Giorgi | Photovoltaic power forecasting using statistical methods: impact of weather data[END_REF]. These forecasting models rely on modeling relationships between influent inputs and the produced output power. Consequently, mathematical model calibration and parameters adjustment process take a long time.

Meanwhile, some intelligent based electrical power generation forecast methods, as expert systems, fuzzy logic, neural networks, are widely used to deal with uncertainties of RES power generation and load demand [START_REF] Espinar | Photovoltaic Forecasting: A state of the art[END_REF][START_REF] Ghods | Different methods of long-term electric load demand forecasting; a comprehensive review[END_REF].

In daily markets, the hourly PV power output for the next day (day D+1) at time step h is represented as the sum of a day ahead hourly forecast PV power ( h v P ~) and the forecast error (ε h Pv ):

Pv h h h v P Pv   ~ = (2)

2)Load Demand Forecasting

For load demand forecast, numerous variables affect directly or indirectly the accuracy. Until now, many methods and models have already been tried out. In [START_REF] Ghods | Different methods of long-term electric load demand forecasting; a comprehensive review[END_REF], several long-term (month or year) load forecasting methods are introduced and are very important for planning and developing future generation, transmission and distribution systems. In [START_REF] Hong | Long term probabilistic load forecasting and normalization with hourly information[END_REF], a long term probabilistic load forecasting method is proposed with three modernized elements: predictive modeling, scenario analysis, and weather normalization. Long-term and short-term load forecast play important roles in the formulation of secure and reliable operating strategies for the electrical power system. The objective is to improve the forecast accuracy in order to optimize power system planning and to reduce costs. ):

The day ahead actual load demand at time step h (

L h h h L L   ~ = (3)
3)Net Demand Forecasting Knowing the PV power forecasting and the load demand forecasting, the net demand forecasting ( h D N ~) for a given time step h is expressed as:

h h h v P L D N   (4)
The real net demand (ND h ) is composed of the forecasted day ahead ND and a forecast error (ε h ND ):

ND h h h D N ND   ~ = (5)

C.Application of Back-Propagation ANN to Forecast

In order to predict the net demand errors, as well as PV and load forecast errors, we have developed several back-propagation (BP) Artificial Neural Networks (ANN) [START_REF] Francois | Orthogonal considerations in the design of neural networks for function approximation[END_REF]. Compared with conventional statistical forecasting schemes, ANN has some additional advantages, such as simplicity in adaptability to online measurements, data error tolerance and lack of any excess information. Since the fundamentals of ANN based predictors can be found in many sources, it will not be recalled again.

III.NET DEMAND UNCERTAINTY ANALYSIS

A.Net Demand Uncertainty In order to simplify the study, the uncertainties coming from conventional generators and network outages are ignored and only load and PV power uncertainties are considered. Then the error of the ND forecasting is representing the ND uncertainty.

Two possible methods are proposed to calculate the forecasted net demand error.

1)First Method: Forecast of the Day-ahead Net Demand Error

The real ND is the difference of the sensed load and the sensed PV power. Based on the historical sensed and forecasted database of the load demand and PV production, the past forecasted ND is calculated as the difference between past forecasted load and past forecasted PV power at time step h. Then by using past real ND ( ). Hence these data are used to calculate the day-ahead forecast of ND errors ND h  ~ (D+1) (Fig. 3). The obtained ND error forecast can be characterized by the mean and variance (respectively, ND h  and ND h  ).

2)Second Method: Calculation from the PV Power and the Load Forecast Errors Estimation

A second method is to define the ND uncertainty as the combination of PV power and load uncertainties. It is generally assumed that PV power and load forecast errors are unrelated random variables. So, firstly the day-ahead PV power and load forecasting errors ( PV h  ~and L h  ~) are estimated independently. Then, last 24-hour load forecast errors and PV power forecast errors are calculated as the difference of the sensed load and PV power, and forecasted load and PV power, respectively (Fig. 4).

The mean values and standard deviations of those forecasting errors can be obtained. Then the ND forecasting error can be attained as a new variable which comes from those two independent variables. The new obtained pdf is also a normal distribution with the following mean and variance [START_REF] Ortega-Vazquez | Estimating the spinning reserve requirements in systems with significant wind power generation penetration[END_REF][START_REF] Bouffard | Stochastic security for operations planning with significant wind power generation[END_REF]:

PV h L h ND h      (6) 2 2 ) ( ) ( PV h L h ND h      (
) ( , ND h ND h   2 ) ( , PV h PV h   2 ) ( , L h L h  

Mean average and Standard deviation

Equations ( 6) and ( 7

B.Assessment of the Forecasting Uncertainty

The predicted errors ( ND h  ~) of the ND forecast ( ND h D N ~) can be obtained with the normal probability density function (Fig. 5).

) , ( 2

1 2 2 ) ( 2 ) ( ND h ND h B ND h ND h B F d e pdf ND h ND h                (8)
The forecasting uncertainty can be represented as upper and lower bound margins around the ND forecast. Bound margins (B) are extracted by a normal inverse cumulative distribution function for a desired probability index x (Fig. 5):

  x B F B x F B ND h ND h ND h ND h    ) , ( : ) , ( 1 -     (9)
Lower bound (Blower)

Normal Inverse

Cumulative Distribution Function (F -1 ) Normal Probability Density Function (F(B|µ ND ,σ ND ))

Upper bound (Bup)

+ + _ + ND h   ND h   Net Demand Forecast

IV.POWER RESERVE QUANTIFICATION

A.Reliability Assessment

Resulting from the uncertainty assessment, the pdf of the forecasted ND errors in a given time step is considered for the calculation of the power reserve [START_REF] Holttinen | Methodologies to determine operating reserves due to increased wind power[END_REF]. To estimate the impact of forecast ND uncertainty, two common reliability assessment parameters are used: the loss of load probability (LOLP) and the expected energy not served (EENS) [START_REF] Li | A multi-state model for the reliability assessment of a distributed generation system via universal generating function[END_REF][START_REF] Wang | Spinning reserve estimation in microgrids[END_REF][START_REF] Liu | Quantifying spinning reserve in systems with significant wind power penetration[END_REF]. LOLP represents the probability that the load demand (L h ) exceeds PV power (P h ) at time step h:

     R h h h pdf P L prob LOLP   )d ( ) 0 ( = ( 10
)
prob(L h -P h >0) is also the probability that the power reserve (R) is insufficient to satisfy the load demand in the time step h.

Meanwhile, EENS measures the magnitude of the load demand not served:

) ( ) 0 ( = h h h h h P L P L prob EENS     (11)
where (L h -P h ) is the missed power in the time step h.

In this situation, the grid operator can either disconnect a part of loads or use the power reserve to increase the power production.

After obtaining each of the next 24 hours forecast ND pdf s , an hourly day ahead reliability assessment can be attained. Electrical system operators can use this reliability to calculate the system security level.

B.Risk-constrained Energy Management

A reserve characteristic according to a risk level and for each time step can be obtained. With a fixed risk index, the operator can then easily quantify the power reserve [START_REF] Vosa | Revision of reserve requirements following wind power integration in island power systems[END_REF]. As shown in Fig. 6, the sum of the shaded areas represents the accepted risk of violation with x % of LOLP. R is the needed power reserve to compensate the remaining power unbalance. So, the reliability assessment can be done with the hourly cumulative distribution function (cdf) obtained from the normal difference distribution of ND errors. Then the cdf represents the probability that the random variable (here the ND error) is less or equal to x.

This assessment has been made under the assumption of a positive hourly forecasted ND. Otherwise, if the forecasted ND is negative, the reserve power for the same reliability level will be unnecessary (the power generation is more than the load demand).

So the reliability has been assessed by considering only positive forecasted ND errors (L h -P h ) for each time step. Then, LOLP is deduced with:

      R h h h pdf P L prob LOLP   )d ( - 1 0) ( - 1 = (12)
When the LOLP equals to the risk index x %, the reserve power (R) covers the remaining probability that the load demand exceeds the PV power generation (blue part in Fig. 6).

V.ILLUSTRATIVE CASE STUDY

A.Presentation and Data Collection

The studied urban microgrid is a 110 kW load peak and is powered with 17 kW PV panels and three micro-gas turbines 30 kW, 30 kW and 60 kW each. Sensed data from our 17 kW PV plant located on the lab roof have been recorded in 2010 and 2013. For the load forecasting, past daily French power consumptions have been scaled to obtain per unit values of locally power consumption with the same characters and dynamics. A part of this database has been used to design the ANN based forecasting tool, a part to assess the estimation quality and a third one to implement the application of the proposed method in a real situation [START_REF] De Rocha | Photovoltaic forecasting with artificial neural networks[END_REF].

The ANN has been trained with past recorded data from the training set to predict hourly PV output power.

   n i k k y y n nRMSE 1 2 ) ( 1 = (13)    n i k k y y n nMAE 1 1 = (14)

B.ANN Based Power Forecast and Net Demand Forecast 1)ANN based PV Power Forecasting

A three-layer ANN has been developed for the PV power generation prediction with:

one input layer including last n hours of measured PV power, of irradiance and of forecasted average temperature (obtained from our local weather information service) (Fig. 7);

one hidden layer with 170 neurons;

-one output layer with the 24 predicted PV power points (for each hour).

Various hidden layer neurons have been tested until getting an nRMSE inferior to 5%. First, 60% of previously sensed data (representing one year of data) have been used for training the ANN based PV power forecasting tool. Next 20% of sensed data are used to create a validation pattern set in order to assess the prediction quality. The test set (with the remaining 20% data) is used to implement the forecast error calculation. Obtained nRMSE and nMAE for next 24 hours PV power predictions are given in Table I. Predicted errors for 120 test days are given in Fig. 8. Absolute values are less than 0.4 p.u. of the PV power output.

The largest errors are in the middle of the day when the PV power production is the highest.

2)ANN based Load Forecasting

Another neural network has been used for load forecast. The load demand prediction model includes: an input layer with last 48 hours load demand measurements and predicted temperatures for next 24 hours, one hidden layer with 70 neurons (in order to get an nRMSE inferior to 4%) and an output layer that predicts next 24 hours load demand. 60% of available data are used for the neural network training, 20% for the validation and 20% for tests. The predicted errors for 120 test days are shown in Fig. 9. As it can be seen, the largest forecast error occurs at 8:00 and 18:00. Yet the total absolute errors are less than 0.2 p.u. of the load demand. Obtained results of nRMSE and nMAE are listed in Table II. Following the method highlighted in Fig. 3, another ANN is applied for ND errors forecast: an input layer with last 24 hours predicted net demand errors, one hidden layer with 70 neurons and an output layer that predicts next 24 hours forecasted net errors.

Application of the first method (Fig. 3) for the time step at 12 am gives: 1.781 0.1282, -

C.Forecasting Uncertainty Assessment

By applying both proposed methods, the uncertainties of PV power forecasting, load forecasting and net forecasting with various probability indices (from 90% to 60%) in a random day are represented as a function of the forecasting data and the predicted errors of forecasting. In order to simplify the explanation, results are given with the second method (corresponding to Fig. 5).

As shown in Fig. 11, the uncertainty of PV power forecasting is higher in the middle of the day, when the PV system generates the highest power. While in the morning (from 6:00 to 10:00) and afternoon (from 17:00 to 21:00), the uncertainty is smaller.

Obviously, PV power forecasting uncertainty increases, and decreases with PV power increase and decrease respectively. Also, the uncertainty is increased when the time horizon is larger. For example, at 10:00 and at 17:00 power outputs are almost at the same level (about 6.5 kW), but uncertainty is larger at 17:00 then at 9:00. The load forecasting has the same variation trend (Fig. 12). Fig. 13 depicts the obtained ND uncertainty with the first method. If the forecasted ND is positive, then additional power sources have to be programmed to cover the difference. Otherwise, if forecasted ND is negative then three actions must be considered to meet the low forecasted demand:

-A part of PV power generators must be switched off (or can work at a sub-optimal level).

-Controllable loads (as electrical vehicles, heating loads.) must be switched on to absorb excess available power.

-Export the available excess energy to the main grid.

D.Power Reserve Calculation with Fixed Risk Indices

The forecasted ND uncertainty assessment has been done with the hourly cumulative distribution function (cdf) obtained from the ND forecast errors. Then, the hourly risk/reserve curve takes into account all the errors from the cdf s . Since the forecast ND errors can be expressed as an x % of the rated power, the PR can be drawn according to the LOLP. Fig. shows the required PR variation according to LOLP and EENS (with the second method in Section III). Therefore, an operating PR under x % of LOLP would cover a part of the forecast ND uncertainty. For example, with 10% of LOLP, the reserve power will be 7 kW and the EENS will be 0.2 kW. In general, this operating reserve is limited not only by the risk indices but also by the availability of micro- On Fig. 15, an assessment of hourly reserve power required with the second method for different LOLP has been deduced.

Much more reserve will be needed when the LOLP rate is very low, which means a high security level. While less reserve power will be needed with a high LOLP rate, but then the risk will be higher. For example with 1% of LOLP, the necessary PR will be 14 kW (EENS is almost zero) at 12:00 am, while the necessary reserve power will be 7 kW with a 10% of LOLP and EENS increases to 0.25 kWh.

If a constant LOLP rate is set, the power reserve for each hour can be obtained. As shown in Fig. 16, with a 1% of LOLP, more power reserve is needed in the middle of the day when larger PV power is generated. Moreover, power reserve with the second method is higher than the method with direct ND forecast. The most likely explanation of this result is because the load forecast uncertainty and PV forecast uncertainty are not totally independent. Sharing a common temperature, integrated PV power uncertainty and load uncertainty is greater than the direct ND forecast uncertainty. This result can be used for power dispatch management.

VI.CONCLUSION

This work proposed a new technique to quantify the power reserve of a microgrid by taking into account the PV power forecasting uncertainty and load forecasting uncertainty. In order to assess these uncertainties, a three-layer BP ANN is used to estimate errors of PV power and load forecastings. Two methods are proposed to obtain the ND forecast uncertainties. With the first method, a probabilistic model is proposed to forecast ND uncertainty distribution by integrating the uncertainties from both PV power and load. The other method is desired to directly forecast the ND errors. The power reserve quantification results demonstrate that with a fixed risk index, the power reserve for next day or next 24 hours can be evaluated to cover the risk.

As the uncertainty from forecasting errors increases with time horizon, future research works are oriented toward the implementation of the intraday adjustment. The dispatch of the calculated power reserve onto micro-gas turbines, controllable loads and also new "PV based active power generators" is also an interesting way to pave.

Fig. 1 .Fig. 2 .

 12 Fig. 1. PV power, solar irradiance and temperature in three continuous days.

  h L ) is assumed to be the sum of the day ahead forecasted load ( h L ~) and an error (ε h L

Fig. 3 .

 3 Fig. 3. Net demand uncertainty calculation from ND error forecast.

Fig. 4 .

 4 Fig. 4. Net demand uncertainty calculation from PV power and load forecasting errors prediction.

Fig. 5

 5 Fig. 5 Net uncertainty calculation at hour h with a given probability.

Fig. 6 .

 6 Fig. 6. Calculation of power reserve requirements (R) based on forecast ND uncertainty ( ND h  ~) with x% of LOLP, at time step h.

24 .Fig. 7 .

 247 Fig. 7. PV power, load forecasting and errors prediction with ANN.

Fig. 8 .

 8 Fig. 8. PV power prediction errors on 120 test days.

Fig. 9 .

 9 Fig. 9. Load prediction on 120 test days.

Fig. 10 .

 10 Fig. 10(c) with obtained parameters:

Fig. 11 .

 11 Fig. 11. PV forecasting with uncertainty (a random day).

Fig. 12 .

 12 Fig. 12. Load forecasting with uncertainty (a random day).

Fig. 13 .

 13 Fig. 13. Next 24 hours NFD with uncertainty (a random day).

Fig. 14 .

 14 Fig. 14. Risk/reserve curve for LOLPh+12 and EENSh+12 at 12:00.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Required power reserve for each hour with x% LOLP. 0 1 2 3 4 5 6 7 8 9 10 0 4 8 12 16 20 24 0 5 10 15 x % of LOLP Time (Hour) Power Reserve (kW)

TABLE I .

 I Errors of the PV Power Forecast with ANN

		nRMSE [%] nMAE [%]
	Training Set	4.67	2.69
	Validation Set	5.58	3.13
	Test Set	5.95	3.12

TABLE II

 II

								nRMSE [%]	nMAE [%]
				Training Set		3.18	2. 45
			Validation Set		3.57	2.76
					Test Set		3.67	2.84
	3)Net Demand Uncertainty
	a) First Method: Direct Net Demand Forecast
			0.4
		Error (p.u.)	0 0.2
			-0.2
			-0.4	0	4	8	12	16	20	24
								Time (hours)
		0.3
		0.2
	Error (p.u.)	-0.1 0 0.1
		-0.2
		-0.3	0	4		8	12	16	20	24
								Time (hours)
					. ERRORS OF THE LOAD DEMAND FORECAST WITH ANN

12

ACKNOWLEDGMENT

The authors would like to thank the China Scholarship Council and Centrale Lille for their cofounding supports


Vineeth S Varma

Irinel-Constantin Morȃrescu

Yezekael Hayel

Continuous time opinion dynamics of agents with multi-leveled opinions and binary actions

Keywords: Opinion dynamics, Social computing and networks, Markov chains, agent based models

This paper proposes and analyzes a stochastic multiagent opinion dynamics model. We are interested in a multileveled opinion of each agent which is randomly influenced by the binary actions of its neighbors. It is shown that, as far as the number of agents in the network is finite, the model asymptotically produces consensus. The consensus value corresponds to one of the absorbing states of the associated Markov system. However, when the number of agents is large, we emphasize that partial agreements are reached and these transient states are metastable, i.e., the expected persistence duration is arbitrarily large. These states are characterized using an N-intertwined mean field approximation (NIMFA) for the Markov system. Numerical simulations validate the proposed analysis.

I. INTRODUCTION

Understanding opinion dynamics is a challenging problem that has received an increasing amount of attention over the past few decades. The main motivation of these studies is to provide reliable tools to fight against different addictions as well as propagation of undesired unsocial behaviors/beliefs. One of the major difficulties related to opinion dynamics is the development of models that can capture many features of a real social network [START_REF] Kerckhove | Modelling influence and opinion evolution in online collective behaviour[END_REF]. Most of the existing models assume that individuals are influenced by the opinion of their neighbors [START_REF] Ising | Contribution to the theory of ferromagnetism[END_REF], [START_REF] Degroot | Reaching a consensus[END_REF], [START_REF] Sznajd-Weron | Opinion evolution in closed community[END_REF], [START_REF] Deffuant | Mixing beliefs among interacting agents[END_REF], [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF], [START_REF] Morȃrescu | Opinion dynamics with decaying confidence: Application to community detection in graphs[END_REF]. Nevertheless, it is very hard to estimate these opinions. In order to relax this constraint of measuring the opinions and to model more realistic behaviors, a mix of continuous opinion with discrete actions (CODA) was proposed in [START_REF] Martins | Continuous opinions and discrete actions in opinion dynamics problems[END_REF]. This model reflects the fact that even if we often face binary choices or actions which are visible by our neighbors, the opinions evolve in a continuous space of values which are not explicitly visible to the neighbors. A multi-agent system with a CODA model was proposed and analyzed in [START_REF] Chowdhury | Continuous opinions and discrete actions in social networks: a multi-agent system approach[END_REF]. It was shown that this deterministic model leads to a variety of asymptotic behaviors including consensus.

Due to the complexity of the opinion dynamics, we believe that stochastic models are more suitable than deterministic ones. Indeed, we can propose a realistic deterministic update rule but many random events will still influence the interaction network and consequently the opinion dynamics. For this reason, we consider that it is important to reformulate the model from [START_REF] Chowdhury | Continuous opinions and discrete actions in social networks: a multi-agent system approach[END_REF] in a stochastic framework, specifically, as an interactive Markov chain. Similar approaches for the Deffuant and Hegselmann-Krause models have been considered in the literature (see for instance [START_REF] Lima | Agent based models and opinion dynamics as markov chains[END_REF], [START_REF] Lorenz | Consensus strikes back in the hegselmann-krause model of continuous opinion dynamics under bounded confidence[END_REF]). Although the asymptotic behavior of the model can be given by characterizing the absorbing states of the Markov chain, the convergence time can be arbitrarily large. Moreover, transient but persistent local agreements, called metastable equilibria, are very interesting because they describe the finite-time behavior of the network. Consequently, we consider in this paper an N-intertwined mean field approximation (NIMFA) based approach in order to characterize the metastable equilibria of the Markov system. It is noteworthy that NIMFA was successfully used to analyze and validate some epidemiological models [START_REF] Van Mieghem | Virus spread in networks[END_REF], [START_REF] Trajanovski | Decentralized protection strategies against sis epidemics in networks[END_REF].

In this work, we model the social network as a multi-agent system in which each agent represents an individual whose state is his opinion. This opinion can be understood as the preference of the agent towards performing a binary action, i.e. action can be 0 or 1. These agents are interconnected through an interaction directed graph, whose edge weights represent the trust given by an agent to his neighbor. We propose continuous time opinion dynamics in which the opinions are discrete and belong to a given set that is fixed a priori. Each agent is influenced randomly by the actions of his neighboring agents and consequently influences its neighbors. Therefore the opinions of agents are an intrinsic variable that is hidden from the other agents, the only visible variable is the action. As an example, consider the opinion of users regarding two products Coca-cola and Pepsi. A user may prefer Coca-cola strongly, while some other users might be more indifferent. However, what the other users see (and is therefore influenced by) is only what the user buys, which is the action taken.

Our goal here is to analyze the behavior of the opinions in the network under the proposed stochastic dynamics. One of the main results states that the opinions always asymptotically reach a consensus defined by one of the extreme opinions. Nevertheless, for large networks, we emphasize that intermediate local agreements are reached and preserved for a long duration of time. The contributions of this paper can be summarized as follows. Firstly, we formulate and analyze a stochastic version of the CODA model proposed in [START_REF] Chowdhury | Continuous opinions and discrete actions in social networks: a multi-agent system approach[END_REF]. Secondly, we characterize the local agreements which are persistent for a long duration by using the NIMFA for the original Markov system. Thirdly, we give a complete characterization of the system behavior under symmetric allto-all connection assumption. Finally, we provide conditions for the preservation of the main action inside one cluster as well as for the propagation of actions.

The rest of the paper is organized as follows. Section II introduces the main notation and concepts and provides a description of the model used throughout the paper. The analysis of the asymptotic behavior of opinions described by this stochastic model is provided in Section III. The presented results are valid for any connected networks with a finite number of agents. Moreover Section III contains the description of the NIMFA model and an method to compute its equilibria. In Section IV, we analyze the particular network in which each agent is connected to all the others. In this case, we show that only three equilibria exist and two of them are stable, and correspond to the absorbing states of the system. Section V presents the theoretical analysis of the system under generic interaction networks. It also emphasizes conditions for the preservation of the main action (corresponding to a metastable state) in some clusters as well as conditions for the action propagation. The results of our work are numerically illustrated in Section VI. The paper ends with some concluding remarks and perspectives for further developments. Preliminaries: We use E for the expectation of a random variable, 1 A (x) the indicator function which takes the value 1 when x ∈ A and 0 otherwise, R + the set of non-negative reals and N = {1, 2, . . . } the set of natural numbers..

II. MODEL

Throughout the paper we consider N ∈ N an even number of possible opinion levels, and the set of agents K = {1, 2, . . . , K} with K ∈ N. Each agent i is characterized at time t ∈ R + by its opinion represented as a scalar X i (t) ∈ Θ where Θ = {θ 1 , θ 2 , . . . , θ N } is the discrete set of possible opinions, such that θ n ∈ (0, 1)\{0.5} and θ n < θ n+1 for all n ∈ {1, 2, . . . , N }. Moreover Θ is constructed such that θ N/2 < 0.5 and θ N/2+1 > 0.5. In the following let us introduce some graph notions allowing us to define the interaction structure in the social network under consideration.

Definition 1 (Directed graph): A weighted directed graph G is a couple (K, A) with K being a finite set denoting the vertices, and A being a K × K matrix, with elements a ij denoting the trust given by i on j. We say that agent j is a neighbor of agent i if a ij > 0. We denote by τ i the total trust in the network for agent i as τ i = K j=1 a ij . Agent i is said to be connected with agent j if G contains a directed path from i to j, i.e. if there exists at least one sequence (i = i 1 , i 2 , . . . , i p+1 = j) such that a i k ,i k+1 > 0, ∀k ∈ {1, 2, . . . , p}.

Definition 2 (Strongly connected):

The graph G is strongly connected if any two distinct agents i, j ∈ K are connected. In the sequel we suppose the following holds true.

Assumption 1: The graph (K, E) modeling the interaction in the network is strongly connected.

(t) = X i (t) ∈ {0, 1} a ij
Trust of agent i in agent j τ i Trust of agent i in the network

τ i = j∈K a ij v i,n (t) Probability that X i (t) = θn R i (t) Influence on i to shift opinion to 1 R i (t) = j∈K a ij Q j (t) L i (t)
Influence on i to shift opinion to 0

L i (t) = j∈K a ij (1 -Q j (t)) r i (t)
Expected influence on i to shift opinion to 1

r i (t) = j∈K a ij N n=N/2+1 v j,n (t) l i (t)
Expected influence on i to shift opinion to 0

l i (t) = j∈K a ij N/2 n=1 v j,n (t) νn(t)
Expected fraction of population with opinion θn, νn

= i∈K v i,n (t) K ν + (t), ν -(t)
Expected fraction of population with action 1, ν

+ (t) = i∈K N n=N/2+1 v i,n (t) K or action 0, ν -= i∈K N/2 n=1 v i,n (t) K ν C n (t) Expected fraction of population in C ⊆ K with opinion θn, ν C n (t) = i∈C v i,n (t) |C| ν C + (t), ν C -(t) Expected fraction of population in C ⊆ K with action 1, ν C + = i∈C N n=N/2+1 v i,n (t) |C| or action 0, ν C -= i∈C N/2 n=1 v i,n (t) |C| TABLE I: Notations used
The action Q i (t) taken by agent i at time t is defined by the opinion X i (t) through the following relation

Q i (t) = X i (t) ,
where . is the nearest integer function. This means that if an agent has an opinion more than 0.5, it will take the action 1 and 0 otherwise. This kind of opinion quantization is suitable for many practical applications. For example, an agent may support the left or right political party, with various opinion levels (opinions close to 0 or 1 represents a stronger preference), however, in the election, the agent's action is to vote with exactly two choices (left or right). Similarly, an agent might have to choose between two cars or other types of merchandise like cola as mentioned in the introduction.

Although its preference for one product is not of the type 0 or 1, its action will be, since it cannot buy fractions of cars, but one of them. For ease of exposition, we provide Table I, a list of notations and their meanings.

A. Opinion dynamics

In this work, we look at the evolution of opinions of the agents based on their mutual influence. We also account for the inertia of opinion, i.e., when the opinion of the agent is closer to 0.5, he is more likely to shift as he is less decisive, whereas someone with a strong opinion (close to 1 or 0) is less likely to shift his opinion as he is more convinced by his opinion. The opinion of agent j may shift towards the actions of its neighbors with a rate β n while X j (t) = θ n . If no action is naturally preferred by the opinion dynamics, then we construct θ n = θ N +1-n and assume that β n = β N +1-n for all n ∈ {1, 2, . . . , N }. At each time t ∈ R + we denote the vector collecting all the opinions in the network by X(t) = (X 1 (t), . . . , X K (t)). Notice that the evolution of X(t) is described by a continuous time Markov process with N K states and its analysis is complicated even for small number opinion levels and relatively small number of agents. The stochastic transition rate of agent i shifting its opinion to the right, i.e. to have opinion θ n+1 when at opinion θ n , with n ∈ {1, 2, . . . , N -1}, is given by

β n N j=1 a ij 1 (0.5,1] (X j (t)) = β n K j=1 a ij Q j (t) = β n R i (t).
Similarly, the transition rate to the left, i.e. to shift from θ n to θ n-1 is given by

β n N j=1 a ij 1 [0,0.5) (X j (t)) = β n N j=1 a ij (1 -Q j (t)) = β n L i (t).
for n ∈ {2, . . . , N }. Therefore, we can write the infinitesimal generator M i,t (a tri-diagonal matrix of size N × N ) for an agent i as:

M i,t =    -β 1 R i (t) β 1 R i (t) 0 . . . β 2 L i (t) -β 2 τ i β 2 R i (t) . . . . . .    (1)
with elements corresponding the n-th rown and m-th column given by M i,t (m, n) and

∀n ∈ {1, . . . , N -1}, M i,t (n, n + 1) = β n R i (t), ∀n ∈ {2, . . . , N }, M i,t (n, n -1) = β n L i (t), ∀|m -n| > 1, M i,t (m, n) = 0 and Mi,t(n, n) =    -β1Ri(t) for n = 1, -βnτi for n ∈ {2, . . . , N -1}, -βN Li(t) for n = N. Let v i,n (t) := E[1 {θn} (X i (t))] = Pr(X i (t) = θ n ) be
the probability for an opinion level θ n for user i at time t. Then, in order to propose an analysis of the stochastic process introduced above, we may consider the mean-field approximation by replacing the transitions by their expectations. Then, the expected transition rate from state n to state n+1 for K → ∞, is given by:

β n K j=1 a ij E 1 (0.5,1] (X j (t)) = β n K j=1 a ij N n=N/2+1 v j,n (t).
We have similar expression for transition between state n and n -1.

III. STEADY STATE ANALYSIS

Define by θ-= (θ 1 , . . . , θ 1 ) and θ+ = (θ N , . . . , θ N ) the states where all the agents in the network have an identical opinion, which correspond to the two extreme opinions.

Proposition 1: Under Assumption 1, the continuous time Markov process X(t), with (1) as the infinitesimal generators corresponding to each agent, has exactly two absorbing states X(t) = θ+ and X(t) = θ-.

Proof: We can verify that θ+ and θare absorbing states by evaluating the transition rates using [START_REF] Kerckhove | Modelling influence and opinion evolution in online collective behaviour[END_REF]. If X(t) = θ-, then X i (t) = θ 1 for all i, and the transition rate is

(1, 0, . . . , 0)M i,t = (-β 1 R i (t), β 1 R i (t), 0, . . . , 0). (2)
but as X i = θ 1 for all i, Q i (t) = 0 for all i and so we have R i (t) = 0 for all i. Therefore the transition rate from this state is 0. We can similarly show that θ+ is also an absorbing state.

Next, we show that no other state can be an absorbing state. Consider any state with at least one agent i such that X i (t) = θ n with 1 < n < N . The transition rate from such a state is never 0, which can is easy to see, as we have

M i,t (n, n) = β n τ i .
This implies that as long as such an agent exists, the global state is not an absorbing state. The only states which are not θ+ , θor such a state are those satisfying the following property. Consider X i (t) = θ 1 for all i ∈ S and X i (t) = θ N for all i ∈ K \ S with S ⊂ K and 1 < |S| < K. As the graph is strongly connected, there is at least one agent k in S, which is directly connected to some agent l in K \ S, i.e. a k,l > 0. The transition rate of this k is given by

(-β 1 R k (t), β 1 R k (t), 0, . . . , 0)
As a k,l > 0 and Q l (t) = 1 as all agents outside S has opinion θ N , we have R k (t) > 0. Therefore such a state is also never an absorbing state, which concludes our proof that θ+ and θare the two absorbing states and no other state is an absorbing state.

Considering the NIMFA approximation, we get that the dynamics of the opinion for an agent i are given by:

vi,1 = -β 1 r i v i,1 + β 2 l i v i,2 vi,n = -β n r i v i,n + β n+1 l i v i,n+1 -β n l i v i,n + β n-1 r i v i,n-1 vi,N = -β N l i v i,N + β N -1 r i v i,N -1 (3)
for all i ∈ K and 1 < n < N where

l i = j∈K a ij E[1 -Q j ] = j∈K N/2 n=1 a ij v j,n , r i = j∈K a ij E[Q j ] = j∈K N n=N/2+1 a ij v j,n . (4)
and n v i,n = 1.

We can easily verify that X i = θ 1 , i.e. v i,1 = 1 for all i is an equilibrium for the above set of equations. When v i,1 = 1 for all i, v i,n = 0 for all n ≥ 2 and as a result, l i = τ i and r i = 0 for all i which gives vi,n = 0 for all i, n. Apart from the extreme solutions θ+ and θ-, the non-linearity of system (3) could give rise to the existence of interior rest points which are locally stable. Such rest points are referred to as metastable states in Physics. Metastability of Markov processes is precisely defined in [START_REF] Huisinga | Phase transitions and metastability in markovian and molecular systems[END_REF], where the exit times from these metastable states are shown to approach infinity.

For a given r i = E[R i (t)], the equilibrium state v * i,n must satisfy the following conditions

0 = -β 1 ri τi v * i,1 + β 2 ( τi-ri τi )v * i,2 0 = -β n v * i,n + β n+1 ( τi-ri τi )v i,n+1 +β n-1 ri τi v * i,n-1 0 = -β N ( τi-ri τi )v * i,N + β N -1 ri τi v * i,N -1 (5)
We can write any v * i,n based on v * i,1 , by simplification as

v * i,n = β 1 β n r i τ i -r i n-1 v * i,1 . (6)
As the sum of v i,n over n must be 1, we can solve for v * i,1

as

v * i,1 = 1 N n=1 β1 βn ri τi-ri n-1 . (7)
We then can use this relationship to construct a fixed-point algorithm that computes a rest-point of the global opinion dynamics for all users.

Algorithm outline: The algorithm involves initializing v to a random value. Then, this v can be used to compute the corresponding r i and l i with (4), which is then used to compute the associated v with [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF]. Repeating this recursively results in a fixed point of the NIMFA which is a potential metastable equilibrium. Other potential equilibria can be found by initializing with another random v.

Additionally, we can obtain some nice properties on the relation between r i and v i,n by studying the following function.

Lemma 1:

Consider the function f : [0, 1] → [0, 1] defined as f (x) := N n=N/2+1 β1 βn x 1-x n-1 N n=1 β1 βn x 1-x n-1 (8)
for all x ∈ [0, 1) and with f (1) = 1. We have that f (x) is a monotonically increasing continuous function and takes the values f (0) = 0, f (0.5) = 0.5 and lim x→1 f (x) = 1.

Proof: We can easily verify that

f (0) = 0 1 + 0 + . . . = 0
As β n is assumed to be symmetric around N/2, i.e. β 1 = β N etc., we have

f (0.5) = N N/2+1 β n N 1 β n = 0.5
In order to simplify differentiations, we use some additional variables, ξ = x 1-x ,

p = β 1 + β 2 ξ + . . . β N/2 ξ N/2-1 and q = β 1 ξ N + β 2 ξ N -1 + . . . β N/2 ξ N/2
which gives f (x) = q p+q . As q = 1 + . . . , q p+q is continuous at all points except when ξ is a pole, i.e., when x = 1.

Therefore if we show that f (x) is continuous at x = 1, f (x) is continious is [0, 1].
For this, we use L'Hôpital's rule to calculate lim x→1 f (x), which is

lim ξ→∞ β 1 ξ N + β 2 ξ N -1 + . . . β N/2 ξ N/2 β 1 + β 2 ξ + . . . β N ξ N (9)
Applying L'Hôpital's rule recursively N times, we get this to be 1. Therefore, we have shown that our function f (x) is continuous in [0, 1]. Next, we show that it is monotonic. We have

(f ) = q (p+q)-(p +q )q (p+q) 2 = q p p 2 (p+q) 2 (10)
where

(•) = d(•)
dξ . We see that dξ dx ≥ 0. Therefore, if q p ≥ 0, then we have the monotonicity and increasing property. This can be verified by simplifying it as

q p = N n=N/2+1 β n ξ n-1 p (11)
Each of these terms in the sum are positive and so we have shown that the first derivative of f (x) w.r.t x is positive.

We can use f ( ri τi-ri ) to calculate the probability that an agent i will take the action 1, i.e., 7) and [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF].

N n=N/2+1 v * i,n = f ( ri τi-ri ) from (

IV. COMPLETE GRAPH WITH IDENTICAL CONNECTIONS

Let us first consider a simple graph structure in which each agent is identically influenced by all the others agents i.e., with a i,j = 1 for all i, j ∈ K and i = j, and a i,i = 0. We use ν n to denote the expected fraction of the population of agents with opinion θ n . This is simply

ν n = E[ i∈K 1(X i = θ n )] K = i∈K v i,n K .
Also recall that we introduced the following notation ν -:= N/2 n=1 ν n and ν + = N n=N/2+1 ν n . Under this specific graph structure, we have

lim K→∞ r i τ i = 1 τ i j∈K a ij N n=N/2+1 v j,n = ν + - N n=N/2+1 v i,n K = ν + (12)
for any i ∈ K. Therefore, using (3) and the fact that νn = i∈K vi,n K we can get the dynamics of νn , ∀n ∈ K. Moreover, following [START_REF] Van Mieghem | Virus spread in networks[END_REF], for large values of K we can approximate this dynamics as:

ν1 K -1 = -β 1 ν 1 ν + + β 2 ν 2 ν -, νn K -1 = -β n ν n + β n+1 ν n+1 ν - +β n-1 ν n-1 ν + , ∀n ∈ {2, 3, . . . , N -1}, νN K -1 = -β N ν N ν -+ β N -1 ν N -1 ν + . ( 13
) where ν + = N n=N/2+1 ν n and ν -= N/2 n=1 ν n .
Theorem 1: When N ≥ 4, β n > 0 and β n = β N -n+1 for all n ∈ {1, 2, . . . , N }, the dynamics described in (13) has exactly three equilibrium points at (1, 0, . . . ), (0, . . . , 0, 1) and at a point with ν + = ν -= 0.5. The first two equilibrium points correspond to the absorbing states and are stable, while the third fixed point is an unstable equilibrium (not metastable).

Proof: Let us notice that:

• if ν 1 = 1 and ν n = 0 for all n > 1 one gets ν + = 0;

• if ν N = 1 with ν n = 0 for all n < N one gets ν -= 0.

Consequently is straightforward to verify that (1, 0, . . . ), (0, . . . , 0, 1) are equilibria of [START_REF] Trajanovski | Decentralized protection strategies against sis epidemics in networks[END_REF]. Another equilibrium point of ( 13) is obtained for ν n = ω βn , with 1 ω := N i=1 1 βi . Indeed, for this point one has ν + = ν -= 0.5 and consequently, the right hand side of (13) becomes ω 2 + ω 2 or -ω + ω 2 + ω 2 which are all 0. Next, we show that no other equilibrium exists. For this, we suppose by contradiction that there exists some 0 < ν + < 1, ν + = 0.5 such that ν + = f (ν + ). If such a ν + exists, then that corresponds to another equilibrium point. However, notice that

ν + 1 -ν + = N i=N/2+1 β1 βi ν+ 1-ν+ i-1 N/2 i=1 β1 βi ν+ 1-ν+ i-1 (14)
replacing ν+ 1-ν+ by ξ and accounting for β n = β N -n+1 (symmetry), we must have

N/2 i=1 β 1 β n (ξ -1)ξ i-1 = 0 (15)
for an equilibrium point. However, ξ = 1 corresponds to ν + = 0.5. As ξ > 0, the above equation is never satisfied unless β 1 = 0 (which also means β N = 0). Hence, by contradiction, we prove that the only three equilibrium points are those described.

In order to study the stability of these equilibrium points, we look at the Jacobian and its eigenvalues. If we denote the Jacobian elements by J i,j , where J i,j = ∂ νi ∂ν j , then for all

1 < i ≤ N 2
, and for all N 2 < j < N , we have:

J 1,1 = β 2 ν 2 -β 1 ν + , J i,i = β i+1 ν i+1 -β i , J j,j = β j-1 ν j-1 -β j , J N,N = β N -1 ν N -1 -β N ν -/
We also have ∀ 2 < i ≤ N/2, J 1,i = β 2 ν 2 , and

∀ N/2 < i ≤ N -2, J 1,i = -β 1 ν 1 , ∀ 2 < i ≤ N/2, J N,i = -β N ν N , ∀ N/2 < i ≤ N -2, J N,i = β N -1 ν N -1 .
For all i, j ∈ {2, 3, . . . , N -1} such that |i-j| > 1 we have J i,j = β i+1 ν i+1 when j ≤ N/2 and J i,j = β i-1 ν i-1 when j > N/2. Next, we have

J 1,2 = β 2 (ν 2 + ν -) and J N -1,N = β N -1 (ν N -1 + ν + ). For all i, j ∈ {2, 3, . . . , N -1} such that |i -j| = 1 we have J i,i+1 = β k ν k + β i+1 ν - where k = i + 1 if i + 1 ≤ N/2 and k = i -1 otherwise; and J i,i-1 = β k ν k + β i-1 ν + where k = i + 1 if i -1 ≤ N/2 and k = i -1 otherwise.
At the absorbing state (1, 0, . . . , 0), we have the Jacobian evaluated to have diagonal elements to be 0, -β 2 , β 3 etc. Additionally, the Jacobian becomes an upper triangular matrix as ν n = 0 for all n > N/2. Therefore, the eigenvalues are the elements of the diagonal, which are non-positive and this corresponds to a stable equilibrium. By symmetry, we have the same result for the other absorbing state.

When N/2 m=1 ν m = 0.5, we have that the fixed point corresponding to this distribution satisfies ν n β n = ω for some K > 0. Thus, the first column of the Jacobian can be written as

(K - β 1 2 , ω + β 1 2 , ω, . . . , K, -K) T
The columns j for 2 ≤ j ≤ N/2 are of the form

(ω, . . . , ω + β j 2 , ω -β j , ω + β j 2 , . . . , ω, -ω) T
where ω + βj 2 is the diagonal term of the Jacobian. The j-th column where N/2 < j ≤ N -1 is of the form

(-ω, ω, . . . , ω + β j 2 , ω -β j , ω + β j 2 , . . . , ω) T
Finally, the N -th column is given by

(-ω, ω, . . . , ω, ω + β N 2 , ω - β N 2 ) T
The above matrix is such that each column has exactly one element which is -ω either at the first row (after column index is more than N/2) or at the last row. If we subtract ω(N -2)I from this matrix, it's determinant becomes 0. This can be verified using the properties of determinants that adding a scalar times a row to another row does not change the determinant. We replace the first row with the sum of all rows and this results in the first row becoming all zeroes. This implies that one of the eigenvalues of the Jacobian is ω(N -2) which is positive. Therefore, the equilibrium at ν + = ν -= 0.5 is unstable. The previous theorem characterizes the behavior of the agents' opinion in all-to-all networks with identical connections. Basically, this result states that beside the two stable equilibria in which all the agents rich consensus we may have a metastable equilibria in which the opinions are symmetrically displaced with respect to 0.5.

V. GENERIC INTERACTION NETWORKS

A way to model generic interaction networks is to consider that they are the union of a number of clusters (see for instance [START_REF] Morȃrescu | Opinion dynamics with decaying confidence: Application to community detection in graphs[END_REF] for a cluster detection algorithm). Basically a cluster C is a group of agents in which the opinion of any agent in C is influenced more by the other agents in C, than agents outside C. When the interactions between cluster are deterministic and very weak, we can use a two time scale-modeling as in [START_REF] Martin | Time scale modeling for consensus in sparse directed networks with time-varying topologies[END_REF] to analyze the overall behavior of the network. In the stochastic framework and knowing only a quantized version of the opinions we propose here a development aiming at characterizing the majority actions in clusters. The notion of cluster can be mathematically formalized as follows.

Definition 3 (Cluster): A subset of agents C ⊂ K defines a cluster when, for all i, j ∈ C and some λ > 0.5 the following inequality holds

a ij ≥ λ τ i |C| . ( 16
)
The maximum λ which satisfies this inequality for all i, j ∈ C is called the cluster coefficient.

One question that can be asked in this context is what are the sufficient conditions for the preservation of actions in a cluster, i.e. regardless of external opinions, agents preserve the majority action inside the cluster C for long time. At the limit, all the agents will have identical opinions corresponding to an absorbing state of the network but, clusters with large enough λ may preserve their action in metastable states (long time).

For any given set of agents C ⊂ K, let us denote that

ν C -= j∈C N/2 n=1 v j,n |C| , and
ν C + = j∈C N n=N/2+1 v j,n |C| .
Those values represent the expected fraction of agents within a set C with action 0 and 1, respectively. We also denote by ν C n , the average probability of agents in a cluster to have opinion θ n , i.e., ν C n = i∈C vi,n

|C|

. Now we can use the definition of a cluster given in ( 16) to obtain the following proposition.

Proposition 2: The dynamics of the average opinion probabilities in a cluster C ⊂ K can be written as:

νC 1 κ = -β 1 ν C 1 λν C + + (1 -λ)δ +β 2 ν C 2 λν C -+ (1 -λ)(1 -δ) νC n κ = -β n ν C n + β n+1 ν C n+1 λν C -+ (1 -λ)(1 -δ) +β n-1 ν C n-1 λν C + + (1 -λ)δ νC N κ = -β N ν C N λν C -+ (1 -λ)(1 -δ) +β N -1 ν C N -1 λν C + + (1 -λ)δ (
≥ j∈C aij τi N n=N/2+1 v j,n ⇒ ri τi ≥ λ |C| j∈C N n=N/2+1 v j,n ⇒ ri τi ≥ λν C + ( 18
)
Applying the same derivation for l i , we get li τi ≥ λν C -. Since

l i + r i = τ i , we always have λν C + ≤ ri τi ≤ λν C + + (1 -λ)
. Thus, we can always rewrite the dynamics of a single agent in the cluster as vi,1

τ i = -β 1 ν 1 λν C + + (1 -λ)δ i +β 2 ν 2 λν C -+ (1 -λ)(1 -δ i ) vi,n τ i = -β n ν n + β n+1 ν n+1 λν C -+ (1 -λ)(1 -δ i ) +β n-1 ν n-1 λν C + + (1 -λ)δ i vi,N τ i = -β N ν N λν C -+ (1 -λ)(1 -δ i ) +β N -1 ν N -1 λν C + + (1 -λ)δ i (19) where δ i ∈ [0, 1]
. By taking the sum of each equation over the cluster and dividing by |C|, we get the averages. Each term on the right hand side has a constant factor of the type

(l i + r i )β m ν n λν C
+ and an additional perturbation term of the type The result above shows that instead of looking at individual opinions of agents inside a cluster, we can provide equation (17) for the dynamics of the expected fraction of agents in a cluster with certain opinions. Using proposition 2 and theorem 1, we can immediately obtain some interesting properties for a cluster.

(l i + r i )β m ν n (1 -λ)δ i .
Corollary 1: Let C be a cluster with coefficient λ → 1. Then the cluster has two metastable equilibria at ν C = (1, 0, . . . , 0), ν C = (0, . . . , 0, 1) and one unstable equilibrium corresponding to ν C + = 0.5. This result holds as if λ → 1, (17) simply becomes [START_REF] Trajanovski | Decentralized protection strategies against sis epidemics in networks[END_REF] and therefore the equilibrium points of (13) must correspond to those of (17), but with the fraction of population being the fraction of agents within the cluster and not the whole graph. For example, consider that

K = C 1 ∪ C 2 with C 1 ∩ C 2 = ∅ and |C 1 |, |C 2 | → ∞. Additionally, a ij = 1 for all i, j ∈ C 1 ,
and for all i, j ∈ C 2 . Finally, we have a set

A 1 ⊂ C 1 and A 2 ⊂ C 2 such that a ij = 0 = a ji = 0 for all i ∈ C 1 \ A 1 and j ∈ C 2 \ A 2 , but a ij = a ji = 1 for all i ∈ A 1 and j ∈ A 2 . If |A 1 | < ∞ and |A 2 | < ∞, then we have that λ 1 → 1, λ 2 → 1
which are the cluster coefficients for C 1 and C 2 . In this example, the graph is connected and so the only two absorbing states must be when all opinions are θ 1 or θ N as shown in proposition 1. Corollary 1 implies that if there are two clusters, each with its λ coefficient close to 1 but not 1, then each cluster can stay with contradicting actions in a metastable state. Such a state is not an absorbing state for the system, but can be held for an arbitrarily long duration.

A. Action preservation

We have shown that if a cluster has its coefficient λ → 1, then it can preserve its action as if it starts with all agents with opinion θ 1 or θ N , this state is a stable state and therefore the local action of the cluster is preserved regardless of the external opinions. However, one problem we want to investigate is if this property hold for other λ < 1. The following proposition provides a necessary condition for a cluster to preserve its action and not collapse to external perturbations.

Proposition 3: A necessary condition for C with coefficient λ to preserve its action in a metastable state is if ∃x ∈ (0.5, 1) such that x = f (λx). If no such x exists, then the only equilibrium when the perturbation term δ = 1 from (17) is at ν C + = 1, and when δ = 0, the equilibrium is at ν C + = 0. Proof: The dynamics of the cluster average opinions follow equation (17). We look for equilibrium points for this dynamics under certain values of δ. To find the equilibrium points, we set the left hand side to 0 in (17). For a given δ, by definition of f (•), we obtain that if ν C + is the fraction of agents with action 1 in the cluster, then it must satisfy

ν C + = f (λν C + + (1 -λ)δ) (20)
Having δ = 0 implies that the agents that interact with the agents in the cluster from outside the cluster have all action 0.

If ν C

+ > 0 is an equilibrium point even with δ = 0, this means that regardless of external opinion, the cluster can preserve an action 1. This is true because f (•) is a monotonic function and therefore f (λν

C + + (1 -λ)δ) ≥ f (λν C + )
for all δ ≥ 0. From the properties of f (•) we know that it is monotonic and takes f (x) = x only when x = 0, .5, 1. Additionally, as f (x) < x when x → 0, and as f (•) is continuous, f (λx) < x for all x < 0.5 except at x = 0. However, x = 0 corresponds to the state ν + = 0 which means that this equilibrium is not preserved. If an x > 0.5 exists such that x = f (λx), then regardless of the actions outside C, i.e. δ, we will have ν C + ≥ x as a possible equilibrium.

By studying the opposite case with the majority action inside the cluster being 0 and external opinion 1, we get 1 -x = f (λ(1 -x)) with 1 -x > 0.5, which means that the same condition holds for preserving the action 0 as well, to get ν C -≥ x.

B. Propagation of actions

In the previous subsection, we have seen that a cluster C can preserve its majority action regardless of external opinion if it has a sufficiently large λ. If there are agents outside C with some connections to agents in C, then this action can be propagated. Let τ C i = j∈C a ij denote the total trust of agent i in the cluster C. Let the cluster C be such that it has λ large enough so that νC + > 0.5 exists where νC + = f (λν C + ). Proposition 4: If the cluster C is preserving an action 1 with at least a fraction νC + of the population in C having action 1, then the probability of any agent i ∈ K \ C to chose action 1 at equilibrium is lower bounded as follows.

Pr(Q i = 1) ≥ f νC + τ C i τ i . (21)
Proof: We know that at equilibrium, we must have:

N n=N/2+1 v i,n = f r i τ i (22)
by definition of f (•), [START_REF] Morȃrescu | Opinion dynamics with decaying confidence: Application to community detection in graphs[END_REF] and [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF]. However, we can lower bound R i as

r i ≥ j∈C a ij N n=N/2+1 v j,n (23)
Since the cluster C has an action 1 by at least νC + of its population, we have

N n=N/2+1 v j,n ≥ f λν C + (24)
for any j ∈ C. As f (λν C + ) = νC + , we have

r i ≥ νC + j∈C a ij (25)
Therefore, r i ≥ νC + τ C i and so we have

N n=N/2+1 v i,n ≥ f νC + τ C i τ i (26)

VI. NUMERICAL RESULTS

For all simulations studied, unless otherwise mentioned, we take Θ = {0.2, 0.4, 0.6, 0.8}, β 1 = β 4 = 0.01 per unit of time and β 2 = β 3 = 0.02 per unit of time.

A. General graph

First, we perform some simulations to validate our NIMFA model for the opinion dynamics. For this purpose, we construct a graph with K = 120 partitioned into two sets B 1 = {1, . . . , 40} and B 2 = {41, . . . , 120}. We take a ij = 1 if a connection exists and 0 otherwise. Connections between agents belonging to the same set (B 1 or B 2 ) are established randomly with a probability of 0.3 and connections from B 1 to B 2 and vice-versa are established with probability of 0.02.

For one such randomly generated graph (K, A), we study the opinion dynamics both by simulating a continuous time Markov chain and by looking at the equilibrium points generated by the recursive search algorithm described in Section III. We always start with the initial state given by

X i (0) = 0.2 for all i ∈ B 1 , X i (0) = 0.8 for all i ∈ B 2 .
Setting this as the initial condition, the algorithm 1 gives us the population fraction of agents at equilibrium with action 1, i.e., ν + to be 0.693. To validate this, we run several simulations of the continuous time Markov process with the same initial state and graph structure, and plot the resulting

K i=1 Q i (t).
This plot is shown in Figure 1. We observe that on average over time, the population with action 1 in all simulations are close to the value obtained from the NIMFA model. Next we focus on the approximation done by NIMFA for individual agent opinions. Table II shows the estimated probability of an agent choosing action 1 for some selected agents. This estimate is computed in each simulation by averaging the action taken over a large time horizon during which the system is in the metastable state. We also observed that the system collapsed to an absorbing state 16 times when we did 1000 simulations. Therefore, we estimate P r(Q i = 1) by averaging over the 984 simulations for which the system stayed in the metastable state.

Notice that, as a result of the starting conditions and the graph structure, most agents in B 1 have a high probability of choosing action 0 at the metastable state (equilibrium of the NIMFA), while most agents in B 2 chose action 1. Some agents like 9 and 10 have trust in some agents of B 1 as well as B 2 . Consequently, these agents constantly shift their opinions resulting in a more random behavior for their actions.

B. Graph with clusters

For the next set of simulations, we take another graph structure, but with the same K as indicated in Figure 2. We first randomly make links between any i, j ∈ K with a probability of 0.05. When such a link exists, a ij = 1 and a ij = 0 otherwise. Then we construct a cluster C 1 , with the agents i = 1, 2, . . . , 40, and C 2 with agents i = 81, 82, . . . , 120. We also label agents 40 < i ≤ 60 as set B 1 and 60 < j ≤ 80 as set B 2 . To provide the relevant cluster structure, we made the edge weights a ij = 1 for all τi ≥ 0.444 for all 60 < i ≤ 80. We find that the largest x satisfying x = f (λ 1 x) is 0.95 and that satisfying x = f (λ 2 x) is 0.94. Therefore, if all agents in C 1 start with opinion 0.2 and all agents in cluster 2 start with opinion 0.8, we predict from proposition 3 that ν C1 -≥ 0.95 and ν C2 + ≥ 0.94 in the metastable state. Additionally, applying proposition 4 yields ν B1 -≥ f (0.95 × 0.714) = 0.85, ν B2 -, ≥ f (0.95×0.444) = 0.324 and ν B2 + ≥ f (0.94×0.444) = 0.315. Simulations of the continuous time Markov chain show that our theoretical results are valid even when the cluster size is 40. Figure 3 plots the population fraction of agents with action 1 within a certain set for one simulation. We look for this value in the clusters C 1 and C 2 as well as the sets B 1 and B 2 . C 1 and C 2 are seen to preserve their actions which are opposite to each other. Since B 1 has a significant trust in C 1 alone, the opinion of C 1 is propagated to B 1 . However, as B 2 trusts both C 1 and C 2 , its opinion is influenced by the two contradicting actions resulting in having some agents with action 1 and the rest with action 0. We repeat this for several simulations with the same graph structure and initial opinions to validate our results.

• i, j ∈ C 1 or i, j ∈ C 2 ,

VII. CONCLUSION

In this paper, we have proposed a stochastic multi-agent opinion dynamics model with binary actions. Agents interact through a network and individual opinions are influenced by neighbors action. Our analysis based on a Markov model of opinion dynamics show a consensus like limiting behavior for a finite number of agents. Whereas, when this number becomes large enough, the stochastic system can enter into a quasistationary regime in which partial agreements are reached. This type of phenomenon has been observed in all-to-all and cluster type topologies. Currently, we have studied the dynamics of opinion without any external control of the network. In the future, we will extend this by accounting for an external entity who has a preferred action (a company selling a product for example) and tries to control the actions of the users in a social network by controlling the opinion of a certain subset of agents. This can be interpreted as a company advertising for its product in order for the other agents in the social network to choose its product over that of a rival company.

  The addition of the first type of terms and division by |C| simply becomes κβ m ν n λν C + , with m being n -1,n or n + 1. The averaging of the perturbation terms results in a value between 0 and κβ m ν n (1 -λ) as all δ i are in [0, 1]. This can be therefore be written as in (17) with a new δ ∈ [0, 1].

Fig. 1 :

 1 Fig. 1: Simulation of

K

  i=1 Q i (t) compared to the NIMFA metastable state. The NIMFA works remarkably well even for K = 120.

Fig. 2 :

 2 Fig. 2: Structure of the graph. Any two agents in K may be connected with a 0.05 probability. All agents within a cluster are connected, and the arrows indicate directed connections.

1 i

 1 making C 1 and C 2 clusters with coefficients λ 1 = 0.833 and λ 2 = 0.816 (for the particular random graph generated for this simulation). • 40 < i ≤ 60 and 1 ≤ j ≤ 20, making agents in B 1 trust C 1 with τ C τi ≥ 0.714 for all 40 < i ≤ 60. • 60 < i ≤ 80 and 1 ≤ j ≤ 20 or 80 < j ≤ 120, making agents in B 2 trust both C 1 and C 2 , with τ

Fig. 3 :

 3 Fig. 3: Simulation of i∈SQ i (t) for S = C 1 , C 2 , B 1 , B 2 .We see that C 1 and C 2 preserve their initial actions as given by proposition 3. We also see that as B 1 follows only C 1 , it's action is close to C 1 . As B 2 follows both C 1 and C 2 who have contradicting actions, it has a very mixed opinion which keeps changing randomly in time.

TABLE II :

 II P r(Q i = 1) at equilibrium for some selected agents. P r(Q i = 1) values mentioned in simulation realization 1, 2 and 3 are calculated by time averaging over t = 1000 to t = 3000 hours and the simulation average is taken over 984 realizations.

	Cluster C 1	Set B1
	|C 1 | = 40	|B1| = 20
	Cluster C 2	Set B2
	|C 1 | = 40	|B2| = 20

  Table III compares the predictions based on propositions 3 and 4 with the values obtained in three of our simulations.

	Set S	ν S + = E[ i∈S Q i ] from
		Propositions	Simulation 1,2,3
	C 1	≤ 0.05	0.002 0.004 0.002
	C 2	≥ 0.94	0.994 0.995 0.995
	B 1	≤ 0.15	0.008 0.006 0.005
	B 2	≥ 0.315, ≤ 0.685	0.476 0.471 0.483

TABLE III :

 III ν S + at equilibrium for the indicated set S. Simulation values are calculated by time averaging over t = 100 to t = 1000 hours.


Bilquis Mohamodhosen
email: bilquis.mohamodhosen@ec-lille.fr

Frédéric Gillon

Abdelmounaïm Tounzi

Loïc Chevallier

Topology Optimisation of a 3D Electromagnetic Device using the SIMP Density-Based Method

Keywords: gradient-based algorithm, relaxation of variables, SIMP density method, weighted objective sum

The presented paper proposes a topology optimisation methodology based on the density-based method SIMP, and applied to a numerical example to validate the former. The approach and methodology are detailed, and the results for a 3D basic electromagnetic example are presented. The non-linear B(H) curve is also taken into account.

I. INTRODUCTION

Topology Optimisation (TO) arouses growing interest in the electromagnetic community as it proposes a new way of tackling the general engineering problem: finding the most optimal design to maximize performance and minimize cost. TO offers the freewill of finding the best design without any layout a priori, given appropriate formulation of optimization objectives and constraints to be consistent with the problem.

Various authors have investigated the topic, presenting different methods to deal with it, but our main interest in this paper will be focused on the SIMP (density-based) method [START_REF] Rozvany | Generalized shape optimization without homogenization[END_REF] due to its ease of application and reproducibility. However, this method presents some weaknesses related to the unavoidable relaxation of variables (discrete to continuous) while using a gradient-based algorithm, giving rise to undesired intermediate variables. This paper presents a methodology devised to solve this problem, and a numerical example using Finite Element (FE) Analysis for validation.

II. PROBLEM FORMULATION

In electromagnetic applications, it is usually desired to maximise an objective function such as force or energy, while using the least Material Quantity (MQ) in the design. The structure is modelled by FE, and calculations are done by a laboratory developed calculation code (code_Carmel) coupled with gradient-based fmincon SQP algorithm via a laboratory developed optimisation platform (Sophemis). This has the advantage of outputting an electromagnetically coherent design at each iteration, and hence increasing the probability of obtaining the best 'feasible' design. In this paper, it is desired to produce the best structure made of iron and air as materials, represented by n number of variables ρ varying from 0 to 1, where 0 is air and 1 is iron. To prevent existence of undefined intermediate materials (e.g. at 0.5) in the final design, we propose to add a Feasibility Factor (FF) to force the values to 0 and 1. The optimisation problem is defined as in equation 1.

(1)

The main objective g(ρ) is maximised, and can also be written as the minimisation of the negative value of g(ρ). The coefficient α is used to weigh the prominence of FF w.r.t the main objective in the problem. If α is less than 0.5, FF will have a reduced, yet substantial importance in problem solving as compared to g(ρ). MQ is constrained to β, where β is the maximum amount of iron allowed in the design.

III. NUMERICAL EXAMPLE

A basic 3D numerical example (fig. 1) is used to validate the methodology. The initial domain used is a cube meshed into 4096 hexahedral elements (fig. 1a) where a Magnetic Potential Difference is applied to the nodes highlighted with white dots. The number of variables used is 64, with therefore 64 hexahedra per variable (green outlined cube). The aim is to find the most optimum design that maximises the energy (hence minimises negative value of energy), with the MQ constrained to 75%. The weightage α is set to 0.25. The final design is given in fig. 1b, where the white zones represent air and black represent iron, and the Magnetic Flux Density B given in fig. 1c. Optimisations are done using a linear, as well as non-linear B(H) behaviour assumption with same initial parameters. Both yield the same topology, with a calculation time of 5 mins and max B of 3.5T for linear case, and 20 mins with max B of 1.3T for non-linear case. In the extended paper version, more details will be given on the application of the methodology to optimise the topology of an electromagnet as in [START_REF] Okamoto | Improvements in Material-Density-Based Topology Optimization for 3-D Magnetic Circuit Design by FEM and Sequential Linear Programming Method[END_REF].

Fig. 1 .

 1 Fig. 1. (a) FE Model, (b) Final design, (c) Magnetic Flux Density, B(T)