Performance Benchmark of Modelica Time-Domain Power System Automated Simulations using Python

In this paper, a benchmark between solvers and Mod-elica tools for time-domain simulations of a power sys-tem model is presented. A Python-based approach is employed to automate Modelica simulations and compute performance metrics. This routine is employed to compare the performance of a commercial (Dymola) against an open-source (OpenModelica) simulation tool with different solver settings. Python scripts are developed to execute a dynamic simulation of a common model for power system studies with 49 states and 420 variables in three different scenarios. This degree of automation makes it easier to change solver settings and tools during execution. The performance of each of the tools is assessed through metrics such as execution time and CPU utilization. The quantitative comparison results provide a clear reference to the performance of the tools and solvers for the execution of time-domain simulations with a signiﬁcant degree of complexity. The commercial tool offers better performance for variable-step solver, but the performance of the open-source software shows signiﬁcantly faster results for ﬁxed-step solvers.


List of Acronyms and Definitions
Definitions ST · Simulation Time: the simulation time is understood as the time it takes for the program to translate, compile and integrate a model.
ET · Execution Time: the execution time is the elapsed time to complete the numerical integration of the compiled model. It also known as integration time. Note that execution time is included as a part of Simulation Time in the context of this paper.
NMT · Normalized Minimum Execution Time: performance metric taken as a function of the execution time of each of the tools. It is defined as where ET D are the execution times for Dymola and OM, respectively, and ET observed is the corresponding integration time of each tool for a given solver obtained from the simulation log.

Acronyms
ET · Execution Time MSE · Mean Square Error NAE · Normalized Absolute Error OM · OpenModelica OMPython · OpenModelica Python Interface OpenIPSL · Open-Instance Power System Library PDI · Python-Dymola Interface ST · Simulation Time

Motivation
Modeling and simulation of power systems have been a habitual practice in the energy industry since the 1960s. The complexity of a power system is steadily increasing to accommodate modern technologies into the existing grid. A more complex system leads to more elaborated models. High-complexity models are directly correlated with computationally expensive tasks (Milano, 2010). In this context, the Modelica language represents an accurate, equation-based, multi-domain solution modeling and simulation alternative. Numerous initiatives such as OpenIPSL have been taken to incorporate into the power system workflow the benefits of the Modelica language (Baudette et al., 2018). On the other hand, the academic, scientific and industrial communities have come to acknowledge the intrinsic benefits of the Modelica language. An outcome of this trend is that the user base has increased significantly during the last years. This has led to the development of many libraries with users coming from a very wide domain spectrum. Nowadays, Modelica stakeholders include students, consulting firms, big laboratories, and industry agents.
Free tools such as OpenModelica are fundamental for learning the language at little to no cost and to set a reference for the Modelica language (Fritzson et al., 2006).
Commercial tools such as Dymola, SystemModeler or SimulationX provide advanced functionalities that satisfy particular requirements from the industry. However, there is no clear guidance for a user on how to select the toolsolver based on its simulation performance exclusively.
This paper intends to provide this guidance. It aims to compare the time-domain simulation performance of the solvers from both Dymola and OpenModelica when subjected to different solver settings (see (Braun et al., 2017) for a detailed analysis of the potential of OpenModelica to solve large-scale models). Since these tools do not have the same features and solvers, we have chosen some of the ones they have in common for benchmarking purposes.

Contribution
This work is relevant to any user of the Modelica language. The tool performance analysis is based on the simulation of a power system model (IEEE 14 bus system), that serves as a representation of a dynamic nonlinear system. We consider three simulation scenarios: an initialization, a line-opening (one discrete event) and two bus faults (two discrete events). The paper reveals the difference in performance within the tools and helps users make an educated choice about the tools to use. The main contributions of the paper are the following: • Quantitative evaluation of Dymola and OpenModelica simulation performance for time-domain simulation of complex dynamic systems (power systems). • Benchmark of different solvers in a dynamic simulation with discrete events. • Implementation of simple Python routines to automate Dymola and OpenModelica time-domain simulations.

Paper Organization
The paper is broken down in the following sections: Section 2 describes the test system and the Modelica library employed to construct it. The experiment setup regarding hardware characteristics and software setup is described in Section 3. In Sections 4 and 5, we discuss performance results of each of the tools with respect to each solver and the corresponding performance metrics. Finally, Section 6 concludes the work.

