G N Kariniotakis

DYNAMIC MODELLING OF LARGE AUTONOMOUS POWER SYSTEMS WITH HIGH PENETRATION FROM RENEWABLE (WIND & HYDRO) POWER SOURCES

Keywords: Dynamic modelling, multimachine power system, wind turbines, hydro turbines; ------

This paper presents a model for the transient stability analysis of isolated power systems including various conventional (thermal, gas, diesel) and renewable (hydro, wind) power units. The objective is to assess the impact of a high integration from renewables and to focus on the interaction between hydro units -which have a special dynamic behaviour, and wind turbines. Detailed models for the power system components are developed. Emphasis is given to the representation of different hydro power plant structures. An algorithm is proposed for the identification of the unknown power system parameters. The examined case study is the one of the power system of New Caledonia, where various types of hydro plants and a wind farm are installed. The developed model permits to define penetration limits for the renewable sources, to define rules for the safe operation of the system and finally it contributes to the dimensioning of new power plant installations.

INTRODUCTION

In recent years there has been a growing interest for electricity generation in autonomous systems using renewable power sources due to the high cost of conventional generation. The aim of the paper is to assess the impact on the quality of service (disturbances on frequency and voltage) that renewable sources may have due to their stochastic nature.

The paper presents an advanced multi-machine model for the simulation of the transient behaviour of an autonomous power system. The time scale of the examined phenomena is between 0-30 seconds or more. Various disturbances can be handled like elements switching, symmetrical short-circuits, operation under turbulent wind, islanding etc. Power systems with steam units, gas turbines, diesel engines, hydro plants and wind parks can be simulated. Detailed models are considered for each one of the above types of units as well as for the grid elements and the load.

Emphasis is given on developing detailed models for the hydroelectric plants. The transient behaviour of these plants depends on structural elements like the tunnel, the surge tank, the conduits, etc. In large systems, where a number of hydro turbines are supplied by a common tunnel, oscillations may appear through the conduits of the turbines. Moreover, the particularity of the dynamic behaviour due to the water hammer effect and the long term oscillations, when a surge tank is present, make it interest to study the interaction between hydro plants and wind parks.

Often, a major problem in transient analysis studies concerns the availability of the power system parameters. The paper presents a methodology to estimate these parameters from measurements. The non-linear simplex optimisation technique is used for this purpose. Emphasis is given to the appropriate definition of the objective function, as well as to the definition of constraints on the parameter values to guarantee their physical meaning.

The case study of the French Island of New Caledonia in the Pacific is examined. This system contains steam, diesel and gas units, a wind farm of 2.7 MW and various types of hydro plants. 30% of the demand is covered by renewable sources (hydraulic and wind). A particularity of this system is the high industrial demand (70% of the total demand) due to the nickel mining industry. The various load components have an important impact on the dynamic behaviour of the system.

DESCRIPTION OF THE MODEL

In [START_REF] Kariniotakis | A general simulation algorithm for the detailed assessment of isolated Diesel-wind turbines interaction[END_REF] a model is presented for the simulation of winddiesel power systems. Models for the mechanical part of a wind turbine, an asynchronous generator, diesel engine, synchronous generator and capacitor banks are presented. A new multimachine power system formulation is developed that permits the representation of the power system in a generic and modular way.

Here, for the simulation of large systems, models for the steam units (with single or double reheat) and their speed governing systems [START_REF]Dynamic models for steam and hydro turbines in power system studies[END_REF] have been considered.

Hydro Turbines Modelling

Two types of hydro plants are of interest here :  of a turbine fed trough a tunnel with or without a surge tank.  of n turbines fed by a common tunnel with or without a surge tank. The modelling of the mechanical part depends on the configuration of the plant. Independently of the type of installation, a synchronous generator [START_REF] Kariniotakis | A general simulation algorithm for the detailed assessment of isolated Diesel-wind turbines interaction[END_REF], a voltage regulator [START_REF] Kariniotakis | A general simulation algorithm for the detailed assessment of isolated Diesel-wind turbines interaction[END_REF] as well as appropriate speed regulator models [START_REF]Dynamic models for steam and hydro turbines in power system studies[END_REF] are associated to each turbine.

Initially, the simplest case of a turbine fed by a simple penstock of length L and cross-section A and without surge tank is considered. By assuming incompressible fluid, the rate of change of flow q in the conduit is given from the laws of momentum in the per unit (p.u.) system as :

  l w h h T dt dq    1 1
where h is the head at the turbine, h l represents penstock head losses and T w is the water time constant defined as :

secs base base w H Q Ag L T 
where Q base and H base are the base flow and head and g is the acceleration of gravity (9.81 m/s 2 ). The p.u. flow rate through the turbine is given by :

h G q 
where G is the per unit gate position. The p.u. mechanical power developed at the turbine is :

       DG q q h A P nl t m
where q nl is the per unit no-load flow and A t is a proportionality constant. In certain installations a number of hydro turbines are fed by a common tunnel. In these cases there is a hydraulic coupling between the turbines. The differential equations describing the change in the flow rate at each conduit are given in a matrix form :

                         2 1 1 2 1 2 1 1 1 h h T T T T T T q q w w w w w w  
where T w is the water time constant of the tunnel and T w1 , T w2 are the water time constants of the two conduits defined in a similar way as for the simple case.

Figure 2 shows the dynamic response of two turbines fed by a common tunnel following a 0.1 p.u. gate opening at the second turbine. Due to this change, a variation of 5% is observed at the mechanical power of the first turbine. The figure shows also the particularity in the dynamic behaviour of the hydro units related to the water hammer effect. The gate opening action corresponds to a load increase at the turbine. However, before load increases there is a first dip as shown in Figure 2. A similar case when the hydro turbines can "see" a load increase as the one simulated here can be during a sharp decrease of the wind power contribution to the load (e.g. wind park disconnection). Then, the power dip of the hydro production can be superposed to the lack of wind generation and have a transient attenuation of the power imbalance. Here, the effect of the speed regulators was neglected in order to simulate the most pessimist case. In reality the magnitude of this dip is partly damped by the regulators' action.

In certain situations when the hydraulic station is in distance from the dam, a surge tank is built near the turbines to damp dynamic phenomena due to flow variations. Figure 3 shows the dynamic response of a system of a single turbine with a surge tank following a 0.1 p.u. load increase. Two models that take into account or neglect the effect of the surge tank are compared. The action of the speed regulator is taken into account. The surge tank does not affect significantly transient behaviour in very short-term. It provokes however oscillations of several quantities that have a period in the order of a minute. Continuous wind power variations may act as disturbances that excite these long-term dynamics.

Figure 4 shows the structure of the largest hydro plant at the island of New Caledonia (multiple units, surge tank).

IDENTIFICATION OF UNKNOWN PARAMETERS

A major problem in the simulation of the dynamic behaviour of a power system is to have accurate values for the various element parameters. Often parameter values are known by the manufacturer but in many cases they do not correspond to effective values due to ageing problems, modifications due to maintenance etc.

Here, a procedure is proposed for the identification of unknown parameters from measurements. The optimisation algorithm is the non-linear simplex method under constraints [START_REF] Kariniotakis | Contribution to the development of an advanced control system for wind-diesel power systems[END_REF]. The general identification scheme is shown in Figure 5 in which : [P] is the vector of the unknown power system parameters.

  

   T m AC AC AC AC y y y t Y , , , 2 1  
is the output of the physical system obtained by the test. i AC y is the measured value of the power system frequency or bus bar voltage at time instant i, i=1,2,...,m.

  

   T m ID ID ID ID y y y P t Y , , , , 2 1  
