Thomas Paris
email: thomas.paris@loria.fr

Alexandre Tan
email: alexandre.tan@inria.fr

Vincent Chevrier
email: vincent.chevrier@loria.fr

Laurent Ciarletta
email: laurent.ciarletta@loria.fr

Study about decomposition and integration of continuous systems in discrete environment

Keywords: Co-simulation, Decomposition, FMI 4

A complex system is one composed of many interacting heterogeneous entities. This kind of system can be dealt with multi-modeling and co-simulation but individual models may also be heterogeneous (continuous, discrete, event-based...). To manage this complexity, we use MECSYCO 1 , a DEVS 2 compliant environment for co-simulation.

MECSYCO handles heterogeneity issues, but the number of models which may interact during a cosimulation of a complex system raises performance issues and it is important to develop performance measurement tools to study these issues with MECSYCO.

In this article we present modular performance measurement tools for MECSYCO. We exemplify the use of these tools on our "Multi-Room Heating" toy model, a scalable continuous system, to assert the tradeoff between accuracy and computational time when integrating continuous system in a discrete modeling environment. We explore the impact of decomposing a continuous system contained in one FMU 3 into several FMUs which interact. Finally we study how, under some conditions, a large model that cannot be solved on one block, can be decomposed into smaller ones, solved and simulated in a co-simulation on MECSYCO without significant loss of accuracy.

Introduction

Modeling a complex system which combines discrete and continuous behaviors is not simple. More generally, it can be hard to model a complex system composed of many heterogeneous entities. One way to deal 1 Multi-agent Environment for Complex-SYstem COsimulation 2 Discrete EVent System specification 3 Functional Mockup Unit 4 Functional Mockup Interface with this complexity is to use a multi-paradigm approach [START_REF] Vangheluwe | An introduction to multi-paradigm modelling and simulation[END_REF]. The complex system is decomposed into subsystems which are modeled separately, and then these models are combined to make a multi-model. The different models interact together to reproduce the behavior of the entire system. In a co-simulation, each model is associated to its own well-adapted simulator and simulators are coupled and exchange data during the simulation (simulators communicate their time and variables to update their inputs and to synchronize their execution). Then it is interesting to be able to reuse existing checked models [START_REF] Gil-Quijano | From biological to urban cells: lessons from three multilevel agent-based models[END_REF]. But reusing existing models raises the heterogeneity issues.

The growing interest for co-simulation in the M&S community leads to the development of the Functional Mockup Interface [START_REF] Blochwitz | The functional mockup interface for tool independent exchange of simulation models[END_REF]. This is an emerging standard which aims to ease model-exchange and co-simulation between different modeling software. FMI is an interface for time dependent model described with differential, discrete and algebraic equations. FMI for co-simulation solves the heterogeneity issues at the execution level, it describes an interface to interact with a model and its simulator. FMUs for co-simulation just need a FMI compliant modeling environment to manage the execution. Several works, such as [START_REF] Schierz | Co-simulation with communication step size control in an fmi compatible master algorithm[END_REF][START_REF] Bertsch | The functional mockup interface-seen from an industrial perspective[END_REF][START_REF] Andersson | Import and export of functional mock-up units in jmodelica[END_REF], have been done to enable the use of the FMI standards in modeling tools, to test their FMI compliancy, or to improve the usage of FMUs.

We are working on MECSYCO [START_REF] Camus | Combining DEVS with multi-agent concepts to design and simulate multi-models of complex systems (WIP)[END_REF], a modeling environment based both on the DEVS formalism [START_REF] Bernard | Theory of modeling and simulation : integration discrete event and continuous complex dynamic systems[END_REF] and Multi-Agent concepts like Agent&Artifact [START_REF] Ricci | Give agents their artifacts: the A&A approach for engineering working environments in MAS[END_REF]. MECSYCO deals with heterogeneity issues and enables the interaction of models from different software and based on different formalisms in a co-simulation, it can handle FMUs for co-simulation. In this paper, we are interested in the simulation of a large number of models with many interactions. Such systems raise performance issues and we need to scale-up. Then it is particularly relevant to give performance measurement tools to MECSYCO to analyze its abilities (computational time, communication time, accuracy of results etc...). The co-simulation of a complex system implies to decompose it into subsystems which interact. Since FMI is an emergent standard for co-simulation, notably for continuous systems, it is relevant to study the impact of decomposing FMUs in MECSYCO.

In this article, we present performance measurement tools on MECSYCO adapted to its Agent&Artifact architecture. After that we use these tools to collect performance indicators when we decompose a continuous system from Modelica into smaller ones contained in FMUs and when we simulate them in the discrete environment of MECSYCO. These simulations aim to exemplify the use of our tools and to evaluate both the use of FMUs with MECSYCO and the impact of the decomposition of a continuous system when we simulate it as a co-simulation of subsystems.

Background

This section presents the different tools used in this article.

MECSYCO

MECSYCO is a software for complex system modeling and simulation. It enables the reuse of models from different modeling tools (with their simulator) to build a multi-model, then the simulation of this multi-model is a co-simulation between models with their simulator. In order to ensure that, MECSYCO is based on the DEVS formalism [START_REF] Bernard | Theory of modeling and simulation : integration discrete event and continuous complex dynamic systems[END_REF] and on Multi-Agent concepts such as A&A (Agent & Artifact) [START_REF] Ricci | Give agents their artifacts: the A&A approach for engineering working environments in MAS[END_REF].

The architecture of MECSYCO is based on the A&A concept, there are four main entities which describe our multi-models and run the simulation:

• M-agent: Model-agents are the dynamic entities of the system, they handle the simulation. Each of them is in charge of one model.

• Model-artifact: Artifacts are tools used by agents to interact with their environment. A modelartifact is used to encapsulate each model into the DEVS formalism and it makes the m-agent able to interact with it. Each m-agent is connected to one model-artifact.

• Coupling-artifact: These artifacts represent the connections between m-agents and enable data exchange between models during the simulation.

Operations on time and data can be performed during the data exchange to resolve representation heterogeneity issues between models.

• Model: In MECSYCO, the models represent preexisting models with their simulator. They come from other modeling software.

