Merge pull request #25 from openEDI/capacitor_fixes

Capacitor fixes

LocalFeeder/FeederSimulator.py

@@ -196,47 +196,95 @@ def get_vnom(self):
         Vnom = np.abs(Vnom)
         return Vnom
 
-    def get_PQs(self):
+    def get_PQs_load(self,static=False):
         num_nodes = len(self._name_index_dict.keys())
 
+        PQ_names= ['' for i in range(num_nodes)]
+        PQ_types= ['' for i in range(num_nodes)]
         PQ_load = np.zeros((num_nodes), dtype=np.complex_)
         for ld in get_loads(dss,self._circuit):
+            self._circuit.SetActiveElement('Load.' + ld["name"])
             for ii in range(len(ld['phases'])):
-                name = ld['bus1'] + '.' + ld['phases'][ii]
-                index = self._name_index_dict[name.upper()]
-                self._circuit.SetActiveElement('Load.' + ld["name"])
-                power = dss.CktElement.Powers()
-                PQ_load[index] += np.complex(power[2 * ii], power[2 * ii + 1])
+                name = ld['bus1'].upper() + '.' + ld['phases'][ii]
+                index = self._name_index_dict[name]
+                if static:
+                    power = complex(ld['kW'],ld['kVar'])
+                    PQ_load[index] += power/len(ld['phases'])
+                else:
+                    power = dss.CktElement.Powers()
+                    PQ_load[index] += np.complex(power[2 * ii], power[2 * ii + 1])
+                PQ_names[index] = name
+                PQ_types[index] = 'Load'
+
+        return PQ_load,PQ_names,PQ_types
+
+    def get_PQs_pv(self,static=False):
+        num_nodes = len(self._name_index_dict.keys())
 
+        PQ_names= ['' for i in range(num_nodes)]
+        PQ_types= ['' for i in range(num_nodes)]
         PQ_PV = np.zeros((num_nodes), dtype=np.complex_)
         for PV in get_pvSystems(dss):
             bus = PV["bus"].split('.')
             if len(bus) == 1:
                 bus = bus + ['1', '2', '3']
             self._circuit.SetActiveElement('PVSystem.'+PV["name"])
-            power = dss.CktElement.Powers()
             for ii in range(len(bus) - 1):
-                index = self._name_index_dict[(bus[0] + '.' + bus[ii + 1]).upper()]
-                PQ_PV[index] += np.complex(power[2 * ii], power[2 * ii + 1])
+                name = bus[0].upper() + '.' + bus[ii + 1]
+                index = self._name_index_dict[name]
+                if static:
+                    power = complex(-1*PV['kW'],-1*PV['kVar']) #-1 because injecting
+                    PQ_PV[index] += power/(len(bus)-1)
+                else:
+                    power = dss.CktElement.Powers()
+                    PQ_PV[index] += np.complex(power[2 * ii], power[2 * ii + 1])
+                PQ_names[index] = name
+                PQ_types[index] = 'PVSystem'
+        return PQ_PV,PQ_names,PQ_types
+
+    def get_PQs_gen(self,static=False):
+        num_nodes = len(self._name_index_dict.keys())
 
+        PQ_names= ['' for i in range(num_nodes)]
+        PQ_types= ['' for i in range(num_nodes)]
         PQ_gen = np.zeros((num_nodes), dtype=np.complex_)
-        PQ_gen_all = np.zeros((num_nodes), dtype=np.complex_)
         for PV in get_Generator(dss):
             bus = PV["bus"]
             self._circuit.SetActiveElement('Generator.'+PV["name"])
-            power = dss.CktElement.Powers()
             for ii in range(len(bus) - 1):
-                index = self._name_index_dict[(bus[0] + '.' + bus[ii + 1]).upper()]
-                PQ_gen_all[index] += np.complex(power[2 * ii], power[2 * ii + 1])
-                PQ_gen[index] += np.complex(power[2 * ii], power[2 * ii + 1])
+                name = bus[0].upper() + '.' + bus[ii + 1]
+                index = self._name_index_dict[name]
+                if static:
+                    power = complex(-1*PV['kW'],-1*PV['kVar']) #-1 because injecting
+                    PQ_gen[index] += power/(len(bus)-1)
+                else:
+                    power = dss.CktElement.Powers()
+                    PQ_gen[index] += np.complex(power[2 * ii], power[2 * ii + 1])
+                PQ_names[index] = name
+                PQ_types[index] = 'Generator'
+        return PQ_gen,PQ_names,PQ_types
+
+    def get_PQs_cap(self,static=False):
+        num_nodes = len(self._name_index_dict.keys())
 