is the output of the simulation model. i ID y is the estimated by the model frequency or voltage value at time instant i, i=1,2,...,m. [U(t)] is the disturbance vector (e.g. load change), which is the input both to the physical system and to the model.

The output of the physical system

    t Y AC
, which is obtained by the test, is compared to the output of the simulation model

    P t Y ID ,
. Below is given the formulation of the optimisation problem for the case that the parameters of a diesel unit are identified. The following parameters are considered as unknown (

  D K R T H P I D , , ,  ) :         m i i ID i AC i f f w P J

Formation of the cost function J([e(t)])

Minimisation of J([e(t)])

Error Of primary importance in the optimisation process is the way the objective function is defined. In the above formulation, this function is the weighted sum of the squared error mismatches between measured and simulated power system frequency. The inclusion of weights is determinant for the success of the identification process. Figure 6 shows a series of frequency measurements performed for a period of 7 seconds. The measured curve can be divided in 3 sectors : between 0-1.1 sec, where the first dip takes place (16 measurements), between 1.1-3.6 sec of the first flat part (30 measur.) and between 3.6-7 sec (45 measur.). If weights are not considered, as is often seen in the literature, the optimisation process risks to converge to a set of parameters that does not properly represent the dynamics of the first dip. The measurements included in the first dip need to be overweighed in the error function.

CASE-STUDY

The case study examined here is the island of New Caledonia located 1500 km east to Australia. Its surface is 18.760 km 2 for a length of 400 km and a width of 40-50 km. The population is 145.000 inhabitants (90.000 for the capital city Nouméa). The island disposes more than 25% on the world's known nickel resources -3 rd producer in the world. The nickel mining industry is the principal base of the economy of the island.

The installed power is 273 MW. The annual production is 1.400 GWh (1994, +9%). The electric demand in the island is between 110-210 MW. 70 % of this demand comes from the industrial sector (peak 126 MW). The peak of the public sector is 79.5 MW (1997). 75% of public demand comes from the capital Nouméa in the southern part of the island. Figure 7 represents a typical day of the power system demand and generation.

Industrial demand is located at a nickel-melting factory and is composed principally by three melting furnaces of 40 MW each. This demand is in general constant in a weekly basis. Important variations are determinant for the form of the daily load curve (see Figure 7). The base load is supplied by diesel units and alternatively by steam or hydro units depending on the season.

The structure of the grid (2 closed loops) is shown in Figure 7. It consists mainly of 150 kV and 33 kV lines. A reactance of 5 MVAr is installed at the level of Népoui plant for the absorption of the excess reactive power generated by the high voltage (150 kV) transmission lines.

In very short-term, transient variations of the industrial load have an important impact on the system frequency. For example, variations of a few MWs may appear in a few seconds. The three furnaces dispose protections for their disconnection which are activated when the system frequency passes the threshold of 48.6 Hz. During the study the following major event was observed. A steam unit (37 MW) was disconnected for some unknown reason. This disconnection provoked a decrease in frequency that activated the protections of the furnaces. The ensemble of the three furnaces (75 MW) was disconnected simultaneously resulting to a frequency deviation to 51.1 Hz. At the same time the underfrequency protection relays of the wind park were activated resulting to the disconnection of the park from the network. Despite the major superposed disturbances observed the system did not collapse.

In order to represent accurately the operating conditions at the grid measurements from the installed SCADA system were considered for the power flows and the voltage at various points of the network. These data permitted to verify assumptions on the transformer settings, the transmission line parameters etc and perform accurate load flow analysis.

Tim e (s) Frequency

In order to analyse the behaviour of the system under higher wind penetration conditions (up to 20%) a wind farm of 20 MW is assumed. Figure 8 shows the frequency deviation following the simultaneous disconnection of all the wind turbines (total 20 MW). This frequency variation is in acceptable levels. The figure shows also the voltage variations at various bus bars of the system.

CONCLUSIONS

This paper extends previous work on the modelling of isolated power systems [START_REF] Kariniotakis | A general simulation algorithm for the detailed assessment of isolated Diesel-wind turbines interaction[END_REF]. Emphasis was given to the representation of hydro units. The particularity of the transient behaviour of these units was discussed in detail through simple examples. The case study of the system of New Caledonia was presented. This system presents great interest due to the various types of units, the high integration from renewables, and its load structure.

The dynamic model proves to be a useful tool for the design of new renewable installations, the definition of maximum penetration limits for wind power, as well as for establishing rules for the safe operation of the power system.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1 : Schematic representation of a plant of two hydro turbines fed by a common tunnel.

Figure 4 :

 4 Figure 4 : Schematic representation of the Yaté plant at the island of New Caledonia ( 4 turbines of 17 MW each ).

  measurement at time instant i and i ID f is the frequency as obtained by the simulation model. H : inertia time constant (sec), TD : diesel engine prime mover and governor time constant (sec) R : regulation parameter (p.u. Hz / p.u. kW) KI : secondary (integral) control (p.u. kW / p.u. Hz) D : load damping coefficient (p.u. kW / p.u. Hz)

Figure 5 :Fig. 6 :

 56 Figure 5 : Identification scheme.

Figure 7 :

 7 Figure 7 : (LEFT) One-line diagram of the system of NewCaledonia. (RIGHT) Generation and demand during a typical day.

Figure 8 :

 8 Figure 8 : Frequency and voltage variations following the disconnection of a wind farm of 20 MW.

ACKNOWLEDGEMENT

The Utility of New Caledonia ENERCAL is greatly thanked for the precious collaboration and the data provided for the purpose of this study.


G Kariniotakis

M Matos

V Miranda

ASSESSMENT OF THE BENEFITS FROM ADVANCED LOAD & WIND POWER FORECASTING IN AUTONOMOUS POWER SYSTEMS

Keywords: Unit Commitment, Wind power forecasting, Load forecasting, Large autonomous power systems;

In this paper, a methodology is developed to assess the benefits from the use of advanced wind power and load forecasting techniques for the scheduling of a medium or large size autonomous power system. The power system scheduling is optimised through a genetic algorithms based unit commitment model that simulates in detail start-up/shut-down procedures of the power units, ramp constraints, generation limits etc. Different types of forecasts are considered as input to the unit commitment model and the operation costs are estimated for each case. Emphasis is given to define appropriate simple forecasting models that can be used on-line. The performance of these models is used as reference to evaluate advanced techniques. The case-study of the Greek island of Crete is examined. A simple model is developed to explain in an intuitive way how the power system structure might attenuate the effect of inaccuracy in forecasts. Finally, the impact of forecasting accuracy on the various power system management functions is discussed.

INTRODUCTION

Although the general operation rules are the same, isolated power systems, when compared to interconnected ones, present additional constraints in their operation.

Generally, in isolated systems, the grid is weaker and spinning reserve needs to be a greater percentage of the installed power (since there is no available help from a neighbouring system) while dynamic security problems are more likely to occur. When a high penetration from renewables, namely wind power, is foreseen, these problems increase, due to the volatility of the power source. In consequence, these systems tend to be managed in a conservative way, leading to an uneconomic operation that increases the already large costs of electricity in islands.

Short-term load forecasting is of primary importance for the efficient management of any power system. In the case of systems with high wind power penetration, as can be the case of isolated systems, forecasts of the wind power production are also required. In large systems containing steam units, scheduling is performed for horizons up to 48 hours ahead. Nowadays, an average prediction error up to ~3% is considered as typical for load forecasting in interconnected systems. However, wind power forecasting models have a significantly lower performance especially for long horizons. If wind power production is considered as "negative load" then, for a system with high penetration (i.e. >20%) it can be assumed that an equivalent high percentage of load is predicted with an error much higher than the typical value of "~3%".