Moreover, in MECSYCO agents use the Chandy-Misra-Bryant algorithm to synchronize their simulator. It is a conservative and decentralized simulation algorithm which enables the distribution of the cosimulation through different computers.

Modelica

Modelica is an object-oriented modeling language adapted to system modeling and simulation [START_REF] Fritzson | Modelica -A unified object-oriented language for system modeling and simulation[END_REF]. Modelica is attractive because it is an object-oriented language which enables a modular and hierarchical construction of models. Moreover, Modelica is an equation-based modeling, hence it is particularly adapted for the design of continuous systems described by differential equation systems. We choose to use Modelica because all these properties enable the hierarchical build of model and then an easy decomposition.

FMI

FMI (Functional Mockup Interface) is an emerging standard which aims to ease model-exchange and cosimulation between different modeling tools [START_REF] Blochwitz | The functional mockup interface for tool independent exchange of simulation models[END_REF]. An FMU is a model exported following the FMI standard. An FMU can be seen as a black box with some inputs and outputs, and an explicit interface to interact. An XML file describes the model (parameters, inputs, outputs etc...) and specifies if some optional C-functions are available. All the functions defined by the FMI interface are stored in a binary. All these elements are stocked in a zip file with extension ".fmu". For co-simulation, FMUs must be managed by a master software which leads the simulation assuming that data exchange between FMUs is restricted to discrete communication points.

Performance measurement tools

This section presents the development and the deployment of our performance measurement tools on MECSYCO.

What do we need to measure?

The first question when we want to analyze a model performance is to wonder what we need to measure. Many parameters must be taken into account to analyze a model performance, they can be classified into four groups [START_REF] Brooks | Choosing the best model: Level of detail, complexity, and model performance[END_REF]:

• Results of the model (easy to understand, accurate, good description of the system behavior)

• The validation of the model (error, accuracy, credibility)

• Resources needed (Construction time and cost, computational time and cost, result analysis time and cost, hardware requirement)

• Future use (portability, reusability)

In our case we focus on measures related to model decomposition and co-simulation. We develop tools for data logging (to compare different execution and compute the accuracy) and for computational time logging. As MECSYCO uses a decentralized simulation algorithm, for distributed simulation we must be able to use the same tools. Consequently we need to log computational time for each model simulated. Additionally, we want to log data about the algorithm execution like the number of events and synchronizations.

Constraints

MECSYCO is an environment for complex system modeling which is still under development. Specific implementations could change and new features will be added in the future. Therefore, to be long-term useful our performance measurement tools must be code-independent and generic. Particularly, many different kinds of models may interact in a MECSYCO co-simulation, then the tools must not be dependent of the nature of model. Naturally, these tools should not disturb the simulation and the measures, especially for time analysis.

Method

Concept

MECSYCO uses the A&A concept, the simulation is conducted by the m-agents which use the artifacts to lead their model and to exchange data. This modular architecture is interesting to collect performance indicators because at each step, agents need to use their artifacts. That means that if we give to our artifact new abilities to log their internal functioning, we are able to get many relevant data about our co-simulation: exchanged events, computational time, the number of internal/external events for each model etc... We choose to use the design pattern decorator which is particularly adapted to our A&A architecture. Indeed, it lets us add some features to our entities without changing their internal functioning. When we decorate an artifact, we give it new features to collect measures each time an agent uses it. The same method can be used on the agents, for instance to measure their computational time.

Implementation

In MECSYCO, the main concepts (agent, artifact...) are directly associated to an implementation. Our main idea is to encapsulate our pre-existing piece of code for agent and artifact into decorators which copy their behavior while adding some features. For example, we have created a model-artifact decorator (Fig. 1), it takes a model-artifact as parameter and copies its behavior by calling its functions. With this architecture we can add some operations before and after the normal model-artifact functions. A model-artifact is used to encapsulate a pre-existing model into DEVS, it implements five main functions. These functions are used for the execution of the simulation algorithm. Adding operations before and after these functions lets us get information about the algorithm such as the number of events exchanged, the computational time to process each step etc... We use the same approach on coupling artifacts to collect data about the events exchanged, and on m-agents to estimate the computational time of the algorithm.

Use

Finally, when we want to get some extra information about a multi-model, we just have to decorate the artifacts (or the agents) used to build it. These new artifacts will log information about events exchanged, computational time, algorithm execution etc... Then we can use them for post-mortem analysis.

To be easy to understand, we choose to create decorators for each specific task (log of times, data, execution etc...). This is not restrictive because decorators are modular and can be composed. We can easily exchange an artifact decorated with another to collect other data or we can decorate several times an artifact to use the features proposed by several tools at the same time. For now, we are able to measure:

• the total computational time for each model

• the computational time to process internal events

• the computational time to process external events

• the number of internal and external events for each model

• the number of synchronizations and exchanged events between two models

Generalization

Our tools let us enhance the behavior of our MECSYCO entities to collect data during the execution. But when we have to test several multi-models, it is very tiresome to change every single MECSYCO artifacts. We want to automate the exchange between different artifacts when we look for performance indicators on a multi-model. The automation of performance indicator collection raise a new issue, on MECSYCO there is no formal structure to build multi-models. Without structure we are unable to create generic functions for the MECSYCO multi-models. Therefore we must add a structure for our multi-models and this structure must be compliant with the actual properties of MECSYCO.

Structuring MECSYCO multi-model

MECSYCO is based on the DEVS formalism, each agent which is in charge of a model using a modelartifact can be considered as an atomic DEVS model. Then a MECSYCO multi-model can be considered as a set of interacting atomic DEVS models. This is helpful because DEVS defines a structure for this [START_REF] Bernard | Theory of modeling and simulation : integration discrete event and continuous complex dynamic systems[END_REF], the

Implementation

We worked on the java implementation of MECSYCO, we have defined a new object multi-model which contains the adapted DEVS coupled model structure. We define an abstract class which contains our new structure and few functions to handle the multi-models (to start a simulation for instance). Now, instead of building a multi-model as a simple process, we define it as a structured object easier to handle and which enables generic analysis. Now with our modular tools we are able to:

• put or remove performance measurement tools

• put several tools on a single entities

• perform measures on subsets of multi-models

Test

