add the possibility to initialize the contingency analysis from the results of the N powerflow

Signed-off-by: DONNOT Benjamin <[email protected]>

Author: BDonnot

CHANGELOG.rst

@@ -37,9 +37,16 @@ TODO: integration test with pandapower (see `pandapower/contingency/contingency.
 - [FIXED] a performance issue for all "XXXSingleSlack" (*eg* KLUSingleSlack) algorithm (filling of the initial
   Jacobian matrix was extremly slow due to the massive 'insert' of data in the eigen sparse matrix instead 
   of relying on the "setFromTriplets" method)
+- [FIXED] an normal "exception" was not catched in the `close()` method of LightSimBackend in case the
+  backend was closed before any grid was loaded.
 - [ADDED] possibility to pickle independantly all part of the grid (*eg* gridmodel.get_lines()
   can be pickled independantly from anything else) **NB** pickling and un-pickling 
   lightsim2grid objects can only be used for the same lightsim2grid version.
+- [ADDED] the `init_from_n_powerflow` property for ContingencyAnalysis (and 
+  ContingencyAnalysisCPP). It allows to chose if the computation of the contingencies
+  are initialized with the complex voltages resulting of the powerflow in N 
+  (`init_from_n_powerflow=True`) or if they are initialized from the 
+  given input vector (`init_from_n_powerflow=False`). Defaults to `False`
 
 [0.12.0] 2026-01-06
 --------------------

benchmarks/benchmark_grid_size.py

@@ -22,6 +22,8 @@
 except ImportError:
     from lightsim2grid import SecurityAnalysis as ContingencyAnalysis
 
+from lightsim2grid.solver import SolverType
+
 from tqdm import tqdm
 import os
 from utils_benchmark import print_configuration, get_env_name_displayed
@@ -46,7 +48,7 @@
               "case300.json",
               "case1354pegase.json",
               "case1888rte.json",
-            #   "GBnetwork.json",  # 2224 buses
+              # "GBnetwork.json",  # 2224 buses
               "case2848rte.json",
               "case2869pegase.json",
               "case3120sp.json",
@@ -235,6 +237,7 @@ def run_grid2op_env(env_lightsim, case, reset_solver,
 
 if __name__ == "__main__":
     prng = np.random.default_rng(42)
+    ls_solver_type = SolverType.KLU
     case_names_displayed = [get_env_name_displayed(el) for el in case_names]
     solver_preproc_solver_time = []
     g2op_speeds = []
@@ -324,6 +327,7 @@ def run_grid2op_env(env_lightsim, case, reset_solver,
                                         load_q,
                                         gen_p_g2op,
                                         sgen_p)
+        env_lightsim.backend.set_solver_type(ls_solver_type)
         # Perform the computation using grid2op
         reset_solver = True  # non default
         nb_step_reset = run_grid2op_env(env_lightsim, case, reset_solver,
@@ -360,7 +364,7 @@ def run_grid2op_env(env_lightsim, case, reset_solver,
                                      env_lightsim.backend.tol)
         time_serie._TimeSerie__computed = True
         a_or = time_serie.compute_A()
-        assert status or computer.nb_solver() == nb_step_pp, f"some powerflow diverge for Time Series for {case_name}: {computer.nb_solved()} "
+        assert status or computer.nb_solved() == nb_step_pp, f"some powerflow diverge for Time Series for {case_name}: {computer.nb_solved()} "
 
         if VERBOSE:
             # print detailed results if needed
@@ -389,6 +393,7 @@ def run_grid2op_env(env_lightsim, case, reset_solver,
             sa.add_single_contingency(i)
             if i >= 1000:
                 break
+        sa.init_from_n_powerflow = True
         p_or, a_or, voltages = sa.get_flows()
         computer_sa = sa.computer
         sa_times.append(computer_sa.total_time() + computer_sa.amps_computation_time())

benchmarks/launch_all_benchmark.sh

@@ -0,0 +1,22 @@
+echo "Solver"
+python benchmark_solvers.py --env_name l2rpn_case14_sandbox --no_test --number 1000
+sleep 1
+echo "\n"
+python benchmark_solvers.py --env_name l2rpn_neurips_2020_track2_small --no_test --number 1000
+sleep 5
+echo "\n\n DC solver"
+python benchmark_dc_solvers.py --env_name l2rpn_case14_sandbox --no_test --number 8000
+echo "\n"
+sleep 1
+python benchmark_dc_solvers.py --env_name l2rpn_neurips_2020_track2_small --no_test --number 8000
+echo "\n\n grid size"
+sleep 5
+python benchmark_grid_size.py
+echo "\n\n Vs pypowsybl"
+sleep 5
+python compare_lightsim2grid_pypowsybl.py --case ieee9
+python compare_lightsim2grid_pypowsybl.py --case ieee14
+python compare_lightsim2grid_pypowsybl.py --case ieee30
+python compare_lightsim2grid_pypowsybl.py --case ieee57
+python compare_lightsim2grid_pypowsybl.py --case ieee118
+python compare_lightsim2grid_pypowsybl.py --case ieee300

benchmarks/test_profile.py

@@ -1,25 +1,140 @@
+#!/home/donnotben/Documents/lightsim2grid/venv_ls/bin/python3
+
 #!/bin/python3
+import pickle
+import time
+
+from kiwisolver import Solver
 import grid2op
 from lightsim2grid import LightSimBackend
+from lightsim2grid.solver import SolverType
 import warnings
+import pandapower as pp
+import numpy as np
 
 # usage:
 # perf record ./test_profile.py
 # perf report
-env_name = "l2rpn_neurips_2020_track2_small"
-with warnings.catch_warnings():
-    warnings.filterwarnings("ignore")
-    env_tmp = grid2op.make(env_name)
+# env_name = "l2rpn_neurips_2020_track2_small"
+
+from  benchmark_grid_size import (
+    get_loads_gens,
+    make_grid2op_env
+)
+prng = np.random.default_rng(42)
+    
+CASE_NAME = "case9241pegase.json"
+NB_TS = 10
 
-    param = env_tmp.parameters
-    param.NO_OVERFLOW_DISCONNECTION = True
-    env = grid2op.make(env_name, backend=LightSimBackend(), param=param)
 
+def my_grid2op_env(case_name, nb_ts, prng):
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore")
+        # env_tmp = grid2op.make(env_name)
+        
+        # load the case file
+        case = pp.from_json(case_name)
+        pp.runpp(case)  # for slack
+        load_p_init = 1.0 * case.load["p_mw"].values
+        load_q_init = 1.0 * case.load["q_mvar"].values
+        gen_p_init = 1.0 * case.gen["p_mw"].values
+        sgen_p_init = 1.0 * case.sgen["p_mw"].values
+        load_p, load_q, gen_p, sgen_p = get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init, prng)
 
-def make_steps(env, nb=1000):
+        slack_gens =  np.zeros((nb_ts, case.ext_grid.shape[0]))
+        if "res_ext_grid" in case:
+            slack_gens += np.tile(case.res_ext_grid["p_mw"].values.reshape(1,-1), (nb_ts, 1))
+        gen_p_g2op = np.concatenate((gen_p, slack_gens), axis=1)  
+        env = make_grid2op_env(case,
+                            case_name,
+                            load_p,
+                            load_q,
+                            gen_p_g2op,
+                            sgen_p)
+            
+        # param = env_tmp.parameters
+        # param.NO_OVERFLOW_DISCONNECTION = True
+        # env = grid2op.make(env_name, backend=LightSimBackend(), param=param)
+        env.reset(seed=0, options={"time serie id": 0})
+        env.backend._timer_preproc = 0
+        env.backend._timer_solver = 0
+        
+        with open(f"gridmodel_{case_name}.pickle", "wb") as f:
+            pickle.dump(obj=env.backend._grid, file=f)
+        return env
+
+
+def make_steps_glop(env, nb=NB_TS, reset_algo=False):
     for i in range(nb):
+        if reset_algo:
+            env.backend._grid.tell_solver_need_reset()
         _ = env.step(env.action_space())
+    print(f"Total time: {env.backend._timer_preproc + env.backend._timer_solver}")
+
 
+def main_glop():
+    env = my_grid2op_env(CASE_NAME, NB_TS, prng)
+    make_steps_glop(env, NB_TS, reset_algo=True)
+    
 
+def main_gridmodel(case_name=CASE_NAME, nb_ts=NB_TS, reset_algo=True, solver_used=SolverType.KLU):
+    time_ls = 0.
+    ls_timer_Fx = 0.
+    ls_timer_solve = 0.
+    ls_timer_initialize = 0.
+    ls_timer_check = 0.
+    ls_timer_dSbus = 0.
+    ls_timer_fillJ = 0.
+    ls_timer_Va_Vm = 0.
+    ls_timer_pre_proc = 0.
+    ls_timer_total_nr = 0.
+    
+    with open(f"gridmodel_{case_name}.pickle", "rb") as f:
+        ls_grid = pickle.load(f)
+    ls_grid.change_solver(solver_used)
+    v_init = ls_grid.dc_pf(np.ones(ls_grid.get_bus_vn_kv().shape[0], dtype=complex) * 1.04, 1, 0.1)
+    for _ in range(nb_ts):
+        if reset_algo:
+            ls_grid.tell_solver_need_reset()
+        beg_ls = time.perf_counter()
+        V = ls_grid.ac_pf(v_init, 10, 1e-6)
+        end_ls = time.perf_counter()
+        time_ls += end_ls - beg_ls
+        if V.shape[0] == 0:
+            raise RuntimeError("Divergence")
+        
+        (timer_Fx_, timer_solve_, timer_initialize_, 
+         timer_check_, timer_dSbus_, timer_fillJ_, 
+         timer_Va_Vm_, timer_pre_proc_, timer_total_nr_
+         ) = ls_grid.get_solver().get_timers_jacobian()
+        ls_grid.unset_changes()  # tell lightsim2grid that the state of the grid is consistent
+        ls_timer_Fx += timer_Fx_
+        ls_timer_solve += timer_solve_
+        ls_timer_initialize += timer_initialize_
+        ls_timer_check += timer_check_
+        ls_timer_dSbus += timer_dSbus_
+        ls_timer_fillJ += timer_fillJ_
+        ls_timer_Va_Vm += timer_Va_Vm_
+        ls_timer_pre_proc += timer_pre_proc_
+        ls_timer_total_nr += timer_total_nr_
+    print(f"Solver used: {solver_used}")
+    print(f"Nb iter: {nb_ts}")
+    print(f"Do reset each step: {reset_algo}")
+    print("--------------------------------------")
+    print("Detailed lightsim2grid timings: ")
+    print(f"Total time {time_ls * 1000.:.2e}ms => {time_ls / nb_ts * 1e3:.2e} ms/pf | {nb_ts / time_ls:.0f} pf /s")
+    print(f"Total time spent in the solver: {1e3 * ls_timer_total_nr:.2e} ms ({100. * ls_timer_total_nr / time_ls:.0f}% of total time spent in lightsim2grid)")
+    print(f"\t Time to pre process Ybus, Sbus etc.: {1e3 * ls_timer_pre_proc:.2e} ms ({100. * ls_timer_pre_proc / ls_timer_total_nr:.0f} % of time in solver)")
+    print(f"\t Time to initialize linear solver {1e3 * ls_timer_initialize:.2e} ms ({100. * ls_timer_initialize / ls_timer_total_nr:.0f} % of time in solver)")
+    print(f"\t Time to compute dS/dV {1e3 * ls_timer_dSbus : .2e} ms ({100. * ls_timer_dSbus / ls_timer_total_nr:.0f} % of time in solver)")
+    print(f"\t Time to fill the Jacobian {1e3 * ls_timer_fillJ:.2e} ms ({100. * ls_timer_fillJ / ls_timer_total_nr:.0f} % of time in solver)")
+    print(f"\t Time to solve the Jacobian linear system: {1e3 * ls_timer_solve:.2e} ms ({100. * ls_timer_solve / ls_timer_total_nr:.0f} % of time in solver)")
+    print(f"\t Time to update Va and Vm {1e3*ls_timer_Va_Vm:.2e} ms ({100. * ls_timer_Va_Vm / ls_timer_total_nr:.0f} % of time in solver)")
+    print(f"\t Time to evaluate p,q mismmatch at each bus {1e3*ls_timer_Fx:.2e} ms ({100. * ls_timer_Fx / ls_timer_total_nr:.0f} % of time in solver)")
+    print(f"\t Time to evaluate cvg criteria {1e3*ls_timer_check:.2e} ms ({100. * ls_timer_check / ls_timer_total_nr:.0f} % of time in solver)")
+    print("--------------------------------------\n")
+    
+    
 if __name__ == "__main__":
-    make_steps(env, nb=1000)
+    # my_grid2op_env(case_name, nb_ts, prng)
+    main_gridmodel(CASE_NAME, NB_TS, reset_algo=True, solver_used=SolverType.KLUSingleSlack)

env_compile_all.sh

@@ -122,6 +122,16 @@ def all_contingencies(self, val):
         raise RuntimeError("Impossible to add new topologies like this. Please use `add_single_contingency` "
                            "or `add_multiple_contingencies`.")
 
+    @property
+    def init_from_n_powerflow(self):
+        return self.computer.init_from_n_powerflow
+    
+    @init_from_n_powerflow.setter
+    def init_from_n_powerflow(self, val):
+        if bool(val) != val:
+            raise ValueError("The `init_from_n_powerflow` attribute must be a boolean.")
+        self.computer.init_from_n_powerflow = bool(val)
+        
     # TODO implement that !
     def __update_grid(self, backend_act):
         raise NotImplementedError("TODO !")

lightsim2grid/contingencyAnalysis.py

@@ -1458,7 +1458,12 @@ def close(self) -> None:
             # self.init_pp_backend.close()  # should not close it, the same init_pp_backend is used when copied
             self._init_pp_backend = None
         self._reset_res_pointers()
-        self._fill_nans()
+        try:
+            self._fill_nans()
+        except AttributeError:
+            # some attributes are not completely filled
+            
+            pass
         self._grid = None
         self.__me_at_init = None
 

lightsim2grid/lightSimBackend.py

@@ -219,7 +219,9 @@ void ContingencyAnalysis::compute(const CplxVect & Vinit, int max_iter, real_typ
         bus_pq_,
         max_iter,
         tol);
-
+    // check if we init the n-1 cases with results from the n cases
+    // or not
+    if(_init_from_n_powerflow) Vinit_solver = _solver.get_V();
     // end of pre processing
     _timer_pre_proc = timer_preproc.duration();
     if(!conv) return;

src/batch_algorithm/ContingencyAnalysis.cpp

@@ -25,7 +25,8 @@ class ContingencyAnalysis final: public BaseBatchSolverSynch
                             _li_coeffs(),
                             _timer_total(0.),
                             _timer_modif_Ybus(0.),
-                            _timer_pre_proc(0.)
+                            _timer_pre_proc(0.),
+                            _init_from_n_powerflow(false)
                             { }
 
         ~ContingencyAnalysis() noexcept = default;
@@ -34,6 +35,9 @@ class ContingencyAnalysis final: public BaseBatchSolverSynch
         ContingencyAnalysis & operator=(ContingencyAnalysis&&) = delete;
         ContingencyAnalysis & operator=(const ContingencyAnalysis&) = delete;
 
+        bool get_init_from_n_powerflow() const noexcept {return _init_from_n_powerflow;}
+        void set_init_from_n_powerflow(bool do_it) noexcept {_init_from_n_powerflow = do_it;}
+
         // utilities to add defaults to simulate
         void add_all_n1(){
             for(int l_id = 0; l_id < n_total_; ++l_id){
@@ -173,6 +177,7 @@ class ContingencyAnalysis final: public BaseBatchSolverSynch
         // sometimes, when i perform some disconnection, I make the graph non connexe
         // in this case, well, i don't use the results of the simulation
         bool check_invertible(const Eigen::SparseMatrix<cplx_type> & Ybus) const;
+
     private:
         // li_default
         std::set<std::set<int> > _li_defaults;  // do not use unordered_set here, we rely on the order for different functions !
@@ -182,5 +187,7 @@ class ContingencyAnalysis final: public BaseBatchSolverSynch
         double _timer_total;  // total time spent in "compute"
         double _timer_modif_Ybus;  // time to update the Ybus between the defaults simulation
         double _timer_pre_proc;  // time to compute the coefficients of the Ybus
+
+        bool _init_from_n_powerflow;
 };
 #endif  //COMPUTERS_H

src/batch_algorithm/ContingencyAnalysis.hpp

@@ -1421,6 +1421,15 @@ between 0 and `n_sub_ * max_nb_bus_per_sub_`
 
     py::class_<ContingencyAnalysis>(m, "ContingencyAnalysisCPP", DocSecurityAnalysis::SecurityAnalysis.c_str())
         .def(py::init<const GridModel &>())
+        .def_property("init_from_n_powerflow",
+                      &ContingencyAnalysis::get_init_from_n_powerflow,
+                      &ContingencyAnalysis::set_init_from_n_powerflow,
+                      R"mydelim(Whether to initialize the complex voltages of "
+                      "each contingencies with the results of a n-powerflow "
+                      "(*ie* a powerflow without any line disconnection) or not. "
+                      "Default: false, meaning each simulation is initialized "
+                      "with the given input vector)mydelim")
+
         // solver control
         .def("change_solver", &ContingencyAnalysis::change_solver, DocGridModel::change_solver.c_str())
         .def("available_solvers", &ContingencyAnalysis::available_solvers, DocGridModel::available_solvers.c_str())