This paper investigates the impact of forecasting accuracy on the various power system management functions. Emphasis is given to the unit commitment (UC) one since it is related to actions like connections or disconnections of power units and hence, it has a major influence on the operation costs of the power system itself.

In order to evaluate the power system operation costs an advanced unit commitment model is used based on the genetic algorithms optimisation technique [START_REF] Miranda | New evolutionary algorithm with chromosome repair and neighborhood digging for unit commitment with fuzzy constraints[END_REF]. The UC model is developed for on-line short-term scheduling of autonomous systems with various types of power units (steam, diesel, combined cycle, hydro, gas turbines, PV plants, wind turbines, etc.). The units characteristics as well as different operating strategies are modelled in detail.

Advanced forecasts are generated using fuzzy logic based models for each case of wind speed and load forecasting. Fuzzy modelling was found to outperform simple as well as other advanced techniques [START_REF] Kariniotakis | Advanced shortterm forecasting of wind power production[END_REF]. Linear ARMA models are also considered as an alternative.

The UC model together with forecasting models have been implemented in a pilot Energy Management System (EMS) for the scheduling of the power system of Crete (project JOULE III: CARE). The island of Crete is taken here as a case-study (peak load: ~300 MW, projected wind power: ~90 MW).

The use of an advanced EMS such as the one under development in the CARE project permits an operation closer to the limits, without jeopardising security. This is achieved mainly by using wind power forecasts based on frequently updated measurements, by on-line monitoring steady-state and dynamic security and by reviewing the unit commitment in a short period cycle.

The paper presents the methodology followed to evaluate the "cost" of forecasts inaccuracy. A detailed and realistic simulation procedure is developed that differentiates this work from similar in the literature [START_REF] Watson | Investigation of the economic benefits to utilities from wind speed forecasting at wind park sites in Crete[END_REF][START_REF] Androutsos | Evaluation of an advanced control system for a medium sized power system with large wind penetration[END_REF]. A simple model is developed to describe in an intuitive way the impact of forecasts accuracy on operation costs.

The UC model is run for the case-study of Crete using simple, advanced as well as "perfect" (actual time-series) forecasts and the operation costs are evaluated for each case. Perfect and simple forecasting define a zone in which the performance of advanced models is situated.

WIND POWER AND LOAD FORECASTING

In this Section the performance of simple and advanced forecasting methods for load and wind power forecasting is analysed for the case study of Crete.

The performance of simple predictors is taken as the lower acceptable performance and is used as reference. The performance of an advanced method is expected to be between that of simple predictors and "perfect" forecasting. Given that it is interesting to evaluate operation costs as a function of different levels of forecasting accuracy one can apply Monte-Carlo techniques to generate pseudo-forecasts of varying accuracy. However, it is of primary importance to use also forecasts directly generated by an advanced model since such forecasts carry characteristics difficult to simulate (e.g. correlated errors according to the time-step).

Concerning load forecasting, hourly time-series of 4 years (1994)(1995)(1996)(1997) have been used. The data up to 1996 were used for the estimation of the advanced models parameters, while the data of 1997 are used for evaluation purposes. The raw load data have been pre-processed to eliminate the effect of power cuts or black-outs.

The simple models considered here to predict load consist in taking the load of the previous day (Day-1) or the load of the same day last week (Day-7) as forecast. The "Day-7" model predicts better the load of Mondays and Sundays, while the "Day-1" model predicts better the rest of the days. A combined predictor according to the type of day has a MAPE (mean abs. percent. error) of 5% for the 1997 data. This performance can be defined as the maximum acceptable error by an advanced model to be worthwhile for on-line use. The "Day-1" model on its own has a 5.7% MAPE error. Advanced models based on ARMA and fuzzy modelling were found to have a performance between 1-5% depending on the time step and the type of model.

Concerning wind speed, a hourly time-series of Crete is also considered. The performance of persistence ("wind speed in the future will be the same as now") as well that of moving averages predictors is evaluated straightforwardly from the datasee Figure 1. An important issue is that from 15 up to 48 hours ahead neither persistence nor a moving average predictor can outperform the use of the mean-value of the time-series

( 3 8. U  m/s) as a predictor [RMSE( U )=4.29 m/s].
The improvement with respect to persistence as obtained by the mean-value predictor is between 1.6% for the 15 th hour ahead and 26.7 % for the 48 th hour ahead.  As a conclusion, for long-term horizons (15-48 hours ahead), an advanced method should outperform persistence by at least 1.5%-26.7 %, according to the time-step, to be worthwhile for on-line use.  Indeed, persistence should be used as a reference for a specific time-step k only if RMSE persist. (k)< data , where  data is the standard deviation of the time series. If not, the performance of an advanced method should be compared to the standard deviation. In the short term (1-20 hours ahead), the advanced timeseries models based on ARMA or fuzzy modelling provide an improvement between 1-13% w.r.t. persistence. For longer horizons ranging between 21-48 hours ahead the mean-value predictor outperforms the advanced time-series models. For such horizons, models based on meteorological information can be used if they provide an improvement w.r.t. persistence higher than 15-26.7% depending on the time step. Otherwise, it is always preferable to use the mean-value predictor.

The above analysis, although restricted to a specific wind speed time series, leads to the useful conclusion that different types of models should be used according to the time-step. When planning horizons are long (24-48 h), then combined forecasting approaches should be followed. A simple combined method is to use persistence for early steps and the mean-value for long steps. Similarly, a combined advanced approach can be to use time-series models in short term and meteorological information based models in long term. The combined approach provides optimal accuracy for the whole horizon and this is beneficial in short term for the planning of the fast units (diesel, gas turbines), as well as in long-term for the slow units (e.g. steam units).

THE UNIT COMMITMENT MODEL

In interconnected systems, unit commitment is usually performed off-line, typically with an horizon of about a week or two (moving window), with hourly time-steps. This gives the basis for performing the economic dispatch every 10 or 15 minutes, most of the times including also reactive power dispatch and perhaps security constraints related to major contingencies. In small isolated systems, on the other hand, a simple unit scheduling is usually necessary, due to the simplicity of the system, even when renewable power sources are present. The latter case was conveniently addressed in the control advice system of [START_REF]Final report of EU-DG XII JOULE II project[END_REF].

In medium-sized or large isolated power systems, however, a different approach is necessary, namely when there are different types of thermal units (steam and gas turbines, diesel units), hydroelectric units and dispersed plants that use wind, solar and other renewable sources. The problem is more stringent if we want to allow for a high penetration of renewables, whose generated power must be forecast and has some degree of uncertainty.

Trying to cope with these needs, economic operation is divided in a unit commitment function and a dispatch function that are performed in sequence, with an optional intermediate decision step that allows the user to take into account information produced by a fast security assessment module. In this scheme, necessary because forecasts have a strong influence in the overall generation schedule, the unit commitment module can no longer be an off-line process as it generally is.

