Merge pull request #2 from yangliu0101/release

Documentation
3 years ago · 5b6625a142
parent 2e648a1bd4 041ae74165
commit 5b6625a142
45 changed files with 1631 additions and 4 deletions
--- a/README.md
+++ b/README.md
@ -1,17 +1,193 @@
 # PowerSAS.m
-**PowerSAS.m** is a robust, efficient and scalable power grid analysis framework based on semi-analytical solutions (SAS) technology. The **PowerSAS.m** is the version for MATLAB/Octave users. It currently provides the following functionalities (more coming soon!):
+## I. Documentation
 [HTML](https://powersasm.readthedocs.io/en/latest/index.html)
 [PDF](https://powersasm.readthedocs.io/_/downloads/en/latest/pdf/)
 ## II. What is PowerSAS.m?
 **PowerSAS.m** is a robust, efficient and scalable power grid analysis framework based on semi-analytical solutions (SAS) technology. [(Click here for more details of the SAS technology)](https://powersasm.readthedocs.io/en/latest/sas_basics.html#). 
 The **PowerSAS.m** is the version for MATLAB/Octave users. It currently provides the following functionalities (more coming soon!):
 * **Steady-state analysis**, including power flow (PF), continuation power flow (CPF), contingency analysis.
 * **Dynamic security analysis**, including voltage stability analysis, transient stability analysis, and flexible user-defined simulation.
 * **Hybrid extended-term simulation** provides adaptive QSS-dynamic hybrid simulation in extended term with high accuracy and efficiency.
-### Key features
+#### Key features
 * **High numerical robustness.** Backed by the SAS approach, the PowerSAS tool provides much better convergence than the tools using traditional Newton-type algebraic equation solvers when solving algebraic equations (AE)/ordinary differential equations (ODE)/differential-algebraic equations(DAE).  
 * **Enhanced computational performance.** Due to the analytical nature, PowerSAS provides model-adaptive high-accuracy approximation, which brings significantly extended effective range and much larger steps for steady-state/dynamic analysis. PowerSAS has been used to solve large-scale system cases with 200,000+ buses.
 * **Customizable and extensible.** PowerSAS supports flexible customization of grid analysis scenarios, including complex event sequences in extended simulation term.  
 ## III. Installation
 #### 1. System requirements
 Matlab (7.1 or later) or GNU Octave (4.0.0 or later).
 #### 2. Installation
 * Extract source code to a directory.
 * Enter the directory in Matlab or GNU Octave.
 * Execute command `setup`. You will see the following sub-directories:
    * `/data`: Stores test system data, simulation settings data, etc.
    * `/example`: Some examples of using PowerSAS.m.
    * `/output`: Stores test result data.
    * `/internal`: Internal functions of PowerSAS.m computation core.
    * `/util`: Auxiliary functions including data I/O, plotting, data conversion, etc.
    * `/logging`: Built-in logging system.
    * `/doc`: Documentation.
 #### 3. Test
 * Execute command `initpowersas` to initialize the environment, then execute `test_powersas` to run tests. You should expect all tests to pass.
 #### 4. Initialization
 To initialize PowerSAS.m, add the main directory of PowerSAS.m to your Matlab or GNU Octave path and run the command `initpowersas`. This will ensure that all the functions of PowerSAS.m are added to the path and thus callable.
 ## IV Basic Usage
 ### 1 Initialization before use
 To initialize PowerSAS.m, add the main directory of PowerSAS to Matlab or GNU Octave path, and execute command `initpowersas`. This will ensure all the functions of PowerSAS be added to the path and thus callable.
 ### 2 Call high-level API -- `runPowerSAS`
 Most grid analysis functionalities can be invoked through the high-level function `runPowerSAS`. `runPowerSAS` is defined as follows:
 ```Matlab
 function res=runPowerSAS(simType,data,options,varargin)
 ``` 
 ##### Input arguments
 * `simType` is a string specifying the type of analysis, which can be one of the following values:
    * `'pf'`: Conventional power flow or extended power flow (for finding steady state of dynamic model).
    * `'cpf'`: Continuation power flow.
    * `'ctg'`: Contingency analysis (line outages).
    * `'n-1'`: N-1 line outage screening.
    * `'tsa'`: Transient stability analysis.
    * `'dyn'`: General dynamic simulation.
 * `data` is the system data to be analyzed. It can be either a string specifying the data file name, or a `SysData` struct. For more information about data format and `SysData` struct, please refer to the "Data format and models" chapter.
 * `options` specifies the options for analysis. If you do not provide `options` argument, or if you simply set the field to empty with `[]`, the corresponding routines will provide default options that will fit most cases. See Advanced Use chapter for more details.
 * `varargin` are the additional input variables depending on the type of analysis. Section 3 Basic analysis funtionalifies will explain more details.
 ##### Output
 Output `res` is a struct containing simulation result, system states, system data, etc.
 * `res.flag`: Flag information returned by the analysis task.
 * `res.msg`: More information as supplemental to the flag information.
 * `res.caseName`: The name of the analyzed case.
 * `res.timestamp`: A string showing the timestamp the analysis started, can be viewed as an unique identifier of the analysis task.
 * `res.stateCurve`: A matrix storing the evolution of system states, where the number of rows equals the number of state variables, and the number of columns equals the number of time points.
 * `res.t`: A vector storing time points corresponding to states in `res.stateCurve`.
 * `res.simSettings`: A struct specifying the simulation settings, including simulation parametersand defined events.
 * `res.eventList`: A matrix showing as the list of events in the system in the analysis task.
 * `res.SysDataBase`: A struct of system data at base state.
 * `res.snapshot`: The snapshot of the system states at the end of anlaysis, which can be used to initilize other analysis tasks.
 To access the system states, we need to further access each kind of state variable in `res.stateCurve`. For example, the commands to extract the voltage from `res.stateCurve` are shown below:
 ```Matlab
 [~,idxs]=getIndexDyn(res.SysDataBase); % Get the indexes of each kind of state variables
 vCurve=res.StateCurve(idxs.vIdx,:); % idxs.vIdx is the row indexes of voltage variables
 ```
 ### 3 Basic analysis functionalities
 #### 3.1 Power flow analysis
 When `simType='pf'`, the `runPowerSAS` function runs power flow analysis. In addition to the conventional power flow model, `runPowerSAS` also integrates an extended power flow to solve the steady state of dynamic models. For example, it will calculate the rotor angles of synchronous generators and slips of induction motors in addition to the network equations.
 To perform power flow analysis, call the `runPowerSAS` function as follows:
 ```Matlab
 res=runPowerSAS('pf',data,options)
 ```
 where the argument `data` can either be a string of file name or a `SysData` struct.  
 Below are some examples:
 ```Matlab
 % Use file name string to specify data
 res1=runPowerSAS('pf','d:/codes/d_003.m'); % Filename can use absolute path
 res2=runPowerSAS('pf','d_003.m'); % If data file is already in the Matlab/Octave path,
                                  % then can directly use file name
 res3=runPowerSAS('pf','d_003'); % Filename can be without '.m'
 res4=runPowerSAS('pf','d_003',setOptions('dataPath','d:/codes'); % Another way to specify data path
 % Use SysData struct to specify data
 SysData=readDataFile('d_003.m','d:/codes'); % Generate SysData struct from data file
 res5=runPowerSAS('pf',SysData); % Run power flow using SysData struct
 ```
 #### 3.2 Continuation Power Flow
 Continuation power flow (CPF) analysis in PowerSAS.m features enhanced efficiency and convergence. To perform continuation power flow analysis, call `runPowerSAS` function as follows:
 ```Matlab
 res=runPowerSAS('cpf',data,options,varargin)
 ```
 where `options` (optional) specifies the options of CPF analysis, and `varargin` are the input arguments:
 * `varargin{1}` (optional) is the ramping direction of load, which is an $\text{N}\times \text{12}$ matrix, the first column is the index of the bus, and the columns 5-10 are the ZIP load increase directions.
 * `varargin{2}` (optional) is the ramping direction of generation power, which is an $\text{N}\times \text{2}$ matrix, the first column is the index of the bus, and the 2nd column is the generation increase directions.
 * `varargin{3}` (optional) is the snapshot of the starting state, with which the computation of starting steady state is skipped.
 Some examples can be found in `example/ex_cpf.m`.
 #### 3.3 Contingency Analysis
 Contingency analysis computes the system states immediately after removing a line/lines. To perform contingency analysis, call `runPowerSAS` as follows:
 ```Matlab
 res=runPowerSAS('ctg',data,options,varargin)
 ```
 where `options` (optional) specifies the options of contingency analysis. When not using customized options, set `options=[]`. And `varargin` are the input arguments:
 * `varargin{1}` (mandatory) is a vector specifying the indexes of lines to be removed simultaneously.
 * `varargin{2}` (optional) is the snapshot of the starting state. With this option, computing the starting steady state is skipped.
 Some examples can be found in `example/ex_ctg.m`.
 #### 3.4 N-1 screening
 N-1 screening is essentially performing a series of contingency analysis, each removing a line from the base state. To perform N-1 screening, call `runPowerSAS` as follows:
 ```Matlab
 res=runPowerSAS('n-1',data,options)
 ```
 The return value `res` is a cell containing each contingency analysis results.
 Some examples can be found in `example/ex_n_minus_1.m`.
 #### 3.5 Transient Stability Analysis
 Transient stability anslysis (TSA) assesses the system dynamic behavior and stability after given disturbance(s). 3-phase balanced fault(s) are the most common disturbances in the TSA. In PowerSAS, the TSA supports the analysis of the combinations of multiple faults. To perform transient Stability Analysis, call `runPowerSAS` in the following way:
 ```Matlab
 res=runPowerSAS('tsa',data,options,varargin)
 ```
 where `options` (optional) specifies the options of TSA. When not using customized options, set `options=[]`. And `varargin` are the input arguments:
 * `varargin{1}` (mandatory) is a $\text{N}\times \text{6}$ matrix specifying the faults:
    * The 1st column is the index of line where the fault happens. 
    * The 2nd column is the relative position of the fault, 0.0 stands for the starting terminal and 1.0 stands for the ending terminal. For example, 0.5 means the fault happens in the middle point of the line.
    * The 3rd and 4th columns are the resistance and reactance of the fault.
    * The 5th and 6th columns specify the fault occurrence and clearing times.
 * `varargin{2}` (optional) is the snapshot of the starting state, with which the computation of starting steady state is skipped.
 By default, the TSA is run for 10 seconds. To change the simulation length, specify in the `options` argument, e.g. `options=setOptions('simlen',30)`.
 Example can be found in `example/ex_tsa.m`.
 ### 4. Plot dynamic analysis results
 PowerSAS provides an integrated and formatted way of plotting the system behavior in the time domain. The function for plotting curves is `plotCurves`. The function is defined as follows:
 ```Matlab
 function plotCurves(figId,t,stateCurve,SysDataBase,variable,...
                    subIdx,lw,dot,fontSize,width,height,sizeRatio)
 ```
 The argument list is explained as follows:
 * `figId`: A positive integer or `[]` specifying the index of figure.
 * `t`: A vector of time instants. If got `res` from `runPowreSAS` function, then input this argument as `res.t`.
 * `stateCurve`: A matrix of system states in time domain, the number of columns should equal to the length of `t`. If got `res` from `runPowreSAS` function, then input this argument as `res.stateCurve`.
 * `SysDataBase`: A SysData struct specifying the base system. If got `res` from `runPowreSAS` function, then input this argument as `res.SysDataBase`.
 * `variable`: A string of variable name to be plotted. Here is a nonexhaustive list:
    * `'v'`: voltage magnitude (pu);
    * `'a'`: voltage angle (rad);
    * `'delta'`: rotor angle of synchronous generators;
    * `'omega'`: deviation of synchronous generator rotor speed;
    * `'s'`: induction motor slips;
    * `'f'`: frequency;
 * `subIdx`: Allows you to pick a portion of the variables to plot e.g., the voltage of some selected buses. Default value is `[]`, which means that all the selected type of variables are plotted.
 * `lw`: Line width. Default value is 1.
 * `dot`: Allows you to choose whether to show data points. 1 means curves mark data dots, and 0 means no data dots are shown on curves. The default value is 0.
 * `fontSize`: Font size of labels. Default value is 12.
 * `width`: Width of figure window in pixels.
 * `height`: Height of figure window in pixels.
 * `sizeRatio`: If `width` or `height` is not specified, the size of the figure is determined by the `sizeRatio` of the screen size. The default value of `sizeRatio` is 0.7.
 * **Customizable and extensible.** PowerSAS supports flexible customization of grid analysis scenarios, including complex event sequences in extended simulation term. 
-### Citing
+## V. Citing
 If you are using **PowerSAS.m** in research work to be published, please include explicit citation of our work in your publication. Please place the following entries in your bibliography:
@ -25,4 +201,4 @@ The corresponding BiBTex citations are given below:
             title={{PowerSAS.m} – An Open-Source Power System Simulation Toolbox Based on Semi-Analytical Solution Technologies}, 
             year={2023},
             doi={10.1109/OAJPE.2023.3245040}
-            }
+            }
--- a/docs/Makefile
+++ b/docs/Makefile
@ -0,0 +1,20 @@
 # Minimal makefile for Sphinx documentation
 #
 # You can set these variables from the command line.
 SPHINXOPTS    = "-W"  # This flag turns warnings into errors.
 SPHINXBUILD   = sphinx-build
 SPHINXPROJ    = PackagingScientificPython
 SOURCEDIR     = source
 BUILDDIR      = build
 # Put it first so that "make" without argument is like "make help".
 help:
 	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
 .PHONY: help Makefile
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
 %: Makefile
 	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
--- a/docs/make.bat
+++ b/docs/make.bat
@ -0,0 +1,36 @@
@ECHO OFF
 pushd %~dp0
 REM Command file for Sphinx documentation
 if "%SPHINXBUILD%" == "" (
 	set SPHINXBUILD=sphinx-build
 )
 set SOURCEDIR=source
 set BUILDDIR=build
 set SPHINXPROJ=PackagingScientificPython
 if "%1" == "" goto help
 %SPHINXBUILD% >NUL 2>NUL
 if errorlevel 9009 (
 	echo.
 	echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
 	echo.installed, then set the SPHINXBUILD environment variable to point
 	echo.to the full path of the 'sphinx-build' executable. Alternatively you
 	echo.may add the Sphinx directory to PATH.
 	echo.
 	echo.If you don't have Sphinx installed, grab it from
 	echo.http://sphinx-doc.org/
 	exit /b 1
 )
 %SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
 goto end
 :help
 %SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
 :end
 popd
--- a/docs/requirements-rtd.txt
+++ b/docs/requirements-rtd.txt
@ -0,0 +1,17 @@
 cvxopt
 numpy
 scipy
 sympy
 pandas
 matplotlib
 openpyxl
 xlsxwriter
 dill
 pathos
 tqdm
 pyyaml
 coloredlogs
 ipython
 numpydoc
 sphinx-copybutton
 sphinx_rtd_theme
--- a/docs/source/_static/.placeholder
+++ b/docs/source/_static/.placeholder
--- a/docs/source/about.rst
+++ b/docs/source/about.rst
@ -0,0 +1,47 @@
 About
 =====
 Contributors
 ~~~~~~~~~~~~
 -  Rui Yao, Argonne National Laboratory <ryao@anl.gov>
 -  Feng Qiu, Argonne National Laboratory <fqiu@anl.gov>
 -  Kai Sun, University of Tennessee, Knoxville <kaisun@utk.edu>
 -  Jianzhe Liu, Argonne National Laboratory <jianzhe.liu@ieee.org>
 Acknowledgement
 ~~~~~~~~~~~~~~~
 This work is supported by the Laboratory Directed Research and
 Development (LDRD) program of Argonne National Laboratory, provided by
 the U.S. Department of Energy Office of Science under Contract
 No. DE-AC02-06CH11357, and the U.S. Department of Energy Office of
 Electricity, Advanced Grid Modeling program under Grant DE-OE0000875.
 Publications
 ~~~~~~~~~~~~
 -  Rui Yao, Kai Sun, Feng Qiu, “Vectorized Efficient Computation of Padé
   Approximation for Semi-Analytical Simulation of Large-Scale Power
   Systems,” IEEE Transactions on Power Systems, 34 (5), 3957-3959,
   2019.
 -  Rui Yao, Yang Liu, Kai Sun, Feng Qiu, Jianhui Wang,“Efficient and
   Robust Dynamic Simulation of Power Systems with Holomorphic
   Embedding”, IEEE Transactions on Power Systems, 35 (2), 938 - 949,
   2020.
 -  Rui Yao, Feng Qiu, “Novel AC Distribution Factor for Efficient Outage
   Analysis”, IEEE Transactions on Power Systems, 35 (6), 4960-4963,
   2020.
 -  Xin Xu, Rui Yao, Kai Sun, Feng Qiu, “A Semi-Analytical Solution
   Approach for Solving Constant-Coefficient First-Order Partial
   Differential Equations”, IEEE Control Systems Letters, 6, 704-709,
   2021.
 -  Rui Yao, Feng Qiu, Kai Sun, “Contingency Analysis Based on
   Partitioned and Parallel Holomorphic Embedding”, IEEE Transactions on
   Power Systems, in press.
 License
 ~~~~~~~
 This software is under BSD license. Please refer to LICENSE.md for
 details.
--- a/docs/source/advanced_usage.rst
+++ b/docs/source/advanced_usage.rst
@ -0,0 +1,499 @@
 Advanced Usage
 ==============
 1. Customize analysis with settings file/struct
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 PowerSAS.m lets you customize your simulation by providing a simulation
 settings interface to specify the events and scenarios to be analyzed.
 To use the customized simulation, call the ``runPowerSAS`` function as
 follows:
 .. code:: matlab
   res=runPowerSAS('dyn',data,options,settings,snapshot)
 Details are explained as follows: #### 1.1 Settings file The input
 argument ``settings`` can be a string specifying the settings file name,
 or a struct containing all the simulation settings.
 When ``settings`` is a string, it should be a valid file name of an .m
 script file containing the settings. Some examples of the settings files
 can be found in the directory ``/data``. The settings file should have
 the following variables:
 -  ``eventList``: A gross list of events. (`more
   details <#variable-eventlist>`__)
 -  ``bsBus``: (for black start simulation only) Black start bus
   information. (`more details <#variable-bsbus>`__)
 -  ``bsSyn``: (for black start simulation only) Generator information on
   black start bus. (`more details <#variable-bssyn>`__)
 -  ``bsInd``: (for black start simulation only) Inductin motor on black
   start bus. (`more details <#variable-bsind>`__)
 -  ``Efstd``: Excitation potential of synchronous generators. (`more
   details <#variable-efstd>`__)
 -  ``evtLine``: List of line addition/outage events. (`more
   details <#variables-evtline-and-evtlinespec>`__)
 -  ``evtLineSpec``: Specifications of line addition/outage events.
   (`more details <#variables-evtline-and-evtlinespec>`__)
 -  ``evtZip``: List of static load addition/shedding events. (`more
   details <#variables-evtzip-evtzipspec-and-evtzipspec2>`__)
 -  ``evtZipSpec``: Specifications of static load addition/shedding
   events. (`more
   details <#variables-evtzip-evtzipspec-and-evtzipspec2>`__)
 -  ``evtZipSpec2``: Alternative specifications of static load
   addition/shedding events. (`more
   details <#variables-evtzip-evtzipspec-and-evtzipspec2>`__)
 -  ``evtInd``: List of induction motor addition/outage events. (`more
   details <#variables-evtind-and-evtindspec>`__)
 -  ``evtIndSpec``: Specifications of induction motor addition/outage
   events. (`more details <#variables-evtind-and-evtindspec>`__)
 -  ``evtSyn``: List of synchronous generator addition/outage events.
   (`more details <#variables-evtsyn-and-evtsynspec>`__)
 -  ``evtSynSpec``: Specifications of synchronous generator
   addition/outage events. (`more
   details <#variables-evtsyn-and-evtsynspec>`__)
 -  ``evtFault``: List of fault occurrence/clearing events. (`more
   details <#variables-evtfault-and-evtfaultspec>`__)
 -  ``evtFaultSpec``: Specifications of fault occurrence/clearning
   events. (`more details <#variables-evtfault-and-evtfaultspec>`__)
 -  ``evtDyn``: List of dynamic ramping events. (`more
   details <#variable-evtdyn>`__)
 -  ``evtDynPQ``: Specifications of PQ bus ramping. (`more
   details <#variable-evtdynpq>`__)
 -  ``evtDynPV``: Specifications of PV bus ramping. (`more
   details <#variable-evtdynpv>`__)
 -  ``evtDynInd``: Specifications of induction motor mechanical load
   ramping. (`more details <#variable-evtdynind>`__)
 -  ``evtDynZip``: Specifications of ZIP load ramping. (`more
   details <#variable-evtdynzip>`__)
 -  ``evtDynSh``: Specifications of shunt compensator ramping. (`more
   details <#variable-evtdynsh>`__)
 -  ``evtDynZipRamp``: Alternative specifications of ZIP load ramping.
   (`more details <#variable-evtdynramp>`__)
 -  ``evtDynTmech``: Specifications of generator mechanical torque
   ramping. (`more details <#variable-evtdyntmech>`__)
 -  ``evtDynPm``: Specifications of generator input active power ramping.
   (`more details <#variable-evtdynpm>`__)
 -  ``evtDynEf``: Specifications of generator excitation potential
   ramping. (`more details <#variable-evtdynef>`__)
 -  ``evtDynVref``: Specifications of exciter reference voltage ramping.
   (`more details <#variable-evtdynvref>`__)
 -  ``evtDynEq1``: Specifications of generator transient excitation
   potential ramping. (`more details <#variable-evtdyneq1>`__)
 1.2 Settings struct
 ^^^^^^^^^^^^^^^^^^^
 Alternatively, the ``settings`` can be a struct containing all the
 previous variables as its fields.
 Example 1: Transient stability analysis (TSA) using settings file
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 In this example, we want to perform transient stability analysis on
 2383-bus system. Here is the scenario of the TSA: \* The total
 simulation period is 10 seconds. \* At 0.5 s, apply faults on the
 starting terminals of lines 74 and 114, the fault resistance is 0.0 and
 the reactance is 0.02. At 0.75 s, clear the faults. \* At 1.5 s, apply
 fault on the starting terminal of line 1674, the fault resistance is 0.0
 and the reactance is 0.1. At 1.95 s, clear the faults.
 The settings file for this simulation is shown below. Other variables
 irrelevant to the fault events are omitted here for the sake of clarity.
 .. code:: matlab
   eventList=[...
      1    0.0000  0.0000    0    1 0.0  0.0000
      1.1  0.5000  0.0000    6    1 0.0  0.0000
      1.2  0.7500  0.0000    6    2 0.0  0.0000
      1.3  1.5000  0.0000    6    3 0.0  0.0000
      1.4  1.9500  0.0000    6    4 0.0  0.0000
      3    10.00   0.0000   99    0 0.0  0.0000
   ];
   % Fault event data
   evtFault=[...
       1   1   2
       2   3   4
       3   5   5
       4   6   6
    ];
   evtFaultSpec=[...
        114,  0.00, 0, 0.02,  0;
         74,  0.00, 0, 0.02,  0;
        114,  0.00, 0, 0.02,  1;
         74,  0.00, 0, 0.02,  1;
       1674,  0.00, 0, 0.1,   0;
       1674,  0.00, 0, 0.1,   1;
    ];
 Assume the settings file is ``settings_polilsh_tsa.m`` and the system
 data file is ``d_dcase2383wp_mod2_ind_zip_syn.m``. We can run the TSA as
 follows:
 .. code:: matlab
   res_2383_st=runPowerSAS('pf','d_dcase2383wp_mod2_ind_zip_syn.m'); % Run steady-state
   res_2383_tsa=runPowerSAS('dyn','d_dcase2383wp_mod2_ind_zip_syn.m',setOptions('hotStart',1),'settings_polilsh_tsa',res_2383_st.snapshot); % Hot start from existing steady-state
   plotCurves(1,res_2383_tsa.t,res_2383_tsa.stateCurve,res_2383_tsa.SysDataBase,'v'); % plot the voltage magnitude curves
 2. Extended-term simulation using hybrid QSS and dynamic engines
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 To accelerate computation — especially for extended-term simulation —
 PowerSAS.m provides an adaptive way to switch between QSS and dynamic
 engines in the course of a simulation. With this feature enabled,
 PowerSAS.m can switch to QSS simulation for better speed on detecting
 the fade-away of transients and switch back to dynamic simulation upon
 detecting transient events.
 For more details on the technical approach, please refer to our paper:
 \* Hybrid QSS and Dynamic Extended-Term Simulation Based on Holomorphic
 Embedding, arXiv:2104.02877
 Example 2 illustrates the use of PowerSAS.m to perform extended-term
 simulation.
 Example 2: Extended-term simulation
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 We want to study the response of a 4-bus system under periodic
 disturbances. The entire simulated process is 500 seconds. Starting at
 60 s and continuing until 270 s, the system undergoes events of
 adding/shedding loads every 30 s.
 The key settings of the simulation are:
 .. code:: matlab
   % settings_d_004_2a_agc.m
   eventList=[...
      1    0.0000  0.0000    0    1 0.0  0.0100
      6   60.0000  0.0000    2    1 0.0  0.0100
      7   90.0000  0.0000    2    2 0.0  0.0100
      8  120.0000  0.0000    2    3 0.0  0.0100
      9  150.0000  0.0000    2    4 0.0  0.0100
     10  180.0000  0.0000    2    1 0.0  0.0100
     11  210.0000  0.0000    2    2 0.0  0.0100
     12  240.0000  0.0000    2    3 0.0  0.0100
     13  270.0000  0.0000    2    4 0.0  0.0100
     18  500.0000  0.0000   99    0 0.0  0.0100
   ];
   % Static load event data
   evtZip=[...
      1    1    1    1
      2    1    2    2
      3    1    3    3
      4    1    4    4
   ];
   evtZipSpec2=[...
      3 100.0000 100.0000 60.0000  0.0648  0.0648  0.0648  0.0359  0.0359  0.0359    0    1
      2 100.0000 100.0000 60.0000  0.0648  0.0648  0.0648  0.0359  0.0359  0.0359    0    1
      3 100.0000 100.0000 60.0000 -0.0648 -0.0648 -0.0648 -0.0359 -0.0359 -0.0359    0    1
      2 100.0000 100.0000 60.0000 -0.0648 -0.0648 -0.0648 -0.0359 -0.0359 -0.0359    0    1
   ];
 First we run the simulation in full-dynamic mode and record time:
 .. code:: matlab
   % Full dynamic simulation 
   tagFullDynStart=tic;
   res_004_fulldyn=runPowerSAS('dyn','d_004_2a_bs_agc.m',[]],'settings_d_004_2a_agc');
   timeFullDyn=toc(tagFullDynStart);
 Then we run the simulation in hybrid QSS & dynamic mode and record time:
 .. code:: matlab
   % Hybrid simulation with dynamic-QSS switching
   tagHybridStart=tic;
   res_004=runPowerSAS('dyn','d_004_2a_bs_agc.m',setOptions('allowSteadyDynSwitch',1),'settings_d_004_2a_agc');
   timeHybrid=toc(tagHybridStart);
 Compare the results:
 .. code:: matlab
   plotCurves(1,res_004_fulldyn.t,res_004_fulldyn.stateCurve,res_004_fulldyn.SysDataBase,'v');
   plotCurves(2,res_004.t,res_004.stateCurve,res_004.SysDataBase,'v');
 And compare the computation time:
 .. code:: matlab
   disp(['Full dynamic simulation computation time:', num2str(timeFullDyn),' s.']);
   disp(['Hybrid simulation computation time:', num2str(timeHybrid),' s.']);
 The complete example can be found in
 ``/example/ex_extended_term_dyn.m``. And the results can also be found
 in our paper: \* Hybrid QSS and Dynamic Extended-Term Simulation Based
 on Holomorphic Embedding, arXiv:2104.02877
 Appendix: Variables in settings
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 variable ``eventList``
 ^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) ##### Table 1. Definition of
 ``eventList`` Column \| Content ——-\| ————- 1 \| Event index (can be an
 integer or a real number) 2 \| Event start time 3 \| Event end time (no
 effect for instant event) 4 \| Type of event (see `Table
 2 <#table-2-event-types>`__) 5 \| Index of event under its type 6 \|
 Simulation method (default 0.0) (see below `Simulation
 methods <#simulation-methods>`__) 7 \| Timestep (default 0.01)
 Table 2. Event types
 ''''''''''''''''''''
 (`back to top <#11-settings-file>`__) Value \| Event type ——-\| ————- 0
 \| Calculate steady-state at start 1 \| Add line 2 \| Add static load 3
 \| Add induction motor load 4 \| Add synchronous generator 6 \|
 Applying/clearing faults 7 \| Cut line 8 \| Cut static load 9 \| Cut
 motor load 10\| Cut synchronous generator 50\| Dynamic process 99\| End
 of simulation
 Simulation methods
 ''''''''''''''''''
 Simulation methods can be specified for each event on the 6th column of
 ``eventList``. It is encoded as a number ``x.yz``, where: \* ``x`` is
 the method for solving differential equation, where 0 - SAS, 1 -
 Modified Euler, 2 - R-K 4, 3 - Trapezoidal rule. \* ``y`` is the method
 for solving algebraic equation, where 0 - SAS, 1 - Newton-Raphson. \*
 ``z`` is whether to use variable time step scheme for numerical
 integration (``x`` is 1, 2 or 3). 0 - Fixed step, 1 - Variable step.
 Note that when ``x=0``, ``y`` and ``z`` are not effective, it
 automatically uses SAS and variable time steps.
 variable ``bsBus``
 ^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Current version only support one
 black start bus and therefore only the first line will be recognized.
 Will expand in the future versions. Column \| Content ——-\| ————- 1 \|
 Bus index 2 \| Active power of Z component load 3 \| Active power of I
 component load 4 \| Active power of P component load 5 \| Reactive power
 of Z component load 6 \| Reactive power of I component load 7 \|
 Reactive power of P component load
 variable ``bsSyn``
 ^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of synchronous generator 2 \| Excitation potential 3 \| Active
 power 4 \| Participation factor for power balancing
 variable ``bsInd``
 ^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of induction motor 2 \| Mechanical load torque
 variable ``Efstd``
 ^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) When there are synchronous
 generators in the system model, ``Efstd`` is needed to compute steady
 state. The ``Efstd`` is a column vector specifying the excitation
 potential of every synchronous generator, or it can also be a scalar
 assigning the excitation potentials of all the generator as the same
 value.
 variables ``evtLine`` and ``evtLineSpec``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) In ``eventList``, when the 4th
 column (event type) equals 1 or 7 (add or cut line, respectively), the
 index of the line events in ``evtLine`` corresponds to the 5th column of
 ``eventList``. ##### variable ``evtLine`` (`back to
 top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \| Index of
 line events (from 5th column of ``eventList``) 2 \| Start index in
 ``evtLineSpec`` 3 \| End index in ``evtLineSpec``
 variable ``evtLineSpec``
 ''''''''''''''''''''''''
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of line 2 \| Add/cut mark, 0 - add line, 1 - cut line 3 \|
 Reserved 4 \| Reserved 5 \| Reserved
 variables ``evtZip``, ``evtZipSpec`` and ``evtZipSpec2``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) In ``eventList``, when the 4th
 column (event type) equals 2 or 8 (add/cut static load), the index of
 the load events in ``evtZip`` corresponds to the 5th column of
 ``eventList``. ##### variable ``evtLine`` (`back to
 top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \| Index of
 load events (from 5th column of ``eventList``) 2 \| Choose
 ``evtZipSpec`` (0) or ``evtZipSpec2`` (1) 3 \| Start index in
 ``evtZipSpec`` or ``evtZipSpec2`` 4 \| End index in ``evtZipSpec`` or
 ``evtZipSpec2``
 variable ``evtZipSpec``
 '''''''''''''''''''''''
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of zip loads in system base state 2 \| Add/cut mark, 0 - add load,
 1 - cut load
 variable ``evtZipSpec2`` (recommended)
 ''''''''''''''''''''''''''''''''''''''
 (`back to top <#11-settings-file>`__) ``evtZipSpec2`` has the same
 format as PSAT ZIP load format, which represents the change of ZIP load.
 Whether the event is specified as add/cut load does not make difference.
 variables ``evtInd`` and ``evtIndSpec``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) In ``eventList``, when the 4th
 column (event type) equals 3 or 9 (add or cut induction motors,
 respectively), the index of the induction motor events in ``evtInd``
 corresponds to the 5th column of ``eventList``. ##### variable
 ``evtInd`` (`back to top <#11-settings-file>`__) Column \| Content ——-\|
 ————- 1 \| Index of induction motor events (from 5th column of
 ``eventList``) 2 \| Start index in ``evtIndSpec`` 3 \| End index in
 ``evtIndSpec``
 variable ``evtIndSpec``
 '''''''''''''''''''''''
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of induction motor 2 \| Event type, 0 - add motor, 1 - change
 state, 2 - cut motor 3 \| Designated mechanical torque 4 \| Designated
 slip
 variables ``evtSyn`` and ``evtSynSpec``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) In ``eventList``, when the 4th
 column (event type) equals 4 or 10 (add or cut synchronous generators,
 respectively), the index of the synchronous generator events in
 ``evtSyn`` corresponds to the 5th column of ``eventList``. #####
 variable ``evtSyn`` Column \| Content ——-\| ————- 1 \| Index of
 synchronous generator events (from 5th column of ``eventList``) 2 \|
 Start index in ``evtSynSpec`` 3 \| End index in ``evtSynSpec``
 variable ``evtSynSpec``
 '''''''''''''''''''''''
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of synchronous generator 2 \| Event type, 0 - add generator, 1 -
 cut generator 3 \| Designated rotor angle (only effective when adding
 generator, NaN means the rotor angle is the same with voltage angle). 4
 \| Designated mechanical power (only effective when adding generator). 5
 \| Designated excitation potential (only effective when adding
 generator).
 variables ``evtFault`` and ``evtFaultSpec``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) In ``eventList``, when the 4th
 column (event type) equals 6 (apply or clear faults, respectively), the
 index of the fault events in ``evtFault`` corresponds to the 5th column
 of ``eventList``. In the current version, we only consider three-phase
 grounding faults. ##### variable ``evtFault`` Column \| Content ——-\|
 ————- 1 \| Index of fault events (from 5th column of ``eventList``) 2 \|
 Start index in ``evtFaultSpec`` 3 \| End index in ``evtFaultSpec``
 variable ``evtFaultSpec``
 '''''''''''''''''''''''''
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of fault line 2 \| Position of fault, 0.0 stands for starting
 terminal and 1.0 stands for ending terminal. 3 \| Resistance of fault. 4
 \| Reactance of fault. 5 \| Event type: 0 - add fault; 1 - clear fault.
 variable ``evtDyn``
 ^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) The ``evtDyn`` variable specifies
 the indexes of ramping events involving various types of components.
 Column \| Content ——-\| ————- 1 \| Index of event 2 \| Start index in
 ``evtDynPQ`` 3 \| End index in ``evtDynPQ`` 4 \| Start index in
 ``evtDynPV`` 5 \| End index in ``evtDynPV`` 6 \| Start index in
 ``evtDynInd`` 7 \| End index in ``evtDynInd`` 8 \| Start index in
 ``evtDynZip`` 9 \| End index in ``evtDynZip`` 10\| Start index in
 ``evtDynSh`` 11\| End index in ``evtDynSh`` 12\| Start index in
 ``evtDynZipRamp`` 13\| End index in ``evtDynZipRamp`` 14\| Start index
 in ``evtDynTmech`` 15\| End index in ``evtDynTmech`` 16\| Start index in
 ``evtDynPm`` 17\| End index in ``evtDynPm`` 18\| Start index in
 ``evtDynEf`` 19\| End index in ``evtDynEf`` 20\| Start index in
 ``evtDynVref`` 21\| End index in ``evtDynVref`` 22\| Start index in
 ``evtDynEq1`` 23\| End index in ``evtDynEq1``
 variable ``evtDynPQ``
 ^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of bus 2 \| Active power ramping rate 3 \| Reactive power ramping
 rate
 variable ``evtDynPV``
 ^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of bus 2 \| Active power ramping rate
 variable ``evtDynInd``
 ^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of induction motor 2 \| Mechanical load torque ramping rate
 variable ``evtDynZip``
 ^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of bus 2 \| ZIP load ramping rate
 variable ``evtDynSh``
 ^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of bus 2 \| Shunt admittance ramping rate
 variable ``evtDynZipRamp``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) ``evtDynZipRamp`` has the same
 format as PSAT ZIP load format, which represents the ramping direction
 of ZIP load.
 variable ``evtDynTmech``
 ^^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of synchronous generator 2 \| Ramping rate of mechanical power
 reference value (TG required)
 variable ``evtDynEf``
 ^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of synchronous generator 2 \| Ramping rate of excitation potential
 variable ``evtDynVref``
 ^^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of synchronous generator 2 \| Ramping rate of exciter reference
 voltage
 variable ``evtDynEq1``
 ^^^^^^^^^^^^^^^^^^^^^^
 (`back to top <#11-settings-file>`__) Column \| Content ——-\| ————- 1 \|
 Index of synchronous generator 2 \| Ramping rate of transient excitation
 potential
--- a/docs/source/basic_usage.rst
+++ b/docs/source/basic_usage.rst
@ -0,0 +1,245 @@
 Using PowerSAS.m: The Basics
 ============================
 1 Initialization before use
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 To initialize PowerSAS.m, add the main directory of PowerSAS to Matlab
 or GNU Octave path, and execute command ``initpowersas``. This will
 ensure all the functions of PowerSAS be added to the path and thus
 callable.
 2 Call high-level API – ``runPowerSAS``
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Most grid analysis functionalities can be invoked through the high-level
 function ``runPowerSAS``. ``runPowerSAS`` is defined as follows:
 .. code:: matlab
   function res=runPowerSAS(simType,data,options,varargin)
 Input arguments
 '''''''''''''''
 -  ``simType`` is a string specifying the type of analysis, which can be
   one of the following values:
   -  ``'pf'``: Conventional power flow or extended power flow (for
      finding steady state of dynamic model).
   -  ``'cpf'``: Continuation power flow.
   -  ``'ctg'``: Contingency analysis (line outages).
   -  ``'n-1'``: N-1 line outage screening.
   -  ``'tsa'``: Transient stability analysis.
   -  ``'dyn'``: General dynamic simulation.
 -  ``data`` is the system data to be analyzed. It can be either a string
   specifying the data file name, or a ``SysData`` struct. For more
   information about data format and ``SysData`` struct, please refer to
   the “Data format and models” chapter.
 -  ``options`` specifies the options for analysis. If you do not provide
   ``options`` argument, or if you simply set the field to empty with
   ``[]``, the corresponding routines will provide default options that
   will fit most cases. See Advanced Use chapter for more details.
 -  ``varargin`` are the additional input variables depending on the type
   of analysis. Section 3 Basic analysis funtionalifies will explain
   more details.
 Output
 ''''''
 Output ``res`` is a struct containing simulation result, system states,
 system data, etc. \* ``res.flag``: Flag information returned by the
 analysis task. \* ``res.msg``: More information as supplemental to the
 flag information. \* ``res.caseName``: The name of the analyzed case. \*
 ``res.timestamp``: A string showing the timestamp the analysis started,
 can be viewed as an unique identifier of the analysis task. \*
 ``res.stateCurve``: A matrix storing the evolution of system states,
 where the number of rows equals the number of state variables, and the
 number of columns equals the number of time points. \* ``res.t``: A
 vector storing time points corresponding to states in
 ``res.stateCurve``. \* ``res.simSettings``: A struct specifying the
 simulation settings, including simulation parametersand defined events.
 \* ``res.eventList``: A matrix showing as the list of events in the
 system in the analysis task. \* ``res.SysDataBase``: A struct of system
 data at base state. \* ``res.snapshot``: The snapshot of the system
 states at the end of anlaysis, which can be used to initilize other
 analysis tasks.
 To access the system states, we need to further access each kind of
 state variable in ``res.stateCurve``. For example, the commands to
 extract the voltage from ``res.stateCurve`` are shown below:
 .. code:: matlab
   [~,idxs]=getIndexDyn(res.SysDataBase); % Get the indexes of each kind of state variables
   vCurve=res.StateCurve(idxs.vIdx,:); % idxs.vIdx is the row indexes of voltage variables
 3 Basic analysis functionalities
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 3.1 Power flow analysis
 ^^^^^^^^^^^^^^^^^^^^^^^
 When ``simType='pf'``, the ``runPowerSAS`` function runs power flow
 analysis. In addition to the conventional power flow model,
 ``runPowerSAS`` also integrates an extended power flow to solve the
 steady state of dynamic models. For example, it will calculate the rotor
 angles of synchronous generators and slips of induction motors in
 addition to the network equations.
 To perform power flow analysis, call the ``runPowerSAS`` function as
 follows:
 .. code:: matlab
   res=runPowerSAS('pf',data,options)
 where the argument ``data`` can either be a string of file name or a
 ``SysData`` struct.
 Below are some examples:
 .. code:: matlab
   % Use file name string to specify data
   res1=runPowerSAS('pf','d:/codes/d_003.m'); % Filename can use absolute path
   res2=runPowerSAS('pf','d_003.m'); % If data file is already in the Matlab/Octave path,
                                     % then can directly use file name
   res3=runPowerSAS('pf','d_003'); % Filename can be without '.m'
   res4=runPowerSAS('pf','d_003',setOptions('dataPath','d:/codes'); % Another way to specify data path
   % Use SysData struct to specify data
   SysData=readDataFile('d_003.m','d:/codes'); % Generate SysData struct from data file
   res5=runPowerSAS('pf',SysData); % Run power flow using SysData struct
 3.2 Continuation Power Flow
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Continuation power flow (CPF) analysis in PowerSAS.m features enhanced
 efficiency and convergence. To perform continuation power flow analysis,
 call ``runPowerSAS`` function as follows:
 .. code:: matlab
   res=runPowerSAS('cpf',data,options,varargin)
 where ``options`` (optional) specifies the options of CPF analysis, and
 ``varargin`` are the input arguments: \* ``varargin{1}`` (optional) is
 the ramping direction of load, which is an
 :math:`\text{N}\times \text{12}` matrix, the first column is the index
 of the bus, and the columns 5-10 are the ZIP load increase directions.
 \* ``varargin{2}`` (optional) is the ramping direction of generation
 power, which is an :math:`\text{N}\times \text{2}` matrix, the first
 column is the index of the bus, and the 2nd column is the generation
 increase directions. \* ``varargin{3}`` (optional) is the snapshot of
 the starting state, with which the computation of starting steady state
 is skipped.
 Some examples can be found in ``example/ex_cpf.m``.
 3.3 Contingency Analysis
 ^^^^^^^^^^^^^^^^^^^^^^^^
 Contingency analysis computes the system states immediately after
 removing a line/lines. To perform contingency analysis, call
 ``runPowerSAS`` as follows:
 .. code:: matlab
   res=runPowerSAS('ctg',data,options,varargin)
 where ``options`` (optional) specifies the options of contingency
 analysis. When not using customized options, set ``options=[]``. And
 ``varargin`` are the input arguments: \* ``varargin{1}`` (mandatory) is
 a vector specifying the indexes of lines to be removed simultaneously.
 \* ``varargin{2}`` (optional) is the snapshot of the starting state.
 With this option, computing the starting steady state is skipped.
 Some examples can be found in ``example/ex_ctg.m``.
 3.4 N-1 screening
 ^^^^^^^^^^^^^^^^^
 N-1 screening is essentially performing a series of contingency
 analysis, each removing a line from the base state. To perform N-1
 screening, call ``runPowerSAS`` as follows:
 .. code:: matlab
   res=runPowerSAS('n-1',data,options)
 The return value ``res`` is a cell containing each contingency analysis
 results.
 Some examples can be found in ``example/ex_n_minus_1.m``.
 3.5 Transient Stability Analysis
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Transient stability anslysis (TSA) assesses the system dynamic behavior
 and stability after given disturbance(s). 3-phase balanced fault(s) are
 the most common disturbances in the TSA. In PowerSAS, the TSA supports
 the analysis of the combinations of multiple faults. To perform
 transient Stability Analysis, call ``runPowerSAS`` in the following way:
 .. code:: matlab
   res=runPowerSAS('tsa',data,options,varargin)
 where ``options`` (optional) specifies the options of TSA. When not
 using customized options, set ``options=[]``. And ``varargin`` are the
 input arguments: \* ``varargin{1}`` (mandatory) is a
 :math:`\text{N}\times \text{6}` matrix specifying the faults: \* The 1st
 column is the index of line where the fault happens. \* The 2nd column
 is the relative position of the fault, 0.0 stands for the starting
 terminal and 1.0 stands for the ending terminal. For example, 0.5 means
 the fault happens in the middle point of the line. \* The 3rd and 4th
 columns are the resistance and reactance of the fault. \* The 5th and
 6th columns specify the fault occurrence and clearing times. \*
 ``varargin{2}`` (optional) is the snapshot of the starting state, with
 which the computation of starting steady state is skipped.
 By default, the TSA is run for 10 seconds. To change the simulation
 length, specify in the ``options`` argument,
 e.g. ``options=setOptions('simlen',30)``.
 Example can be found in ``example/ex_tsa.m``.
 4. Plot dynamic analysis results
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 PowerSAS provides an integrated and formatted way of plotting the system
 behavior in the time domain. The function for plotting curves is
 ``plotCurves``. The function is defined as follows:
 .. code:: matlab
   function plotCurves(figId,t,stateCurve,SysDataBase,variable,...
                       subIdx,lw,dot,fontSize,width,height,sizeRatio)
 The argument list is explained as follows: \* ``figId``: A positive
 integer or ``[]`` specifying the index of figure. \* ``t``: A vector of
 time instants. If got ``res`` from ``runPowreSAS`` function, then input
 this argument as ``res.t``. \* ``stateCurve``: A matrix of system states
 in time domain, the number of columns should equal to the length of
 ``t``. If got ``res`` from ``runPowreSAS`` function, then input this
 argument as ``res.stateCurve``. \* ``SysDataBase``: A SysData struct
 specifying the base system. If got ``res`` from ``runPowreSAS``
 function, then input this argument as ``res.SysDataBase``. \*
 ``variable``: A string of variable name to be plotted. Here is a
 nonexhaustive list: \* ``'v'``: voltage magnitude (pu); \* ``'a'``:
 voltage angle (rad); \* ``'delta'``: rotor angle of synchronous
 generators; \* ``'omega'``: deviation of synchronous generator rotor
 speed; \* ``'s'``: induction motor slips; \* ``'f'``: frequency; \*
 ``subIdx``: Allows you to pick a portion of the variables to plot e.g.,
 the voltage of some selected buses. Default value is ``[]``, which means
 that all the selected type of variables are plotted. \* ``lw``: Line
 width. Default value is 1. \* ``dot``: Allows you to choose whether to
 show data points. 1 means curves mark data dots, and 0 means no data
 dots are shown on curves. The default value is 0. \* ``fontSize``: Font
 size of labels. Default value is 12. \* ``width``: Width of figure
 window in pixels. \* ``height``: Height of figure window in pixels. \*
 ``sizeRatio``: If ``width`` or ``height`` is not specified, the size of
 the figure is determined by the ``sizeRatio`` of the screen size. The
 default value of ``sizeRatio`` is 0.7.
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@ -0,0 +1,203 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 # ANDES documentation build configuration file, created by
 # sphinx-quickstart on Thu Jun 28 12:35:56 2018.
 #
 # This file is execfile()d with the current directory set to its
 # containing dir.
 #
 # Note that not all possible configuration values are present in this
 # autogenerated file.
 #
 # All configuration values have a default; values that are commented out
 # serve to show the default.
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
 #
 # import os
 # import sys
 # sys.path.insert(0, os.path.abspath('.'))
 # -- General configuration ------------------------------------------------
 # If your documentation needs a minimal Sphinx version, state it here.
 #
 # needs_sphinx = '1.0'
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = [
    'sphinx.ext.autodoc',
    'sphinx.ext.autosummary',
    'sphinx.ext.githubpages',
    'sphinx.ext.intersphinx',
    'sphinx.ext.mathjax',
    'sphinx.ext.viewcode',
 ]
 # Configuration options for plot_directive. See:
 # https://github.com/matplotlib/matplotlib/blob/f3ed922d935751e08494e5fb5311d3050a3b637b/lib/matplotlib/sphinxext/plot_directive.py#L81
 plot_html_show_source_link = False
 plot_html_show_formats = False
 # Generate the API documentation when building
 autosummary_generate = True
 numpydoc_show_class_members = False
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
 # The suffix(es) of source filenames.
 # You can specify multiple suffix as a list of string:
 #
 # source_suffix = ['.rst', '.md']
 source_suffix = '.rst'
 # The master toctree document.
 master_doc = 'index'
 # General information about the project.
 project = 'powerSAS.m'
 copyright = ''
 author = ''
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
 # built documents.
 # The short X.Y version.
 version = "0.0"
 # The full version, including alpha/beta/rc tags.
 release = "0.0.0"
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
 #
 # This is also used if you do content translation via gettext catalogs.
 # Usually you set "language" from the command line for these cases.
 language = None
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
 # This patterns also effect to html_static_path and html_extra_path
 exclude_patterns = []
 # The name of the Pygments (syntax highlighting) style to use.
 pygments_style = 'sphinx'
 # If true, `todo` and `todoList` produce output, else they produce nothing.
 todo_include_todos = False
 # -- Options for HTML output ----------------------------------------------
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
 html_theme = 'sphinx_rtd_theme'
 import sphinx_rtd_theme
 html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
 # documentation.
 #
 # html_theme_options = {}
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
 html_static_path = ['_static']
 # Custom sidebar templates, must be a dictionary that maps document names
 # to template names.
 #
 # This is required for the alabaster theme
 # refs: http://alabaster.readthedocs.io/en/latest/installation.html#sidebars
 html_sidebars = {
    '**': [
        'relations.html',  # needs 'show_related': True theme option to display
        'searchbox.html',
    ]
 }
 # -- Options for HTMLHelp output ------------------------------------------
 # Output file base name for HTML help builder.
 htmlhelp_basename = 'powerSAS.m'
 # -- Options for LaTeX output ---------------------------------------------
 latex_elements = {
    # The paper size ('letterpaper' or 'a4paper').
    #
    'preamble': r'\DeclareUnicodeCharacter{2588}{-}',
    'papersize': 'letterpaper',
    # The font size ('10pt', '11pt' or '12pt').
    #
    'pointsize': '11pt',
    # Additional stuff for the LaTeX preamble.
    #
    # 'preamble': '',
    # Latex figure (float) alignment
    #
    # 'figure_align': 'htbp',
 }
 # Grouping the document tree into LaTeX files. List of tuples
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).
 # latex_documents = [
 #     (master_doc, 'andes.tex', 'ANDES Manual',
 #      'Hantao Cui', 'manual'),
 # ]
 # -- Options for manual page output ---------------------------------------
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
 # man_pages = [
 #     (master_doc, 'andes', 'ANDES Manual',
 #      [author], 1)
 # ]
 # -- Options for Texinfo output -------------------------------------------
 # Grouping the document tree into Texinfo files. List of tuples
 # (source start file, target name, title, author,
 #  dir menu entry, description, category)
 # texinfo_documents = [
 #     (master_doc, 'andes', 'ANDES Manual',
 #      author, 'andes', 'Python Software for Symbolic Power System Modeling and Numerical Analysis',
 #      'Miscellaneous'),
 # ]
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
    'python': ('https://docs.python.org/3/', None),
    'numpy': ('https://docs.scipy.org/doc/numpy/', None),
    'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None),
    'pandas': ('https://pandas.pydata.org/pandas-docs/stable', None),
    'matplotlib': ('https://matplotlib.org', None),
 }
 # Favorite icon
 html_favicon = 'images/curent.ico'
 # Disable smartquotes to display double dashes correctly
 smartquotes = False
 # import and execute model reference generation script