-        Qcap = [0] * num_nodes
+        PQ_names= ['' for i in range(num_nodes)]
+        PQ_types= ['' for i in range(num_nodes)]
+        PQ_cap = np.zeros((num_nodes), dtype=np.complex_)
         for cap in get_capacitors(dss):
             for ii in range(cap["numPhases"]):
-                index = self._name_index_dict[cap["busname"].upper() + '.' + cap["busphase"][ii]]
-                Qcap[index] = -cap["power"][2 * ii - 1]
+                name = cap["busname"].upper() + '.' + cap["busphase"][ii]
+                index = self._name_index_dict[name]
+                if static:
+                    power = complex(0,-1*cap['kVar']) #-1 because it's injected into the grid
+                    PQ_cap[index] += power/cap["numPhases"]
+                else:
+                    PQ_cap[index] = np.complex(0,cap["power"][2 * ii + 1])
+                PQ_names[index] = name
+                PQ_types[index] = 'Capacitor'
+
+        return PQ_cap,PQ_names,PQ_types
+
 
-        return PQ_load + PQ_PV + PQ_gen #+ 1j * np.array(Qcap)  # power injection
 
 
     def get_loads(self):

LocalFeeder/dss_functions.py

@@ -39,7 +39,7 @@ def get_y_matrix_file(dss):
     dss.run_command('batchedit load..* enabled=false')
     dss.run_command('batchedit generator..* enabled=false')
     dss.run_command('batchedit pvsystem..* enabled=false')
-    dss.run_command('Batchedit Capacitor..* enabled=yes')
+    dss.run_command('Batchedit Capacitor..* enabled=false')
     dss.run_command('batchedit storage..* enabled=false')
     dss.run_command('CalcVoltageBases')
     dss.run_command('set maxiterations=20')
@@ -277,7 +277,7 @@ def get_capacitors(dss):
         capname = dss.CktElement.Name()
         NumPhase = dss.CktElement.NumPhases()
         bus = dss.CktElement.BusNames()[0]
-        kvar = dss.run_command('? ' + capname + '.kVar')
+        kvar = dss.Capacitors.kvar()
         datum["name"] = capname
         temp = bus.split('.')
         datum["busname"] = temp[0]

LocalFeeder/sender_cosim.py

@@ -9,11 +9,11 @@
 from typing import List, Tuple
 import numpy as np
 from datetime import datetime, timedelta
-from gadal.gadal_types.data_types import Complex,Topology,VoltagesReal,VoltagesImaginary,PowersReal,PowersImaginary, AdmittanceMatrix, VoltagesMagnitude, VoltagesAngle
+from gadal.gadal_types.data_types import Complex,Topology,VoltagesReal,VoltagesImaginary,PowersReal,PowersImaginary, AdmittanceMatrix, VoltagesMagnitude, VoltagesAngle, Injection
 
 logger = logging.getLogger(__name__)
 logger.addHandler(logging.StreamHandler())
-logger.setLevel(logging.INFO)
+logger.setLevel(logging.DEBUG)
 
 test_se = False
 
@@ -99,6 +99,7 @@ def get_phase(name):
     ]
 
     unique_ids = sim._AllNodeNames
+    snapshot_run(sim)
 
     logger.debug("y-matrix")
     logger.debug(y_matrix)
@@ -125,11 +126,31 @@ def get_phase(name):
            values = phases,
            ids = unique_ids
     )
-
+    all_PQs = {}
+    all_PQs['load'] = sim.get_PQs_load(static=True) # Return type is PQ_values, PQ_names, PQ_types all with same size
+    all_PQs['pv'] = sim.get_PQs_pv(static=True)
+    all_PQs['gen'] = sim.get_PQs_gen(static=True)
+    all_PQs['cap'] = sim.get_PQs_cap(static=True)
+
+    PQ_real = []
+    PQ_imaginary = []
+    PQ_names = []
+    PQ_types = []
+    for i in range(len(all_PQs['load'][0])):
+        for key in all_PQs:
+            if all_PQs[key][1][i] !='':
+                PQ_real.append(-1*all_PQs[key][0][i].real) # injections are negative singe PQ values are positive
+                PQ_imaginary.append(-1*all_PQs[key][0][i].imag) # injections are negative singe PQ values are positive
+                PQ_names.append(all_PQs[key][1][i])
+                PQ_types.append(all_PQs[key][2][i])
+    power_real = PowersReal(ids = PQ_names, values = PQ_real, equipment_type = PQ_types)
+    power_imaginary = PowersImaginary(ids = PQ_names, values = PQ_imaginary, equipment_type = PQ_types)
+    injections = Injection(power_real=power_real, power_imaginary = power_imaginary)
+
     topology = Topology(
         admittance=admittancematrix,
         base_voltage_angles=base_voltageangle,
-        injections={},
+        injections=injections,
         base_voltage_magnitudes=base_voltagemagnitude,
         slack_bus=slack_bus,
     )