The unit commitment module itself determines what generators will be on the grid in the next intervals, trying to optimise all the costs involved (running costs related to fuel consumption, ramping costs in thermal and diesel generators, shut down and start-up costs), while taking into consideration all the technical constraints (power balance, minimum down and start-up times, technical limits of generators, ramping limits of thermal units, maximum wind penetration, spinning reserve requirements, etc.). The module also produces a pre-dispatch, that is, an approximation of the set points of each generator or wind park. The unit commitment is run in two cycles; an "external" cycle with a period of four hours that aims to generate guidelines (end-of-period constraints) for an "internal" basic cycle run each 20 minutes. The external cycle is so designed to account for the slow units like steam generators. Thus, it considers an horizon of 48 hours ahead. The internal cycle considers an horizon of 8 hours ahead. The output of the internal cycle is a pre-dispatch operation scenario for all the power units in the next 8 hours. Economic dispatch is then performed for the next 20 min. period, in order to propose set points to the generators. The input to the two UC cycles and the economic dispatch are forecasts for load and wind power. The block diagram in figure 2 shows the overall scheme.

ESTIMATION OF OPERATION COSTS

This Section describes the steps followed to simulate the power system operation in a realistic way in order to accurately estimate operation costs.

1.a) Definition of the power units operation schedules. The unit commitment model runs each 20 minutes using forecasted load and wind power values as input. At each time step t it produces a schedule UC forec (t) for the next 8 hours having an operation cost F' forec (t). The estimation of this cost is based on the forecasted load and wind power. 1.b) In real operation conditions the operators apply the proposed UC forec (t) schedule in order to cover the real load of the system. The real operation cost F forec (t) is thus calculated as a function of the real and not the forecasted load and wind power. 1.c) Over a long period T of time however, the total operation cost is not the sum of the F forec (t)'s costs. This is due to the sliding window scheme considered here, which has the following particularity: each 20 minutes the operators receive advice on the actions to be taken in the next 20 minutes up to 8 hours ahead. It is obvious that consecutive schedules propose actions for periods that overlap. The operators will execute only the actions corresponding to the next 20 minutes, while for 40, 60, etc. minutes ahead they will wait for the updated schedules. It is thus assumed that the real operation cost will be the cost of the actions proposed at the first step of each schedule. The sequence of the "first step ahead" actions is denoted as UC T forec . The corresponding cost F T forec is estimated for the real load and wind power conditions (by a simulation of the power system operation in time). 2.a) The unit commitment program is run using "perfect" forecasts as input, that is, using the real load and wind power values instead of forecasts. As in step (1.a) a series of actions UC perf (t) (connections/disconnections of power units) are obtained for the next 8 hours. The operation cost is estimated as F perf (t) . 2.b) As in step (1.c) the sequence of the "first step ahead" actions is denoted as UC T perf . The power system operation is simulated under the UC T perf series of decisions to estimate the operation cost F T perf of the system over a long period T. 3.) The economic impact over T is computed as the difference :

Forecasting accuracy impact(T)= F T forec -F T perf . The step 1 can be repeated for a number of different forecasting techniques having different accuracy. When simple techniques are used for both wind speed and load forecasting, then the operation cost is denoted as F T simple . The values F T simple and F T perf define an envelop where the cost of advanced methods should be located. The size of this envelope indicates the size of investment worthwhile to do on advanced forecasting techniques.

Operation costs vs forecasting accuracy.

This paragraph presents a simple model developed to explain the effect of the forecasting accuracy on operation costs. More specifically, the cost of spinning reserve is considered. The reasoning can be easily extended to include other types of costs like start-up/shut-down costs of units, loss of load costs in case of contigencies, etc. It is assumed that forecasting "inaccuracy" k takes values between 0 for accurate forecasts and 1 for erroneous ones. Let's suppose that the forecasting model provides a load profile F for the next planning horizon (e.g. 48 hours).

Taking into account accuracy, the real load of the system is F(1k). If the UC program considers forecasts as they are and schedules units to cover F, then there will be situations of excessive spinning reserve (when real load is near F(1-k)) or situations with lack of reserve (when load is greater than F). A different strategy is that the UC program schedules units to cover the worst case of a load L + =F(1+k). Then, there is no risk of lack of reserve but costs will be greater most of times without justification since load will be less than estimated F(1+k). The fraction kF acts indeed as a required spinning reserve to the unit commitment optimisation procedure. The nominal power of the units scheduled to operate in order to cover this load will be P N =L + + where is the spinning reserve due to the fact that P N is a step function with values higher or equal than L + . The value of depends on the system structure, that is on the size of the power units with respect to the system load. Now, if the real load in the system is L + =F(1+k), then the observed spinning reserve is SR=. In this case the forecasts inaccuracy does not involve a cost.

If the real load in the system is L -=F(1-k), then the observed spinning reserve is SR=2Fk+. This last expression tells indeed that : SR=f(accuracy)+f(system structure).

The term depending on the system structure is related to the fact that the unit commitment function 'fits' a step process (connection/disconnection of units) to a continuous one, which is the load profile. It is evident that inaccuracy in forecasts is partly damped by a filtering effect.

Figure 3 depicts the above concepts. For simplicity, only the power of the peak units to be connected within the examined period is shown. From the figure it is evident that forecasting accuracy is critical in some regions of the load curve : at 06:00 and at 24:00 when the curve slope is steep. At that points the required reserve constraint may trigger the connection of additional units. The figure shows also that the spinning reserve term  is high at certain periods and overweighs on the spinning reserve term that depends on the forecasts uncertainty. Certainly the situation can be different if one considers higher wind power penetration levels and higher uncertainty for wind power especially for longer time steps. However, the continuous update of forecasts and UC schedule (sliding window scheme) tend to correct decisions. For example, according to the figure, the UC advises the operators at 24:00 that they should connect five units at 06:00. As time will approach 06:00 these decisions will be updated using forecasts with higher accuracy due to the smaller lead time. This however would not be the situation if the examined power system had more slow (steam) than fast units (gas, diesel).

CONCLUSIONS

This paper presented a methodology for the evaluation of the benefits from the use of sophisticated forecasting techniques in the case of autonomous systems. An advanced unit commitment model is used to generate realistic schedules for the power system operation and estimate properly related costs. For this purpose, it receives load and wind power forecasts in the form of intervals of uncertainty, instead of crisp numbers. Moreover, it adopts fuzzy constraints to accommodate uncertainty and risk. Specifically, fuzzy constraints for wind penetration permit to reach a compromise between potential risk and the reduced costs of increasing wind generation.

This methodology has been applied for the case study of Crete. The difference in operation costs for the cases of perfect and simple forecasting was found to be less than 1.5% for a projected wind penetration of 20%. The small value of this difference can be explained if one takes into account the arguments developed in the previous sections, especially on the system structure. Future stages of this work include a more detailed assessment by focusing on special load situations (extreme weather conditions, special days, etc.) and on important wind power variations (e.g. due to wind fronts). The benefits from the use of advanced forecasting and scheduling techniques will be also evaluated on-line following the installation of the CARE EMS system in Crete during 1999. This evaluation will focus also on management functions like economic dispatch and fast security assessment.

The economic dispatch function is performed for the first time step of the planning horizon. It produces set points for the generators based on the forecasted load and wind power values. In case of inaccurate forecasts the Automatic Generation Control device will try to establish optimal dispatch having as target bad forecasts.

The role of the dynamic security assessment function is to evaluate the degree of security of the proposed by the UC schedules. Security is checked against a number of preselected disturbances like a major wind power variation. Schedules that might result to large frequency excursions are rejected. When the UC output is evaluated as unsafe, the values of control parameters like the level of spinning reserve, can be properly adjusted and the unit commitment can be re-run to produce a new schedule. Inaccurate forecasts might lead to situations where excessive spinning reserve is allocated because the operation schedule is erroneously characterised as unsafe and inversely.