In this section, we exemplify the use of our performance measurement tools on a basic continuous system example using FMUs. A continuous system is a good starting point for various reasons. Continuous systems are very common in the modeling and simulation field, and there are a lot of dedicated tools (like Dymola) to simulate them. These specific tools are quite accurate and allow to easily compute a reference solution to compare with the results of the different. We choose to run our tests with a toy example (a basic thermic system) easy to decompose and to make more complex.

Presentation of "Multi-Room Heating"

"Multi-Room Heating" is a model which represents temperature evolution in four rooms under the influence of the outside temperature, this model is closed to the one in [START_REF] Gilpin | Cosimulation solutions using aa4mm-fmi applied to smart space heating models[END_REF]. It is a continuous system defined by two main equations. The equation

C * dT (t) dt = Q(t)
represents the behavior of a room. T (t) is the temperature inside the room, C the heat capacity of the volume of air and Q(t) is the incoming heat flow. This incoming heat flow is determined with the equation

Q i (t) = G * (T a (t) -T b (t))
which represents the behavior of a wall i. Q(t) is the heat flow through the wall, G is the thermal conductance of the wall, T a (t) and T b (t) represents the temperatures at the two faces of the wall.

If we consider that the outside temperature is T out (t) = A * sin(t * f ) + B where A is a amplitude, B an offset temperature and f a frequency adapted to have the evolution of the outside temperature in a day (86400 seconds). The problem with a single room becomes:

C * dT (t) dt = G * (A * sin(t * f ) + B -T (t))
This equation has an analytical solution:

T (t) = T (0) * e -G C * t + A * G 2 * (sin(t * f )-C * f G * cos(t * f )+ C * f G * e -G C * t ) G 2 +C 2 * f 2 -B * e -G C * t + B
This analytical solution could have been interesting but the purpose of this kind of models is not the modeling of one room but the modeling of many rooms interacting (to represents a building for instance). In most cases the analytical solution does not exist or is too hard to compute and we must compute approximate solutions. That is why we use the results of Dymola as a reference for our experiments.

This continuous system is interesting because it can be easily complexified by adding rooms and connections between these rooms. Analytical solutions become quickly complicated to find and we have to use approximate solutions. Our main example is constructed with four rooms (Fig. 2).

Simulations

We simulate our "Multi-Room Heating" using Modelica first, on Dymola, this simulation is used as a reference. Then, we export the entire model in a single FMU, and finally we decompose it into several parts (rooms, walls and the outside temperature), each exported in FMUs. These FMUs are connected to rebuild the "Multi-Room Heating" and are used in a cosimulation on MECSYCO. This allows us to compare the accuracy of MECSYCO, depending on the time step size, with our reference results. Performing the same tests with a single FMU containing the entire model allow us to evaluate the impact of decomposing a continuous system into several FMUs instead of using one single FMU. We choose to export our FMUs with JModelica because it is an open-source Modelica platform. We import the FMUs in the java version of MECSYCO with a FMI model-artifact based on JavaFMI, an open-source java library. In our simulations our systems represents the evolution of the temperature of rooms during 10 days. We choose to use the same constant time step size for each FMU. We follow our multi-model structure to build our test models, and then we run several simulations with our new tools encapsulating our artifacts. These simulations let us test our performance measurement tools.

Results

We compare our results to a set of results from Dymola to compute the accuracy. Previous tests show that the evolution of accuracy is approximately the same for the 4 rooms, so we only display the evolution of error in the room 1. The Figure 3 shows the behavior of our system for three days.

The computational time on Dymola is not displayed because it is very closed to the computational time obtained with one FMU on MECSYCO.

The Figure 4 shows our results concerning the computational time. In both cases (one FMU or several FMUs), the computational time is inversely proportional to the time step size. Moreover, the computational time is increased when several FMUs are used. Indeed, in our particular example, the FMUs are very simple and most of the global computing time is spent for the communications.

The Figure 5 represents the evolution of the accuracy. As expected, when using a single FMU the difference with the results of Dymola is independent of the time step size and are very similar to the results of Dymola. Our test using several FMUs shows that the difference with Dymola varies almost linearly depending on the time step size.

These results could be used to figure out a tradeoff between computing time and accuracy, but they depend on the system under study and the objectives. The tradeoff between computing time and accuracy must be determined for each system. This first experiment mainly allows to verify that the data collected by the measurement tools are consistent and validate their implementation. Then our tools can be used for further performance tests.

Further tests

Adding discrete events

Our first model is a pure continuous system which evolves quite slowly. Now we want to study the impact of adding discrete events which increase our system variation. For that, we add a heater to our model on Modelica. When the temperature inside a room is lower than a value, a heater is put on and heats the room to reach quickly a maximal temperature, i.e. the temperature inside the room stays inside an interval. For this test, we choose to add an heater to rooms 1, 2 and 3 with respectively 293.15 K, 293.15 K and 288.15 K for the setpoint temperatures (see Fig. 6).

We run our test with this new feature and, as expected, we find that MECSYCO is still accurate even if the error increases more rapidly with the time step size. That shows that the tradeoff between accuracy and computational time depends closely on the model.

Making rooms more complex

We find that we can decompose a continuous system and simulate it using several FMUs without a significant loss of accuracy. Now we want to test that with more complex versions of our FMUs. By more complex, we mean that these new FMUs contain more equations or compute more calculation at each step. We try that according two ways: adding useless equations in our Modelica models and adding an algorithm (a matrix product) at each step. With both of these methods we run tests with an increasing number of extra computations.

In our case, the computational time increases linearly with the number of calculations per FMU, but this increase is quite low, we stay in the range of seconds. There is no impact on the accuracy, but at the end we find a limit in term of memory. We try to generate our single FMU of the entire system with JModelica first but its FMU export does not handle too many equations (about 30000) due to memory limitations. Dymola is more resilient so we use it to export our FMUs for this test. With MECSYCO as with JModelica, we find a memory limit due to the size of our FMU. In this case, the decentralized algorithm of MECSYCO combines with the possibility of decomposing our sys-tem let us overcame this limitation by distributing our FMUs co-simulation on several computers.

Simulating large sets of rooms

