An Open Source Power System Simulator in Python for Efficient Prototyping of WAMPAC Applications
AAn Open Source Power System Simulator in Pythonfor Efficient Prototyping of WAMPAC Applications
Hallvar Haugdal, Kjetil Uhlen
Department of Electrical Power EngineeringNorwegian University of Science and TechnologyTrondheim, [email protected], [email protected]
Abstract —An open source software package for performingdynamic RMS simulation of small to medium-sized power sys-tems is presented, written entirely in the Python programminglanguage. The main objective is to facilitate fast prototyping ofnew wide area monitoring, control and protection applicationsfor the future power system by enabling seamless integration withother tools available for Python in the open source community,e.g. for signal processing, artificial intelligence, communicationprotocols etc. The focus is thus transparency and expandabilityrather than computational efficiency and performance.The main purpose of this paper, besides presenting the codeand some results, is to share interesting experiences with thepower system community, and thus stimulate wider use andfurther development. Two interesting conclusions at the currentstage of development are as follows:First, the simulation code is fast enough to emulate real-timesimulation for small and medium-size grids with a time step of 5ms, and allows for interactive feedback from the user during thesimulation. Second, the simulation code can be uploaded to anonline Python interpreter, edited, run and shared with anyonewith a compatible internet browser. Based on this, we believe thatthe presented simulation code could be a valuable tool, both forresearchers in early stages of prototyping real-time applications,and in the educational setting, for students developing intuitionfor concepts and phenomena through real-time interaction witha running power system model.
Index Terms —Dynamic RMS simulation, wide area monitor-ing, control and protection, real-time simulation, Python
I. I
NTRODUCTION
Python ranks high among the worlds most popular program-ming languages. Although very slow compared to many otheralternatives (e.g. C++), it has the advantage of being easyto learn, read and write. Development of new applications istherefore generally considered to be much easier and faster inPython. Furthermore, a rich library of open source packagesis available, making it straightforward to make use of artificialintelligence, communication protocols, powerful visualizationsand more.In this paper, we present an open source package for per-forming dynamic RMS-simulations, written entirely in Python.With this, we seek to contribute with a highly transparent,easily modifiable and expandable power system simulator thatallows tight integration with powerful tools that already existin the open source community. We believe that open solutionsare essential for developing the wide area monitoring, controland protection applications of the future power system. The initial motivation for development of such a packagearose during research towards a controller for damping ofpower oscillations. The particular study required a non-lineardynamic power system simulator with detailed generator mod-els, capable of incorporating signal processing techniques likeKalman Filters in the simulation loop, as well as the abilityto perform eigenvalue analysis on linearizations of the model.This was found to be surprisingly difficult.Using Python as a scripting language is standard func-tionality in many commercial tools, for instance DIgSILENTPowerFactory, and is well described in the documentationof the software. Including Python in the dynamic simulationloop is much more difficult, and is not standard functionality.Interfacing Python with dynamic simulations in PowerFactorycan be achieved with the repository described in [1]. Thisfunctions by making PowerFactory call a dll extension anumber of times during each simulation step. The dll file,compiled from C code, calls Python. This potentially expandsthe functionality of PowerFactory a great deal, but might tosome appear complicated to use, for instance requiring someknowledge on compilation C code.PSS/E is very well interfaced with Python through thepsspy module. However, small signal stability analysis requiresseparate modules to be installed. These modules are not thatwell interfaced with Python, and potentially complicates theworkflow. Also, the full, linearized system matrix, which mightbe required in some applications, is not directly available.Among open source alternatives for dynamic power systemsimulations, we have DPsim [2], which is a simulator writtenin highly efficient C++ code, specifically developed for real-time- and co-simulation. ANDES [3] is a Python software forsymbolic modelling and numerical analysis of power systems.The Open-Instance Power System Library (OpenIPSL) [4],[5], written in the Modellica language, specifically targetsunambiguous model sharing among utilities and researchers.The above mentioned tools are powerful enough to simulatesystems with thousands of buses. This is necessary for thesetools to be able to perform simulations on detailed models ofreal transmission grids. However, for research or educationalpurposes, smaller aggregated systems with tens to hundredsof buses are often used. Frequently used test systems includethe Kundur Two-Area System, the IEEE 9, 14, 39 and 68-bussystems and the Nordic 32 and 44-test systems. Systems of a r X i v : . [ ee ss . S Y ] J a n his size can be simulated entirely in Python, without strictlyrequiring sophisticated performance-boosting techniques thatin some cases compromise transparency and expandability.Having developed the simulation code, the encouragingexperience obtained this far regarding performance and easeof use motivates us to present the code in this paper.The entire code is developed in Python, where the corefunctionality relies only on standard packages like NumPy [6]and SciPy [7]. This comes with the advantage that the softwareis cross platform-compatible and easily deployable, and doesnot require compilation. Further, the software can easily beuploaded to an online Python interpreter and shared, editedand run in the cloud, which could represent a major benefitfor reproducibilty of research.The code is available in a repository called DynPSSimPy[8] (Dynamic Power System Simulator in Python) on GitHub,released under the GPLv3 license. By following links pro-vided on the GitHub-page, simple examples can be launchedonline (using Binder [9]) to demonstrate some of the basicfunctionality of the software.Finally, it should be mentioned that the package is inits very early stages of development at this point, with alimited number of dynamic models and somewhat limitedfunctionality in some areas.The rest of this paper is structured as follows: SectionII goes into detail on the core functionality; in Section IIIsimulation results are validated against commercial simula-tion software; in Section IV we elaborate on the real-timesimulation functionality; finally, Sections V and VI containdiscussion and conclusions.II. C ORE FUNCTIONALITY
In dynamic analysis of large scale power systems, theac electrical variables (voltages and currents) are usuallyrepresented as phasors when the focus is on electromechanicaldynamics and primary control. The relevant models appliedare referred to as RMS-models. The fundamental problemthat needs to be solved by the RMS-model simulator is theintegration of the differential equations of all the dynamicmodels in the system, while at the same time making sure thatthe algebraic equations representing the network are satisfied.
A. Dynamic RMS-simulation
The Differential Algebraic Equations (DAE) describing thedynamics of the system can be written on the form ˙ x = f ( x , y )0 = g ( x , y ) (1)The function f describes the differential equations for allstates, while g describes the algebraic equations representingthe network equations. x is the vector of state variables, while y is the vector of algebraic variables.With the dynamic models implemented this far, the alge-braic variables are constituted solely by the bus voltages of the simulated system. The set of algebraic equations is linear,and can be written on the form YV = I inj ( x ) (2)where Y is the admittance matrix, I inj ( x ) is the vector ofcurrent injections and V = y the vector of bus voltages.Solving the algebraic equations is very efficient with thesystem sizes considered this far (up to 45 buses), allowingus to convert the DAE system into a system of OrdinaryDifferential Equations (ODE) without significant deteriorationof performance.First, we rewrite the algebraic equations as y = V = Y − I inj ( x ) = h ( x ) (3)(In practice, of course, the admittance matrix is not inverted foreach time step, but rather it is solved by an efficient algorithmfor sparse systems of linear equations, for instance using thefunction scipy.sparse.linalg.spsolve.)Further, we eliminate the algebraic variables: ˙ x = f ( x , h ( x )) = f ( x ) (4)This system of ODEs can now be integrated with anysuitable integration method, for instance the Runge Kutta-implementation found in scipy.integrate.RK45 [7].For real-time simulation, it makes sense to use integrationmethods with fewer evaluations of the ODE function per steptaken by the solver. This allows the step size to be decreased,which in turn results in more accurate representation ofcontinuously changing inputs. The Euler- or Modified Eulermethods, requiring one or more evaluations, are thus betteralternatives than the fourth order Runge Kutta-method, whichrequires four evaluations per time step. B. Initialization of dynamic simulation
To make sure that the system stays on the initial operatingpoint after the simulation started (assuming no disturbancesare applied), the system must be initialized such that alltime-derivatives are equal to zero. Essentially, this meansdetermining the state vector x that gives f ( x ) (5)Initialization of the simulation can be summarized as fol-lows:1) Perform a Newton-Rhapson power flow calculation todetermine the bus voltages and active and reactive powerinjections at each bus.2) From the power flow result, determine the terminalvoltage and the active and reactive power produced byeach generator.3) Perform Kron-reduction of the dynamic admittance ma-trix. All other buses than the generator buses can beeliminated. However, other buses that should not beeliminated can also be specified. For instance, lineoutages require less computation effort if the buses ofhe line in question are present in the reduced system,which might be interesting when performing real-timesimulation.4) The number of states is determined, and the ordering ofthe state vector is determined. According to the dynamicmodel definitions in the model library and the numberof instances of each model, the appropriate number ofstates is added. (For instance, for 10 generators with 6states each, 60 states is added in the state vector.)5) The state vector that makes all the differential equationsequal to zero is found. This is done by solving thealgebraic equations resulting from setting the differentialequations equal to zero. C. Modal Analysis
Modal analysis is performed by numerical linearization ofthe system of ODEs. We find the linearized system matrix asthe Jacobian of the ODE function, as follows: A = ∂f ( x ) ∂ x (6)Finally, we perform an eigendecomposition of the systemmatrix to get the eigenvalues and eigenvectors of the sys-tem at this operating point (for instance using the functionnumpy.linalg.eig). D. Dynamic models
At present, one generator model, one AVR, one GOV andone PSS model is implemented. Generator controls can be ap-plied on arbitrary units with varying parameters. Furthermore,the software is written such that implementation of additionalmodels is possible.The formulation of the differential equations representingthe synchronous machine is based on the formulation in [10]: M ∆ ˙ ω = P m − P e ˙ δ =∆ ωT (cid:48) d ˙ E (cid:48) q = E f − E (cid:48) q − I d ( X d − X (cid:48) d ) T (cid:48) q ˙ E (cid:48) d = − E (cid:48) d + I q ( X q − X (cid:48) q ) T (cid:48)(cid:48) d ˙ E (cid:48)(cid:48) q = E (cid:48) q − E (cid:48)(cid:48) q − I d ( X (cid:48) d − X (cid:48)(cid:48) d ) T (cid:48)(cid:48) q ˙ E (cid:48)(cid:48) d = E (cid:48) d − E (cid:48)(cid:48) d + I q ( X (cid:48) q − X (cid:48)(cid:48) q ) (7)Subtransient saliency is neglected, i.e. we require X (cid:48)(cid:48) d = X (cid:48)(cid:48) q . Interfacing each generator with the grid is achieved usingthe following equations: (cid:20) E (cid:48)(cid:48) d E (cid:48)(cid:48) q (cid:21) = (cid:20) V d V q (cid:21) + (cid:20) R X (cid:48)(cid:48) q X (cid:48)(cid:48) q R (cid:21) (cid:20) I d I q (cid:21) (8)Using that X (cid:48)(cid:48) d = X (cid:48)(cid:48) q = X (cid:48)(cid:48) , this can be simplified to E (cid:48)(cid:48) = E (cid:48)(cid:48) d + jE (cid:48)(cid:48) q = ( V d + jV q ) + ( R + jX (cid:48)(cid:48) )( I d + jI q )= V + ( R + jX (cid:48)(cid:48) ) I (9) Further, we establish the Northon Equivalent of the syn-chronous machine, where the current source is given by I no = E (cid:48)(cid:48) /X (cid:48)(cid:48) (10)and the shunt impedance by Z no = R + jX (cid:48)(cid:48) (11)The impedances Z no of each of the synchronous machinesare included in the admittance matrix used during dynamicsimulation, Y in (2), appearing as contributions on the diag-onal entries. The currents I no appear as contributions to thecurrent injection vector I inj ( x ) .Finally, the electrical power is given by P e = E (cid:48)(cid:48) d I d + E (cid:48)(cid:48) q I q (12)For generator controls, the standard models SEXS forAVR, TGOV1 for turbine/governor and STAB1 for PSS areimplemented. E. Implementation
Each dynamic model is defined as a class in Python. Withineach class, the function for calculating the state derivatives ofthe dynamic model is defined. Functions for initialization andcomputing current injections are also implemented for modelswhere this is required.Evaluating the complete state derivative vector in (4) isachieved by looping through the state derivative functionsof all the individual dynamic models included in the powersystem model.For the implemented models described above, the statederivative functions make use of vectorization, which sig-nificantly improves the computational efficiency. A furtherincrease in computational efficiency can be achieved by lever-aging Just-In-Time compilation of the derivative functionsof each model, available in Python using Numba [11]. Thisintroduces some overhead, but speeds up the state derivativefunctions significantly once the compilation is done.III. V
ALIDATION
The simulator is validated by comparing results fromthe Kundur Two-Area System, obtained using the presentedPython package, with results obtained when simulating amodel with identical parameters in DigSILENT PowerFac-tory. The ”Standard/Detailed model 2.2” synchronous machinemodel is chosen for all synchronous machines in PowerFac-tory. This is also a sixth order model, but the differentialequations are formulated using fluxes instead of voltages. Sat-uration is neglected, and all the loads are constant impedanceloads.Figure 1 shows results from simulating a short circuit onGenerator 1 in the Kundur Two-Area System, where the resultsfrom the Python simulation are shown alongside referenceresults from PowerFactory. The angle, speed and subtransientd- and q-axis voltages are shown. The results are very similarin this case, substantiating the validity of the model. . . . . . A n g l e [ r a d ] . . . . Sp ee d [ pu ] . . . . . Sub t r a n s i e n t q - a x i s v o l t ag e [ pu ] . . . Sub t r a n s i e n t d - a x i s v o l t ag e [ pu ] PythonPowerFactory G1 G2 G3 G4
Fig. 1. The results from simulating the Kundur Two-Area System in Pythonare shown alongside corresponding results from a model with identical param-eters simulated in PowerFactory. The transient d- and q-axis voltages are notshown, since these are not directly available from the PowerFactory simulation(since in PowerFactory the dynamic generator equations are formulated interms of fluxes instead of voltages). The results are not identical, but verysimilar.
IV. R
EAL -T IME S IMULATION
For smaller grids, the simulator is fast enough to performreal-time simulation with a time step of around 5 ms. Toperform real-time simulation, a timer is started at the beginningof the simulation, which is used to synchronize the simulationtime with the wall clock time. After each step taken by thesolver, the simulator is set on hold until the next time step isdue. Although Python is by no means a language tailored forreal-time simulation, experience shows that the impression ofcontinuous real-time interaction with the simulated model can T i m ee l a p s e dp r . s t e p [ m s ] Ideal (5 ms)ActualCalc. only
Fig. 2. The figure indicates how the synchronization between exact time andsimulation time is performed during real-time simulation of the IEEE 39-bussystem, with a total of 123 states. The simulation step in this case is 5 ms(red, dashed line). The actual loop time varies somewhat around the ideal looptime (blue line). It can also be seen that the time spent only on the calculationof each step (green line) is on average around 1 to 2 ms, i.e. well below 5ms. be achieved.Fig. 2 shows the amount of time spent on each time stepby the solver during real-time simulation of the IEEE 39-bus system with a time step of 5 ms. All generators aremodelled using the sixth order model, and all generatorsexcept the largest synchronous machine (which represents therest of USA and Canada) are equipped with the AVR, PSSand turbine-governor models presented in section II-D. Themodel has a total of 123 states in this case, and is simulatedusing the Modified Euler method, with constant step sizeand one correction iteration. Numba was used for Just-In-Time compilation of the state derivative functions to increasecomputational efficiency. The simulation was performed on alaptop with an Intel i7 CPU with two 2.8 GHz-cores and 16GB RAM.It can be seen that the calculation time (shown in green) onaverage is around 1 to 2 ms, which is fast enough for a timestep of 5 ms (shown in red). The calculation time surpasses5 ms in some occasions. Also, the actual time spent on eachtime step (shown in blue) varies approximately ± ig. 3. The GUI during real-time simulation of the Kundur Two-Area System is shown. The GUI is composed of a phasor plot displaying the rotor angle andexcitation voltage magnitude ( E f ∠ δ ) for each generator (upper left), and curves for generator power angles and speed (right). The controls (lower left) allowsthe user to interact with the simulation through a console (allowing commands to be issued to modify the admittance matrix, generator control parametersetc.), activate/deactivate AVR, governor or PSS of any of the generators, or disconnect or connect one or more of the eight lines. If the AVR is deactivated,the excitation of the individual generators can be controlled manually by sliders. GUIs to be implemented easily, with customized live plots andcontrols. A screen shot of the GUI is shown in Fig. 3 duringreal-time simulation of the Kundur Two-Area system. TheGUI in this case shows live phasors representing the excitationvoltage and rotor angle for each generator, as well as speedand angle plots. Additionally, the controls in the lower left partof the figure allows lines to be disconnected, and generatorcontrols to be activated/deactivated.V. D
ISCUSSION
Based on the presented results, we elaborate on the follow-ing ideas:
A. Desktop Real-Time Simulator
It is found that real-time simulations of small to mediumsized systems can be performed on an average laptop com-puter, without requiring specialized operative systems orhardware. For researchers, this potentially simplifies testingreal-time capability and performing live demonstrations ofWAMPAC applications significantly. In the educational set-ting, interactive real-time simulation and visualization allowsstudents to get hands on experience with a running simulationmodel, which is valuable for developing intuitive understand-ing of different stability problems, control systems etc.
B. Reproducibilty of research
For researchers, it is often more time consuming thannecessary to reproduce results from other simulation studies.Even small power system models have hundreds of parametersthat need to be correct for researchers to be able to accuratelyreproduce results. Furthermore, the specific implementation ofdynamic models might differ across simulation tools. If theentire code for reproducing results for a particular study couldbe made available for the community, complete with the powersystem model, dynamic models and the simulation code, thiscould significantly reduce the time spent on establishing thereference case.The presented simulation code being developed entirely inPython (relying only on standard packages for the core func-tionality) comes with the advantage that it is straightforwardto upload the code to an online Python interpreter. This makesit easy to share, edit and run simulations in the cloud. Thiscould make reproducing results even easier, as the simulationrequired to produce the published result could be performedin the internet browser.VI. C
ONCLUSION
Dynamic RMS simulations of small to medium sized gridscan be performed efficiently in Python using the presentedsimulation software. This enables tight integration of powerystem simulations with other open source packages availablefor Python, which is the main objective. Although transparencyand expandability is prioritized over computational efficiencyand performance, it is found that real-time simulation witha reasonably short time step is possible. The possibility ofuploading the entire code required for reproducing results froma publication is also emphasized.Finally, adding that the package is free, easily deployableand cross-platform compatible, we believe that it could bea valuable tool, both for researchers prototyping WAMPACapplications, and in the educational setting.R
EFERENCES[1] C. D. Lopez, M. Cvetkovic, and P. Palensky, “Enhancing PowerFactoryDynamic Models with Python for Rapid Prototyping,” in ,vol. 2019-June. IEEE, 6 2019, pp. 93–99. [Online]. Available:https://ieeexplore.ieee.org/document/8781432/[2] M. Mirz, S. Vogel, G. Reinke, and A. Monti, “DPsim A dynamic phasorreal-time simulator for power systems,”
SoftwareX , vol. 10, p. 100253,2019. [Online]. Available: https://doi.org/10.1016/j.softx.2019.100253[3] H. Cui, F. Li, and K. Tomsovic, “Hybrid Symbolic-NumericFramework for Power System Modeling and Analysis,”
IEEETransactions on Power Systems , 2020. [Online]. Available:https://ieeexplore.ieee.org/document/9169830/ [4] L. Vanfretti, T. Rabuzin, M. Baudette, and M. Murad,“iTesla Power Systems Library (iPSL): A Modelica library forphasor time-domain simulations,”
SoftwareX , vol. 5, pp. 84–88,2016. [Online]. Available: http://dx.doi.org/10.1016/j.softx.2016.05.001https://linkinghub.elsevier.com/retrieve/pii/S2352711016300097[5] M. Baudette, M. Castro, T. Rabuzin et al. , “OpenIPSL: Open-Instance Power System Library Update 1.5 to iTesla PowerSystems Library (iPSL): A Modelica library for phasor time-domainsimulations,”
SoftwareX , vol. 7, pp. 34–36, 2018. [Online]. Available:https://doi.org/10.1016/j.softx.2018.01.002[6] C. R. Harris, K. J. Millman, S. J. van der Walt et al. , “Arrayprogramming with NumPy,”
Nature et al. , “SciPy 1.0: fundamentalalgorithms for scientific computing in Python,”
Nature Methods , vol. 17,no. 3, pp. 261–272, 2020.[8] H. Haugdal, “DynPSSimPy,” 11 2020. [Online]. Available:http://doi.org/10.5281/zenodo.4290126[9] M. Bussonnier, J. Forde, J. Freeman et al. , “Binder 2.0 - Reproducible,interactive, sharable environments for science at scale,” 2018, pp. 113–120.[10] J. Machowski, J. Bumby, and J. Bialek,
Power system dynamics: Stabilityand control , 2nd ed. Wiley, 2008.[11] S. K. Lam, A. Pitrou, and S. Seibert, “Numba,” in