The impact of accurate forecasting has been examined here from an Energy Management System point of view. When an EMS is used as a decision support tool for the operators, forecasts are displayed through the man machine interface. The operator, using his experience, can assess visually if forecasts are in acceptable levels or not. If he is confident that they are acceptable he will be ready to follow the EMS recommendations. This would not be the case for "visually inaccurate" forecasts (especially for the peaks). Then, even if the implemented management functions are able to handle or "filter" inaccuracy and still produce economic solutions, the system operators would tend to reject them. Forecasts accuracy influences thus the acceptability itself of an EMS.

Figure 1 :

 1 Figure 1 : RMSE performance of persistence and moving average predictors for 48 hours ahead.

Figure 1

 1 Figure 1 shows the RMSE (root mean square error) performance of persistence and of simple predictors based on an average of n past values. Persistence outperforms any moving average predictor if it is used for an horizon up to 20 hours ahead. Between 21 and 48 hours ahead its use is either indifferent or not recommended since moving averages are outperforming. An important issue is that from 15 up to 48 hours ahead neither persistence nor a moving average predictor can outperform the use of the mean-value of the time-series ( 3 8. U  m/s) as a predictor [RMSE( U )=4.29 m/s]. The improvement with respect to persistence as obtained by the mean-value predictor is between 1.6% for the 15 th hour ahead and 26.7 % for the 48 th hour ahead.  As a conclusion, for long-term horizons (15-48 hours ahead), an advanced method should outperform persistence by at least 1.5%-26.7 %, according to the time-step, to be worthwhile for on-line use.  Indeed, persistence should be used as a reference for a

Figure 2 :

 2 Figure 2: General scheme of the Economic Operation Procedure

Figure 3 :

 3 Figure 3 : Unit commitment scheme for 24 hours ahead.

ACKNOWLEDGEMENT

This work is performed under the JOULE III project "CARE"-JOR-CT92-0119, funded in-part by the European Commission.


Cécile Hardouin

Adoption dynamics: sequential or synchronous modelling

Keywords: standard adoption, sequential dynamics, synchronous dynamics, partial parallelism, Markov chain, ergodicity

This paper deals with the choice of dynamics in spatial simulation and modelling. In economical context, N agents choose between two technological standards according to a local assignment rule. The adoption dynamics is sequential if the choices are made one after the other; it is synchronous or partially synchronous if all or some part of the agents choose simultaneously. This paper points out differences among the three dynamics, especially in their evolution and limit configurations.

Introduction

In many applications, we are interested in the study of the evolution of a system, in space and in time. For instance it may concern competition between different species in ecology, diffusion of technological innovation involving social behaviour, particles system in physics, fluid spread model, disease propagation, image sequences. Some models are deterministic such as cellular automata ( [START_REF] Chopard | Cellular automata and lattice Boltzmann techniques : an approach to model and simulate complex systems[END_REF]), differential equations for instance, while others are based on a stochastic framework, like Gibbs dynamics, probabilistic automata , or particle systems.

Our framework is standard's adoption but the results may apply to other applications where the entourage plays a major role in one's decision making process (see for instance [START_REF] Banerjee | A simple model of herd behaviour[END_REF], [START_REF] Bikhchandani | Learning from the behaviour of others: conformity, fads, and informational cascades[END_REF]). We consider a finite set of sites S = {1, 2, • • • , N }. Each site i is associated to an agent who makes a choice X i in a state space E. The state space E can be finite or not. In all the following, we will assume for simplicity E = {-1, +1}, which is associated to a choice between two competitive technologies (see for instance [START_REF] Arthur | Competing technologies, increasing returns, and lockin by historical events[END_REF]). When this choice X i depends of the local context, we say that there is spatial coordination, the spatial dependency being positive if there is cooperation between the agents, and negative in case of competition. Guyon and Hardouin proposed in [START_REF] Hardouin | Standards adoption dynamics and test for non spatial coordination[END_REF] tests for spatial coordination allowing to distinguish between "independent" and influenced choices. We consider in this work probabilistic assignment rules depending on the neighbourhood, and compare different adoption dynamics: sequential, synchronous or partially synchronous. We do not provide new theoretical results but we gather and clarify some results, hoping to shed light on different modellings, offering information tools about the ins and outs. Indeed, the Gibbs sampler theory is well known, synchronous parallelism a little bit less, and partial parallelism is not much investigated in the literature, though the three dynamics are considered in standard's adoption. The main result gives the limit behaviour of these dynamics; first, whatever the dynamics, sequential, synchronous or partially synchronous, it is ergodic; the proof involves standard Markov chains properties; Second, we show that the limit distribution differs with the dynamics; consequences are important, directly linked with the choice of the model for standards adoption or other economic application; one has to figure out that this primary choice will lead to different economic situations. The invariant distribution is the generating distribution for the sequential dynamics, but it remains generally unknown for the other dynamics. However, its characterization is possible for the synchronous dynamics in the particular case of the Ising model; and we can apply some simulated annealing properties to the partially synchronous dynamics.

A scan or a sweep of S is a tour of all the sites. When the scans are sequential and indefinitely repeated, the agents make their decision one by one; then we get the well known Gibbs sampler and it is possible to characterize the probability distribution of the limit configuration. When the dynamics is synchronous, all the agents make their decision simultaneously, there is still ergodicity but it is difficult to explicit the limit distribution (See [START_REF] Geman | Random fields and inverse problem in imaging[END_REF] for a full description). Finally, partial synchronous dynamics run step by step, a significant part of the changes happening simultaneously at each step.

Our purpose is not to discuss about the choice of the dynamics; this depends only on the economic situation to be modeled. We just want to point out that, for a same local assignment rule, the configurations of the systems can differ widely according to a synchronous or sequential course.

In section 2, we briefly describe the difference between deterministic and probabilistic assignment rules, through standard examples. Then we present the sequential, synchronous and partially synchronous dynamics in section 3, followed by their ergodic properties in section 4. Some illustrating examples are given in section 5.

Assignment rules

Let us specify the model and give some notations. S is equipped with a symmetric graph G and i, j denotes a pair of neighbouring sites i and j. If A is a subset of S, we denote ∂A = {i ∈ S, i / ∈ A and ∃j ∈ A s.t. i and j are neighbouring sites} the neighbourhood of A, and ∂i = ∂{i}. Let us note

x = (x 1 , x 2 , • • • , x N ) a realization of X = (X 1 , X 2 , • • • , X N ) in Ω = E S ; for a subset A ⊂ S, x A (resp.
x A ) is the configuration x on A (resp. outside of A), and x i = x {i} . Finally, |A| denotes the cardinal number of A.

The agent i makes his choice according to a local assignment rule π i (. | x ∂i ) depending on x ∂i . We give below two commonly used examples of deterministic and probabilistic rules.

Example 1

Deterministic Majority choice Let S = {1, 2, ..., n} 2 be a square lattice of size n × n, with the 4 nearest neighbours system. The agent i chooses the state +1 (resp. -1) if +1 (resp.

-1) is majority among his neighbours, and makes a choice at random in case of equality. If we add the assumption that the agent also takes into account his own advice, or private information, then there is always a majority state among the 5 sites of ∂i∪{i} and the rule is deterministic; the system is a kind of cellular automata ( [START_REF] Cox | Coalescing random walks and voter model consensus times on the torus Z d[END_REF]). For those two rules, general consensus (same state everywhere) is an absorbing state. That means that if the number of scans is large, one of the technological standards will emerge and dominate, ending unique. The only point is to know which one of the standards will disappear, and the necessary number of scans to determine the winner. The answer depends on the rule, that is the kind of majority, and on the initial rates of the two standards.