Another way to make our system more complex is to add rooms and walls. This adds more equations and more connections in our decomposed system. We simulate our system with a varying number of rooms interconnected by walls. We just want to test our computation limits. This raises interesting questions about the build of such multi-models because it is too long to connect each model by hand. To construct these multi-models on MECSYCO we use our multi-model structure to define a parametric multi-model where we can choose the number of rooms and the way they are connected. We define a rectangular set of rooms where they are connected like a grid. We use this model to test our tools on a large set of models, thus we are able to simulate a set of about 615 FMUs and to get their computational time. We verify also that the use of our tools does not limit the number of models we can load.

Study of different decompositions of a system

Until now, we have compared two versions of our "Multi-room Heating" example, one where we considered it as a whole with a single FMU (Fig. 7) and the others where it is totally decomposed into 14 components (Fig. 8). Now we study different ways to decompose our system and to figure out the best in terms of balance between accuracy and computational time. To do so, we propose five different ways to generate FMUs for the co-simulation of our system (from the whole version to the totally decomposed one). As the tools we have implemented are modular, they work whatever the configuration of the FMUs. We use them to collect the computational time and to compute the accuracy of each configuration.

First decomposition, 2 FMUs

For this first decomposition, we get out the outside temperature of our FMU (Fig. 9) and we obtain two FMUs interacting. Figure 11: "Multi-Room Heating" as three FMUs, one for the outside temperature, one for rooms and the last one for walls.

Second decomposition, 3 FMUs

For the second decomposition we split our set of rooms and walls in two almost equal parts (Fig. 10). We get three FMUs.

Third decomposition, 3 FMUs

Another way to get three FMUs is to gather the elements by nature: Rooms, Wall, Outside Temperature (Fig. 11).

Study

Now we have 5 versions of the same multi-model (one FMU, 14 FMUs and the three previous possible decompositions), for each configuration we generate the appropriate FMUs and we build a MECSYCO multimodel. To compare these different configurations we choose to observe the subset of each configuration which contains the "Room 1" to log the temperature of this room at each step and to get the computing time of this FMU. We choose to collect also the total computing time of each multi-model. For each version we observe the accuracy and the computing time of a subset of the multi-model, and the computing time of the entire model. With our new multi-model structure and our decorated artifacts we build the generic tools we need and we use them on our multi-models.

Collected results

The Figure 12 displays computing time measures, the more there are FMUs the more we need time to compute the simulation. It is natural since all these FMUs are simple and most of the computing time is the time needed to exchange data and to synchronize the FMUs. The Figure 13 represents an evaluation of the communication time as the difference between the time spent to compute events in one FMU and the global time of the agent managing this FMU. This communication Finally the Figure 14 shows the evolution of the max difference of temperature between each configuration of the multi-model and the results of Dymola, depending on the time step size. This kind of result enables to determine a tradeoff between accuracy and computational time. For example, the third decomposition is not useful since it has the same evolution than the version with 14 FMUs. The second decomposition is more interesting, there is probably a compensation of errors since it is more closed to the results on Dymola with a time step of 900s than with a time step of 600s. Finally, if we consider that we are accurate enough with an error below 0.4K, the second decomposition is the best because we can keep a time step of 1200s.

Conclusion

We presented a set of modular performance measurement tools for MECSYCO, a co-simulation platform based both on DEVS and on Multi-Agent concepts. We applied these tools to study the impact of decomposition of a continuous system toy example to illustrate their use in the context of decomposition and integration of complex systems.

These first preliminary results are mainly focusing on solving continuous systems with multiple FMUs and on comparing their accuracy to a reference provided by the Dymola tool. Our experiments show that the use of one FMU provided results very closed to the ones of Dymola. When decomposing a complex system into several FMUs, our experiments highlight that our tools enable to study the tradeoff between accuracy and computational time when simulating several FMUs on MECSYCO. This result can be especially interesting when a single FMU can not be executed on a single machine because of memory consumption and imposes to assess different decomposition strategies.

The tools proposed in this article are not specific to continuous systems and we plan to use them on hybrid systems (with continuous and discrete models). They enable more advanced work on MECSYCO performances for potential improvements, for example a time step optimization like in [START_REF] Schierz | Co-simulation with communication step size control in an fmi compatible master algorithm[END_REF]. An approach like in [START_REF] Glinsky | Devstone: a benchmarking technique for studying performance of devs modeling and simulation environments[END_REF] could be interesting to evaluate MECSYCO on various kinds of models. This can lead also to the study of better deployment strategies for our multi-models. The structuration of MECSYCO multi-model is a first step for the integration of the concept of composition in MECSYCO, to enhance the reusability of multi-models and to enable a hierarchical build of them.

Figure 1 :

 1 Figure 1: UML diagram of our model-artifact decorator.

Figure 2 :

 2 Figure 2: "Multi-Room Heating" model.

Figure 3 :

 3 Figure 3: Evolution of the temperature of 4 rooms under the influence of an outdoor temperature during 3 days.

Figure 4 :

 4 Figure 4: Evolution of the mean computational time (with min and max bounds) depending on the time step size, for our model with one FMU or several FMUs (logarithmic scale).

Figure 5 :

 5 Figure 5: Evolution of the mean and max difference between our models (one FMU and several FMUs) and our reference result, depending on the time step size.

Figure 6 :

 6 Figure 6: Graph showing the evolution of the temperature of 4 heated rooms under the influence of the outdoor temperature, view for 3 days.

Figure 7 :

 7 Figure 7: "Multi-Room Heating" as a single FMU.

Figure 8 :

 8 Figure 8: "Multi-Room Heating" decomposed as 14 FMUs.

Figure 9 :

 9 Figure 9: "Multi-Room Heating" as two FMUs.

Figure 10 :

 10 Figure 10: "Multi-Room Heating" as three FMUs.

Figure 12 :

 12 Figure 12: Evolution of the global computing time for each configuration depending on the time step size.

Figure 13 :

 13 Figure 13: Evolution of the communication time for the agent which "contains" the room 1 for each configuration depending on the time step size.

Figure 14 :

 14 Figure 14: Evolution of the max difference of temperature with Dymola depending on the time step size.

ACKNOWLEDGMENTS

We would like to thank Jean-Philippe Tavella from EDF R&D for his contribution on FMU export with Dymola.