--- a/docs/source/data_and_models.rst
+++ b/docs/source/data_and_models.rst
@ -0,0 +1,21 @@
 Data and Models
 ===============
 1. Supported data formats
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 Currently PowerSAS.m supports extended PSAT (Matlab) data format.
 Support for other formats and data format conversion features will be
 added in future versions.
 2. Extension of PSAT (Matlab) format
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 2.1 Automatic generation control (AGC) model
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Here PowerSAS provides a simple AGC model. It is named as ``Agc.con`` in
 data files and it is a :math:`\text{N}\times \text{4}` matrix. Each
 column is defined as below: Column \| Content ——-\| ————- 1 \| Bus index
 2 \| Reciprocal of turbine governor gain on bus 3 \| Effective damping
 ratio on bus 4 \| Reciprocal of AGC control time constant
--- a/docs/source/images/curent.ico
+++ b/docs/source/images/curent.ico
--- a/docs/source/images/diagrams/ieeest.png
+++ b/docs/source/images/diagrams/ieeest.png
--- a/docs/source/images/example-kundur/efd.png
+++ b/docs/source/images/example-kundur/efd.png
--- a/docs/source/images/example-kundur/omega.png
+++ b/docs/source/images/example-kundur/omega.png
--- a/docs/source/images/example-npcc/omega.png
+++ b/docs/source/images/example-npcc/omega.png
--- a/docs/source/images/example-tgov1/tgov1.png
+++ b/docs/source/images/example-tgov1/tgov1.png
--- a/docs/source/images/example-tgov1/tgov1_class.png
+++ b/docs/source/images/example-tgov1/tgov1_class.png
--- a/docs/source/images/example-tgov1/tgov1_eqns.png
+++ b/docs/source/images/example-tgov1/tgov1_eqns.png
--- a/docs/source/images/example-wecc/omega.png
+++ b/docs/source/images/example-wecc/omega.png
--- a/docs/source/images/misc/doc-screenshot.png
+++ b/docs/source/images/misc/doc-screenshot.png
--- a/docs/source/images/misc/ieeeg1-screenshot.png
+++ b/docs/source/images/misc/ieeeg1-screenshot.png
--- a/docs/source/images/sponsors/curent.jpg
+++ b/docs/source/images/sponsors/curent.jpg
--- a/docs/source/images/sponsors/doe.png
+++ b/docs/source/images/sponsors/doe.png
--- a/docs/source/images/sponsors/inl.jpg
+++ b/docs/source/images/sponsors/inl.jpg
--- a/docs/source/images/sponsors/llnl.jpg
+++ b/docs/source/images/sponsors/llnl.jpg
--- a/docs/source/images/sponsors/nsf.jpg
+++ b/docs/source/images/sponsors/nsf.jpg
--- a/docs/source/images/tutorial/xlsx-bus.png
+++ b/docs/source/images/tutorial/xlsx-bus.png
--- a/docs/source/images/tutorial/xlsx-pq.png
+++ b/docs/source/images/tutorial/xlsx-pq.png
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@ -0,0 +1,134 @@
 .. powerSAS.m documentation master file, created by
   sphinx-quickstart on 02/21/2023.
   You can adapt this file completely to your liking, but it should at least
   contain the root `toctree` directive.
 .. raw:: html
   <embed>
   <h1 style="letter-spacing: 0.4em; font-size: 2.5em !important;
   margin-bottom: 0; padding-bottom: 0"> powerSAS.m </h1>
   <p style="color: #00746F; font-variant: small-caps; font-weight: bold;
   margin-bottom: 2em">
   Rrobust, Efficient and Scalable Power Grid Analysis Framework based on Semi-Analytical Solutions (SAS) Technology</p>
   </embed>
 ****
 Home
 ****
 .. PowerSAS.m is a robust, efficient and scalable power grid analysis framework based on **semi-analytical solutions (SAS)** technology. 
 .. The PowerSAS.m is the version for MATLAB/Octave users. It currently provides the following functionalities (more coming soon!)
 .. * Steady-state analysis, including power flow (PF), continuation power flow (CPF), contingency analysis.
 .. * Dynamic security analysis, including voltage stability analysis, transient stability analysis, and flexible user-defined simulation.
 .. * Hybrid extended-term simulation provides adaptive QSS-dynamic hybrid simulation in extended term with high accuracy and efficiency.
 PowerSAS.m
 ==========
 **PowerSAS.m** is a robust, efficient and scalable power grid analysis
 framework based on semi-analytical solutions (SAS) technology. The
 **PowerSAS.m** is the version for MATLAB/Octave users. It currently
 provides the following functionalities (more coming soon!):
 -  **Steady-state analysis**, including power flow (PF), continuation
   power flow (CPF), contingency analysis.
 -  **Dynamic security analysis**, including voltage stability analysis,
   transient stability analysis, and flexible user-defined simulation.
 -  **Hybrid extended-term simulation** provides adaptive QSS-dynamic
   hybrid simulation in extended term with high accuracy and efficiency.
 Key features
 ~~~~~~~~~~~~
 -  **High numerical robustness.** Backed by the SAS approach, the
   PowerSAS tool provides much better convergence than the tools using
   traditional Newton-type algebraic equation solvers when solving
   algebraic equations (AE)/ordinary differential equations
   (ODE)/differential-algebraic equations(DAE).
 -  **Enhanced computational performance.** Due to the analytical nature,
   PowerSAS provides model-adaptive high-accuracy approximation, which
   brings significantly extended effective range and much larger steps
   for steady-state/dynamic analysis. PowerSAS has been used to solve
   large-scale system cases with 200,000+ buses.
 -  **Customizable and extensible.** PowerSAS supports flexible
   customization of grid analysis scenarios, including complex event
   sequences in extended simulation term.
 .. ANDES is a Python-based free software package for power system simulation, control and analysis.
 .. It establishes a unique **hybrid symbolic-numeric framework** for modeling differential algebraic
 .. equations (DAEs) for numerical analysis. Main features of ANDES include
 .. * a unique hybrid symbolic-numeric approach to modeling and simulation that enables descriptive DAE modeling and
 ..   automatic numerical code generation
 .. * a rich library of transfer functions and discontinuous components (including limiters, dead-bands, and
 ..   saturation) available for prototyping models, which can be readily instantiated as multiple devices for
 ..   system analysis
 .. * industry-grade second-generation renewable models (solar PV, type 3 and type 4 wind),
 ..   distributed PV and energy storage model
 .. * comes with the Newton method for power flow calculation, the implicit trapezoidal method for time-domain
 ..   simulation, and full eigenvalue calculation
 .. * strictly verified models with commercial software. ANDES obtains identical time-domain simulation results for
 ..   IEEE 14-bus and NPCC system with GENROU and multiple controller models. See the verification link for details.
 .. * developed with performance in mind. While written in Python, ANDES comes with a performance package and can
 ..   finish a 20-second transient simulation of a 2000-bus system in a few seconds on a typical desktop computer
 .. * out-of-the-box PSS/E raw and dyr file support for available models. Once a model is developed, inputs from a
 ..   dyr file can be readily supported
 .. * an always up-to-date equation documentation of implemented models
 .. ANDES is currently under active development. To get involved,
 .. * Follow the tutorial at
 ..   `https://andes.readthedocs.io <https://andes.readthedocs.io/en/stable/tutorial.html>`_
 .. * Checkout the Notebook examples in the
 ..   `examples folder <https://github.com/cuihantao/andes/tree/master/examples>`_
 .. * Try ANDES in Jupyter Notebook
 ..   `with Binder <https://mybinder.org/v2/gh/cuihantao/andes/master>`_
 .. * Download the PDF manual at
 ..   `download <https://andes.readthedocs.io/_/downloads/en/stable/pdf/>`_
 .. * Report issues in the
 ..   `GitHub issues page <https://github.com/cuihantao/andes/issues>`_
 .. * Learn version control with
 ..   `the command-line git <https://git-scm.com/docs/gittutorial>`_ or
 ..   `GitHub Desktop <https://help.github.com/en/desktop/getting-started-with-github-desktop>`_
 .. * If you are looking to develop models, read the
 ..   `Modeling Cookbook <https://andes.readthedocs.io/en/stable/modeling.html>`_
 .. This work was supported in part by the Engineering Research Center Program of
 .. the National Science Foundation and the Department of Energy under NSF Award
 .. Number EEC-1041877 and the `CURENT <https://curent.utk.edu>`_ Industry Partnership Program.
 .. **ANDES is made open source as part of the CURENT Large Scale Testbed project.**
 .. ANDES is developed and actively maintained by `Hantao Cui <https://cui.eecps.com>`_.
 .. See the GitHub repository for a full list of contributors.
 .. toctree::
   :caption: powerSAS.m Manual
   :maxdepth: 3
   :hidden:
   about.rst
   basic_usage.rst
   advanced_usage.rst
   data_and_models.rst
   troubleshooting.rst
   installation.rst
   sas_basics.rst
   copyright.rst
 .. toctree::
   :hidden:
   :caption: API References
   :maxdepth: 3
 Indices and tables
 ==================
 * :ref:`genindex`
 * :ref:`modindex`
 * :ref:`search`
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@ -0,0 +1,41 @@
 Installation
 ============
 1. System requirements
 ~~~~~~~~~~~~~~~~~~~~~~
 Matlab (7.1 or later) or GNU Octave (4.0.0 or later).
 .. _installation-1:
 2. Installation
 ~~~~~~~~~~~~~~~
 -  Extract source code to a directory.
 -  Enter the directory in Matlab or GNU Octave.
 -  Execute command ``setup``. You will see the following
   sub-directories:
   -  ``/data``: Stores test system data, simulation settings data, etc.
   -  ``/example``: Some examples of using PowerSAS.m.
   -  ``/output``: Stores test result data.
   -  ``/internal``: Internal functions of PowerSAS.m computation core.
   -  ``/util``: Auxiliary functions including data I/O, plotting, data
      conversion, etc.
   -  ``/logging``: Built-in logging system.
   -  ``/doc``: Documentation.
 3. Test
 ~~~~~~~
 -  Execute command ``initpowersas`` to initialize the environment, then
   execute ``test_powersas`` to run tests. You should expect all tests
   to pass.
 4. Initialization
 ~~~~~~~~~~~~~~~~~
 To initialize PowerSAS.m, add the main directory of PowerSAS.m to your
 Matlab or GNU Octave path and run the command ``initpowersas``. This
 will ensure that all the functions of PowerSAS.m are added to the path
 and thus callable.