Example 2

Probabilistic Ising rule Let us consider an Ising type model; the Ising model was introduced by physicists in ferromagnetism to represent spins configurations. It is frequently used for standards' adoption or voter models (see for instance [START_REF] Galam | Rational group decision making: a random field Ising model at T=0[END_REF]). We note N i (x) = N (x ∂i ) = j∈∂i x j ; then agent i chooses state +1 with probability:

π i (x i | x ∂i ) = exp x i (α + βN i (x)) exp(α + βN i (x)) + exp -(α + βN i (x)) . ( 1
)
This probability is nothing but the conditional distribution probability of a Gibbs field on Ω with a joint distribution given by:

π(x) = Z -1 exp{α i∈S x i + β i,j x i x j } (2)
The normalization constant Z of the joint distribution, which is often computationally intractable, does not occur in the expression of conditional local distributions (1); that's the reason why we work with these ones.

In this example, the choice of the agent depends on two parameters α and β setting the marginal distributions and spatial correlation.The parameter α is a measure of the global frequency of +1 and -1; α > 0 strengthens states +1 while α < 0 increases the number of states -1, and α = 0 balances the two standards. The parameter β (the inverse temperature) determines the resemblance or dissimilarity between neighbouring sites. There is cooperation between neighbouring sites if β > 0, while β < 0 ensures competition. If β = 0, the assignment is independent of the neighbourhood.

If we set α = 0 and β > 0, β rather large, the rule meets the previous majority choice. Then, a deterministic rule can be approximated by a probabilistic rule.

In all the following, we consider the case of a probabilistic rule in terms of conditional distributions; this implies to define properly the underlying model. Let us consider the general positive distribution π on the configuration set Ω:

π = {π(x), x ∈ Ω}, with π(x) > 0 for all x and Ω π(x) = 1 (3)
The positivity condition allows us to define, for all A and x A , the conditional probabilities π A (. | x A ), particularly the conditional distributions {π i (. | x i ), i ∈ S}. When π A (. | x A ) depends only on x ∂A , π refers to a Markov random field, as in the example above (2) (see for instance [START_REF] Guyon | Random fields on a network: modeling, statistics and applications[END_REF]).

Sequential, synchronous dynamics

Let us describe the three dynamics involved in standards adoption. We start with an initial state denoted x = (x 1 , x 2 , • • • , x N ) which represents the initial standards adopted by agents 1, ..., N . Then, one agent (sequential dynamics), or all agents (synchronous dynamics), or some of the agents (partially synchronous dynamics) choose a new standard, each agent i making his choice according to a local probabilistic rule π i (x i | x ∂i ).

The sequential dynamics

A sequential dynamics is defined by a sequence of scans of the set of sites S. For instance, we browse the sites 1 to N , sequentially, in this order.

• step 1: the initial state is

x = (x 1 , x 2 , • • • , x N ).
• step 2: we browse the sites 1 to N ; first, agent 1 turns from

x 1 to y 1 according to π 1 (y 1 | x 2 , • • • , x N ); then agent 2 turns from x 2 to y 2 according to π 2 (y 2 | y 1 , x 3 , • • • , x N )
; and so on; at the k-th site, agents 1 to k -1 have made a choice and their state is (y 1 , y 2 , ..., y k-1 ) while agents k + 1 to N are still in the state (x k+1 , ..., x N ); we (uniquely) relax (modify) the k-th value according to the local conditional assignment rule, and conditionally to the previous configuration;

x k → y k according to π k (y k | y 1 , y 2 , • • • , y k-1 , x k+1 , • • • , x N ).
• step 3: we come back to step 1 with the new initial state y = (y 1 , y 2 , • • • , y N ). Then a scan of S changes the configuration x to the new one y

= (y 1 , y 2 , • • • , y N ) in N steps.
Some variants are possible:

1. the route to visit all sites can be different from one scan to the other;

2. an individual site can be visited several times during the scan, the important point being to visit all sites;

3. the order of the visits can be chosen at random; 4. the release can also be done by groups of sites, one group followed by another one, with S being the union of the groups, and each group evolving internally sequentially.

As we will see, all those sequential procedures are asymptotically equivalent, leading to the same stationary distribution ( [START_REF] Geman | Random fields and inverse problem in imaging[END_REF], [START_REF] Guyon | Random fields on a network: modeling, statistics and applications[END_REF]).

The synchronous dynamics

Synchronous dynamics is also called total parallelism; in fact, all the sites are relaxed simultaneously and a scan of S is realized in 1 step.

• step 1 : the initial state is

x = (x 1 , x 2 , • • • , x N ).
• step 2: we release simultaneously all the states, getting x to y = (y 1 , y 2 , • • • , y N ) with the simultaneous rule on each site k,

x k → y k according to π k (y k | x k ) = π k (y k | x 1 , x 2 , • • • , x k-1 , x k+1 , • • • , x N ), k ∈ S.
• step 3 : we come back to step 1 with the new initial state y. Partial parallelism. Between sequential and synchronous dynamics, we can define partially synchronous dynamics where M of the N agents choose simultaneously.

Let M be an integer, 1 ≤ M ≤ N .

• step 1 : the initial state is

x = (x 1 , x 2 , • • • , x N ).
• step 2: we choose a subset A of S with M elements (|A| = M ) and we simultaneously modify the values of the sites in A, while the other sites remain unchanged, getting x = (x A , x A ) to y = (y A , x A ) with the simultaneous rule on each site of A,

x k → y k according to π k (y k | x k ) = π k (y k | x 1 , x 2 , • • • , x k-1 , x k+1 , • • • , x N ), k ∈ A, x k → x k , k / ∈ A.
• step 3 : we come back to step 1 with the new initial state y = (y A , x A ).

The ratio τ = M N is called the parallelism rate; M = 1 corresponds to the sequential dynamics, while M = N defines the synchronous one. Let us precise that iterating the dynamics, we choose a new subset A at each step 2. This hybrid dynamics depends on the way of choosing A. It can be chosen at random, for instance with a uniform distribution giving the same weight to the N M subsets of S with M elements, or not, for instance we fix a covering of S with subsets of cardinal M .

More generally, we can consider several rates of active sites; let A 1 , ..., A n be some subsets of S such that ∪ i A i = S; at each step, we choose a subset A i of S with probability γ(A i ) > 0 and we update the sites of A i .

We distinguish a particular case of partial parallelism; let us assume that there is a neighbourhood graph on S; a coding subset C is a subset of S such that any two sites of C are not neighbours with respect to this graph. Then, let us consider a partition of coding subsets {C 1 , C 2 , ...C K } of S; if we run the previous partial algorithm with these coding subsets, then it meets the sequential dynamics.

Ergodicity

General results

Let us consider the same generating distribution π for each dynamics; we assume that, for each dynamics, we repeat the scans a large number of times; the following result shows that the final configurations differ from each other; specifically, in the case of sequential dynamics, the generating distribution π is stationary, whereas it is not the case for other dynamics. It underlines the importance of the kind of dynamics when choosing such an agent-based model. The result is obtained writing the dynamics in terms of Markov chains.

Let us note

σ k the k-th scan, x = X(k) = (X 1 (k), X 2 (k), • • • , X N (k)
) and y = X(k + 1) the configurations before and after the k th scan; let us write P = (P (x, y)) x,y∈Ω , the dynamics' transition matrix for the scan σ k defined by:

P σ k (x, y) = P (X(k + 1) = y | X(k) = x), x, y ∈ Ω
The following properties hold for X = (X(k), k ≥ 0), the evolution of these configurations.