Siddhartha Kumar Khaitan
email: skhaitan@iastate.edu

MASTER: A JAVA Based Work-Stealing Technique For Parallel Contingency Analysis

Keywords: Concurrency, multithreading, work-stealing, power systems, parallel contingency analysis

In this paper, we present MASTER, a Java based multithreaded work-stealing technique for parallel contingency analysis in power systems. MASTER analyzes contingencies using time domain simulation and scales contingency analysis task to multiple cores using multithreading in Java. To achieve load balancing, MASTER uses efficient implementation of workstealing algorithm. We discuss several implementation issues and design time choices which are crucial for achieving efficient implementation. Experiments performed with contingencies of a 13029 bus power system shows that MASTER provides high computational gains and provides much better load-balancing than the conventional scheduling techniques.

I. INTRODUCTION

As power systems are growing in their size and complexity, the probability of critical failures is also increasing. To take preventive and corrective actions towards such disturbances, power system operators rely on obtaining information from analysis of a large number of contingencies. Conventionally, due to the small size of power systems, the number of contingencies to be analyzed has also remained small and hence, sequential computing platforms have fulfilled the computational needs of the contingency analysis task. However, with increasing system size and requirements of analyzing higher order (N -x) contingencies, sequential computing platforms are proving to be insufficient in catering to the computational demands of the contingency analysis. Recently, researchers have proposed use of high-performance computing (HPC) resources , however, parallelization also brings the issue of achieving load balancing 1 and maximizing resource usage efficiency [START_REF] Khaitan | Dynamic load balancing and scheduling for parallel power system dynamic contingency analysis[END_REF]. The conventional scheduling approaches (e.g. [START_REF] Santos | Distributed contingency analysis: practical issues[END_REF][START_REF] Morante | Pervasive grid for large-scale power systems contingency analysis[END_REF][START_REF] Mittal | Real time contingency analysis for power grids[END_REF]) either provide poor load-balancing or do not scale well to large number of processors. Thus, the state-of-the-art in power systems calls for efficient parallelization of contingency analysis task, in a manner that provides load balancing and resource usage efficiency maximization.

In this paper, we present MASTER, a Java based multithreaded work-stealing technique for parallel contingency analysis in power systems. MASTER conducts time domain simulation of power system for analyzing contingencies. For parallelization of contingency analysis task, MASTER employs multithreading using Java (Section III). To achieve load balancing on the available cores, MASTER uses the workstealing algorithm (Section IV). While work-stealing algorithm has been used in other domains (e.g. [START_REF] Zhou | Renderants: interactive reyes rendering on GPUs[END_REF][START_REF] Van Nieuwpoort | Efficient load balancing for wide-area divide-and-conquer applications[END_REF]), to the best of our knowledge, our work is among the first to implement a multithreaded JAVA implementation of work-stealing to achieve load-balancing in power systems.

While work-stealing is known to be efficient in theory, achieving the expected efficiency in practice requires careful tuning of the parameters depending on the problem domain and the hardware platform used. An inefficient implementation is likely to significantly slow-down the algorithm [START_REF] Saha | Enabling scalability and performance in a large scale cmp environment[END_REF][START_REF] Vrba | Evaluating the run-time performance of kahn process network implementation techniques on shared-memory multiprocessors[END_REF]. For this reason, we explore several algorithmic factors which impact the efficiency of work-stealing implementation and also discuss the actual choices made in MASTER for efficient implementation. We also present the details of MASTER scheduler along with its implementation in Java.

MASTER is envisioned as a security assessment tool for efficiently analyzing a large number of contingencies and assisting the system personnel in power system control centers. By virtue of using JAVA implementation, MASTER is portable and architecture independent. We have evaluated MASTER on a large power system with 13029 buses. Dynamic contingency analysis simulations have been performed on a multicore processor using 2, 4, 8 and 16 threads. The results show that MASTER offers high computational gains and provides much better load-balanced allocation than a conventional scheduling technique, namely centralized task-queue scheduling.

The rest of the paper is organized as follows. Section II discusses the related work in power system contingency analysis and scheduling techniques. Section III presents the design and architecture of MASTER. Section IV discusses the implementation of work-stealing algorithm. Section V presents the experimental results on evaluation of MASTER and comparison with a conventional scheduling approach. The conclusion and the future work are discussed in Section VI.

II. BACKGROUND

High performance computing (HPC) is a promising approach for accelerating computation intensive and performance-critical tasks using parallelization technique [START_REF] Raju | Domain decomposition based high performance parallel computing[END_REF][START_REF] Raju | High performance computing using out-of-core sparse direct solvers[END_REF][START_REF] Khaitan | A Survey Of High-performance Computing Approaches in Power Systems[END_REF]. However, parallelization also brings the requirement of efficient scheduling for achieving load-balancing [START_REF] Khaitan | Dynamic load balancing and scheduling for parallel power system dynamic contingency analysis[END_REF][START_REF] Khaitan | A High-Performance Parallelization and Load-Balancing Approach for Modern Power-Systems[END_REF]. For accelerating contingency analysis using parallelization, researchers have used different platforms such as MPI (message passing interface) [START_REF] Mittal | Real time contingency analysis for power grids[END_REF] and GPU (graphics processing unit). For parallelization, we use multithreading and discuss the reasoning behind the choice in Section III.

Scheduling techniques for achieving load-balancing in HPC, can be broadly divided into two categories, namely "static" and "dynamic" scheduling. Static scheduling works by using a fixed assignment of tasks on the available processors. The limitation of static scheduling is that when tasks have heterogeneous execution times, static scheduling fails to achieve load balancing. Work-stealing [START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF] is a dynamic scheduling method, which works by allocating tasks to worker threads to finish and letting a worker with no pending task (called thief) steal a task from another worker (called victim).

While Mittal et al. [START_REF] Mittal | Real time contingency analysis for power grids[END_REF] perform load-balancing for steady state contingency analysis, we conduct dynamic contingency analysis which requires higher computation time and resources and hence, also imposes larger penalty of load imbalance.

