Merge branch 'dev_reduced_3_phase'

README.md

@@ -33,3 +33,31 @@ python identify_network.py --help
 python plot_results.py --help
 ```
 
+The files in the [ROOT_REPO]/conf folder also provide variables that can be modified. 
+Comments in their contents explain what they do.
+
+## Simulation for CDC21 manuscript
+Simulations for our manuscript *Bayesian error-in-variable models for the identification of distribution grids* 
+([available here](https://infoscience.epfl.ch/record/290186?&ln=en)) have been generated using the following commands.
+``` shell script
+python simulate_network.py -vle --network cigre_mv_feeder1
+python identify_network.py -vlw --network cigre_mv_feeder1 --max-iterations 50
+python plot_results.py -v --network cigre_mv_feeder1 --color-scale 1.0
+```
+The parameters in [ROOT_REPO]/conf need however to be adjusted as follows.
+``` python script
+conf.simulation.days to 1
+conf.simulation.time_steps to 500
+conf.identification.lambda_eiv = 4e7
+conf.identification.lambdaprime = 20
+conf.identification.use_tls_diag = True
+```
+DISCLAIMER: Due to newly introduced seed resets for reproductibility and potential differences in load profiles, 
+the results may not be exactly the same as in the paper.
+
+## Simulation for other manuscript
+``` shell script
+python simulate_network.py -vqte --network ieee123center
+python identify_network.py -vqe --max-iterations 20 --network ieee123center --phases 012
+python plot_results.py -v --network ieee123center
+```

conf/conf.py

@@ -2,20 +2,24 @@
 import numpy as np
 import pandas as pd
 import pkgutil
+import cupy
 
+# Simulation parameters
 seed = 131
+noise_level = 1
 
+# Technical settings
 ROOT_DIR = Path(__file__).parent.parent.absolute()
 DATA_DIR = ROOT_DIR / "data"
 SIM_DIR = DATA_DIR / "simulations_output"
 
 TEST_SIM_DIR = ROOT_DIR / "test" / "data" / "simulations_output"
 
-np.set_printoptions(threshold=np.inf)
+np.set_printoptions(threshold=np.inf, linewidth=np.inf)
 pd.set_option('display.max_rows', 500)
 pd.set_option('display.max_columns', 500)
 
 GPU_AVAILABLE = pkgutil.find_loader('cupy')
 # Use this to switch GPU if one is in use.
 # If all GPUs are used you are in bad luck, wait for your turn
-CUDA_DEVICE_USED = 3
+cupy.cuda.Device(3).use()

conf/identification.py

@@ -1,10 +1,15 @@
-lambda_lasso = 1e-7
-lambda_eiv = 3e-5#2e13 #4e11 for non-constant loads reduced net
+from conf.conf import noise_level
+
+# Regularization parameters
+lambda_lasso = 4e-10
+lambda_eiv = 5e-7
 lambdaprime = 200
 
+# Tolerance for stopping the iterative algorithm
 abs_tol = 1e-13
 rel_tol = 1e-16
 
+# What should be regularized with the Bayesian prior
 use_tls_diag = False
 contrast_each_row = True
 regularize_diag = False

conf/simulation.py

@@ -1,23 +1,32 @@
 import numpy as np
+from conf.conf import noise_level
 
+# Definition of hidden nodes
 observed_nodes = [1, 3, 4, 6, 8, 9, 10, 12, 15, 16, 17, 18, 19, 22, 24, 26, 28,
                   32, 36, 37, 39, 40, 43, 44, 46, 47, 49, 50, 51, 52, 53, 55]  # About 60% of the nodes
 hidden_nodes = list(set(range(56-1)) - set(np.array(observed_nodes)-1))
-
 hidden_nodes = []
-constant_load_hidden_nodes = False
 
+# Assumption of constant Z loads on hidden nodes
+constant_load_hidden_nodes = True
+
+# Week during which the sample size starts
 selected_weeks = np.array([12])
-days = len(selected_weeks)
-time_steps = 24 * 60# 15000
 
-current_magnitude_sd = 1e-4
-voltage_magnitude_sd = 1e-4
-current_phase_sd = 1e-4#0.01*np.pi/180
-voltage_phase_sd = 1e-4#0.01*np.pi/180
+# Total sample size on which block averaging is performed
+days = 7*len(selected_weeks)
+
+# Number of averaged blocks
+time_steps = 24 * 60 * 7  # 15000 is max capable for GPU RTX 3090
+
+# Sensor sampling frequency
 measurement_frequency = 50 # [Hz] # 50
+
+# Safety factor for sensor dimension (compared to nominal node power)
 safety_factor = 4
 
-#for 3ph
-#days = len(selected_weeks)
-#hidden_nodes = []
+# Do not touch that, rather change noise_level in conf.conf
+current_magnitude_sd = 1e-4 * noise_level
+voltage_magnitude_sd = 1e-4 * noise_level
+current_phase_sd = 0.01*np.pi/180 * noise_level
+voltage_phase_sd = 0.01*np.pi/180 * noise_level

identify_network.py

@@ -4,6 +4,7 @@
 from src.identification import run, run3ph
 import conf.identification
 from conf import simulation
+from src.simulation.lines import admittance_phase_to_sequence, measurement_phase_to_sequence
 
 
 @click.command()
@@ -12,18 +13,28 @@
 @click.option('--standard', "-s", is_flag=True, help='Redo only standard methods')
 @click.option('--bayesian-eiv', "-b", is_flag=True, help='Redo only Bayesian eiv methods')
 @click.option('--continue-id', "-c", is_flag=True, help='Is the matrix laplacian')
-@click.option('--three-phased', "-t", is_flag=True, help='Identify asymmetric network')
+@click.option('--phases', "-p", type=str, default="1", help='Phases or sequences to identify')
+@click.option('--sequence', "-q", is_flag=True, help='Use zero/positive/negative sequence values')
+@click.option('--exact', "-e", is_flag=True, help='Use exact values for prior')
 @click.option('--laplacian', "-l", is_flag=True, help='Is the matrix laplacian')
 @click.option('--weighted', "-w", is_flag=True, help='Use covariance matrices as weights')
+@click.option('--uncertainty', "-u", is_flag=True, help='Analyse error covariance')
 @click.option('--verbose', "-v", is_flag=True, help='Activates verbosity')
 
-def identify(network, max_iterations, standard, bayesian_eiv, continue_id, three_phased, laplacian, weighted, verbose):
+def identify(network, max_iterations, standard, bayesian_eiv, continue_id,
+             phases, sequence, exact, laplacian, weighted, uncertainty, verbose):
     if verbose:
         def pprint(a):
             print(a)
     else:
         pprint = lambda a: None
 
+    if (phases != "012" and phases != "0" and phases != "1" and phases != "2" and sequence) or \
+            (phases != "123" and phases != "1" and phases != "2" and phases != "3" and not sequence):
+        raise NotImplementedError("Error: two phases identification not implemented, " +
+                                  "try to identify phases one by one separately")
+    three_phased = phases == "012" or phases == "123"
+
     redo_standard_methods = standard or (not standard and not bayesian_eiv)
     redo_STLS = bayesian_eiv or (not standard and not bayesian_eiv)
 
@@ -33,13 +44,28 @@ def pprint(a):
 
     pprint("Loading network simulation...")
     sim_STLS = np.load(conf.conf.DATA_DIR / ("simulations_output/" + name + ".npz"))
-    noisy_current = sim_STLS['i']
-    noisy_voltage = sim_STLS['v']
-    current = sim_STLS['j']
-    voltage = sim_STLS['w']
-    y_bus = sim_STLS['y']
-    pmu_ratings = sim_STLS['p']
-    fparam = sim_STLS['f']
+    noisy_current, noisy_voltage = sim_STLS['i'], sim_STLS['v']
+    current, voltage, y_bus = sim_STLS['j'], sim_STLS['w'], sim_STLS['y']
+    pmu_ratings, fparam, phases_ids = sim_STLS['p'], sim_STLS['f'], sim_STLS['h']
+
+    if not three_phased and np.any(phases_ids != 1):
+        if not np.any(phases_ids == int(phases)):
+            raise IndexError("Error: Trying to identify phase or sequence " + phases +
+                             " but no data is available for it")
+        else:
+            pprint("Identifying " + ("sequence" if sequence else "phase") + " " + phases + "...")
+
+        # Keep chosen phase/sequence
+        idx_todel = (phases_ids != int(phases)).nonzero()
+        y_bus = np.delete(np.delete(y_bus, idx_todel, axis=0), idx_todel, axis=1)
+        pmu_ratings = np.delete(pmu_ratings, idx_todel)
+        [noisy_voltage, noisy_current, voltage, current] = [np.delete(m, idx_todel, axis=1) for m in
+                                                            [noisy_voltage, noisy_current, voltage, current]]
+    pprint("Done!")
+
+    pprint("Saving reference network result...")
+    sim_STLS = {'y': y_bus, 'p': phases, 'l': laplacian}
+    np.savez(conf.conf.DATA_DIR / ("simulations_output/reference_results_" + name + ".npz"), **sim_STLS)
     pprint("Done!")
 
     # Updating variance
@@ -49,21 +75,8 @@ def pprint(a):
     current_phase_sd = simulation.current_phase_sd/np.sqrt(fparam)
 
     y_ols, y_tls, y_lasso = run_fcns.standard_methods(name, noisy_voltage if redo_standard_methods else None,
-                                                      noisy_current, laplacian, max_iterations, verbose)
-
-    fim_wtls, cov_wtls, expected_max_rrms = run_fcns.eiv_fim(name, voltage, current,
-                                                             voltage_magnitude_sd + 1j * voltage_phase_sd,
-                                                             current_magnitude_sd + 1j * current_phase_sd,
-                                                             pmu_ratings, y_bus, laplacian, verbose)
-    print(expected_max_rrms*100)
-    meaned_volts = voltage-np.tile(np.mean(voltage, axis=0), (voltage.shape[0], 1))
-    #print(np.linalg.eigvals(meaned_volts.T @ meaned_volts))
-    print(np.std(noisy_voltage - voltage, axis=0))
-    print(np.std(noisy_current - current, axis=0))
-
-    # TODO: implement 3-phase Bayesian eiv identification
-    if three_phased:
-        return
+                                                      noisy_current, phases_ids if phases == "012" else None,
+                                                      laplacian, max_iterations, verbose)
 
     if continue_id:
         pprint("Loading previous bayesian eiv identification result...")
@@ -75,11 +88,15 @@ def pprint(a):
     else:
         y_init = (y_tls + y_tls.T).copy()/2
 
-    y_sparse_tls_cov, sparse_tls_cov_iterations = run_fcns.bayesian_eiv(name, noisy_voltage, noisy_current,
+    y_exact = y_bus if exact and not laplacian else None
+
+    y_sparse_tls_cov, sparse_tls_cov_iterations = run_fcns.bayesian_eiv(name, noisy_voltage, noisy_current, phases_ids,
                                                                         voltage_magnitude_sd + 1j*voltage_phase_sd,
                                                                         current_magnitude_sd + 1j*current_phase_sd,
-                                                                        pmu_ratings, y_init, laplacian, weighted,
-                                                                        max_iterations if redo_STLS else 0, verbose)
+                                                                        pmu_ratings, y_init, y_exact, laplacian,
+                                                                        weighted, max_iterations if redo_STLS else 0,
+                                                                        verbose)
+
     if continue_id:
         sparse_tls_cov_iterations = sparse_tls_cov_old_iterations.extend(sparse_tls_cov_iterations)
 
@@ -88,5 +105,26 @@ def pprint(a):
         np.savez(conf.conf.DATA_DIR / ("simulations_output/bayesian_results_" + name + ".npz"), **sim_STLS)
         pprint("Done!")
 
+    # Error covariance analysis
+    if uncertainty:
+        fim_wtls, cov_wtls, expected_max_rrms = run_fcns.eiv_fim(name, voltage, current,
+                                                                 voltage_magnitude_sd + 1j * voltage_phase_sd,
+                                                                 current_magnitude_sd + 1j * current_phase_sd,
+                                                                 pmu_ratings, y_bus, None, laplacian, verbose)
+
+        fim_est, cov_est, expected_max_rrms_uncertain = run_fcns.eiv_fim(name, noisy_voltage, noisy_current,
+                                                                         voltage_magnitude_sd + 1j * voltage_phase_sd,
+                                                                         current_magnitude_sd + 1j * current_phase_sd,
+                                                                         pmu_ratings, y_sparse_tls_cov, None, laplacian,
+                                                                         verbose)
+        print("Cramer rao bound (exact+approx): " + str(expected_max_rrms*100)
+              + ", " + str(expected_max_rrms_uncertain*100))
+
+        if not three_phased:
+            pprint("Saving uncertainty results...")
+            sim_STLS = {'e': cov_wtls, 'b': cov_est}
+            np.savez(conf.conf.DATA_DIR / ("simulations_output/uncertainty_results_" + name + ".npz"), **sim_STLS)
+            pprint("Done!")
+
 if __name__ == '__main__':
     identify()