Proposition 1 Let X = (X(k), k ≥ 0) be the dynamics generated by probability π.

(1) X is an ergodic Markov chain on Ω.

(2) For a sequential dynamics, the invariant distribution is π.

(3) For a synchronous dynamics, the invariant distribution is ν, and ν differs from π.

(4) For a partially synchronous dynamics with M ≥ 2, τ = M N , the invariant distribution is λ τ , and λ τ differs from π.

We give hereafter the main lines of the proof and refer the reader for instance to [START_REF] Bremaud | Markov Chains. Gibbs fields, Monte Carlo simulation and queues[END_REF] or [START_REF] Norris | Markov chains. Cambridge series in statistical and probabilistic mathematics[END_REF], [START_REF] Winkler | Image analysis, random fields and Markov chain Monte Carlo methods[END_REF] for general results on Markov chains. We complete the results with properties on the limit distributions characterization.

Sequential dynamics

The transition is for one scan

P σ (x, y) = N i=1 π i (y i | y 1 , y 2 , • • • , y i-1 , x i+1 , • • • , x N ),
and P σ (x, y) > 0 for every x, y. If we have different ways of scanning, we note P k = P σ k and µ the distribution of X(0); X is an inhomogeneous chain ( [START_REF] Bremaud | Markov Chains. Gibbs fields, Monte Carlo simulation and queues[END_REF]) and the distribution of X(k) is

X(k) ∼ µP 1 P 2 P 3 • • • P k .
On the other hand, if the scanning order is always the same, σ k ≡ σ for all k ≥ 1, the chain is homogeneous with transition probabilities P = P σ , and

X(k) ∼ µP k .
It is easy to see that π is invariant for each P σ which is strictly positive; therefore π is the stationary distribution. For instance in the homogeneous case, we write ∀x ∈ Ω, P k (x, y) -→ k→∞ π(y).

Hence, if we repeat the scans a large number of times, the "final" layout of the standards depends on π and its parameters.

Application: this result enables one to simulate any law π; it suffices to use it as generating distribution in the sequential dynamics.This procedure is the well-known Gibbs sampler ( [START_REF] Winkler | Image analysis, random fields and Markov chain Monte Carlo methods[END_REF]).

Synchronous dynamics

The general case

Let us write the transition

Q(x, y) = N i=1 π i (y i | x 1 , x 2 , • • • , x i-1 , x i+1 , • • • , x N ) = N i=1 π i (y i | x i ) > 0.
Again, this expression is strictly positive for every x, y, which ensures the ergodicity of the dynamics:

∀x ∈ Ω, Q k (x, y) -→ k→∞ ν(y)
But the stationary distribution ν is different from π (the generating distribution which is also the invariant distribution for the sequential case), and moreover in most cases, not explicit. Indeed, π is no more invariant for Q.

ν verifies νQ = ν; then ν is a left eigenvector associated to Q and the eigenvalue 1. The search of this eigenvector is difficult because of the high dimension of the matrix Q (originated by the large cardinal number of Ω) while in the sequential case, this search of the eigenvector is trivially solved since πP = π.

However, if the π i are the conditional distributions of the nearest neighbours Ising model, it is possible to write ν; we detail this in the following example.

The Ising model

We consider the torus S = {1, 2, • • • , n} 2 , n being even, equipped with the 4 nearest neighbours system. The generating distribution π is the trimmed Gibbs law (1) :

π i (x i | x ∂i ) = exp x i (α + βN i (x)) 2 cosh(α + βN i (x)) , with N i (x) = j:|i-j|=1
x j .

Then the transition for one synchronous scan is