Modelica Power System Model
The IEEE14 bus system 1 represents a part of the Midwestern USA American Electric Power System as of February of 1962. The single-line diagram of the system can be seen in Figure 1. This model was chosen because it is a widely used testing system for an initial assessment in power system dynamical studies since it has a significant number of variables and states (420 and 49, respectively) which makes it a common factor in such simulation-based 1 https://icseg.iti.illinois.edu/ieee-14-bus-system/ studies (Milano, 2010). For this reason, its dynamic simulation a challenge to the tools and the CPU.

Experiment Specifications
To make sure that the results are reproducible, this section details the conditions under which the experiments were performed regarding hardware setup and software characteristics.

Hardware and Software Setup
The characteristics of the computer used to run the simulations are shown in Table 1 To assess solver performance correctly, numerical integration must run in only one processor. While this is a default option in Dymola, we need to specify this option explicitly in OpenModelica before starting any simulation since it defaults to multicore execution. This is done thanks to the flag setCommandLineOptions("-n=1").

Simulation Scenarios
To properly measure solver performance in diverse dynamic conditions, we will consider the following three scenarios of the IEEE14 bus model: system initialization, time-domain simulation with one line opening, and system response with two faults.

IEEE14 System Initialization (S 1 )
This scenario corresponds to a system with no disturbing events. The power flow of the model is modified so that the numerical simulation has problems during initialization. The provided initial conditions are such that the dynamic system is not initially at an equilibrium point, thus forcing the system to look for an acceptable steady-state condition at the beginning of the integration process. This increases the computational task and challenges the solver since the integration does not start with all state derivatives equal to zero.

Line Opening (S 2 )
Besides the aforementioned bad initialization condition, we introduce a line opening to disturb the system from steady-state and excite nonlinear dynamics. This kind of scenario is used to study system-wide stability when two sub-areas are disconnected from each other. The line opening corresponds to the connection between buses 2 and 4 (B2 and B4). The line will open from both ends at time t = 60 s and will re close at t = 61.5 s.

Bus Faults (S 3 )
In this case, the system will face two three-phase to ground faults at different times. This configuration is used to test the resiliency and stability of the system. By having two faults, the numerical complexity of the simulation increases, creating a more adverse scenario for the solvers to come up with a solution. Fault 1 occurs at bus 4 (B4) starting at t = 20 s and being removed at t = 21.2 s. Fault

Solver Selection
The performance of the time-domain simulation depends not only on the dynamic condition to be analyzed but also on the solver selection. In this regard, OpenModelica and Dymola contain a wide variety of different integration methods and three of them are going to be used and thus briefly described in this study. The Differential Algebraic System Solver (dassl) is an implicit, high-order, variablestep solver with time-step control. This solver is set as default solver in both OpenModelica and Dymola. The Euler method is another solver available in both software packages and it is an explicit (Forward Euler), first-order, fixed time-step solver. Finally, the last solver used in this study is the runge kutta. Dymola allows the user to chose between second, third and fourth order Runge-Kutta methods but in this work, only the fourth order is used since it is also available in OpenModelica. This solver is an explicit, fourth-order, fixed time-step solver. This paper will benchmark the performance of the tools with each of the mentioned solvers for the different scenarios of the test power system.

Time-step Selection
Since dassl is a variable-step solver with step-size control, there is no need to select a specific time-step for the solver. The selection of an adequate number of intervals is necessary to plot and analyze the results. For both tools, 5000 was found to be a reasonable number of simulation intervals. Moreover, to use the capabilities of a DAE solver to their full extent, we enable the newly incorporated DAEmode in Dymola by enabling the flag Advanced.Define.DAEsolver = true (Henningsson et al., 2019). In OpenModelica, to set similar settings we use the command setCommandLineOptions("daeMode=true").
On the other hand, it is important to select an adequate step size T s for fixed-step solvers in order to guarantee that the algorithm is operating in its region of convergence. To get an upper bound for T s , we performed a linear analysis of the system in Dymola employing the library Modelica_LinearSystems2. After determining the time constant of the fastest mode (τ ≈ 1 ms), we found that T s = 0.5 ms was a reasonable value to capture the effects of the fastest mode, guaranteeing numerical convergence for both solvers, Euler and Runge Kutta. The selected time-step size implies that 240,000 simulation intervals are going to be needed for a simulation time of 120 s.