In this paper, we evaluate single-order contingencies, although our framework permits analysis of cascading contingencies also. Previous work [START_REF] Khaitan | On-line cascading event tracking and avoidance decision support tool[END_REF][START_REF] Mccalley | Dynamic decision-event trees for rapid response to unfolding events in bulk transmission systems[END_REF] has proposed methods for identifying higher order contingencies and capturing their time delayed cascading effects on the network. Different events of contingencies can be seen as taking place as function of time along the dynamic event tree (DET), which accounts for both slow and fast dynamics. Each of these points along the DET can serve as the initiating event for the contingency analysis in the MASTER framework. These would include low frequency, high impact higher order contingencies. We also allow dynamic updating of system conditions which could make a high order contingency more probable due to changes in system configuration or maintenance. In this way, dependent events can be properly captured in the MASTER framework.

III. MASTER: DESIGN AND ARCHITECTURE

MASTER uses a multithreaded scheduler for implementing work-stealing algorithm. In comparison to multiprocessing, our choice of multithreading is guided by the fact that threads are lightweight compared to processes and hence, threadswitching generally incurs less overhead than switching between the processes. Further, threads share the same address space and hence intercommunication between them is also less expensive than that between processes.

MASTER scheduler is designed using Java. The multithreading built into the Java language and runtime platform provides support for running multiple concurrent threads. Compared to languages such as C and C++, the benefit of using Java is that it includes multithreading primitives as part of the language itself and as part of its libraries. Thus, a portable implementation of multithreading can be easily achieved. In contrast, languages such as C do not have built-in multithreading capabilities and hence, to implement multithreading, they must make calls to operating system multithreading primitives, which harms efficiency and portability.

The work-flow of MASTER is explained in Figure 1. The MASTER code is compiled into platform-independent byte-code, which runs on the Java virtual machine (JVM). The JVM runs as a single process which internally spawns many threads. When the MASTER code running inside the JVM asks for another thread, JVM starts another thread which can run on a different CPU core. In this manner, multithreading is achieved on a multicore processor. In hardware, the management of threads is done by the operating system (OS), and the JVM uses the facility provided by the OS. More implementation details are provided in Section IV.

GPUs and MPI are other widely-used parallelization platforms. However, the limitation of GPUs is that they do not share memory with CPUs and are generally slower than CPUs in executing serial tasks. Hence, a significant fraction of execution time is spent in transferring data from GPU to CPU and vice-versa. Similarly, use of MPI requires reserving the nodes which leads to delays and resource wastage. Also interprocess communication used in MPI is expected to be slower than the inter-thread communication used in multithreading. Further, a common limitation of these platforms is that modern ISOs have simulators with thousands of lines of legacy code and porting them to these platforms will incur significant economic and time overhead. In contrast, MASTER scheduler runs separate to the main application program and provides ability to parallelize. Moreover, MASTER uses multithreading with Java and hence, it can be easily implemented on the commonly available multicore processors.

For contingency analysis, MASTER performs time domain simulation using a high-speed cross-platform power system simulator [START_REF] Khaitan | High performance computing for power system dynamic simulation[END_REF][START_REF] Khaitan | TDPSS: a scalable time domain power system simulator for dynamic security assessment[END_REF]. The simulator uses efficient numerical algorithms (linear, non-linear and integration solvers) to solve the DAEs (differential and algebraic equations) appearing in power-system modeling [START_REF] Khaitan | Numerical methods for on-line power system load flow analysis[END_REF]. The simulator provides different models of power system components such as governors, exciters and generators. The overall contingency analysis steps are explained in Figure 2.

From an architectural perspective, the simulator operates in two modes, namely anticipatory mode and emergency mode. The simulator analyzes most of the critical disturbances in anticipatory mode and storing the results for future use. This helps in providing the power system operator the ability to take timely corrective response at the crucial time after the fault. It provides a cross-platform simulation framework, and has been developed using object-oriented programming (OOP) approach, which facilitates easy development and extension of simulator features.

IV. WORK-STEALING BASED SCHEDULING

For achieving load-balancing on available threads, MAS-TER uses work-stealing algorithm. It has been shown in [START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF], that using work-stealing, the time taken in executing a fullystrict computation 2 with P workers is given by:

T = T 1 /P + O(T ∞ ) (1)
Here T 1 denotes the serial execution time and T ∞ denotes the estimated execution time with infinite number of workers. Also, the storage space required scales only linearly with the number of processors [START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF]. Thus, work-stealing algorithm is efficient in terms of time and space requirement. Workstealing has been shown to impose low overheads and has good scalability [START_REF] Khaitan | Achieving load-balancing in power system parallel contingency analysis using X10 programming language[END_REF][START_REF] Khaitan | PARAGON: An Approach for Parallelization of Power System Contingency Analysis Using Go Programming Language[END_REF][START_REF] Hoffmann | Performance evaluation of task pools based on hardware synchronization[END_REF][START_REF] Khaitan | Parallelizing power system contingency analysis using D programming language[END_REF] and hence, it has been used in several parallel programming models.

However, due to practical considerations, several algorithmic factors influence the efficiency of work-stealing implementation [START_REF] Khaitan | EmPower: An efficient load balancing approach for massive dynamic contingency analysis in power systems[END_REF][START_REF] Khaitan | SCALE: A hybrid MPI and multithreading based work stealing approach for massive contingency analysis in power systems[END_REF]. In the following, we discuss them in more detail and also mention the actual choices made in MASTER.

A. Choice of Algorithmic Options

1. Task dependency and ordering: In the application domains where tasks have interdependence or orderingconstraints (e.g. [START_REF] Min | Hierarchical work stealing on manycore clusters[END_REF]), scheduling algorithm needs to fulfill these constraints for correctness. In our work, each task represents analysis of a single contingency and hence, tasks are mutually-independent, and maintaining ordering is not required. Further, the policy used for enqueuing tasks in the task-queues also becomes unimportant and hence we safely ignore this factor in the design of MASTER.