Q(x, y) = i∈S exp y i (α + βN i (x) 2 cosh(α + βN i (x))
One can show that the invariant distribution ν for the transition Q is

ν(x) = Γ -1 exp{α i∈S x i } i∈S cosh(α + βN i (x))
where Γ is a normalization constant. Indeed, let's write

ν(x)Q(x, y) = Γ -1 exp{α i∈S x i } exp{α i∈S y i } i∈S exp βy i N i (x) ;
since for instance for i = (s, t), y s,t (x s-1,t + x s+1,t ) = x s,t (y s-1,t + y s+1,t ), the previous equality is also

ν(x)Q(x, y) = Γ -1 exp{α i∈S x i } exp{α i∈S y i } i∈S exp βx i N i (y), that is ν(x)Q(x, y) = ν(y)Q(y, x). The transition matrix Q is then ν-reversible, which implies that ν is Q-invariant;
From this explicit expression for ν, let us underline two important differences between the sequential and synchronous dynamics:

(i) We have seen that the π i are the conditional distribution probability of a Gibbs field on Ω with a joint distribution given by

π(x) = Z -1 exp{ i∈S αx i + β i∈S j:|i-j|=1 x i x j }.
From this form, we see that π is a Markov distribution and the cliques of the associated neighbourhood graph (the four nearest neighbours) are of order 1 and 2, made up with singletons and pairs of sites at distance 1, which are of two types, horizontal or vertical.

Similarly, we write

ν(x) = Γ -1 exp{ i∈S αx i + i∈S log{cosh(α + β j:|i-j|=1
x j )}}.

Then ν characterizes a Markov field, like π. But the neighbourhood system is quite different: the cliques are the singletons and the four nearest neighbours for π, while they are the singletons and the squares of sites at distance √ 2 for ν.

(ii) Let us denote S + the subset of sites i = (u, v) with u + v even (the black fields on a chequer board), S -the complementary subset (the white fields), and x + (resp. x -) the configuration on S + (resp. S -).

We define

ν + (x + ) = Γ -1 2 exp{a i∈S + x + i } i∈S - cosh(a + bN i (x +
)), and

ν -(x -) = Γ -1 2 exp{a i∈S - x - i } i∈S + cosh(a + bN i (x -)).
We have ν(x) = ν + (x + )ν -(x -): contrary to the sequential dynamics, the synchronous evolutions on S + and S -are independent from each other.

Figures 1 and2 in section 5 illustrate the difference between these dynamics.

Partially synchronous dynamics

Let us denote R the transition matrix of this dynamics. We also denote 1 the indicator function.

The general case

We choose the subset A at random, for instance uniformly in the set of the subsets of cardinal number M (with τ = M N ). Therefore,

R(x, y) = N M -1 A⊂S:|A|=M 1(x A = y A ) k∈A π k (y k | x A ) .
This transition is positive for every x, y. Hence the partially synchronous dynamics is ergodic, with the stationary distribution λ τ ,

∀x ∈ Ω, R k (x, y) k→∞ -→ λ τ (y).
If M ≥ 2, we verify that π is not invariant for R, and then λ τ = π.

Again, λ τ is not explicit but we some properties established for simulated annealing apply.

Proposition 2 ([15] Theorem 2.7 and [START_REF] Trouvé | Massive parallelization of simulated annealing: a mathematical study[END_REF]) (i) If the parallelism rate tends to zero, and for a fix value of the interaction parameter β > 0, we get a continuity property with

lim τ →0 λ τ = π . ( 4
)
(ii) On the other hand, if the interaction parameter β tends to infinity, and for a fix value of the parallelism rate τ, 0 < τ < 1, we get a limit distribution that does not depend anymore on the parallelism rate:

lim β→+∞ λ τ (β) = λ τ = λ(∞) if 0 < τ < 1 (5)
This is no more true for τ = 1. Therefore we can have discontinuity for large β.

These properties are illustrated by Figure 4 in the next section.

More generally, let us consider A 1 , ..., A n some subsets of S; at each step, we choose a subset A i of S with probability γ(A i ) > 0 and we update the sites of A i ; the parallelism rate is no longer fixed but determined by γ; the dynamics is ergodic with limit distribution λ γ if and only if ∪ i A i = S ( [START_REF] Trouvé | Noyaux partiellement synchrones et parallèlisation du recuit simulé[END_REF]).

The coding case

Let us assume that π has a Markov property with respect to a neighbourhood graph. A subset C of S is called a coding subset if for all i, j of C, i = j, i and j are not neighbours for the markovian structure. Let {C 1 , C 2 , ...C K } be a partition of coding subsets of S; at each step we choose a coding subset C k of S and we update its |C k | sites. We define for k = 1, ..., K,

R C k (x, y) = 1(x C k = y C k ) s∈C k π s (y C k s | x).
Then the transition matrix is R(x, y) = R C1 ...R C K (x, y); but for each k and each C k = {s 1 , ...s |C k | } there is no interaction between the sites of C k that all simultaneously change, and all conditional probabilities π s (y

s | x C k ) = π s (y s | x ∂s ), s ∈ C k , depend only of the values of C k ; this leads to write R C k (x, y) = π s1 (y s1 | x ∂s1 )π s2 (y s2 | x ∂s2 )...π s |C k | (y s |C k | | x ∂s |C k | )
; finally R coincides with the transition probability P for a sequential sweep of S. As for a simple example, consider the square lattice for S with the four nearest neighbours system, C is the subset of the "black" nodes, and C is the white ones. Obviously, C and C are coding subsets and S = C ∪ C; changing simultaneously the black sites then the white ones leads to the same result as changing all the sites one by one.

Some illustrating examples

We propose a simulation experiment to illustrate the differences of the previous dynamics. As for the generating distribution we consider the Ising model defined in the previous section with the four or eight nearest neighbours system, with conditional distributions

π i (x i | x ∂i ) = exp x i (α + βV i (x) + γW i (x)) 2 cosh(a + βV i (x) + γW i (x))
where V i (x) is the sum of the four nearest neighbours (at distance 1) of site i and W i (x) is the sum of the four diagonal neighbours (at distance √ 2) of site i.If γ = 0, then we come back to the four nearest neighbours. We will consider α = 0, that is +1 and -1 occur with the same probability.

We initialize at random and for one initialization we simulate the three dynamics described above: the sequential, synchronous and partially synchronous dynamics. Each simulation of the sequential and synchronous dynamics is obtained running 600 scans on a square toric lattice of size 64 × 64. Of course, similar results occur for non toric lattices but we avoid here some disturbances which may be caused by edge effects.

In the case of partial parallelism, we iterate the scans until each site has been visited at least 600 times. We also compute the empirical spatial correlation ρ1 at distance 1 (based on the four nearest neighbours) and ρ8 based on the eight nearest neighbours.

The parameters are β and γ, and the parallelism rate τ. First we present below some examples of realizations obtained from the same generating distribution π with the three algorithms, for different values of β and γ (the parallelism rate is fixed to τ = 50%). Figure 1 present resulting configurations which strongly visually differ, illustrating Proposition 1 (π = ν = λ 0.5 ). We get horizontal stripes for the sequential algorithm; we still meet these stripes together with vertical ones in the partially synchronous case; but the totally synchronous layout is completely different and features a chessboard. Most often in the context of standards' adoption, interaction parameters are positive, agents inclining to imitate their "neighbours"; in this case, we obtain similar configurations ensuing from the sequential and partially synchronous dynamics, unlike the synchronous case; Figure 3 illustrates this case with the set of parameters β = 1, γ = 0; similar phenomena occur for β = γ = 0.5 or β = 0.5, γ = 0, or β = 0.3, γ = 0; in each case, we get similar images for sequential and partially synchronous dynamics, the clustering being more or less emphasised with respect to the interaction parameter value; but the synchronous realization is different. It is not so clear for small values of β; for instance, Figure 3 (last row) presents realizations issued for β = 0.2, γ = 0; it seems now that the three configurations visually look like each other; and the final frequencies of sites +1 are very close, between 49.44% and 50.49%; however, the empirical correlations highlights the differences and isolates the synchronous case. Then we consider only the partially synchronous dynamics. Coming back to the previous section, both limit distributions π and ν are known for the Ising model (sequential and synchronous dynamics), but we can't characterize λ τ for the partially synchronous dynamics. We repeat hundred times each simulation and compute the mean of the empirical correlations ρ1 and ρ8 We draw the evolution of those two correlations when τ increases from 0 to 1 by steps of 5%. We know from property (5) that if β tends to infinity, the stationary distribution doesn't depend of τ, 0 < τ < 1; in practice, this means that we should observe some constancy for large values of β; in fact, we observe that the correlations are constant from 5% to 95% in the cases β > 0 (β ≥ 0.5 is large enough) and γ = 0, but also if β < 0 (β ≤ -0.5) and γ = 0. In both cases we observe constancy or a light gap between the sequential dynamics with τ = 0 and the partial dynamics with τ = 5%, which illustrates property (4); most of all, we observe discontinuity between τ = 0.95 and τ = 1, as foreseen again by [START_REF] Chopard | Cellular automata and lattice Boltzmann techniques : an approach to model and simulate complex systems[END_REF]. The same thing occurs if we permute parameters β and γ (for instance β = 0, γ = 1 or β = 0, γ = -1), see Figure 4. Finally, if the parameters have different signs, we observe different behaviours; property (5) doesn't hold anymore. We observe constant correlations in the case β = -0.5, γ = 0.5; in fact the observed images at different parallelism rates are all similar to the sequential and synchronous dynamics cases, leading to chessboard like configurations (see Figure 5). On the other hand, if β > 0, and γ < 0, the behaviour of the correlations is quite different; the correlations are quite equal for small parallelism rate and then slowly meet the total synchronous correlations values. It seems there is a threshold rate from which the configurations change, from the correlations point of view (see Figure 5). In conclusion, our simulations allow us to illustrate the theoretical results; we have shown that the choice of the dynamics is very important in standards adoption context, as well as in other application fields. Moreover, we point out that except in specific examples, the limit ergodic distributions remain unknown for the synchronous and partially synchronous choices; they coincide or strongly differ for specific values of the dynamics parameters, as presented in the previous examples. We may suppose that the partial synchronous case coincides with the sequential dynamics for positive interaction parameters, whenever the parallelism rate remains much less than 1 (less than more or less 90% to 95% in our simulations).

Fig. 1 : 1 sequentiel

 11 Fig. 1: Simulations of fields with parameters β = 1, γ = -1 for different dynamics Sequential, τ = 0 Partially synchronous τ = 50% Synchronous, τ = 1 sequentiel beta=1 gamma=-1

Fig. 2 : 1 sequentiel- 1 ,

 211 Fig. 2: Simulations of fields with parameters β = -0.5, γ = 0.5 for different dynamics Sequential, τ = 0 Partially synchronous τ = 50% Synchronous, τ = 1 sequentiel beta=-0.5 gamma=0.5

Fig. 3 : 1 sequential

 31 Fig. 3: Simulations of fields with positive interaction parameters for different dynamics β = 1, γ = 0 Sequential, τ = 0 Partially synchronous τ = 50% Synchronous, τ = 1 sequentiel beta=1 gamma=0

Fig. 4 :

 4 Fig. 4: Correlations ρ1 - * -and ρ8 -o-

Fig. 5 :

 5 Fig. 5: Correlations ρ1 - * -and ρ8 -o-