Benchmark Metrics
In order to understand and accurately compare the two tools the paper focuses on two simulation features to compare: • Simulation Time (ST) corresponds to the time it takes for a program to complete all of the routines for each scenario comprising model translation, compilation and execution. The discussion of the results of the simulation time are found in Section 4.1, with special remark on Execution Time (ET). • CPU Utilization is the percentage of central processing unit (CPU) that is being used at any time during the execution. Results for CPU utilization can be found in Section 4.2.

Code Structure
The complete code to perform the experiments and analyze the resulting data can be found in GitHub 3 . The execution of the simulations is automated through Python using the Dymola API (Python-Dymola Interface) and the OpenModelica Python Interface (OMPython) (Lie et al., 2018). The details of the Dymola routine can be seen in the file dymola_simulation.py. Likewise, the OpenModelica commands are included in the file om_simulation.py.
To measure performance we execute the routine in the script measurement_performance.py. It measures each of the performance metrics every 0.2 s while the code is running in a different parallel process.
The main program is contained in the file 01_modelica_tool_performance_benchmark.py.

Performance Results
Before presenting the performance results, we validate the simulation outputs of the three scenarios for Dymola and OpenModelica for all solvers. We employed the Normalized Absolute Error (NAE) and the Mean Square Error (MSE) defined in Equation (1) to quantify the numerical difference between the outcomes of each tool.
NAE shows how different the Dymola and OpenModelica results are throughout the simulation. MSE outputs a quantitative validation of the results of both tools (Devore and Berk, 2012). Full details can be seen in Table 3.
The numerical behavior of the simulation during initialization (runge kutta solver) can be observed in Figure 3 for the voltage magnitude signal at Buses 2 and 4. An initial transient behavior can be seen at the beginning of the integration time. This is not desired in a dynamic simulation since numerical convergence to a steady-state solution is not guaranteed given the fact that the solver starts from a guessing point with non-zero derivatives. The non-steady state behavior at the on-set of the simulation is due to the fact the initial guess used in the model (the so-called power flow) is not close enough to an equilibrium for the initialization routine to solve for a more precise set of initial values. A more complex initialization problem will better benchmark the capabilities of the tools. Despite this, Dymola and OM produce almost the same results, with an NAE in the order of 10 −3 .
Likewise, for the runge kutta solver, Figures 4 and 5 show the simulation results for the line opening (voltage magnitude at buses 2 and 4) and the double bus fault (voltage at affected buses 4 and 14) scenarios, respectively. Both Figures reveal how there is a minimal error between the results of both tools. Based on these results, it is concluded that fixed-step solvers can be applied to reduce discrepancies between different Modelica tools.  The complete collection of plots for all solvers and simulation scenarios can be found online in the GitHub repository in the Notebook 02_Data_PostProcessing_ SimulationResultPlotting.ipynb.

Simulation Time
The information regarding simulation time is presented for all scenarios and solvers in Table 2. We must underline that simulation time includes compilation, translation and actual integration (execution time).
A clear conclusion from this information is that the variable-step solver is the most convenient for an initial analysis of the conditions of the system with an important amount of detail. Nevertheless, considering the information about MSE, a fixed-step solver shows advantages to reduce the numerical discrepancy between tools running the same model. The cost is a considerable increase in simulation time.

CPU Utilization
Since each instance of Dymola/OpenModelica was constrained to run only on one core, we expect exactly one processor to be responsible for numerical integration while a simulation is being carried out. The CPU usage of the assigned execution core is 100% due to the heavy numerical task of the simulation.  An interesting outcome of our experiments is that several CPUs are involved in the execution process but just one is performing the simulation tasks at a given time. We can detail this behavior in Figure 6 for a Dymola simulation using the runge kutta method of the bus fault sce-nario. Simulation starts in Core 1 where the CPU usage is at a 100% at the beginning of the running time. Afterward, it is delegated to Core 5. Finally, Core 10 completes the execution of the program. This behavior is due to a task scheduling routine in the processor level that dispatches to different cores the compilation, translation, and integration sub tasks.
Similar behavior happens with another solver and OpenModelica (Figure 7) in which the simulation started in Core 2, then was briefly assigned to Core 5 and was finished in Core 6. All the graphics can be detailed in the GitHub repository inside the Jupyter Notebook called 03_DataPostprocessing_CPU_Usage.ipynb.