--- a/docs/source/old/about.md
+++ b/docs/source/old/about.md
--- a/docs/source/old/advanced_usage.md
+++ b/docs/source/old/advanced_usage.md
--- a/docs/source/old/basic_usage.md
+++ b/docs/source/old/basic_usage.md
--- a/docs/source/old/data_and_models.md
+++ b/docs/source/old/data_and_models.md
--- a/docs/source/old/img/accuracy_039_1.png
+++ b/docs/source/old/img/accuracy_039_1.png
--- a/docs/source/old/img/accuracy_039_2.png
+++ b/docs/source/old/img/accuracy_039_2.png
--- a/docs/source/old/img/comp_speed_458.png
+++ b/docs/source/old/img/comp_speed_458.png
--- a/docs/source/old/img/comp_time_039.png
+++ b/docs/source/old/img/comp_time_039.png
--- a/docs/source/old/img/contingency_458.png
+++ b/docs/source/old/img/contingency_458.png
--- a/docs/source/old/img/ssa_benchmarking.png
+++ b/docs/source/old/img/ssa_benchmarking.png
--- a/docs/source/old/index.md
+++ b/docs/source/old/index.md
--- a/docs/source/old/index.rst
+++ b/docs/source/old/index.rst
@ -0,0 +1,31 @@
 PowerSAS.m
 ==========
 **PowerSAS.m** is a robust, efficient and scalable power grid analysis
 framework based on semi-analytical solutions (SAS) technology. The
 **PowerSAS.m** is the version for MATLAB/Octave users. It currently
 provides the following functionalities (more coming soon!):
 -  **Steady-state analysis**, including power flow (PF), continuation
   power flow (CPF), contingency analysis.
 -  **Dynamic security analysis**, including voltage stability analysis,
   transient stability analysis, and flexible user-defined simulation.
 -  **Hybrid extended-term simulation** provides adaptive QSS-dynamic
   hybrid simulation in extended term with high accuracy and efficiency.
 Key features
 ~~~~~~~~~~~~
 -  **High numerical robustness.** Backed by the SAS approach, the
   PowerSAS tool provides much better convergence than the tools using
   traditional Newton-type algebraic equation solvers when solving
   algebraic equations (AE)/ordinary differential equations
   (ODE)/differential-algebraic equations(DAE).
 -  **Enhanced computational performance.** Due to the analytical nature,
   PowerSAS provides model-adaptive high-accuracy approximation, which
   brings significantly extended effective range and much larger steps
   for steady-state/dynamic analysis. PowerSAS has been used to solve
   large-scale system cases with 200,000+ buses.
 -  **Customizable and extensible.** PowerSAS supports flexible
   customization of grid analysis scenarios, including complex event
   sequences in extended simulation term.