@@ -139,7 +160,6 @@ def get_phase(name):
         f.write(topology.json())
     pub_topology.publish(topology.json())
 
-    snapshot_run(sim)
 
     granted_time = -1
     current_hour = 0
@@ -161,25 +181,35 @@ def get_phase(name):
         sim.solve(current_hour,current_second)
 
         feeder_voltages = sim.get_voltages_actual()
-        PQ_node = sim.get_PQs()
+        all_PQs = {}
+        all_PQs['load'] = sim.get_PQs_load(static=False) # Return type is PQ_values, PQ_names, PQ_types all with same size
+        all_PQs['pv'] = sim.get_PQs_pv(static=False)
+        all_PQs['gen'] = sim.get_PQs_gen(static=False)
+        all_PQs['cap'] = sim.get_PQs_cap(static=False)
+
+        PQ_injections_all = all_PQs['load'][0]+ all_PQs['pv'][0] + all_PQs['gen'][0]+ all_PQs['cap'][0]
+
         logger.debug("Feeder Voltages")
         logger.debug(feeder_voltages)
         logger.debug("PQ")
-        logger.debug(PQ_node)
+        logger.debug(PQ_injections_all)
         logger.debug("Calculated Power")
         Cal_power = feeder_voltages * (Y.conjugate() @ feeder_voltages.conjugate()) / 1000
-        errors = PQ_node + Cal_power
-        PQ_node[:3] = -Cal_power[:3]
+        errors = PQ_injections_all + Cal_power
+        sort_errors = np.sort(np.abs(errors))
         logger.debug("errors")
         logger.debug(errors)
+        if np.any(sort_errors[:-3] > 1):
+            raise ValueError('Power balance does not hold')
+        PQ_injections_all[sim._source_indexes] = -Cal_power[sim._source_indexes]
         power_balance = (feeder_voltages * (Y.conjugate() @ feeder_voltages.conjugate()) / 1000)
         logger.debug(power_balance)
         indices, = np.nonzero(np.abs(errors) > 1)
         logger.debug("Indices with error > 1")
         logger.debug(indices)
         logger.debug([sim._AllNodeNames[i] for i in indices])
         logger.debug("Power, Voltages, and Calculated Power at Indices")
-        logger.debug(PQ_node[indices])
+        logger.debug(PQ_injections_all[indices])
         logger.debug(feeder_voltages[indices])
         logger.debug(power_balance[indices])
 
@@ -190,22 +220,14 @@ def get_phase(name):
            values = phases,
            ids = unique_ids
         )
-        topology = Topology(
-            admittance=admittancematrix,
-            base_voltage_angles=base_voltageangle,
-            injections={},
-            base_voltage_magnitudes=base_voltagemagnitude,
-            slack_bus=slack_bus,
-        )
-        pub_topology.publish(topology.json())
 
         logger.debug('Publish load ' + str(feeder_voltages.real[0]))
         voltage_magnitudes = np.abs(feeder_voltages.real + 1j* feeder_voltages.imag)
         pub_voltages_magnitude.publish(VoltagesMagnitude(values=list(voltage_magnitudes), ids=sim._AllNodeNames, time = current_timestamp).json())
         pub_voltages_real.publish(VoltagesReal(values=list(feeder_voltages.real), ids=sim._AllNodeNames, time = current_timestamp).json())
         pub_voltages_imag.publish(VoltagesImaginary(values=list(feeder_voltages.imag), ids=sim._AllNodeNames, time = current_timestamp).json())
-        pub_powers_real.publish(PowersReal(values=list(PQ_node.real), ids=sim._AllNodeNames, time = current_timestamp).json())
-        pub_powers_imag.publish(PowersImaginary(values=list(PQ_node.imag), ids=sim._AllNodeNames, time = current_timestamp).json())
+        pub_powers_real.publish(PowersReal(values=list(PQ_injections_all.real), ids=sim._AllNodeNames, time = current_timestamp).json())
+        pub_powers_imag.publish(PowersImaginary(values=list(PQ_injections_all.imag), ids=sim._AllNodeNames, time = current_timestamp).json())
 
         logger.info('end time: '+str(datetime.now()))
 

test_system.json

@@ -47,7 +47,7 @@
                 "smartds_year": "2017",
                 "start_date": "2017-05-01 00:00:00",
                 "increment_value": 900,