Performance Evaluation Metrics
A score was proposed to quantify the performance differences between the tools and the solvers. The score is obtained from the data generated for all simulations and solvers. This single metric makes it simpler to directly compare the performance of Dymola versus OpenModelica. From Table 2 the Execution Time (ET) for each scenario and solver were employed. These metrics were obtained directly from the program logs and measured in Python. Notice that the time registered using OMPython is slightly larger than the reported by the simulation log due to the communication interface between Python and OM. The translation and compilation time were not taken into account since this information is only available in the Dymola developer version, not in the release version.
The Normalized Minimum Execution Time score (NMT) of each scenario per solver is computed as where ET observed is the ET for a particular solver in Dymola or OpenModelica), and min(ET D , ET OM ) is the minimum execution time between both tools for a specific solver. Clearly, NMT [solver] lies between 0 and 1. The higher the NMT is, the faster the simulation will run for a particular selected solver. At a first glance, this metric might be counter-intuitive since a better solver/tool combination would reduce execution time. However, we propose an increasing score metric due to the fact that users are more familiar to higher scores for better performance. Therefore, the larger the NMT is, the faster a particular solver will run. The NMT metric results are presented in Table 4. The performance of Dymola is remarkably better using dassl. Nevertheless, OM shows a smaller execution time than Dymola for fixed-step solvers (as can be seen from Table  2, the NMT scores and 5). This conclusion can be further detailed in Table 5 where a direct comparison between the execution time for the tools with the different solvers for each scenario is presented.
The NMT scores highlight that the performance of Dymola in terms of execution time is remarkably better for variable-step solvers. The relative advantage of selecting one tool with respect to the other can be computed from the NMT directly. For instance, Dymola runs 47.3x faster than OM for the first scenario using dassl which can be computed by a direct comparison of the ET listed in Table 5. The NMT score of OM for S 1 is 0.0211 which is 1/0.0211 = 47.3 times smaller than the corresponding Dymola metric reflecting the relative difference in execution time.
For the variable-step solver, the discrepancy between the tools can be attributed to the performance of the dassl solver in all simulations thanks to the aforementioned improvements for DAEMode inside Dymola (Henningsson et al., 2019).
The execution time of OM is faster than the one of Dymola for all scenarios when a fixed-step solver is used. The NMT scores show a relative advantage between 3.4 and 6.8 times favoring OM. We have contacted Dassault Systèmes about the performance of the simulations of the IEEE 14 Bus System using runge kutta methods (including euler) as integrator and GCC for compilation. Dassault reports that bug fixes have been made in the GCC runtime libraries, leading to CPU times are about 3 -4 times faster, on par with the run times given when compiling with Visual Studio under Windows 10. Dassault informs that the updated libraries will be part of Dymola 2021.
We should point out that the scope on the potential optimization features has been limited to the use of the flag Evaluate = true in Dymola and -d=evaluateAllParameters in OM, which is standard practice when attempting to improve simulation performance.
The detailed step-by-step computations of the scores can be found in GitHub in the Notebook 05_BenchmarkMetrics.ipynb.

Conclusions and Future Work
The paper presents a concrete analysis of the time-domain simulation performance of Modelica-based tools for different solvers in the context of large-scale nonlinear dynamic systems. The presented results can help a user to choose a tool depending on the final application, and lead to improvements in Modelica tools. The methodology of this benchmark can be extended to virtually any platform or Modelica tool.
We benchmarked the time-domain simulation performance of two popular Modelica tools, Dymola and Open-Modelica, for a dynamic power system simulation using the IEEE 14 bus system. We considered several scenarios that challenge numerical solvers differently. Thanks to Python scripting, we were able to change automatically the simulation settings while directly measuring the performance of the computer instead of relying on simulation logs. Python functions also made it quicker to analyze straightforwardly the big set of data regarding simulation results and computer performance.
For the proposed heuristic score, we found out that OpenModelica performs better than Dymola in terms of execution time for fixed-step solvers while Dymola shows faster results when using a variable-step solver (see Table 5). Despite this, we must warn the reader that this conclusion is based upon only a particular system. Further research has to be done to include more test systems. Moreover, the use of a fixed-step solver has the main advantage of The tool and solver benchmark results are expected to be reproduced in a larger system such as the Nordic 44 with ≈ 1300 states and 6300 variables. This system requires a considerable amount of RAM given its large number of states. Therefore, future work is related to the performance analysis in a 64-core machine with 512 GB of RAM for the N44 system considering different types of simulations and various initialization parameters.
In this appendix, all performance results of the different simulation experiments are presented. In Table 5