diff --git a/README.md b/README.md index 4159197..e9574fb 100644 --- a/README.md +++ b/README.md @@ -105,11 +105,17 @@ obs, info = env.reset(seed=2026) next_obs, reward, terminated, truncated, info = env.step(action) ``` +## Supported Examples + +The supported example entrypoints are the Python scripts described in `examples/README.md`. + +The notebooks are retained only as supplementary reference material and are no longer treated as the primary supported workflow. + ## Repository Structure - `rl_adn/`: package source code - `tests/`: smoke and domain validation tests -- `examples/`: script-first quickstart plus notebooks +- `examples/`: supported scripts plus supplementary notebooks - `docs/`: Sphinx documentation sources ## Highlights @@ -130,7 +136,7 @@ The library was originally released alongside the RL-ADN research paper on optim 1. Run `examples/quickstart_env.py` for the minimal package-backed environment flow. 2. Read the typed config surface through `rl_adn.make_env_config(...)`. 3. Try fixed and pooled topology scenarios before moving on to GNN-based experiments. -4. Use notebooks only as supplementary material after the script workflow is clear. +4. Treat notebooks as archival/supplementary material after the script workflow is clear. ## Current Limits diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000..630dd87 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,18 @@ +# Examples + +This directory is now script-first. + +## Supported scripts + +- `quickstart_env.py`: minimal environment reset/step walkthrough +- `custom_env_config.py`: typed config customization example +- `topology_scenarios.py`: Phase A topology scenario sampling example +- `training_smoke.py`: minimal algorithm/environment interaction smoke + +## Supplementary notebooks + +The notebooks in this directory are kept as supplementary reference material only. + +- They are not part of the supported automated verification surface. +- They may lag behind the current script-first APIs. +- Prefer the Python scripts above for reproducible and maintained workflows. diff --git a/pyproject.toml b/pyproject.toml index bd538f2..ccd7d64 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -98,3 +98,4 @@ indent-style = "space" [tool.ruff.lint.per-file-ignores] "tests/*_node_network_powerflow_test.py" = ["E402"] +"rl_adn/algorithms/utility.py" = ["F822"] diff --git a/rl_adn/algorithms/__init__.py b/rl_adn/algorithms/__init__.py index 39f3ec5..d345d0e 100644 --- a/rl_adn/algorithms/__init__.py +++ b/rl_adn/algorithms/__init__.py @@ -1,9 +1,6 @@ """RL algorithm exports for RL-ADN.""" -from rl_adn.algorithms.DDPG import AgentDDPG -from rl_adn.algorithms.PPO import AgentPPO -from rl_adn.algorithms.SAC import AgentSAC -from rl_adn.algorithms.TD3 import AgentTD3 +from importlib import import_module __all__ = [ "AgentDDPG", @@ -11,3 +8,19 @@ "AgentSAC", "AgentTD3", ] + +_LAZY_EXPORTS = { + "AgentDDPG": ("rl_adn.algorithms.DDPG", "AgentDDPG"), + "AgentPPO": ("rl_adn.algorithms.PPO", "AgentPPO"), + "AgentSAC": ("rl_adn.algorithms.SAC", "AgentSAC"), + "AgentTD3": ("rl_adn.algorithms.TD3", "AgentTD3"), +} + + +def __getattr__(name: str): + if name not in _LAZY_EXPORTS: + raise AttributeError(f"module 'rl_adn.algorithms' has no attribute {name!r}") + + module_name, attr_name = _LAZY_EXPORTS[name] + module = import_module(module_name) + return getattr(module, attr_name) diff --git a/rl_adn/algorithms/evaluation.py b/rl_adn/algorithms/evaluation.py new file mode 100644 index 0000000..465b220 --- /dev/null +++ b/rl_adn/algorithms/evaluation.py @@ -0,0 +1,52 @@ +from __future__ import annotations + +from typing import Any + +import torch + +from rl_adn.algorithms.env_api import reset_env, step_env + + +def get_episode_return(env, act, device) -> tuple[float, int, float, float, float, float, list[Any]]: + """Evaluate one episode with the current policy and return aggregate metrics.""" + env.train = False + episode_return = 0.0 + violation_time = 0 + reward_for_power = 0.0 + reward_for_good_action = 0.0 + reward_for_penalty = 0.0 + violation_value = 0.0 + state_list = [] + + state = reset_env(env) + for _ in range(env.episode_length): + state_tensor = torch.as_tensor((state,), device=device, dtype=torch.float32) + action_tensor = act(state_tensor) + action = action_tensor.detach().cpu().numpy()[0] + next_state, reward, done, info = step_env(env, action) + state_list.append(state) + + post_control_voltage = info["post_control_voltage_pu"] + for node_index in env.battery_nodes: + violation = min(0.0, 0.05 - abs(1.0 - post_control_voltage[node_index])) + if violation < 0: + violation_time += 1 + violation_value += violation + + reward_breakdown = info["reward_breakdown"] + reward_for_power += reward_breakdown["economic"] + reward_for_penalty += reward_breakdown["voltage_penalty"] + episode_return += reward + state = next_state + if done: + break + + return ( + episode_return, + violation_time, + violation_value, + reward_for_power, + reward_for_good_action, + reward_for_penalty, + state_list, + ) diff --git a/rl_adn/algorithms/replay.py b/rl_adn/algorithms/replay.py new file mode 100644 index 0000000..7808d1d --- /dev/null +++ b/rl_adn/algorithms/replay.py @@ -0,0 +1,218 @@ +from __future__ import annotations + +import math +import os +from typing import Tuple + +import torch +from torch import Tensor + +from rl_adn.algorithms.training_config import Config + + +class ReplayBuffer: + """Replay buffer for off-policy algorithms.""" + + def __init__( + self, + max_size: int, + state_dim: int, + action_dim: int, + gpu_id: int = 0, + num_seqs: int = 1, + if_use_per: bool = False, + args: Config | None = None, + ) -> None: + self.args = args or Config() + self.p = 0 + self.if_full = False + self.cur_size = 0 + self.add_size = 0 + self.add_item = None + self.max_size = max_size + self.num_seqs = num_seqs + self.device = torch.device(f"cuda:{gpu_id}" if (torch.cuda.is_available() and gpu_id >= 0) else "cpu") + + self.states = torch.empty((max_size, num_seqs, state_dim), dtype=torch.float32, device=self.device) + self.actions = torch.empty((max_size, num_seqs, action_dim), dtype=torch.float32, device=self.device) + self.rewards = torch.empty((max_size, num_seqs), dtype=torch.float32, device=self.device) + self.undones = torch.empty((max_size, num_seqs), dtype=torch.float32, device=self.device) + + self.if_use_per = if_use_per + if if_use_per: + self.sum_trees = [SumTree(buf_len=max_size) for _ in range(num_seqs)] + self.per_alpha = getattr(self.args, "per_alpha", 0.6) + self.per_beta = getattr(self.args, "per_beta", 0.4) + else: + self.sum_trees = None + self.per_alpha = None + self.per_beta = None + + def update(self, items: Tuple[Tensor, ...]) -> None: + self.add_item = items + states, actions, rewards, undones = items + assert states.shape[1:] == (self.args.num_envs, self.args.state_dim) + assert actions.shape[1:] == (self.args.num_envs, self.args.action_dim) + assert rewards.shape[1:] == (self.args.num_envs,) + assert undones.shape[1:] == (self.args.num_envs,) + self.add_size = rewards.shape[0] + + new_pointer = self.p + self.add_size + if new_pointer > self.max_size: + self.if_full = True + split_index = self.max_size - self.p + new_pointer -= self.max_size + + self.states[self.p : self.max_size], self.states[0:new_pointer] = states[:split_index], states[-new_pointer:] + self.actions[self.p : self.max_size], self.actions[0:new_pointer] = actions[:split_index], actions[-new_pointer:] + self.rewards[self.p : self.max_size], self.rewards[0:new_pointer] = rewards[:split_index], rewards[-new_pointer:] + self.undones[self.p : self.max_size], self.undones[0:new_pointer] = undones[:split_index], undones[-new_pointer:] + else: + self.states[self.p : new_pointer] = states + self.actions[self.p : new_pointer] = actions + self.rewards[self.p : new_pointer] = rewards + self.undones[self.p : new_pointer] = undones + + if self.if_use_per and self.sum_trees is not None: + data_ids = torch.arange(self.p, new_pointer, dtype=torch.long, device=self.device) + if new_pointer > self.max_size: + data_ids = torch.fmod(data_ids, self.max_size) + for sum_tree in self.sum_trees: + sum_tree.update_ids(data_ids=data_ids.cpu(), prob=10.0) + + self.p = new_pointer + self.cur_size = self.max_size if self.if_full else self.p + + def sample(self, batch_size: int) -> Tuple[Tensor, Tensor, Tensor, Tensor, Tensor]: + sample_len = self.cur_size - 1 + indices = torch.randint(sample_len * self.num_seqs, size=(batch_size,), requires_grad=False) + time_indices = torch.fmod(indices, sample_len) + seq_indices = torch.div(indices, sample_len, rounding_mode="floor") + return ( + self.states[time_indices, seq_indices], + self.actions[time_indices, seq_indices], + self.rewards[time_indices, seq_indices], + self.undones[time_indices, seq_indices], + self.states[time_indices + 1, seq_indices], + ) + + def sample_for_per(self, batch_size: int) -> Tuple[Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor]: + beg = -self.max_size + end = (self.cur_size - self.max_size) if self.cur_size < self.max_size else -1 + + assert batch_size % self.num_seqs == 0 + sub_batch_size = batch_size // self.num_seqs + is_indices = [] + is_weights = [] + for env_index, sum_tree in enumerate(self.sum_trees or []): + sampled_indices, sampled_weights = sum_tree.important_sampling(batch_size, beg, end, self.per_beta) + is_indices.append(sampled_indices + sub_batch_size * env_index) + is_weights.append(sampled_weights) + + index_tensor = torch.hstack(is_indices).to(self.device) + weight_tensor = torch.hstack(is_weights).to(self.device) + + time_indices = torch.fmod(index_tensor, self.cur_size) + seq_indices = torch.div(index_tensor, self.cur_size, rounding_mode="floor") + return ( + self.states[time_indices, seq_indices], + self.actions[time_indices, seq_indices], + self.rewards[time_indices, seq_indices], + self.undones[time_indices, seq_indices], + self.states[time_indices + 1, seq_indices], + weight_tensor, + index_tensor, + ) + + def td_error_update_for_per(self, is_indices: Tensor, td_error: Tensor) -> None: + prob = td_error.clamp(1e-8, 10).pow(self.per_alpha).squeeze(-1) + batch_size = td_error.shape[0] + sub_batch_size = batch_size // self.num_seqs + for env_index, sum_tree in enumerate(self.sum_trees or []): + start = env_index * sub_batch_size + end = start + sub_batch_size + sum_tree.update_ids(is_indices[start:end].cpu(), prob[start:end].cpu()) + + def save_or_load_history(self, cwd: str, if_save: bool) -> None: + item_names = ( + (self.states, "states"), + (self.actions, "actions"), + (self.rewards, "rewards"), + (self.undones, "undones"), + ) + + if if_save: + for item, name in item_names: + if self.cur_size == self.p: + buffer_item = item[: self.cur_size] + else: + buffer_item = torch.vstack((item[self.p : self.cur_size], item[0 : self.p])) + torch.save(buffer_item, f"{cwd}/replay_buffer_{name}.pth") + return + + expected_files = [f"{cwd}/replay_buffer_{name}.pth" for _, name in item_names] + if not all(os.path.isfile(path) for path in expected_files): + return + + max_sizes = [] + for item, name in item_names: + buffer_item = torch.load(f"{cwd}/replay_buffer_{name}.pth") + max_size = buffer_item.shape[0] + item[:max_size] = buffer_item + max_sizes.append(max_size) + assert all(size == max_sizes[0] for size in max_sizes) + self.cur_size = self.p = max_sizes[0] + self.if_full = self.cur_size == self.max_size + + +class SumTree: + """Binary tree used for prioritized experience replay.""" + + def __init__(self, buf_len: int) -> None: + self.buf_len = buf_len + self.max_len = (buf_len - 1) + buf_len + self.depth = math.ceil(math.log2(self.max_len)) + self.tree = torch.zeros(self.max_len, dtype=torch.float32) + + def update_id(self, data_id: int, prob: float = 10.0) -> None: + tree_id = data_id + self.buf_len - 1 + delta = prob - self.tree[tree_id] + self.tree[tree_id] = prob + for _ in range(self.depth - 2): + tree_id = (tree_id - 1) // 2 + self.tree[tree_id] += delta + + def update_ids(self, data_ids: Tensor, prob: Tensor = 10.0) -> None: + leaf_ids = data_ids + self.buf_len - 1 + self.tree[leaf_ids] = prob + for _ in range(self.depth - 2): + parent_ids = torch.div(leaf_ids - 1, 2, rounding_mode="floor").unique() + left_ids = parent_ids * 2 + 1 + right_ids = left_ids + 1 + self.tree[parent_ids] = self.tree[left_ids] + self.tree[right_ids] + leaf_ids = parent_ids + + def get_leaf_id_and_value(self, value: float) -> Tuple[int, float]: + parent_id = 0 + for _ in range(self.depth - 2): + left_id = min(2 * parent_id + 1, self.max_len - 1) + right_id = left_id + 1 + if value <= self.tree[left_id]: + parent_id = left_id + else: + value -= self.tree[left_id] + parent_id = right_id + return parent_id, float(self.tree[parent_id]) + + def important_sampling(self, batch_size: int, beg: int, end: int, per_beta: float) -> Tuple[Tensor, Tensor]: + values = (torch.arange(batch_size) + torch.rand(batch_size)) * (self.tree[0] / batch_size) + leaf_ids, leaf_values = list(zip(*[self.get_leaf_id_and_value(v) for v in values])) + leaf_ids = torch.tensor(leaf_ids, dtype=torch.long) + leaf_values = torch.tensor(leaf_values, dtype=torch.float32) + + indices = leaf_ids - (self.buf_len - 1) + assert indices.max() < self.buf_len + + prob_ary = leaf_values / self.tree[beg:end].min() + weights = torch.pow(prob_ary, -per_beta) + return indices, weights diff --git a/rl_adn/algorithms/torch_utils.py b/rl_adn/algorithms/torch_utils.py new file mode 100644 index 0000000..607fec8 --- /dev/null +++ b/rl_adn/algorithms/torch_utils.py @@ -0,0 +1,26 @@ +from __future__ import annotations + +from typing import Iterable + +import torch +from torch import nn + + +def get_optim_param(optimizer: torch.optim.Optimizer) -> list[torch.Tensor]: + params_list: list[torch.Tensor] = [] + for params_dict in optimizer.state_dict()["state"].values(): + params_list.extend(value for value in params_dict.values() if isinstance(value, torch.Tensor)) + return params_list + + +def build_mlp(dims: Iterable[int]) -> nn.Sequential: + dims = list(dims) + if len(dims) < 2: + raise ValueError("build_mlp expects at least an input and output dimension") + + layers: list[nn.Module] = [] + for index in range(len(dims) - 1): + layers.append(nn.Linear(dims[index], dims[index + 1])) + if index < len(dims) - 2: + layers.append(nn.ReLU()) + return nn.Sequential(*layers) diff --git a/rl_adn/algorithms/training_config.py b/rl_adn/algorithms/training_config.py new file mode 100644 index 0000000..3f5df28 --- /dev/null +++ b/rl_adn/algorithms/training_config.py @@ -0,0 +1,97 @@ +from __future__ import annotations + +import os +from pprint import pprint +from typing import Any + +import numpy as np + + +class Config: + """Configuration container for RL-ADN training workflows.""" + + def __init__(self, agent_class=None, env_class=None, env_args: dict[str, Any] | None = None): + self.agent_class = agent_class + self.env_class = env_class + self.env_args = dict(env_args or {}) + self.if_off_policy = self.get_if_off_policy() + + default_env_args = { + "env_name": None, + "num_envs": 1, + "max_step": 96, + "state_dim": None, + "action_dim": None, + "if_discrete": None, + } + default_env_args.update(self.env_args) + self.env_args = default_env_args + self.env_name = self.env_args["env_name"] + self.num_envs = self.env_args["num_envs"] + self.max_step = self.env_args["max_step"] + self.state_dim = self.env_args["state_dim"] + self.action_dim = self.env_args["action_dim"] + self.if_discrete = self.env_args["if_discrete"] + + self.gamma = 0.99 + self.reward_scale = 1.0 + self.net_dims = (64, 32) + self.learning_rate = 6e-5 + self.clip_grad_norm = 3.0 + self.state_value_tau = 0.0 + self.soft_update_tau = 5e-3 + + if self.if_off_policy: + self.batch_size = 64 + self.target_step = 512 + self.buffer_size = int(1e6) + self.repeat_times = 1.0 + self.if_use_per = False + else: + self.batch_size = 128 + self.target_step = 2048 + self.buffer_size = None + self.repeat_times = 8.0 + self.if_use_vtrace = False + + self.random_seed = 0 + self.num_episode = 2000 + self.gpu_id = 0 + self.num_workers = 2 + self.num_threads = 8 + self.learner_gpus = 0 + + self.run_name = None + self.cwd = None + self.if_remove = True + self.train = True + + def init_before_training(self) -> None: + import torch + + np.random.seed(self.random_seed) + torch.manual_seed(self.random_seed) + torch.set_num_threads(self.num_threads) + torch.set_default_dtype(torch.float32) + + if self.cwd is None: + agent_name = self.agent_class.__name__[5:] if self.agent_class else "agent" + run_name = self.run_name or "default" + self.cwd = f"./{agent_name}/{run_name}" + + if self.if_remove: + import shutil + + shutil.rmtree(self.cwd, ignore_errors=True) + os.makedirs(self.cwd, exist_ok=True) + + def get_if_off_policy(self) -> bool: + agent_name = self.agent_class.__name__ if self.agent_class else "" + on_policy_names = ("SARSA", "VPG", "A2C", "A3C", "TRPO", "PPO", "MPO") + return all(agent_name.find(name) == -1 for name in on_policy_names) + + def print(self) -> None: + pprint(vars(self)) + + def to_dict(self) -> dict[str, Any]: + return vars(self) diff --git a/rl_adn/algorithms/utility.py b/rl_adn/algorithms/utility.py index 36709f8..26765e1 100644 --- a/rl_adn/algorithms/utility.py +++ b/rl_adn/algorithms/utility.py @@ -1,662 +1,30 @@ -import math -import os -from typing import Tuple - -import numpy as np -import torch -from torch import Tensor, nn - -from rl_adn.algorithms.env_api import reset_env, step_env - - -class Config: - """ - Configuration class for setting up and managing parameters for the agent and environment. - - Attributes: - num_envs (int): Number of environments. - agent_class (class): Class of the agent. - if_off_policy (bool): Indicates whether the DRL algorithm is off-policy or on-policy. - env_class (class): Class of the environment. - env_args (dict): Arguments for the environment. - env_name (str): Name of the environment. - max_step (int): Maximum number of steps in an episode. - state_dim (int): Dimension of the state vector. - action_dim (int): Dimension of the action vector. - if_discrete (bool): Indicates if the action space is discrete. - gamma (float): Discount factor for future rewards. - reward_scale (float): Scale of the reward. - net_dims (tuple): Dimensions of the MLP layers. - learning_rate (float): Learning rate for network updates. - clip_grad_norm (float): Gradient clipping norm. - state_value_tau (float): Tau for normalizing state and value. - soft_update_tau (float): Tau for soft target update. - batch_size (int): Batch size for training. - target_step (int): Number of steps for target update. - buffer_size (int): Size of the replay buffer. - repeat_times (float): Number of times to update the network with the replay buffer. - if_use_per (bool): Indicates if PER (Prioritized Experience Replay) is used. - if_use_vtrace (bool): Indicates if V-trace is used. - random_seed (int): Random seed for reproducibility. - num_episode (int): Number of episodes for training. - gpu_id (int): GPU ID for training. - num_workers (int): Number of workers for data collection. - num_threads (int): Number of threads for PyTorch. - learner_gpus (int): GPU ID for the learner. - run_name (str): Name of the run for data storage. - cwd (str): Current working directory. - if_remove (bool): Flag to remove the current working directory. - train (bool): Flag to indicate training mode. - - Methods: - init_before_training(): Initializes settings before training starts. - get_if_off_policy(): Determines if the agent is off-policy based on its name. - print(): Prints the configuration in a readable format. - to_dict(): Converts the configuration to a dictionary. - """ - - def __init__(self, agent_class=None, env_class=None, env_args=None): - self.num_envs = None - self.agent_class = agent_class # agent = agent_class(...) - self.if_off_policy = self.get_if_off_policy() # whether off-policy or on-policy of DRL algorithm - - """Argument of environment""" - self.env_class = env_class # env = env_class(**env_args) - self.env_args = env_args # env = env_class(**env_args) - if env_args is None: # dummy env_args - env_args = { - "env_name": None, - "num_envs": 1, - "max_step": 96, - "state_dim": None, - "action_dim": None, - "if_discrete": None, - } - env_args.setdefault("num_envs", 1) # `num_envs=1` in default in single env. - env_args.setdefault("max_step", 96) # `max_step=12345` in default, which is a large enough value. - self.env_name = env_args["env_name"] # the name of environment. Be used to set 'cwd'. - self.num_envs = env_args["num_envs"] # the number of sub envs in vectorized env. `num_envs=1` in single env. - self.max_step = env_args["max_step"] # the max step number of an episode. 'set as 12345 in default. - self.state_dim = env_args["state_dim"] # vector dimension (feature number) of state - self.action_dim = env_args["action_dim"] # vector dimension (feature number) of action - self.if_discrete = env_args["if_discrete"] # discrete or continuous action space - """Arguments for reward shaping""" - self.gamma = 0.99 # discount factor of future rewards - self.reward_scale = 2**0 # an approximate target reward usually be closed to 256 - - """Arguments for training""" - self.net_dims = (64, 32) # the middle layer dimension of MLP (MultiLayer Perceptron) - self.learning_rate = 6e-5 # the learning rate for network updating - self.clip_grad_norm = 3.0 # 0.1 ~ 4.0, clip the gradient after normalization - self.state_value_tau = 0 # the tau of normalize for value and state `std = (1-std)*std + tau*std` - self.soft_update_tau = 5e-3 # 2 ** -8 ~= 5e-3. the tau of soft target update `net = (1-tau)*net + tau*net1` - if self.if_off_policy: # off-policy - self.batch_size = int(64) # num of transitions sampled from replay buffer. - self.target_step = int(512) # collect horizon_len step while exploring, then update networks - self.buffer_size = int(1e6) # ReplayBuffer size. First in first out for off-policy. - self.repeat_times = 1.0 # repeatedly update network using ReplayBuffer to keep critic's loss small - self.if_use_per = False # use PER (Prioritized Experience Replay) for sparse reward - else: # on-policy - self.batch_size = int(128) # num of transitions sampled from replay buffer. - self.target_step = int(2048) # collect horizon_len step while exploring, then update network - self.buffer_size = None # ReplayBuffer size. Empty the ReplayBuffer for on-policy. - self.repeat_times = 8.0 # repeatedly update network using ReplayBuffer to keep critic's loss small - self.if_use_vtrace = False # use V-trace + GAE (Generalized Advantage Estimation) for sparse reward - self.random_seed = 521 - self.num_episode = 2000 - self.buffer_size = 500000 # capacity of replay buffer - """Arguments for device""" - self.gpu_id = int(0) # `int` means the ID of single GPU, -1 means CPU - self.num_workers = 2 # rollout workers number pre GPU (adjust it to get high GPU usage) - self.num_threads = 8 # cpu_num for pytorch, `torch.set_num_threads(self.num_threads)` - self.random_seed = 0 # initialize random seed in self.init_before_training() - self.learner_gpus = 0 # `int` means the ID of single GPU, -1 means CPU - - """arguments for creating data storage directory""" - - self.run_name = None - """Arguments for save and plot issues""" - self.cwd = None # current work directory. None means set automatically - self.if_remove = True # remove the cwd folder? (True, False, None:ask me) - self.train = True - - def init_before_training(self): - np.random.seed(self.random_seed) - torch.manual_seed(self.random_seed) - torch.set_num_threads(self.num_threads) - torch.set_default_dtype(torch.float32) - if self.cwd is None: - agent_name = self.agent_class.__name__[5:] - self.cwd = f"./{agent_name}/{self.run_name}" - - """remove history""" - if self.if_remove is None: - self.if_remove = bool(input(f"| Arguments PRESS 'y' to REMOVE: {self.cwd}? ") == "y") - if self.if_remove: - import shutil - - shutil.rmtree(self.cwd, ignore_errors=True) - print(f"| Arguments Remove cwd: {self.cwd}") - else: - print(f"| Arguments Keep cwd: {self.cwd}") - os.makedirs(self.cwd, exist_ok=True) - - def get_if_off_policy(self) -> bool: - agent_name = self.agent_class.__name__ if self.agent_class else "" - on_policy_names = ("SARSA", "VPG", "A2C", "A3C", "TRPO", "PPO", "MPO") - return all([agent_name.find(s) == -1 for s in on_policy_names]) - - def print(self): - from pprint import pprint - - pprint(vars(self)) # prints out args in a neat, readable format - - def to_dict(self): - return vars(self) - - -def get_optim_param(optimizer: torch.optim) -> list: # backup - """ - Extracts parameters from the optimizer state. - - Args: - optimizer (torch.optim): The optimizer from which to extract parameters. - - Returns: - list: A list of parameters extracted from the optimizer state. - """ - params_list = [] - for params_dict in optimizer.state_dict()["state"].values(): - params_list.extend([t for t in params_dict.values() if isinstance(t, torch.Tensor)]) - return params_list - - -def build_mlp(dims: [int]) -> nn.Sequential: # MLP (MultiLayer Perceptron) - """ - Builds a Multi-Layer Perceptron (MLP) network. - - Args: - dims (list of int): A list containing the dimensions of each layer in the MLP. - - Returns: - nn.Sequential: The constructed MLP network. - """ - net_list = list() - for i in range(len(dims) - 1): - net_list.extend([nn.Linear(dims[i], dims[i + 1]), nn.ReLU()]) - del net_list[-1] # remove the activation of output layer - return nn.Sequential(*net_list) - - -class ReplayBuffer: # for off-policy - """ - Replay Buffer for storing and sampling experiences for off-policy reinforcement learning algorithms. - - Attributes: - max_size (int): Maximum size of the buffer. - state_dim (int): Dimension of the state space. - action_dim (int): Dimension of the action space. - gpu_id (int): GPU ID for storing the buffer. - num_seqs (int): Number of sequences in the buffer. - if_use_per (bool): Flag to use Prioritized Experience Replay. - args (Config): Configuration object with additional parameters. - - Methods: - update(items): Updates the buffer with new experiences. - sample(batch_size): Samples a batch of experiences from the buffer. - sample_for_per(batch_size): Samples a batch with prioritization. - td_error_update_for_per(is_indices, td_error): Updates the priorities based on TD error. - save_or_load_history(cwd, if_save): Saves or loads the buffer history. - """ - - def __init__( - self, - max_size: int, - state_dim: int, - action_dim: int, - gpu_id: int = 0, - num_seqs: int = 1, - if_use_per: bool = False, - args: Config = Config(), - ): - self.p = 0 # pointer - self.if_full = False - self.cur_size = 0 - self.add_size = 0 - self.add_item = None - self.max_size = max_size - self.num_seqs = num_seqs - self.device = torch.device(f"cuda:{gpu_id}" if (torch.cuda.is_available() and (gpu_id >= 0)) else "cpu") - self.args = args - """The struction of ReplayBuffer (for examples, num_seqs = num_workers * num_envs == 2*4 = 8 - ReplayBuffer: - worker0 for env0: sequence of sub_env0.0 self.states = Tensor[s, s, ..., s, ..., s] - self.actions = Tensor[a, a, ..., a, ..., a] - self.rewards = Tensor[r, r, ..., r, ..., r] - self.undones = Tensor[d, d, ..., d, ..., d] - <-----max_size-----> - <-cur_size-> - ↑ pointer - sequence of sub_env0.1 s, s, ..., s a, a, ..., a r, r, ..., r d, d, ..., d - sequence of sub_env0.2 s, s, ..., s a, a, ..., a r, r, ..., r d, d, ..., d - sequence of sub_env0.3 s, s, ..., s a, a, ..., a r, r, ..., r d, d, ..., d - worker1 for env1: sequence of sub_env1.0 s, s, ..., s a, a, ..., a r, r, ..., r d, d, ..., d - sequence of sub_env1.1 s, s, ..., s a, a, ..., a r, r, ..., r d, d, ..., d - sequence of sub_env1.2 s, s, ..., s a, a, ..., a r, r, ..., r d, d, ..., d - sequence of sub_env1.3 s, s, ..., s a, a, ..., a r, r, ..., r d, d, ..., d - - D: done=True - d: done=False - sequence of transition: s-a-r-d, s-a-r-d, s-a-r-D s-a-r-d, s-a-r-d, s-a-r-d, s-a-r-d, s-a-r-D s-a-r-d, ... - <------trajectory-------> <----------trajectory---------------------> <----------- - """ - self.states = torch.empty((max_size, num_seqs, state_dim), dtype=torch.float32, device=self.device) - self.actions = torch.empty((max_size, num_seqs, action_dim), dtype=torch.float32, device=self.device) - self.rewards = torch.empty((max_size, num_seqs), dtype=torch.float32, device=self.device) - self.undones = torch.empty((max_size, num_seqs), dtype=torch.float32, device=self.device) - - self.if_use_per = if_use_per - if if_use_per: - self.sum_trees = [SumTree(buf_len=max_size) for _ in range(num_seqs)] - self.per_alpha = getattr(args, "per_alpha", 0.6) # alpha = (Uniform:0, Greedy:1) - self.per_beta = getattr(args, "per_beta", 0.4) # alpha = (Uniform:0, Greedy:1) - """PER. Prioritized Experience Replay. Section 4 - alpha, beta = 0.7, 0.5 for rank-based variant - alpha, beta = 0.6, 0.4 for proportional variant - """ - else: - self.sum_trees = None - self.per_alpha = None - self.per_beta = None - - def update(self, items: Tuple[Tensor, ...]): - """ - Updates the replay buffer with new experience tuples. - - Args: - items (Tuple[Tensor, ...]): A tuple containing tensors of states, actions, rewards, and undones. - Each tensor should have a shape that matches the expected dimensions - for states, actions, rewards, and undones respectively. - - Description: - This method updates the replay buffer with new experiences. It handles the buffer's internal - pointers and ensures that new data is added correctly, even when the buffer is full. If the buffer - is full, it starts overwriting the oldest data. In case of using Prioritized Experience Replay (PER), - it updates the sum trees with new priorities. - """ - self.add_item = items - states, actions, rewards, undones = items - assert states.shape[1:] == (self.args.num_envs, self.args.state_dim) - assert actions.shape[1:] == (self.args.num_envs, self.args.action_dim) - assert rewards.shape[1:] == (self.args.num_envs,) - assert undones.shape[1:] == (self.args.num_envs,) - self.add_size = rewards.shape[0] - - p = self.p + self.add_size # pointer - if p > self.max_size: - self.if_full = True - p0 = self.p - p1 = self.max_size - p2 = self.max_size - self.p - p = p - self.max_size - - self.states[p0:p1], self.states[0:p] = states[:p2], states[-p:] - self.actions[p0:p1], self.actions[0:p] = actions[:p2], actions[-p:] - self.rewards[p0:p1], self.rewards[0:p] = rewards[:p2], rewards[-p:] - self.undones[p0:p1], self.undones[0:p] = undones[:p2], undones[-p:] - else: - self.states[self.p : p] = states - self.actions[self.p : p] = actions - self.rewards[self.p : p] = rewards - self.undones[self.p : p] = undones - - if self.if_use_per: - """data_ids for single env""" - data_ids = torch.arange(self.p, p, dtype=torch.long, device=self.device) - if p > self.max_size: - data_ids = torch.fmod(data_ids, self.max_size) - - """apply data_ids for vectorized env""" - for sum_tree in self.sum_trees: - sum_tree.update_ids(data_ids=data_ids.cpu(), prob=10.0) - - self.p = p - self.cur_size = self.max_size if self.if_full else self.p - - def sample(self, batch_size: int) -> Tuple[Tensor, Tensor, Tensor, Tensor, Tensor]: - """ - Samples a batch of experiences from the replay buffer. - - Args: - batch_size (int): The size of the batch to sample. - - Returns: - Tuple[Tensor, Tensor, Tensor, Tensor, Tensor]: A tuple containing batches of states, actions, rewards, - undones, and next_states. Each tensor in the tuple has - dimensions corresponding to the batch size. - - Description: - This method randomly samples a batch of experiences from the replay buffer. It is typically used in - off-policy algorithms where random sampling of experiences is required for training the agent. - """ - - sample_len = self.cur_size - 1 - - ids = torch.randint(sample_len * self.num_seqs, size=(batch_size,), requires_grad=False) - ids0 = torch.fmod(ids, sample_len) # ids % sample_len - ids1 = torch.div(ids, sample_len, rounding_mode="floor") # ids // sample_len - - return ( - self.states[ids0, ids1], - self.actions[ids0, ids1], - self.rewards[ids0, ids1], - self.undones[ids0, ids1], - self.states[ids0 + 1, ids1], - ) # next_state - - def sample_for_per(self, batch_size: int) -> Tuple[Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor]: - """ - Samples a batch of experiences using Prioritized Experience Replay. - - Args: - batch_size (int): The size of the batch to sample. - - Returns: - Tuple[Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor]: A tuple containing batches of states, - actions, rewards, undones, next_states, - importance sampling weights, and indices. - Each tensor in the tuple has dimensions - corresponding to the batch size. - - Description: - This method samples experiences using Prioritized Experience Replay (PER). It uses importance sampling - to give more priority to experiences with higher expected learning value. This is particularly useful - in scenarios where some experiences may be more significant than others for learning. - """ - beg = -self.max_size - end = (self.cur_size - self.max_size) if (self.cur_size < self.max_size) else -1 - - """get is_indices, is_weights""" - is_indices: list = [] - is_weights: list = [] - - assert batch_size % self.num_seqs == 0 - sub_batch_size = batch_size // self.num_seqs - for env_i in range(self.num_seqs): - sum_tree = self.sum_trees[env_i] - _is_indices, _is_weights = sum_tree.important_sampling(batch_size, beg, end, self.per_beta) - is_indices.append(_is_indices + sub_batch_size * env_i) - is_weights.append(_is_weights) - - is_indices: Tensor = torch.hstack(is_indices).to(self.device) - is_weights: Tensor = torch.hstack(is_weights).to(self.device) - - ids0 = torch.fmod(is_indices, self.cur_size) # is_indices % sample_len - ids1 = torch.div(is_indices, self.cur_size, rounding_mode="floor") # is_indices // sample_len - return ( - self.states[ids0, ids1], - self.actions[ids0, ids1], - self.rewards[ids0, ids1], - self.undones[ids0, ids1], - self.states[ids0 + 1, ids1], # next_state - is_weights, # important sampling weights - is_indices, # important sampling indices - ) - - def td_error_update_for_per(self, is_indices: Tensor, td_error: Tensor): # td_error = (q-q).detach_().abs() - """ - Updates the priorities in the sum trees based on the TD error. - - Args: - is_indices (Tensor): Tensor containing indices of sampled experiences. - td_error (Tensor): Tensor containing the Temporal Difference (TD) error for each sampled experience. - - Description: - This method updates the priorities in the sum trees for each experience based on the provided TD error. - It is an essential part of the Prioritized Experience Replay mechanism, ensuring that experiences - with higher TD error (and thus potentially higher learning value) have a higher chance of being sampled. - """ - prob = td_error.clamp(1e-8, 10).pow(self.per_alpha).squeeze(-1) - - # self.sum_tree.update_ids(is_indices.cpu(), prob.cpu()) - batch_size = td_error.shape[0] - sub_batch_size = batch_size // self.num_seqs - for env_i in range(self.num_seqs): - sum_tree = self.sum_trees[env_i] - slice_i = env_i * sub_batch_size - slice_j = slice_i + sub_batch_size - - sum_tree.update_ids(is_indices[slice_i:slice_j].cpu(), prob[slice_i:slice_j].cpu()) - - def save_or_load_history(self, cwd: str, if_save: bool): - """ - Saves or loads the replay buffer history to/from disk. - - Args: - cwd (str): The current working directory where the buffer history will be saved or loaded from. - if_save (bool): A flag indicating whether to save (True) or load (False) the buffer history. - - Description: - This method either saves the current state of the replay buffer to disk or loads it from disk. - This is useful for persisting the replay buffer across different training sessions or for - transferring the buffer state between different instances. - """ - item_names = ( - (self.states, "states"), - (self.actions, "actions"), - (self.rewards, "rewards"), - (self.undones, "undones"), - ) - - if if_save: - for item, name in item_names: - if self.cur_size == self.p: - buf_item = item[: self.cur_size] - else: - buf_item = torch.vstack((item[self.p : self.cur_size], item[0 : self.p])) - file_path = f"{cwd}/replay_buffer_{name}.pth" - print(f"| buffer.save_or_load_history(): Save {file_path}") - torch.save(buf_item, file_path) - - elif all([os.path.isfile(f"{cwd}/replay_buffer_{name}.pth") for item, name in item_names]): - max_sizes = [] - for item, name in item_names: - file_path = f"{cwd}/replay_buffer_{name}.pth" - print(f"| buffer.save_or_load_history(): Load {file_path}") - buf_item = torch.load(file_path) - - max_size = buf_item.shape[0] - item[:max_size] = buf_item - max_sizes.append(max_size) - assert all([max_size == max_sizes[0] for max_size in max_sizes]) - self.cur_size = self.p = max_sizes[0] - self.if_full = self.cur_size == self.max_size - - -class SumTree: - """ - Binary Search Tree for efficient sampling in Prioritized Experience Replay. - - Attributes: - buf_len (int): Length of the buffer. - max_len (int): Maximum length of the tree. - depth (int): Depth of the tree. - tree (Tensor): Tensor representing the tree structure. - - Methods: - update_id(data_id, prob): Updates a single node in the tree. - update_ids(data_ids, prob): Updates multiple nodes in the tree. - get_leaf_id_and_value(v): Retrieves the leaf ID and value for a given value. - important_sampling(batch_size, beg, end, per_beta): Performs important sampling for a batch. - """ - - def __init__(self, buf_len: int): - """ - Initializes the SumTree object. - - Args: - buf_len (int): The length of the buffer for which this SumTree is being used. - - Description: - This method initializes a SumTree data structure. The SumTree is a binary tree where each node's - value is the sum of its children's values. This structure is particularly useful for efficiently - implementing Prioritized Experience Replay (PER) in reinforcement learning. - """ - self.buf_len = buf_len # replay buffer len - self.max_len = (buf_len - 1) + buf_len # parent_nodes_num + leaf_nodes_num - self.depth = math.ceil(math.log2(self.max_len)) - - self.tree = torch.zeros(self.max_len, dtype=torch.float32) - - def update_id(self, data_id: int, prob=10): # 10 is max_prob - """ - Updates the priority of a single data point in the SumTree. - - Args: - data_id (int): The index of the data point in the buffer. - prob (float, optional): The new priority value for the data point. Defaults to 10, which is considered the maximum priority. - - Description: - This method updates the priority of a single data point in the SumTree. It adjusts the values in the tree - to maintain the sum property after the update. This is used in PER to adjust the sampling probability of experiences. - """ - tree_id = data_id + self.buf_len - 1 - - delta = prob - self.tree[tree_id] - self.tree[tree_id] = prob - - for depth in range(self.depth - 2): # propagate the change through tree - tree_id = (tree_id - 1) // 2 # faster than the recursive loop - self.tree[tree_id] += delta - - def update_ids(self, data_ids: Tensor, prob: Tensor = 10.0): # 10 is max_prob - """ - Updates the priorities of multiple data points in the SumTree. - - Args: - data_ids (Tensor): A tensor of indices of the data points in the buffer. - prob (Tensor, optional): A tensor of new priority values for the data points. Defaults to 10 for each, which is considered the maximum priority. - - Description: - This method updates the priorities of multiple data points in the SumTree simultaneously. It ensures that - the sum property of the tree is maintained after the updates. This method is typically used in batch updates - in PER. - """ - l_ids = data_ids + self.buf_len - 1 - - self.tree[l_ids] = prob - for depth in range(self.depth - 2): # propagate the change through tree - p_ids = torch.div(l_ids - 1, 2, rounding_mode="floor").unique() # parent indices - l_ids = p_ids * 2 + 1 # left children indices - r_ids = l_ids + 1 # right children indices - self.tree[p_ids] = self.tree[l_ids] + self.tree[r_ids] - - l_ids = p_ids - - def get_leaf_id_and_value(self, v) -> Tuple[int, float]: - """Retrieve the leaf node index and priority value for a cumulative-sum lookup.""" - - p_id = 0 # the leaf's parent node - - for depth in range(self.depth - 2): # propagate the change through tree - l_id = min(2 * p_id + 1, self.max_len - 1) # the leaf's left node - r_id = l_id + 1 # the leaf's right node - if v <= self.tree[l_id]: - p_id = l_id - else: - v -= self.tree[l_id] - p_id = r_id - return p_id, self.tree[p_id] # leaf_id and leaf_value - - def important_sampling(self, batch_size: int, beg: int, end: int, per_beta: float) -> Tuple[Tensor, Tensor]: - """ - Performs important sampling to select indices and compute weights for experiences. - - Args: - batch_size (int): The number of samples to draw. - beg (int): The beginning index for sampling. - end (int): The ending index for sampling. - per_beta (float): The beta parameter for PER, controlling the degree of importance sampling. - - Returns: - Tuple[Tensor, Tensor]: A tuple containing tensors of indices and corresponding weights for the sampled experiences. - - Description: - This method performs important sampling based on the priorities in the SumTree. It is used in PER to - select experiences non-uniformly, giving more priority to experiences with higher expected learning value. - """ - # get random values for searching indices with proportional prioritization - values = (torch.arange(batch_size) + torch.rand(batch_size)) * (self.tree[0] / batch_size) - - # get proportional prioritization - leaf_ids, leaf_values = list(zip(*[self.get_leaf_id_and_value(v) for v in values])) - leaf_ids = torch.tensor(leaf_ids, dtype=torch.long) - leaf_values = torch.tensor(leaf_values, dtype=torch.float32) - - indices = leaf_ids - (self.buf_len - 1) - if indices.min() < 0: - print(f"the wrong indice is{indices.min()}") - print(f"the whole indices is {indices}") - # assert 0 <= indices.min() - assert indices.max() < self.buf_len - - prob_ary = leaf_values / self.tree[beg:end].min() - weights = torch.pow(prob_ary, -per_beta) - return indices, weights - - -def get_episode_return(env, act, device): - """ - Calculates the return of an episode. - - Args: - env: The environment to interact with. - act: The action function to use. - device: The device to perform computations on. - - Returns: - Tuple containing episode return, violation time, violation value, rewards for power, good actions, and penalties, and the list of states. - """ - env.train = False - episode_return = 0.0 - violation_time = 0 - reward_for_power = 0 - reward_for_good_action = 0 - reward_for_penalty = 0 - state_list = [] - state = reset_env(env) - print(f"the year:{env.year},month:{env.month},day:{env.day} is used for testing this episode") - - violation_value = 0.0 - - for i in range(env.episode_length): - s_tensor = torch.as_tensor((state,), device=device, dtype=torch.float) - a_tensor = act(s_tensor) - - action = a_tensor.detach().cpu().numpy()[0] - next_state, reward, done, info = step_env(env, action) - state_list.append(state) - - post_control_voltage = info["post_control_voltage_pu"] - for node_index in env.battery_nodes: - violation = min(0, 0.05 - abs(1.0 - post_control_voltage[node_index])) - if violation < 0: - violation_time += 1 - violation_value += violation - - reward_for_power += info["reward_breakdown"]["economic"] - reward_for_good_action += 0 - reward_for_penalty += info["reward_breakdown"]["voltage_penalty"] - episode_return += reward - state = next_state - if done: - break - return ( - episode_return, - violation_time, - violation_value, - reward_for_power, - reward_for_good_action, - reward_for_penalty, - state_list, - ) +"""Legacy compatibility facade for algorithm utilities.""" + +from importlib import import_module + +__all__ = [ + "Config", + "ReplayBuffer", + "SumTree", + "build_mlp", + "get_episode_return", + "get_optim_param", +] + +_LAZY_EXPORTS = { + "Config": ("rl_adn.algorithms.training_config", "Config"), + "ReplayBuffer": ("rl_adn.algorithms.replay", "ReplayBuffer"), + "SumTree": ("rl_adn.algorithms.replay", "SumTree"), + "build_mlp": ("rl_adn.algorithms.torch_utils", "build_mlp"), + "get_episode_return": ("rl_adn.algorithms.evaluation", "get_episode_return"), + "get_optim_param": ("rl_adn.algorithms.torch_utils", "get_optim_param"), +} + + +def __getattr__(name: str): + if name not in _LAZY_EXPORTS: + raise AttributeError(f"module 'rl_adn.algorithms.utility' has no attribute {name!r}") + + module_name, attr_name = _LAZY_EXPORTS[name] + module = import_module(module_name) + return getattr(module, attr_name) diff --git a/rl_adn/benchmarks/__init__.py b/rl_adn/benchmarks/__init__.py index d9350b3..be0546f 100644 --- a/rl_adn/benchmarks/__init__.py +++ b/rl_adn/benchmarks/__init__.py @@ -1,5 +1,17 @@ """Benchmark exports for RL-ADN.""" -from rl_adn.benchmarks.pyomo_timeseries_pandapower import construct_opf_model +from rl_adn.benchmarks.pyomo_timeseries_pandapower import ( + BatterySpec, + DispatchBenchmarkData, + construct_opf_model, + convert_dict_to_pd, + convert_indexed_values_to_frame, +) -__all__ = ["construct_opf_model"] +__all__ = [ + "BatterySpec", + "DispatchBenchmarkData", + "construct_opf_model", + "convert_dict_to_pd", + "convert_indexed_values_to_frame", +] diff --git a/rl_adn/benchmarks/pyomo_timeseries_pandapower.py b/rl_adn/benchmarks/pyomo_timeseries_pandapower.py index 56f3784..035d7e7 100644 --- a/rl_adn/benchmarks/pyomo_timeseries_pandapower.py +++ b/rl_adn/benchmarks/pyomo_timeseries_pandapower.py @@ -1,81 +1,122 @@ -"""in this script, we consider constraints corresponding to the current multi battery environment, -that is to say the voltage regulation only considers that constrainting the limitation connected to the nodes connected to batteries, -Using 2020 11-6 as the example""" +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Any, Mapping import numpy as np import pandas as pd -from pyomo.environ import * - -from rl_adn import PowerNetEnv, make_env_config - - -def construct_opf_model(Vnom, Vmin, Vmax, Data_Network): - # Data Processing - battery_parameters = { - "capacity": 1.0, # MW.h - "max_charge": 0.3, # MW - "max_discharge": 0.3, # MW - "efficiency": 1, - "degradation": 0, # euro/kw - "max_soc": 0.8, - "min_soc": 0.2, - "initial_soc": 0.4, - } - TIMES = Data_Network["TIMES"] - NODES = Data_Network["NODES"] - LINES = Data_Network["LINES"] - Tb = Data_Network["Tb"] - PD = Data_Network["PD"] - QD = Data_Network["QD"] - R = Data_Network["R"] - X = Data_Network["X"] - BATTERY_NODES = Data_Network["BATTERY_NODES"] - # Type of Model - model = ConcreteModel() - # Define Sets - model.NODES = Set(initialize=NODES) - model.LINES = Set(initialize=LINES) - model.TIMES = Set(initialize=TIMES) - - # Define Parameters - model.Vnom = Param(initialize=Vnom, mutable=False) - model.Vmin = Param(initialize=Vmin, mutable=False) - model.Vmax = Param(initialize=Vmax, mutable=False) - model.Tb = Param(model.NODES, initialize=Tb, mutable=True) - # model.PD = Param(model.TIMES,model.NODES, initialize=0, mutable=True) # Node demand - model.QD = Param(model.TIMES, model.NODES, initialize=0, mutable=True) # Node demand - model.R = Param(model.LINES, initialize=R, mutable=False) # Line resistance - model.X = Param(model.LINES, initialize=X, mutable=False) # Line resistance - ## define parameters for battery - model.battery_initial_soc = Param(default=battery_parameters["initial_soc"]) - model.battery_capacity = Param(default=battery_parameters["capacity"]) - model.battery_soc_max = Param(default=battery_parameters["max_soc"]) - model.battery_soc_min = Param(default=battery_parameters["min_soc"]) - model.battery_max_change = Param(default=battery_parameters["max_charge"]) - - # define initialize PD - def PD_init_rule(model, time, node): - model.PD[time, node] = PD[time, node] - return model.PD[time, node] - - model.PD = Param(model.TIMES, model.NODES, initialize=PD_init_rule) - - def R_init_rule(model, i, j): - return model.R[i, j] - - model.RM = Param(model.LINES, initialize=R_init_rule) # Line resistance - - def X_init_rule(model, i, j): - return model.X[i, j] - - model.XM = Param(model.LINES, initialize=X_init_rule) # Line resistance - - # Define Variables - model.P = Var(model.TIMES, model.LINES, initialize=0) # Acive power flowing in lines - model.Q = Var(model.TIMES, model.LINES, initialize=0) # Reacive power flowing in lines - model.I = Var(model.TIMES, model.LINES, initialize=0) # Current of lines +def _require_pyomo(): + try: + from pyomo.environ import ConcreteModel, Constraint, Objective, Param, Set, Var, minimize + except ImportError as exc: + raise ImportError("Pyomo benchmark support requires the optional dependency 'pyomo'.") from exc + + return ConcreteModel, Constraint, Objective, Param, Set, Var, minimize + + +@dataclass(frozen=True) +class BatterySpec: + capacity_mwh: float = 1.0 + max_charge_mw: float = 0.3 + max_discharge_mw: float = 0.3 + efficiency: float = 1.0 + degradation_eur_per_kw: float = 0.0 + max_soc: float = 0.8 + min_soc: float = 0.2 + initial_soc: float = 0.4 + time_interval_minutes: float = 15.0 + + +@dataclass(frozen=True) +class DispatchBenchmarkData: + times: tuple[int, ...] + nodes: tuple[int, ...] + lines: tuple[tuple[int, int], ...] + tb: dict[int, int] + pd: np.ndarray + qd: np.ndarray + r: dict[tuple[int, int], float] + x: dict[tuple[int, int], float] + battery_nodes: frozenset[int] + price: np.ndarray + battery: BatterySpec = field(default_factory=BatterySpec) + + @classmethod + def from_mapping( + cls, + data_network: Mapping[str, Any], + *, + battery: BatterySpec | None = None, + ) -> "DispatchBenchmarkData": + times = tuple(data_network["TIMES"]) + nodes = tuple(data_network["NODES"]) + lines = tuple(data_network["LINES"]) + tb = dict(data_network["Tb"]) + pd_array = np.asarray(data_network["PD"], dtype=float) + qd_input = data_network.get("QD") + qd_array = np.zeros_like(pd_array) if qd_input is None else np.asarray(qd_input, dtype=float) + price = np.asarray(data_network["PRICE"], dtype=float) + if pd_array.shape[0] != len(times): + raise ValueError("PD must have one row per time step") + if qd_array.shape != pd_array.shape: + raise ValueError("QD must match PD shape") + if price.shape[0] != len(times): + raise ValueError("PRICE must have one value per time step") + return cls( + times=times, + nodes=nodes, + lines=lines, + tb=tb, + pd=pd_array, + qd=qd_array, + r=dict(data_network["R"]), + x=dict(data_network["X"]), + battery_nodes=frozenset(data_network["BATTERY_NODES"]), + price=price, + battery=battery or BatterySpec(), + ) + + +def construct_opf_model(v_nom: float, v_min: float, v_max: float, data_network: Mapping[str, Any] | DispatchBenchmarkData): + """Build the Pyomo dispatch model for ESS scheduling on a radial feeder.""" + ConcreteModel, Constraint, Objective, Param, Set, Var, minimize = _require_pyomo() + data = data_network if isinstance(data_network, DispatchBenchmarkData) else DispatchBenchmarkData.from_mapping(data_network) + + model = ConcreteModel() + model.NODES = Set(initialize=data.nodes) + model.LINES = Set(initialize=data.lines) + model.TIMES = Set(initialize=data.times) + + model.Vnom = Param(initialize=v_nom, mutable=False) + model.Vmin = Param(initialize=v_min, mutable=False) + model.Vmax = Param(initialize=v_max, mutable=False) + model.Tb = Param(model.NODES, initialize=data.tb, mutable=True) + model.QD = Param( + model.TIMES, + model.NODES, + initialize=lambda _, time, node: float(data.qd[time, node]), + mutable=False, + ) + model.R = Param(model.LINES, initialize=data.r, mutable=False) + model.X = Param(model.LINES, initialize=data.x, mutable=False) + model.battery_initial_soc = Param(default=data.battery.initial_soc) + model.battery_capacity = Param(default=data.battery.capacity_mwh) + model.battery_soc_max = Param(default=data.battery.max_soc) + model.battery_soc_min = Param(default=data.battery.min_soc) + model.battery_max_change = Param(default=data.battery.max_charge_mw) + model.PD = Param( + model.TIMES, + model.NODES, + initialize=lambda _, time, node: float(data.pd[time, node]), + ) + model.RM = Param(model.LINES, initialize=lambda _, i, j: model.R[i, j]) + model.XM = Param(model.LINES, initialize=lambda _, i, j: model.X[i, j]) + + model.P = Var(model.TIMES, model.LINES, initialize=0) + model.Q = Var(model.TIMES, model.LINES, initialize=0) + model.I = Var(model.TIMES, model.LINES, initialize=0) model.SOC = Var( model.TIMES, model.NODES, @@ -83,15 +124,9 @@ def X_init_rule(model, i, j): bounds=(model.battery_soc_min, model.battery_soc_max), ) - # we set energy>0 is discharge, also only when no slack bus we put battery - def energy_change_rule(model, time, i): - if i not in BATTERY_NODES: - tem = 0.0 - model.energy_change[time, i].fixed = True - else: - tem = 0.0 - model.energy_change[time, i].fixed = False - return tem + def energy_change_rule(model, time, node): + model.energy_change[time, node].fixed = node not in data.battery_nodes + return 0.0 model.energy_change = Var( model.TIMES, @@ -100,123 +135,74 @@ def energy_change_rule(model, time, i): bounds=(-model.battery_max_change, model.battery_max_change), ) - def PS_init_rule(model, time, i): - # for time in model.TIMES: - if model.Tb[i].value == 0: - temp = 0.0 - model.PS[time, i].fixed = True - else: - temp = 0.0 - return temp - - model.PS = Var(model.TIMES, model.NODES, initialize=PS_init_rule) # Active power of the SS - - def QS_init_rule(model, time, i): - # for time in model.TIMES: - if model.Tb[i].value == 0: - temp = 0.0 - model.QS[time, i].fixed = True - else: - temp = 0.0 - return temp - - model.QS = Var(model.TIMES, model.NODES, initialize=QS_init_rule) # Reactive power of the SS - - # price init rule - def PRICE_init_rule(model, time): - return PRICE[time] - - model.PRICE = Param(model.TIMES, initialize=PRICE_init_rule, mutable=False) - - # Voltage of nodes - def Voltage_init(model, time, i): - # for time in model.TIMES: - if model.Tb[i].value == 1: - temp = model.Vnom - model.V[time, i].fixed = True - else: - temp = model.Vnom - model.V[time, i].fixed = False - return temp - - model.V = Var(model.TIMES, model.NODES, initialize=Voltage_init) - - # Define Objective Function,minimize the optimal power loss. Actually, when we only have one source from the grid, - """Since each element of model.LINES is a tuple of two integers, - you will need to use two indices to access the variable indexed by model.LINES. - For example, if you define a variable P indexed by both model.LINES and model.TIMES, - you would access the value of P for the line (1,2) at time t=1 using model.P[1, (1,2)].""" - # def act_loss(model): - # return (sum(sum(model.RM[i, j] * (model.I[time,(i, j)] ** 2) for i, j in model.LINES)for time in model.TIMES)) - - # here we create another objective: minimizing the imported power from external grid - # def min_power_ext_grid(model): - # - # return (sum(sum(model.PS[time,node]for node in model.NODES)for time in model.TIMES)) - - # Update the objective function to minimize the cost of buying energy from the external grid - def min_cost_ext_grid(model): - return sum(sum(model.PS[time, node] * model.PRICE[time] for node in model.NODES) for time in model.TIMES) - - model.obj = Objective(rule=min_cost_ext_grid, sense=minimize) - - # Update the objective function to earn money from battery dispatch - # def max_benefits_dispatch_battery(model): - # return sum(sum(model.energy_change[time, node] * model.PRICE[time] for node in model.NODES) for time in model.TIMES) - # - # model.obj = Objective(rule=max_benefits_dispatch_battery,sense=maximize) - - # model.obj = Objective(rule=min_power_ext_grid) - # model.obj = Objective(rule=act_loss) - # we need to revise this part for adding time constraint into here. - # %% Define Constraints - # define soc update constraint + def substation_active_rule(model, time, node): + if model.Tb[node].value == 0: + model.PS[time, node].fixed = True + return 0.0 + + def substation_reactive_rule(model, time, node): + if model.Tb[node].value == 0: + model.QS[time, node].fixed = True + return 0.0 + + model.PS = Var(model.TIMES, model.NODES, initialize=substation_active_rule) + model.QS = Var(model.TIMES, model.NODES, initialize=substation_reactive_rule) + model.PRICE = Param(model.TIMES, initialize=lambda _, time: float(data.price[time]), mutable=False) + + def voltage_init(model, time, node): + if model.Tb[node].value == 1: + model.V[time, node].fixed = True + return model.Vnom + + model.V = Var(model.TIMES, model.NODES, initialize=voltage_init) + model.obj = Objective( + rule=lambda model: sum( + sum(model.PS[time, node] * model.PRICE[time] for node in model.NODES) + for time in model.TIMES + ), + sense=minimize, + ) + + interval_hours = data.battery.time_interval_minutes / 60.0 def soc_update_rule(model, time, node): - if node not in BATTERY_NODES: + if node not in data.battery_nodes: return Constraint.Skip if time == model.TIMES.first(): - return ( - model.SOC[time, node] - == model.battery_initial_soc - (model.energy_change[time, node] * 15.0 / 60.0) / model.battery_capacity - ) - else: - return ( - model.SOC[time, node] - == model.SOC[model.TIMES.prev(time), node] - - (model.energy_change[time, node] * 15.0 / 60.0) / model.battery_capacity - ) + return model.SOC[time, node] == model.battery_initial_soc - ( + model.energy_change[time, node] * interval_hours + ) / model.battery_capacity + return model.SOC[time, node] == model.SOC[model.TIMES.prev(time), node] - ( + model.energy_change[time, node] * interval_hours + ) / model.battery_capacity model.constaint_soc_update = Constraint(model.TIMES, model.NODES, rule=soc_update_rule) - # for line k consumption == injection - def active_power_flow_rule(model, time, k): - + def active_power_flow_rule(model, time, node): return ( - sum(model.P[time, (j, i)] for j, i in model.LINES if i == k) + sum(model.P[time, (j, i)] for j, i in model.LINES if i == node) - sum( - model.P[time, (i, j)] + model.RM[i, j] * (model.I[time, (i, j)] ** 2) for i, j in model.LINES if k == i + model.P[time, (i, j)] + model.RM[i, j] * (model.I[time, (i, j)] ** 2) + for i, j in model.LINES + if node == i ) - + model.PS[time, k] - + model.energy_change[time, k] - == model.PD[time, k] + + model.PS[time, node] + + model.energy_change[time, node] + == model.PD[time, node] ) - model.active_power_flow = Constraint(model.TIMES, model.NODES, rule=active_power_flow_rule) - - def reactive_power_flow_rule(model, time, k): + def reactive_power_flow_rule(model, time, node): return ( - sum(model.Q[time, (j, i)] for j, i in model.LINES if i == k) + sum(model.Q[time, (j, i)] for j, i in model.LINES if i == node) - sum( - model.Q[time, (i, j)] + model.XM[i, j] * (model.I[time, (i, j)] ** 2) for i, j in model.LINES if k == i + model.Q[time, (i, j)] + model.XM[i, j] * (model.I[time, (i, j)] ** 2) + for i, j in model.LINES + if node == i ) - + model.QS[time, k] - == model.QD[time, k] + + model.QS[time, node] + == model.QD[time, node] ) - model.reactive_power_flow = Constraint(model.TIMES, model.NODES, rule=reactive_power_flow_rule) - - ## role of voltage drop def voltage_drop_rule(model, time, i, j): return ( model.V[time, i] ** 2 @@ -225,38 +211,34 @@ def voltage_drop_rule(model, time, i, j): - model.V[time, j] ** 2 ) == 0 - model.voltage_drop = Constraint(model.TIMES, model.LINES, rule=voltage_drop_rule) - - def define_current_rule(model, time, i, j): - return (model.I[time, (i, j)] ** 2) * (model.V[time, j] ** 2) == model.P[time, (i, j)] ** 2 + model.Q[ - time, (i, j) - ] ** 2 - - model.define_current = Constraint(model.TIMES, model.LINES, rule=define_current_rule) - - # here the current limit is over 0, representing that current can only from i to j, instead of versa. - # we change this step and try to calculate it according to the result three phase one - - def current_limit_rule(model, time, i, j): - return (0, model.I[time, (i, j)], None) - - # if we cancel this, then things to error - model.current_limit = Constraint(model.TIMES, model.LINES, rule=current_limit_rule) + def current_definition_rule(model, time, i, j): + return (model.I[time, (i, j)] ** 2) * (model.V[time, j] ** 2) == ( + model.P[time, (i, j)] ** 2 + model.Q[time, (i, j)] ** 2 + ) - def voltage_limit_rule(model, time, i): - if i in BATTERY_NODES: - return (model.Vmin, model.V[time, i], model.Vmax) - return Constraint.Skip + model.active_power_flow = Constraint(model.TIMES, model.NODES, rule=active_power_flow_rule) + model.reactive_power_flow = Constraint(model.TIMES, model.NODES, rule=reactive_power_flow_rule) + model.voltage_drop = Constraint(model.TIMES, model.LINES, rule=voltage_drop_rule) + model.define_current = Constraint(model.TIMES, model.LINES, rule=current_definition_rule) + model.current_limit = Constraint(model.TIMES, model.LINES, rule=lambda model, time, i, j: (0, model.I[time, (i, j)], None)) + model.voltage_limit = Constraint( + model.TIMES, + model.NODES, + rule=lambda model, time, node: (model.Vmin, model.V[time, node], model.Vmax) + if node in data.battery_nodes + else Constraint.Skip, + ) + return model - model.voltage_limit = Constraint(model.TIMES, model.NODES, rule=voltage_limit_rule) - return model +def convert_indexed_values_to_frame(data: Mapping[tuple[int, int], Any]) -> pd.DataFrame: + columns = sorted({key[1] for key in data}) + frame = pd.DataFrame(columns=columns) + for (row_index, column_index), value in data.items(): + frame.loc[row_index, column_index] = value + return frame.sort_index().sort_index(axis=1) -def convert_dict_to_pd(data: dict): - df = pd.DataFrame(columns=list(set([k[1] for k in data.keys()]))) - for key, value in data.items(): - df.loc[key[0], key[1]] = value - # df=df.iloc[:,1:] - # df=df.drop(df.columns[0],axis=1) - return df +def convert_dict_to_pd(data: Mapping[tuple[int, int], Any]) -> pd.DataFrame: + """Legacy alias for ``convert_indexed_values_to_frame``.""" + return convert_indexed_values_to_frame(data) diff --git a/rl_adn/data_augment/__init__.py b/rl_adn/data_augment/__init__.py index e69de29..b010cda 100644 --- a/rl_adn/data_augment/__init__.py +++ b/rl_adn/data_augment/__init__.py @@ -0,0 +1,5 @@ +"""Data augmentation helpers for RL-ADN.""" + +from rl_adn.data_augment.data_augment import ActivePowerDataManager, TimeSeriesDataAugmentor + +__all__ = ["ActivePowerDataManager", "TimeSeriesDataAugmentor"] diff --git a/rl_adn/data_augment/data_augment.py b/rl_adn/data_augment/data_augment.py index 2329609..23a45a9 100644 --- a/rl_adn/data_augment/data_augment.py +++ b/rl_adn/data_augment/data_augment.py @@ -1,263 +1,190 @@ +from __future__ import annotations + import re -from datetime import timedelta +from datetime import datetime, timedelta +from pathlib import Path +from typing import Literal import numpy as np import pandas as pd -from copulas.multivariate import GaussianMultivariate -from multicopula import EllipticalCopula from scipy.optimize import brentq from scipy.stats import norm -from sklearn.mixture import GaussianMixture from rl_adn.data import GeneralPowerDataManager +AugmentationMethod = Literal["GMC", "GMM", "TC"] -class ActivePowerDataManager(GeneralPowerDataManager): - """ - A subclass of GeneralPowerDataManager that adds a method for retrieving - active power data specifically. - """ - - def __init__(self, datapath: str) -> None: - """ - Initialize the ActivePowerDataManager object with the path to the data file. - """ - if not datapath: - raise ValueError("Please input the correct datapath") - self.df = pd.read_csv(datapath, index_col="date_time") - self.df.index = pd.to_datetime(self.df.index) +def _require_gmm(): + try: + from sklearn.mixture import GaussianMixture + except ImportError as exc: + raise ImportError("GMM-based data augmentation requires the optional dependency 'scikit-learn'.") from exc + return GaussianMixture - # Assuming the data interval is consistent throughout the dataset - self.time_interval = int((self.df.index[1] - self.df.index[0]).seconds / 60) - print(f"Data time interval: {self.time_interval} minutes") - def get_active_power_data(self) -> np.ndarray: - """ - Retrieve and preprocess active power data from the dataset. - """ - self.df.interpolate(method="linear", inplace=True) - self.df["day"] = self.df.index.date - self.df["time"] = self.df.index.time +def _require_gmc(): + try: + from copulas.multivariate import GaussianMultivariate + except ImportError as exc: + raise ImportError("GMC-based data augmentation requires the optional dependency 'copulas'.") from exc + return GaussianMultivariate - count_per_day = self.df.groupby("day").size() - expected_time_steps = 24 * 60 / self.time_interval - days_with_extra_steps = count_per_day[count_per_day > expected_time_steps] - if not days_with_extra_steps.empty: - self.df = self.df[~self.df["day"].isin(days_with_extra_steps.index)] +def _require_tc(): + try: + from multicopula import EllipticalCopula + except ImportError as exc: + raise ImportError("TC-based data augmentation requires the optional dependency 'multicopula'.") from exc + return EllipticalCopula - active_power_columns = [col for col in self.df.columns if re.fullmatch(r"active_power(_\w+)?", col)] - active_power_df = self.df[active_power_columns].copy() - active_power_df["day"] = self.df["day"] - active_power_df["time"] = self.df["time"] - reshaped_active_power_df = active_power_df.set_index(["day", "time"]).stack().reset_index().rename(columns={"level_2": "node", 0: "value"}) - - grouped_active_power_df = reshaped_active_power_df.groupby("time")["value"].apply(list).reset_index() - - reshaped_df = pd.DataFrame(grouped_active_power_df["value"].tolist(), index=grouped_active_power_df["time"]) +class ActivePowerDataManager(GeneralPowerDataManager): + """Specialized data manager for active-power augmentation workflows.""" - active_power_array = reshaped_df.to_numpy().T + def get_active_power_data(self) -> np.ndarray: + if not self.active_power_cols: + raise ValueError("No active power columns were found in the dataset") + + active_power_df = self.df[self.active_power_cols].copy() + active_power_df["day"] = self.df.index.normalize() + active_power_df["time"] = self.df.index.time + + expected_time_steps = int(24 * 60 / self.time_interval) + complete_days = active_power_df.groupby("day").size() + complete_days = complete_days[complete_days == expected_time_steps].index + active_power_df = active_power_df[active_power_df["day"].isin(complete_days)] + if active_power_df.empty: + raise ValueError("No complete days were found for active-power augmentation") + + reshaped = active_power_df.set_index(["day", "time"]).stack().unstack("time") + active_power_array = reshaped.to_numpy(dtype=float) active_power_array = active_power_array[~np.isnan(active_power_array).any(axis=1)] - + if active_power_array.size == 0: + raise ValueError("Active power data contains only NaN values after reshaping") return active_power_array class TimeSeriesDataAugmentor: - def __init__(self, data, augmentation_model_name="GMC"): - """ - Initialize the data augmentor with a data manager instance and the selected augmentation model. - Additional parameters can be set here if required. - """ + """Generate synthetic active-power time series using copula or GMM methods.""" + + def __init__(self, data: ActivePowerDataManager, augmentation_model_name: AugmentationMethod = "GMC"): self.data = data self.augmentation_model_name = augmentation_model_name self.augmentation_model = None + self.n_models = int(24.0 * 60.0 / self.data.time_interval) + self.gmm_models = [] + self._create_augmentation_model() - self._create_augmentation_model() # Create the augmentation model upon initialization + def _create_augmentation_model(self) -> None: + active_power_array = self.data.get_active_power_data() + GaussianMixture = _require_gmm() + best_components = [self._bic_value(active_power_array[:, i].reshape(-1, 1), 20, GaussianMixture) for i in range(self.n_models)] + self.gmm_models = [GaussianMixture(n_components=count).fit(active_power_array[:, i].reshape(-1, 1)) for i, count in enumerate(best_components)] - def _create_augmentation_model(self): - """ - Private method to create the augmentation model based on the chosen method, e.g., GMC. - GMC: Data augmentation using Gaussian Mixture Copulas - GMM: Data augmentaiton using Gaussian Mixture models - TC: Data augmentaiton using T Copulas - """ if self.augmentation_model_name == "GMC": - # Extract data from the data manager - active_power_array = self.data.get_active_power_data() - - # Determine the best number of components for each GMM model, for one day now it is 96 time steps - self.n_models = int(24.0 * 60.0 / self.data.time_interval) - - # print(f'data manager time interval {self.data.time_interval}') - best_components = [self._bic_value(active_power_array[:, i].reshape(-1, 1), 20) for i in range(self.n_models)] - - # Fit the GMM models - self.gmm_models = [GaussianMixture(n_components=bc).fit(active_power_array[:, i].reshape(-1, 1)) for i, bc in enumerate(best_components)] - - # Transform the data to standard format for copula fitting - std_input_data = np.empty((active_power_array.shape[0], self.n_models)) - for i in range(self.n_models): - std_input_data[:, i] = np.array([self._gmm_cdf(self.gmm_models[i], x) for x in active_power_array[:, i]]).reshape(1, -1) - self.copula = GaussianMultivariate() - self.copula.fit(std_input_data) - - # Assign the copula as the augmentation model - self.augmentation_model = self.copula - - if self.augmentation_model_name == "GMM": - active_power_array = self.data.get_active_power_data() - - # Determine the best number of components for each GMM model, for one day now it is 96 time steps - self.n_models = int(24.0 * 60.0 / self.data.time_interval) - - # print(f'data manager time interval {self.data.time_interval}') - best_components = [self._bic_value(active_power_array[:, i].reshape(-1, 1), 20) for i in range(self.n_models)] - - # Fit the GMM models - self.gmm_models = [GaussianMixture(n_components=bc).fit(active_power_array[:, i].reshape(-1, 1)) for i, bc in enumerate(best_components)] - + GaussianMultivariate = _require_gmc() + standard_input_data = np.empty((active_power_array.shape[0], self.n_models)) + for index in range(self.n_models): + standard_input_data[:, index] = np.array([self._gmm_cdf(self.gmm_models[index], value) for value in active_power_array[:, index]]) + copula = GaussianMultivariate() + copula.fit(standard_input_data) + self.augmentation_model = copula + elif self.augmentation_model_name == "GMM": self.augmentation_model = self.gmm_models - - if self.augmentation_model_name == "TC": - active_power_array = self.data.get_active_power_data() - self.n_models = int(24.0 * 60.0 / self.data.time_interval) - - self.tc_model = EllipticalCopula(active_power_array.T) - self.tc_model.fit() - self.augmentation_model = self.tc_model - - def _gmm_cdf(self, gmm, x): - """ - Convert CDF to pseudo-observations in the interval [0, 1]. - """ - cdf = 0 - for n in range(gmm.n_components): - cdf += gmm.weights_[n] * norm.cdf(x, gmm.means_[n, 0], np.sqrt(gmm.covariances_[n, 0])) - return cdf - - def _inverse_gmm_cdf(self, gmm, percentile): - """ - Find the inverse of the CDF for a given percentile using a GMM model. - """ - - def f(x): + elif self.augmentation_model_name == "TC": + EllipticalCopula = _require_tc() + tc_model = EllipticalCopula(active_power_array.T) + tc_model.fit() + self.augmentation_model = tc_model + else: + raise ValueError(f"Unsupported augmentation_model_name: {self.augmentation_model_name}") + + def _gmm_cdf(self, gmm, value: float) -> float: + cdf = 0.0 + for component in range(gmm.n_components): + cdf += gmm.weights_[component] * norm.cdf( + value, + gmm.means_[component, 0], + np.sqrt(gmm.covariances_[component, 0]), + ) + return float(cdf) + + def _inverse_gmm_cdf(self, gmm, percentile: float) -> float: + def root_fn(x: float) -> float: return self._gmm_cdf(gmm, x) - percentile - return brentq(f, -3000, 3000) + return float(brentq(root_fn, -3000, 3000)) - def _bic_value(self, data, n_components_range): - """ - Compute BIC value for different numbers of components to determine the best model. - """ + @staticmethod + def _bic_value(data: np.ndarray, n_components_range: int, gaussian_mixture_cls) -> int: bic_values = [] for n_components in range(1, n_components_range): - gmm = GaussianMixture(n_components=n_components).fit(data) + gmm = gaussian_mixture_cls(n_components=n_components).fit(data) bic_values.append(gmm.bic(data)) - best_component = np.argmin(bic_values) + 1 - return best_component - - def check_data_format(self): - """ - Verify that the data matches the expected format for augmentation. - """ - # Implementation details to verify data format - - def augment_data(self, num_nodes, num_days, start_date): - """ - Perform data augmentation using the specified model and parameters. - """ + return int(np.argmin(bic_values) + 1) + + def augment_data(self, num_nodes: int, num_days: int, start_date: datetime) -> pd.DataFrame: + num_samples = num_days * num_nodes if self.augmentation_model_name == "GMC": - num_samples = num_days * num_nodes - print("The number of samples is", num_samples) - - generated_pesudo_obs = np.empty((0, self.n_models)) - count = 0 - while True: - # print('days is generated in sample', count) - gen_one_sample = np.array(self.copula.sample(1)) - if gen_one_sample.min() > 0 and gen_one_sample.max() < 1: - count += 1 - generated_pesudo_obs = np.vstack((gen_one_sample, generated_pesudo_obs)) - if count == num_samples: - break - # print(' the pesudo data is now sampled and next process is to transfer it to the realistic data') - tran_samples = np.empty((generated_pesudo_obs.shape[0], generated_pesudo_obs.shape[1])) - for i in range(self.n_models): - tran_samples[:, i] = np.array([self._inverse_gmm_cdf(self.gmm_models[i], u) for u in generated_pesudo_obs[:, i]]) - print(f"the {i} model columns now is calculated") - tran_samples = tran_samples.flatten() - - if self.augmentation_model_name == "GMM": - num_samples = num_days * num_nodes - print("The number of samples is", num_samples) - - # generating the data + generated_pseudo_obs = self._sample_gmc(num_samples) + transformed_samples = np.empty_like(generated_pseudo_obs) + for index in range(self.n_models): + transformed_samples[:, index] = np.array([self._inverse_gmm_cdf(self.gmm_models[index], u) for u in generated_pseudo_obs[:, index]]) + flattened_samples = transformed_samples.flatten() + elif self.augmentation_model_name == "GMM": gmm_samples = np.empty((num_samples, self.n_models)) - for i in range(self.n_models): - gmm_samples[:, i] = self.gmm_models[i].sample(num_samples)[0].reshape(-1) - - tran_samples = gmm_samples.flatten() - - if self.augmentation_model_name == "TC": - num_samples = num_days * num_nodes - print("The number of samples is", num_samples) - - # generating the data - TC_samples = np.empty((0, self.n_models)) - count = 0 - while True: - gen_one_sample = np.array(self.tc_model.sample(1)).reshape(1, -1) - - # cancel inf - if not np.isinf(gen_one_sample).any(): - count += 1 - # print(count,num_samples) - TC_samples = np.vstack((gen_one_sample, TC_samples)) - if count == num_samples: - break - tran_samples = TC_samples.flatten() - - # Initialize lists to hold timestamps and node indices + for index in range(self.n_models): + gmm_samples[:, index] = self.gmm_models[index].sample(num_samples)[0].reshape(-1) + flattened_samples = gmm_samples.flatten() + elif self.augmentation_model_name == "TC": + tc_samples = self._sample_tc(num_samples) + flattened_samples = tc_samples.flatten() + else: + raise ValueError(f"Unsupported augmentation_model_name: {self.augmentation_model_name}") + timestamps = [] node_index = [] + time_step = timedelta(minutes=self.data.time_interval) + for day_offset in range(num_days): + for node_id in range(1, num_nodes + 1): + timestamps.extend([start_date + timedelta(days=day_offset) + slot * time_step for slot in range(self.n_models)]) + node_index.extend([f"active_power_node_{node_id}" for _ in range(self.n_models)]) - # Populate timestamps and node_index for each day and each node - - for day in range(num_days): - for node in range(1, num_nodes + 1): - time_step = timedelta(minutes=self.data.time_interval) - timestamps.extend([start_date + timedelta(days=day) + i * time_step for i in range(self.n_models)]) - node_index.extend([f"active_power_node_{node}" for _ in range(self.n_models)]) - - # Create DataFrame - synthetic_data_df = pd.DataFrame({"date_time": timestamps, "node": node_index, "value": tran_samples}) - - # Pivot the DataFrame to get it into the desired format + synthetic_data_df = pd.DataFrame({"date_time": timestamps, "node": node_index, "value": flattened_samples}) augmented_df = synthetic_data_df.pivot(index="date_time", columns="node", values="value").reset_index() - # Reorder the columns based on we need - active_power_cols = self.sort_columns(augmented_df.columns, r"active_power(_\w+)?") - reactive_power_cols = self.sort_columns(augmented_df.columns, r"reactive_power(_\w+)?") - renewable_active_power_cols = self.sort_columns(augmented_df.columns, r"renewable_active_power(_\w+)?") - renewable_reactive_power_cols = self.sort_columns(augmented_df.columns, r"renewable_reactive_power(_\w+)?") - price_cols = self.sort_columns(augmented_df.columns, r"price(_\w+)?") - # Combine columns in the specified order - ordered_columns = ["date_time"] + active_power_cols + reactive_power_cols + renewable_active_power_cols + renewable_reactive_power_cols + price_cols - ordered_augmented_df = augmented_df[ordered_columns] - - return ordered_augmented_df - - def save_augmented_data(self, augmented_df, file_name): - augmented_df.to_csv(file_name, index=False) - print("The data file is stored:", file_name) - - def sort_columns(self, columns, pattern): - def sort_key(col_name): - parts = col_name.split("_") + ordered_columns = ["date_time"] + self.sort_columns(augmented_df.columns, r"active_power(_\w+)?") + return augmented_df[ordered_columns] + + def _sample_gmc(self, sample_count: int) -> np.ndarray: + generated_pseudo_obs = np.empty((0, self.n_models)) + while generated_pseudo_obs.shape[0] < sample_count: + sampled = np.array(self.augmentation_model.sample(1)) + if sampled.min() > 0 and sampled.max() < 1: + generated_pseudo_obs = np.vstack((sampled, generated_pseudo_obs)) + return generated_pseudo_obs[:sample_count] + + def _sample_tc(self, sample_count: int) -> np.ndarray: + tc_samples = np.empty((0, self.n_models)) + while tc_samples.shape[0] < sample_count: + sampled = np.array(self.augmentation_model.sample(1)).reshape(1, -1) + if not np.isinf(sampled).any(): + tc_samples = np.vstack((sampled, tc_samples)) + return tc_samples[:sample_count] + + @staticmethod + def save_augmented_data(augmented_df: pd.DataFrame, file_name: str | Path) -> Path: + output_path = Path(file_name) + augmented_df.to_csv(output_path, index=False) + return output_path + + @staticmethod + def sort_columns(columns, pattern: str): + def sort_key(column_name: str) -> int: + parts = column_name.split("_") if parts[-1].isdigit(): return int(parts[-1]) - return 0 # Default sort value for non-numeric endings + return 0 - filtered_cols = [col for col in columns if re.fullmatch(pattern, col)] + filtered_cols = [column for column in columns if re.fullmatch(pattern, column)] return sorted(filtered_cols, key=sort_key) diff --git a/tests/test_legacy_cleanup_smoke.py b/tests/test_legacy_cleanup_smoke.py new file mode 100644 index 0000000..612f607 --- /dev/null +++ b/tests/test_legacy_cleanup_smoke.py @@ -0,0 +1,79 @@ +from __future__ import annotations + +import importlib +import importlib.util +from pathlib import Path + +import pandas as pd +import pytest + + +def test_algorithm_utility_facade_exports_expected_symbols(): + from rl_adn.algorithms import utility + + assert hasattr(utility, "Config") + assert set(utility.__all__) == {"Config", "ReplayBuffer", "SumTree", "build_mlp", "get_episode_return", "get_optim_param"} + + if importlib.util.find_spec("torch") is not None: + for name in ("ReplayBuffer", "SumTree", "build_mlp", "get_episode_return", "get_optim_param"): + assert hasattr(utility, name) + + +def test_benchmark_module_imports_without_pyomo_side_effects(): + module = importlib.import_module("rl_adn.benchmarks.pyomo_timeseries_pandapower") + assert hasattr(module, "DispatchBenchmarkData") + assert hasattr(module, "construct_opf_model") + + +def test_benchmark_frame_conversion_orders_rows_and_columns(): + from rl_adn.benchmarks.pyomo_timeseries_pandapower import convert_indexed_values_to_frame + + frame = convert_indexed_values_to_frame({(1, 2): 5, (0, 1): 3, (0, 2): 4}) + assert list(frame.index) == [0, 1] + assert list(frame.columns) == [1, 2] + assert frame.loc[0, 1] == 3 + + +def test_active_power_data_manager_extracts_day_node_matrix(tmp_path: Path): + from rl_adn.data_augment.data_augment import ActivePowerDataManager + + periods = 96 + frame = pd.DataFrame( + { + "date_time": pd.date_range("2021-01-01", periods=periods, freq="15min", tz="UTC"), + "active_power_node_1": range(periods), + "active_power_node_2": range(periods, periods * 2), + } + ) + csv_path = tmp_path / "active_power.csv" + frame.to_csv(csv_path, index=False) + + manager = ActivePowerDataManager(str(csv_path)) + active_power_data = manager.get_active_power_data() + + assert active_power_data.shape == (2, periods) + + +def test_construct_opf_model_requires_pyomo_when_dependency_missing(): + from rl_adn.benchmarks.pyomo_timeseries_pandapower import BatterySpec, DispatchBenchmarkData, construct_opf_model + + data = DispatchBenchmarkData( + times=(0,), + nodes=(0, 1), + lines=((0, 1),), + tb={0: 1, 1: 0}, + pd=[[0.0, 0.0]], + qd=[[0.0, 0.0]], + r={(0, 1): 0.1}, + x={(0, 1): 0.1}, + battery_nodes=frozenset({1}), + price=[1.0], + battery=BatterySpec(), + ) + + if importlib.util.find_spec("pyomo") is None: + with pytest.raises(ImportError): + construct_opf_model(1.0, 0.95, 1.05, data) + else: + model = construct_opf_model(1.0, 0.95, 1.05, data) + assert model is not None