-                "number_of_timesteps": 36
+                "number_of_timesteps": 96
             }
         },
         {
@@ -57,7 +57,7 @@
                 "feeder_file": "gadal-ieee123/qsts/master.dss",
                 "start_date": "2017-01-01 00:00:00",
                 "run_freq_sec": 900,
-                "number_of_timesteps": 36,
+                "number_of_timesteps": 96,
                 "start_time_index": 0,
                 "topology_output": "../../outputs/topology.json"
             }

wls_federate/state_estimator_federate.py

@@ -155,15 +155,17 @@ def state_estimator(parameters: AlgorithmParameters, topology, P, Q, V, initial_
     logging.debug(delta)
     X0 = np.concatenate((delta, Vabs))
     logging.debug(X0)
+    ang_low =  np.concatenate(([-1e-5], np.ones(num_node -1)* (- np.inf)))
+    ang_up =  np.concatenate(([1e-5], np.ones(num_node -1)* ( np.inf)))
+    mag_low =  np.ones(num_node)* (- np.inf)
+    mag_up =  np.ones(num_node)* (np.inf)
+    low_limit = np.concatenate((ang_low, mag_low))
+    up_limit = np.concatenate((ang_up, mag_up))
     # Weights are ignored since errors are sampled from Gaussian
     # Real dimension of solutions is
     # 2 * num_node - len(knownP) - len(knownV) - len(knownQ)
     if len(knownP) + len(knownV) + len(knownQ) < num_node * 2:
         #If not observable 
-        low_limit = np.concatenate((np.ones(num_node)* (- np.pi - np.pi/6),
-                                    np.ones(num_node)*0.90))
-        up_limit = np.concatenate((np.ones(num_node)* (np.pi + np.pi/6),
-                                    np.ones(num_node)*1.05))
         res_1 = least_squares(
             residual,
             X0,
@@ -181,7 +183,7 @@ def state_estimator(parameters: AlgorithmParameters, topology, P, Q, V, initial_
             residual,
             X0,
             jac=cal_H,
-            # bounds = (low_limit, up_limit),
+            bounds = (low_limit, up_limit),
             #method = 'lm',
             verbose=2,
             ftol=tol,
@@ -254,39 +256,40 @@ def run(self):
 
         self.initial_ang = None
         self.initial_V = None
+        topology = Topology.parse_obj(self.sub_topology.json)
+        slack_index =  None
+        if not isinstance(topology.admittance, AdmittanceMatrix):
+            raise "Weighted Least Squares algorithm expects AdmittanceMatrix as input"
+
+        for i in range(len(topology.admittance.ids)):
+            if topology.admittance.ids[i] == topology.slack_bus[0]:
+                slack_index = i
+
         while granted_time < h.HELICS_TIME_MAXTIME:
 
-            topology = Topology.parse_obj(self.sub_topology.json)
             if not self.sub_voltages_magnitude.is_updated():
                 granted_time = h.helicsFederateRequestTime(self.vfed, h.HELICS_TIME_MAXTIME)
                 continue
 
             logger.info('start time: '+str(datetime.now()))
 
-            slack_index =  None
-            if not isinstance(topology.admittance, AdmittanceMatrix):
-                raise "Weighted Least Squares algorithm expects AdmittanceMatrix as input"
-
-            for i in range(len(topology.admittance.ids)):
-                if topology.admittance.ids[i] == topology.slack_bus[0]:
-                    slack_index = i
 
             voltages = VoltagesMagnitude.parse_obj(self.sub_voltages_magnitude.json)
+            power_Q = PowersImaginary.parse_obj(self.sub_power_Q.json)
+            power_P = PowersReal.parse_obj(self.sub_power_P.json)
+            knownP = get_indices(topology, power_P)
+            knownQ = get_indices(topology, power_Q)
             knownV = get_indices(topology, voltages)
+
             if self.initial_V is None:
-                self.initial_V = np.mean(
-                    np.array(voltages.values) / np.array(topology.base_voltage_magnitudes.values)[knownV])
-
-            #if self.initial_V is None:
-               # self.initial_V = 1.025 #*np.array(topology.base_voltages)
+                #Flat start or using average measurements
+                if len(knownP) + len(knownV) + len(knownQ) > len(topology.admittance.ids) * 2:
+                    self.initial_V = 1.0
+                else:
+                    self.initial_V = np.mean(np.array(voltages.values) / np.array(topology.base_voltage_magnitudes.values)[knownV])
             if self.initial_ang is None:
                 self.initial_ang = np.array(topology.base_voltage_angles.values)
 
-
-
-            power_P = PowersReal.parse_obj(self.sub_power_P.json)
-            power_Q = PowersImaginary.parse_obj(self.sub_power_Q.json)
-
             voltage_magnitudes, voltage_angles = state_estimator(
                 self.algorithm_parameters,
                 topology, power_P, power_Q, voltages, initial_V=self.initial_V,