diff --git a/LocalFeeder/FeederSimulator.py b/LocalFeeder/FeederSimulator.py index 5d4e8f0..8ed4503 100644 --- a/LocalFeeder/FeederSimulator.py +++ b/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): diff --git a/LocalFeeder/dss_functions.py b/LocalFeeder/dss_functions.py index 9a7e595..4967700 100644 --- a/LocalFeeder/dss_functions.py +++ b/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] diff --git a/LocalFeeder/sender_cosim.py b/LocalFeeder/sender_cosim.py index acb527e..e42f464 100644 --- a/LocalFeeder/sender_cosim.py +++ b/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,17 +181,27 @@ 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) @@ -179,7 +209,7 @@ def get_phase(name): 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())) diff --git a/errors.png b/errors.png index 5d57e7f..d1d4017 100644 Binary files a/errors.png and b/errors.png differ diff --git a/test_system.json b/test_system.json index 20c5470..e887c7a 100644 --- a/test_system.json +++ b/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" } diff --git a/voltage_angles_0.png b/voltage_angles_0.png index 43afb5c..52c3651 100644 Binary files a/voltage_angles_0.png and b/voltage_angles_0.png differ diff --git a/voltage_angles_95.png b/voltage_angles_95.png index 268fd6f..4523c1e 100644 Binary files a/voltage_angles_95.png and b/voltage_angles_95.png differ diff --git a/voltage_magnitudes_0.png b/voltage_magnitudes_0.png index 6d85f78..cff3b5a 100644 Binary files a/voltage_magnitudes_0.png and b/voltage_magnitudes_0.png differ diff --git a/voltage_magnitudes_95.png b/voltage_magnitudes_95.png index 03606c2..dc479ef 100644 Binary files a/voltage_magnitudes_95.png and b/voltage_magnitudes_95.png differ diff --git a/wls_federate/state_estimator_federate.py b/wls_federate/state_estimator_federate.py index 23f77dc..62ea8ac 100644 --- a/wls_federate/state_estimator_federate.py +++ b/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,