2. Organization and management of task-queues: The choice of data-structure used for storing task-queues affects the efficiency of task enqueuing and dequeuing. A naïve approach is to keep a single-ended queue where both owner and thief dequeue from the same end. However, this is likely to lead to contention. Another choice is to keep double-ended queue, where the owner dequeues from one end and the thief dequeues from another [START_REF] Dinan | Scalable work stealing[END_REF]. A variant of this is partitioning the per-thread task queue into a private and a public region [START_REF] Min | Hierarchical work stealing on manycore clusters[END_REF]. The thief only accesses the public region of the queue, while the private region is accessed by the owner thread only. This alleviates the need of providing locks in the private region of the queue. However, this also requires separate management of the private and public region and suitably transferring the task from one region to another, when an imbalance is created. In general, double-ended queues (known as deque) are expected to be simpler than the queue with private-public region. For this reason, we have used double-ended queues in MASTER.

Victim selection:

The choice of a victim-selection strategy (also called polling strategy) impacts the communication overhead incurred in each stealing action. Further, when tasks have ordering-constraints, completing the tasks in a desired order may greatly reduce the overall completion time. Some commonly used victim selection policies are:

1) Fixed: The victim is selected in fixed order (e.g.

1, 2, 3, . . . etc.). 2) Random: Any of the available workers is selected randomly as the victim [START_REF] Fluet | A scheduling framework for generalpurpose parallel languages[END_REF]. 3) Local: The thief tries to steal from the nearest-neighbor or one of the close neighbors. Guo et al. [START_REF] Guo | Work-first and help-first scheduling policies for async-finish task parallelism[END_REF] provision stealing within a pre-defined locality domain. Min et al. [START_REF] Min | Hierarchical work stealing on manycore clusters[END_REF] propose a "hardware topology-aware hierarchical victim selection" strategy where a victim is selected in hierarchical manner, i.e., victim processor is first searched in the same socket [START_REF] Matsuba | Single IP address cluster for internet servers[END_REF] and only if no victim is found, the algorithm searches for victims in other sockets or nodes [START_REF] Matsuba | Single IP address cluster for internet servers[END_REF][START_REF] Rabenseifner | Hybrid MPI/OpenMP parallel programming on clusters of multi-core SMP nodes[END_REF]. 4) Circular: The victim is selected in circular manner, i.e., a thief i searches in the order i+1, i+2, . . . , N, 1, . . . , i-1.

Of these, MASTER uses circular polling since unlike fixed polling, it does not cause contention on initial threads, viz. 1, 2,3 etc. Further, unlike random polling, it does not require generating random numbers and is likely to spread the stealing requests more uniformly to the available threads, thus reducing contention. Also, since the tasks do not have ordering constraints and threads share memory, choosing victim in topology-aware manner is not required.

4. Task-stealing granularity: Task-stealing granularity refers to the number of tasks which are allowed to be stolen in each stealing request. While some implementations provision stealing a fixed number of tasks [START_REF] Guo | Work-first and help-first scheduling policies for async-finish task parallelism[END_REF][START_REF] Frigo | The implementation of the Cilk-5 multithreaded language[END_REF][START_REF] Narang | Dynamic distributed scheduling algorithm for state space search[END_REF], other implementations allow stealing a variable number of tasks in each stealing action [START_REF] Min | Hierarchical work stealing on manycore clusters[END_REF][START_REF] Dinan | Scalable work stealing[END_REF][START_REF] Cong | Solving large, irregular graph problems using adaptive work-stealing[END_REF], which can be as large as half the number of tasks available on the victim process (termed as "steal-half"). For sake of simplicity, in this paper, we provision stealing only one task at a time. This enables fine-grain load balancing and also keeps the quantum of work stolen small. For higher number of threads, steal-half policy may provide better load-balancing [START_REF] Dinan | Scalable work stealing[END_REF]. Exploration of this policy is planned as part of future work with larger number of threads.

Proactive task-pushing:

A variant of work-stealing algorithm provisions actively giving tasks to the starving threads, when the amount of local tasks exceeds a given threshold so that the idle-wait time of the starving threads is reduced [START_REF] Wu | Task-pushing: a scalable parallel GC marking algorithm without synchronization operations[END_REF][START_REF] Khaitan | Proactive task scheduling and stealing in master-slave based load balancing for parallel contingency analysis[END_REF]. The disadvantage of task-pushing, however, is that it increases communication overhead and causes congestion in the worst case. For this reason, MASTER does not employ proactive task-pushing.

B. Algorithm Pseudo-code

Algorithm 1 depicts the pseudo-code for MASTER scheduling algorithm. For sake of simplicity, in this paper, we assume that all the tasks have already arrived, although MASTER works equally well for the case when tasks arrive dynamically. At the time of initialization, MASTER creates a number of threads using the Thread class. All the threads carry the same priority. Each object of the Thread class, encapsulates the data and functions (methods) associated with separate threads of execution. Each of these threads has task queues to which it enqueues the newly created tasks and from which it dequeues the tasks for execution. Since the contingencies can be analyzed independently, the tasks can be executed in any order. Each task is run using exec() command, which executes a given command in a separate process. Through this command, the MASTER scheduler interfaces with the power system simulator for consistency analysis. This command returns a Process object, which is used by the scheduler for managing the process created, which includes waiting for the process to finish and handling errors (if any) occurred in the process execution.

The task queues are implemented as deques (i.e. double ended queues, which support element insertion and removal at both ends). The worker executes the task from the front side and stealing requests are responded to from the rear side. For this, Java's BlockingDeque class is used.

A main thread handles initial task distribution and periodically checks whether all the task-queues are empty. For implementation of periodic wait, the sleep() method is used which causes a thread to enter the not runnable state until the specified time has expired. Use of sleep() provides an efficient means of making processor resources available to the other threads. Either for polling or stealing, a worker thread needs exclusive and consistent access to the task-queues of the victim thread, while the victim thread itself needs to dequeue the task for execution. To ensure consistency, MASTER uses the built-in mechanism of Java to provide synchronization. The thread which first accesses the deque gets the lock on it and meanwhile, no other thread can acquire the lock.

To avoid the overhead of creating threads for every task and destroying them at the completion of the task, MASTER uses thread pools. The tasks are passed to the thread pool and when a thread becomes idle, the task is assigned to that thread. If no new task arrives, the idle threads wait for a new task. By reusing threads for multiple tasks, the overhead of threadcreation is amortized over many tasks. Further, use of threadpool enables limiting the number of running threads which helps in avoiding resource thrashing caused by excessive number of threads.