--- a/docs/source/old/installation.md
+++ b/docs/source/old/installation.md
--- a/docs/source/old/sas_basics.md
+++ b/docs/source/old/sas_basics.md
--- a/docs/source/sas_basics.rst
+++ b/docs/source/sas_basics.rst
@ -0,0 +1,157 @@
 SAS and PowerSAS.m: The Story
 =============================
 1. What are Semi-Analysical Solutions (SAS)?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Semi-analytical solutions (SAS) is a family of computational methods
 that uses certain analytical formulations (e.g., power series, fraction
 of power series, continued fractions) to approximate the solutions of
 mathematical problems. In terms of formulation, they are quite different
 from the commonly used numerical approaches e.g., Newton-Raphson method
 for solving algebraic equations, Runge-Kutta and Trapezoidal methods for
 solving differential equations. The parameters of SAS still need to be
 determined through some (easier and more robustness-guaranteed)
 numerical computation, and thus these methods are called
 semi-analytical.
 2. What are the advantages of SAS?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 In power system modeling and analysis, SAS has been proven to have the
 following features:
 -  **High numerical robustness.** Steady-state analysis usually requires
   solving nonlinear algebraic equations. Traditional tools usually use
   Newton-Raphson method or its variants, whose results can be highly
   dependent on the selection of starting point and they suffer from
   non-convergence problem. In contrast, SAS provides much better
   convergence thanks to the high-level analytical nature.
 -  **Enhanced computational performance.** In dynamic analysis, the
   traditional numerical integration approaches are essentially
   lower-order methods, which are confined to small time steps to avoid
   too-rapid error accumulation. These tiny time steps severely restrict
   the computation speed. In contrast, SAS provides high-order
   approximation, enabling much larger effective time steps and faster
   computation speed.
 -  **More accurate event-driven simulation.** For complex system
   simulation, it is common to simulate discrete events. Traditional
   numerical integration methods only provide solution values on
   discrete time steps and thus may incur substantial errors predicting
   events. In contrast, SAS provides an analytical form of solution as a
   continuous function, and thus can significantly reduce event
   prediction errors.
 3. How is the performance of PowerSAS.m?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 3.1 Benchmarking with traditional methods on Matlab
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 PowerSAS.m shows advantages in both computational robustness and
 efficiency over the traditional approaches.
 On **steady-state analysis**, we have done several benchmarking with
 traditional methods. For example, we test the steady-state contingency
 analysis on PowerSAS.m and Newton-Raphson (NR) method and its variants
 on Matlab. The test is performed on a reduced Eastern-Interconnection
 (EI) system and we tested on 30,000 contingency scenarios. The results
 suggest that the traditional methods have about 1% chance of failing to
 deliver correct results, while SAS has delivered all the correct
 results.
 For more details, please refer to our recent paper:
 -  Rui Yao, Feng Qiu, Kai Sun, “Contingency Analysis Based on
   Partitioned and Parallel Holomorphic Embedding”, IEEE Transactions on
   Power Systems, in press.
 On **dynamic analysis**, we have compared with serveral most commonly
 used traditional numerical approaches for solving ODE/DAEs, including
 modified Euler, Runge-Kutta, and trapezoidal methods. Tests of
 transient-stability analysis on IEEE 39-bus system model and large-scale
 mdodified Polish 2383-bus system model have verified that SAS has
 significant advantages over the traditional methods in both accuracy and
 efficiency.
 **Accuracy comparison on IEEE 39-bus system (1) – Comparison with
 fixed-time-step traditional methods** |accuracy_039_1|
 **Accuracy comparison on IEEE 39-bus system (2) – Comparison with
 variable-time-step traditional method** |accuracy_039_2|
 **Computation time comparison on IEEE 39-bus system**
 .. figure:: https://user-images.githubusercontent.com/96191387/184000437-6aa9150e-d4b1-4297-b982-61e3e68bc2b8.png
   :alt: comp_time_039
   comp_time_039
 For more details, please refer to our recent paper:
 -  Rui Yao, Yang Liu, Kai Sun, Feng Qiu, Jianhui Wang,“Efficient and
   Robust Dynamic Simulation of Power Systems with Holomorphic
   Embedding”, IEEE Transactions on Power Systems, 35 (2), 938 - 949,
   2020.
 3.2 Benchmarking with PSS/E
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 3.2.1 Static Security Region (SSR)
 ''''''''''''''''''''''''''''''''''
 Static Security Region (SSR) is an important decision-support tool
 showing region of stable operating points. However, there are often
 challenges on convergence when computing SSRs, especially near the
 boundaries. So SSR can be used for benchmarking the numerical robustness
 of computational methods.
 We test SSR on IEEE 39-bus system by varying active power of buses 3&4.
 The active power of buses 3&4 are sampled uniformly over the interval of
 [-4000, 4000] MW. The figure below shows the SSR derived by PSS/E and
 PowerSAS.m. It shows that PSS/E result have some irregular outliers
 (about 0.1% of the samples) outside of the SSR and actually are not
 correct solutions of power flow equations. In contrast, PowerSAS.m
 correctly identifies the SSR.
 .. figure:: https://user-images.githubusercontent.com/96191387/184000532-d838e7c4-7dc3-4fd6-98ad-486a596ef33d.png
   :alt: ssa_benchmarking
   ssa_benchmarking
 3.2.2 N-k Contingency analysis
 ''''''''''''''''''''''''''''''
 Contingency ananlysis also has convergence challenges due to large
 disturbances. Here we perform benchmarking between PSS/E (with and
 without non-divergence options) and PowerSAS.m on the N-25 contingency
 analysis on a reduced eastern-interconnection (EI) system with 458
 buses. We increase the load & generation level by 15%, 20%, and 20.7%,
 respectively, as 3 different loading scenarios (loading margin is
 20.791%). In each scenario, we randomly choose 5000 N-25 contingency
 samples.
 .. figure:: https://user-images.githubusercontent.com/96191387/184000600-6ac3101f-d8bc-49bb-b85d-4cea43ab3549.png
   :alt: contingency_458
   contingency_458
 The figure shows the percentage of correct results using different
 tools. It can be seen that PSS/E has some chance to deliver incorrect
 results, and the chance increases with loading level. In contrast,
 PowerSAS.m still returns results all correctly.
 We also compared the computation speeds of PowerSAS.m and PSS/E. The
 figure below shows the average contingency analysis computation time of
 on the 458-bus system. The results show that SAS’s speed is comparable
 to and even faster than PSS/E’s.
 .. figure:: /img/comp_speed_458.png
   :alt: x
   x
 .. |accuracy_039_1| image:: https://user-images.githubusercontent.com/96191387/183999952-362734f7-d40c-4d27-aa79-eb48bdebcebf.png
 .. |accuracy_039_2| image:: https://user-images.githubusercontent.com/96191387/184000210-90382d81-06bb-4cf6-a423-b8588579e0fd.png