Skip to content

Commit

Permalink
Fix initial set points and adjust timing
Browse files Browse the repository at this point in the history
  • Loading branch information
Joseph McKinsey committed Feb 21, 2024
1 parent 3bb3273 commit 5fdf758
Show file tree
Hide file tree
Showing 3 changed files with 93 additions and 78 deletions.
43 changes: 26 additions & 17 deletions LocalFeeder/FeederSimulator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,14 @@
import xarray as xr
from botocore import UNSIGNED
from botocore.config import Config
from dss_functions import (get_capacitors, get_generators, get_loads,
get_pvsystems, get_voltages)
from oedisi.types.data_types import (Command, InverterControl,
InverterControlMode)
from dss_functions import (
get_capacitors,
get_generators,
get_loads,
get_pvsystems,
get_voltages,
)
from oedisi.types.data_types import Command, InverterControl, InverterControlMode
from pydantic import BaseModel
from scipy.sparse import coo_matrix, csc_matrix

Expand Down Expand Up @@ -53,7 +57,7 @@ class FeederConfig(BaseModel):
existing_feeder_file: Optional[str] = None
sensor_location: Optional[str] = None
start_date: str
number_of_timesteps: float
number_of_timesteps: int
run_freq_sec: float = 15 * 60
start_time_index: int = 0
topology_output: str = "topology.json"
Expand Down Expand Up @@ -701,7 +705,8 @@ def set_properties_to_inverter(self, inverter: str, inv_control: InverterControl
)
if inv_control.vwcontrol is not None:
vw_curve = self.create_xy_curve(
inv_control.vwcontrol.voltage, inv_control.vwcontrol.power_response,
inv_control.vwcontrol.voltage,
inv_control.vwcontrol.power_response,
)
dss.Text.Command(f"{inverter}.voltwatt_curve={vw_curve.split('.')[1]}")
dss.Text.Command(
Expand All @@ -713,26 +718,30 @@ def set_properties_to_inverter(self, inverter: str, inv_control: InverterControl
dss.Text.Command(f"{inverter}.Mode = {inv_control.mode.value}")

def set_pv_output(self, pv_system, p, q):
"""Sets the P and Q values for a PV system in OpenDSS
"""
"""Sets the P and Q values for a PV system in OpenDSS"""
max_pv = self.get_max_pv_available(pv_system)
#pf = q / ((p**2 + q **2)**0.5)
# pf = q / ((p**2 + q **2)**0.5)

obj_name = f"PVSystem.{pv_system}"
if max_pv <=0 or p == 0:
if max_pv <= 0 or p == 0:
Warning("Maximum PV Value is 0")
obj_val = 100
q=0
q = 0
elif p < max_pv:
obj_val = p/float(max_pv) *100
obj_val = p / float(max_pv) * 100
else:
obj_val = 100
ratio = float(max_pv)/p
q = q*ratio #adjust q value to that it matches the kw output
command = [Command(obj_name=obj_name,obj_property="%Pmpp",val=str(obj_val)), Command(obj_name=obj_name,obj_property="kvar",val=str(q)), Command(obj_name=obj_name,obj_property="%Cutout", val="0"), Command(obj_name=obj_name,obj_property="%Cutin", val="0")]
ratio = float(max_pv) / p
q = q * ratio # adjust q value to that it matches the kw output
command = [
Command(obj_name=obj_name, obj_property="%Pmpp", val=str(obj_val)),
Command(obj_name=obj_name, obj_property="kvar", val=str(q)),
Command(obj_name=obj_name, obj_property="%Cutout", val="0"),
Command(obj_name=obj_name, obj_property="%Cutin", val="0"),
]
self.change_obj(command)

def get_max_pv_available(self,pv_system):
def get_max_pv_available(self, pv_system):
dss.PVsystems.First()
irradiance = None
pmpp = None
Expand All @@ -744,7 +753,7 @@ def get_max_pv_available(self,pv_system):
break
if irradiance is None or pmpp is None:
raise ValueError(f"Irradiance or PMPP not found for {pv_system}")
return irradiance*pmpp
return irradiance * pmpp

def get_available_pv(self):
pv_names = []
Expand Down
40 changes: 27 additions & 13 deletions LocalFeeder/sender_cosim.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,22 @@
import xarray as xr
from FeederSimulator import FeederConfig, FeederSimulator
from oedisi.types.common import BrokerConfig
from oedisi.types.data_types import (AdmittanceMatrix, AdmittanceSparse,
CommandList, EquipmentNodeArray,
Injection, InverterControlList,
MeasurementArray, PowersImaginary,
PowersReal, Topology, VoltagesAngle,
VoltagesImaginary, VoltagesMagnitude,
VoltagesReal)
from oedisi.types.data_types import (
AdmittanceMatrix,
AdmittanceSparse,
CommandList,
EquipmentNodeArray,
Injection,
InverterControlList,
MeasurementArray,
PowersImaginary,
PowersReal,
Topology,
VoltagesAngle,
VoltagesImaginary,
VoltagesMagnitude,
VoltagesReal,
)
from scipy.sparse import coo_matrix

logger = logging.getLogger(__name__)
Expand Down Expand Up @@ -322,9 +331,7 @@ def go_cosim(
sub_invcontrol.option["CONNECTION_OPTIONAL"] = 1

pv_set_key = (
"unused/pv_set"
if "pv_set" not in input_mapping
else input_mapping["pv_set"]
"unused/pv_set" if "pv_set" not in input_mapping else input_mapping["pv_set"]
)

sub_pv_set = vfed.register_subscription(pv_set_key, "")
Expand All @@ -340,8 +347,14 @@ def go_cosim(
pub_topology.publish(initial_data.topology.json())

granted_time = -1
for request_time in range(0, int(config.number_of_timesteps)):
request_time = 0
while request_time < int(config.number_of_timesteps):
granted_time = h.helicsFederateRequestTime(vfed, request_time)
assert (
granted_time <= request_time + deltat
), f"granted_time: {granted_time} past {request_time}"
if granted_time >= request_time + deltat:
request_time += 1

current_index = int(granted_time) # floors
current_timestamp = datetime.strptime(
Expand Down Expand Up @@ -398,7 +411,8 @@ def go_cosim(
voltage_magnitudes = np.abs(current_data.feeder_voltages)
pub_voltages_magnitude.publish(
VoltagesMagnitude(
**xarray_to_dict(voltage_magnitudes), time=current_timestamp,
**xarray_to_dict(voltage_magnitudes),
time=current_timestamp,
).json()
)
pub_voltages_real.publish(
Expand Down Expand Up @@ -430,7 +444,7 @@ def go_cosim(
MeasurementArray(
**xarray_to_dict(sim.get_available_pv()),
time=current_timestamp,
units="kWA"
units="kWA",
).json()
)

Expand Down
88 changes: 40 additions & 48 deletions omoo_federate/OMOO.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -385,9 +385,10 @@ def opf_run(self, V, P_wopv, Q_wopv):
np.delete(P_wopv, self.slack_bus),
np.delete(Q_wopv, self.slack_bus),
)
# Initial P and Q setpoint
P_0, Q_0 = (
self.pv_frame["avai"].values / self.base_power,
self.pv_frame["avaiQ"].values / self.base_power,
np.zeros(len(self.pv_frame)),
)
Ppv, Qpv = np.zeros(self.num_node), np.zeros(self.num_node)
for pp, pv_ii in enumerate(self.pv_index):
Expand All @@ -405,11 +406,11 @@ def opf_run(self, V, P_wopv, Q_wopv):
)
if len(ind) == 0:
P_0, Q_0 = P_0 * self.base_power, Q_0 * self.base_power
Flag = False
set_power = False
logger.debug("Skip this step since no violation")
logger.debug(f"minimum V_hat is {np.min(V_hat)}")
logger.debug(f"maximum V_hat is {np.max(V_hat)}")
return P_0, Q_0, Flag, V_hat
return P_0, Q_0, set_power, V_hat
# logger.debug('ind')
# logger.debug(ind)
else:
Expand Down Expand Up @@ -482,11 +483,10 @@ def opf_run(self, V, P_wopv, Q_wopv):
logger.debug(
f"Target bounds are [{self.parameters.Vmax}, {self.parameters.Vmin}]"
)
Flag = True
return (
Pk_last * self.base_power,
Qk_last * self.base_power,
Flag,
True,
V_hat_final,
)

Expand Down Expand Up @@ -615,10 +615,9 @@ def run(self):
voltages_imag = VoltagesImaginary.parse_obj(
self.sub_voltages_imaginary.json
)
if voltages is None:
voltages = measurement_to_xarray(
voltages_real
) + 1j * measurement_to_xarray(voltages_imag)
voltages = measurement_to_xarray(
voltages_real
) + 1j * measurement_to_xarray(voltages_imag)
logger.debug(np.max(np.abs(voltages) / v))
assert topology.base_voltage_magnitudes.ids == list(voltages.ids.data)

Expand All @@ -634,9 +633,12 @@ def run(self):
MeasurementArray.parse_obj(self.sub_available_power.json)
)

split_power = available_power / pv_injections.ids.groupby("equipment_ids").count().rename({"equipment_ids": "ids"})
available_power = split_power.loc[pv_injections.equipment_ids].assign_coords(ids=pv_injections.ids)

split_power = available_power / pv_injections.ids.groupby(
"equipment_ids"
).count().rename({"equipment_ids": "ids"})
available_power = split_power.loc[
pv_injections.equipment_ids
].assign_coords(ids=pv_injections.ids)

pv = pd.DataFrame()
pv["name"] = pv_ratings.equipment_ids.data
Expand All @@ -645,8 +647,6 @@ def run(self):

# This needs to be fixed.
pv["avai"] = available_power
pv["avaiQ"] = np.sqrt(pv_ratings**2 - available_power**2) # TODO: Get a better pf limit

bus_to_index = {
v: i for i, v in enumerate(topology.base_voltage_magnitudes.ids)
}
Expand All @@ -661,10 +661,8 @@ def run(self):
logger.debug("PVframe")
logger.debug(pv)

if power_P is None:
power_P = PowersReal.parse_obj(self.sub_power_P.json)
if power_Q is None:
power_Q = PowersImaginary.parse_obj(self.sub_power_Q.json)
power_P = PowersReal.parse_obj(self.sub_power_P.json)
power_Q = PowersImaginary.parse_obj(self.sub_power_Q.json)
assert topology.base_voltage_magnitudes.ids == power_P.ids
assert topology.base_voltage_magnitudes.ids == power_Q.ids
ts = time.time()
Expand All @@ -681,7 +679,9 @@ def run(self):
self.w_mag,
)

P_set, Q_set, flag, V_hat = opf.opf_run(np.abs(voltages), power_P, power_Q)
P_set, Q_set, set_power, V_hat = opf.opf_run(
np.abs(voltages), power_P, power_Q
)
power_set = P_set + 1j * Q_set
power_factor = power_set.real / (np.abs(power_set) + 1e-7)
pmpp = power_set.real / pv["kVarRated"]
Expand All @@ -691,50 +691,42 @@ def run(self):
te = time.time()
logger.debug(f"OMOO takes {(te-ts)/60} (min)")

power_set_xr = xr.DataArray(
power_set,
coords={"equipment_ids": pv.loc[:, "name"]}
).groupby("equipment_ids").sum()
power_set_xr = (
xr.DataArray(power_set, coords={"equipment_ids": pv.loc[:, "name"]})
.groupby("equipment_ids")
.sum()
)

available_total_xr = available_power.groupby("equipment_ids").sum()

pv_settings = []
command_list = []
# We should test against the new interface
for i in range(len(power_set_xr)):
assert available_total_xr.equipment_ids[i] == power_set_xr.equipment_ids[i]
assert (
available_total_xr.equipment_ids[i] == power_set_xr.equipment_ids[i]
)
if np.isclose(available_total_xr[i], 0):
continue
#if np.isclose(pv["avai"][i], 0) and np.isclose(pv["avaiQ"][i], 0):
# continue
pv_settings.append((power_set_xr.equipment_ids.data[i], power_set_xr.values[i].real, power_set_xr.values[i].imag))
# command_list.append(
# Command(
# obj_name=pv.loc[i, "name"],
# obj_property="PF",
# val=str(np.sign(Q_set[i]) * power_factor[i]),
# # increases Q until it runs out of rating,
# # and then decreases P
# )
# )
# command_list.append(
# Command(
# obj_name=pv.loc[i, "name"],
# obj_property="%Pmpp",
# val=str(100 * pmpp[i]),
# # % of rating!
# )
# )
pv_settings.append(
(
power_set_xr.equipment_ids.data[i],
power_set_xr.values[i].real,
power_set_xr.values[i].imag,
)
)
command_list_obj = CommandList(__root__=command_list)
logger.debug(command_list_obj)
# Turn P_set and Q_set into commands
#self.pub_P_set.publish(command_list_obj.json())
self.pub_P_set.publish(json.dumps(pv_settings))
if set_power:
self.pub_P_set.publish(json.dumps(pv_settings))

logger.info("end time: " + str(datetime.now()))

import pdb; pdb.set_trace()
granted_time = h.helicsFederateRequestTime(self.vfed, 1000)
# There should be a HELICS way to do this? Set resolution?
previous_time = granted_time
while granted_time <= previous_time + 1:
granted_time = h.helicsFederateRequestTime(self.vfed, 1000)

self.destroy()

Expand Down

0 comments on commit 5fdf758

Please sign in to comment.