In a multithreaded scheduling algorithm, termination detection requires reading the status of task-queues of all threads. In work-stealing, detecting termination is even more challenging, since threads can steal tasks from other threads, and hence, finishing of the task-queues of a thread does not necessarily imply algorithm termination. To detect termination, MASTER works as follows. The main thread periodically reads the status of all task-queues and when all of queues are empty, it sets allQueuesEmpty to true. When this is true, all threads stop trying to steal more tasks. Each thread i, which has finished its tasks, sets isFinished[i] to true. Meanwhile, the main thread periodically checks whether all the threads have finished and when that is true, the algorithm is terminated.

C. Features and Efficiency

From Amdahl's law [START_REF] Amdahl | Validity of the single processor approach to achieving large scale computing capabilities[END_REF], the speedup obtained from parallelizing a program is limited by the fraction of time spent in the sequential part. For this reason, use of a single (global) locked queue might present a scalability bottleneck since the accesses to the global queue are required to be serialized. Since the MASTER scheduler uses work-stealing with distributed task-queues, it automatically avoids the need of a global taskqueue. Further, since the tasks are allocated from the local task-queue, the latency of access is reduced.

In power systems, contingencies can be analyzed independent of each other and maintaining an ordering between them is not required. MASTER leverages this fact to further speedup its work-stealing implementation, since regardless of whether other tasks have been completed, MASTER can schedule a task to the first free thread. Also, maintaining a table of completed tasks and reading from it to ensure task-ordering is also avoided, which leads to more optimized implementation.

MASTER leverages several features of Java, such as easy memory management, platform-independence and portability. Thus, MASTER can be used on a wide variety of platforms. MASTER has been designed to be thread-safe and thus, the access to shared data (e.g. task-queues) does not cause conflict or race conditions. Also, due to the object-oriented design features of Java, MASTER scheduler can be easily extended.

D. Limitations and Possible Extensions

Multithreading introduces the overhead of thread-switching. When the processor switches from executing one thread to another, the context (e.g. local data and program pointer) of the current thread needs to be saved and the context of next thread needs to be loaded. With increasing number of worker threads, the contention for the shared resources also increases.

A possible improvement in load balancing can be done by provisioning stealing of the tasks from the thread with the largest number of remaining tasks. The overhead in this approach is the need of finding lengths of task-queues of all tasks and comparing them to find the largest sized queue. This extension is planned as a part of our future work.

In case different contingencies have different priorities (e.g., due to their severity or probability of occurrence), MASTER can be easily extended to set thread priorities and assigning tasks with high priorities to the workers (threads) with high priorities. Also, by allocating tasks such that the tasks with higher priorities are executed before those with lower priorities, task-priorities can be handled.

V. RESULTS

Experimental Platform: We simulate a large power system with 13029 buses, 431 generators, 12488 branches and 5950 loads. We show wall-clock time, since it is a metric of relevance to end-user. We simulate 128, 256 and 512 contingencies using 2, 4, 8 and 16 threads and compare the results with serial execution. Each contingency models fault occurrence and clearing in single or multiple system components such as bus and branch. We run the simulations on a multicore processor and different threads run on different cores.

Performance Results: We compare the work-stealing implementation to centralized task-queue scheduling where the main thread holds the queue of ready tasks and all worker threads request task from it on finishing their tasks.

Table I summarizes the time taken in simulation of different contingencies with different number of threads. Here C denotes the number of tasks (contingencies) and N denotes the number of threads. For C contingencies, let T (C, N ) denote the time taken with N threads and T (C, 1) be the time taken using serial execution. To gain insights into the working of different scheduling methods, we define speedup values S(C, N ), for any scheduling method as follows.

S(C, N ) =

T (C, 1)

T (C, N ) (2)
Figure 3 shows these speedup values for different methods. Clearly, MASTER provides highest computational gains and outperforms the conventional techniques. A crucial limitation of centralized task-queue scheduling is that one thread (main thread) becomes busy with the task of allocating the works and performs no useful work.

With higher number of threads, contention increases and resultant synchronization serializes the execution. For this reason, the speedup does not scale linearly with the number of threads. Still, MASTER achieves high computation gains (of the order of 13.5× speedup with 16-threads) and much better load balancing than conventional scheduling techniques.

VI. CONCLUSION

In this paper, we presented MASTER, a multithreaded work-stealing technique for efficient parallelization and loadbalancing of contingency analysis task in power systems.

Fig. 3 :

 3 Fig. 3: Speedup Results with MASTER and conventional scheduling techniques for different number of worker threads and contingencies

Initialization 1. System Model 2. Dynamic and steady state data Initial state Final state Input DAE Kernel 1. Integrator 2. Nonlinear Solver 3. Linear Solver 1. Stability Analysis 2. Failure Detection Output Time domain Simulation

	MASTER Scheduler
	Code
	Java compiler
	Compiled Code
	in Bytecode
	Main Thread			Java Virtual Machine (JVM)
	(started by
	JVM)
	Initialization
	Threads started	T 1	T 2	T N	Each thread performs its task
		T 1	T 2	T N
	Time
	Finalization	T 1	T 2	T N
	Program End
	Fig. 1: Implementation of multithreading in MASTER using Java
	Fig. 2: Single Contingency Analysis

TABLE I :

 I Simulation time (seconds). N=1 refers to serial execution, Centr. = centralized task-queue scheduling

			C=128		C=256		C=512
	N	Centr.	MASTER	Centr.	MASTER	Cent.	MASTER
	1	2657	2657	5328	5328	10780	10780
	2	2670	1334	5321	2698	10682	5353
	4	894	672	1813	1367	3575	2703
	8	413	358	817	717	1610	1404
	16	228	199	436	393	893	790

In this paper, we use 'load-balancing' as it is used in parallel processing field (distribution of tasks to threads in a manner that the completion time of different workers is as close as possible) and not as it is used in powersystems.

Fully-strict computation is one where the data dependencies of a thread go to its parent only and not to any other thread. Since contingencies are independent and thus, there are no data dependencies. Hence, the theoretical guarantee provided in[START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF] holds for our experiments also.

We discussed several design-time choices made for efficient implementation. The experiments show that MASTER offers large computation gains.