import pandas as pd
import networkx as nx
import numpy as np
import itertools
# --- Constants and Assumptions ---
# These should be clearly stated and can be modified.
VOLTAGE_KV = 10.0 # Line voltage in kV
ROOT_3 = np.sqrt(3)
BASE_DG_CAPACITY_KW = 300.0 # Initial capacity for each DG
N_DG = 8
# Failure rates from problem description
FAILURE_RATE_DG_PERCENT = 0.5 / 100.0
# FAILURE_RATE_USER_PERCENT = 0.5 / 100.0 # Not directly used in this simplified line-fault model for widespread outages
# FAILURE_RATE_SWITCH_PERCENT = 0.2 / 100.0 # Assuming switch failures manifest as line failures or inability to operate tie
FAILURE_RATE_LINE_PER_KM = 0.002 # Per km per year (assuming rates are annual)
# Costs (placeholders - these are critical for actual risk values)
# Value of Lost Load (VoLL) in monetary units per kW per hour.
# For risk = P * C, if P is annual probability, C should be impact of one event.
# Let's define C_loss as total kW unserved * a severity factor.
# Or, if we want an annual risk cost: P_annual_fault * kW_unserved * hours_outage * cost_per_kWh
# For simplicity, using $/kW of unserved load for the consequence C.
COST_VOLL_PER_KW = 10.0 # Example: $10 per kW of unserved load
AVG_OUTAGE_DURATION_H = 4 # Example: average hours for an outage, if converting to energy
# Cost of Overload (Consequence C_over)
# This can be complex: accelerated aging, tripping, damage.
# Simplified: A penalty if any line is overloaded in a given state.
COST_PENALTY_FOR_ANY_OVERLOAD = 1000.0 # Example: $1000 penalty if system is in an overloaded state
# Or, a cost per MWh of overloaded energy, or per overloaded line.
# Line and Feeder Capacities
# Main feeder rated current from problem: 220A.
# P_rated_feeder_kW = ROOT_3 * VOLTAGE_KV * FEEDER_RATED_CURRENT_A * 1.0 (pf=1)
# = 1.732 * 10 * 220 = 3810.4 kW (approx 3.8 MW, problem says 2.2MW for 220A, implies lower pf or different basis)
# Let's use current as the primary limit.
FEEDER_RATED_CURRENT_A = 220.0
# Assumption for individual line segments: For this model, we'll assume all lines
# have a rated current equal to the main feeder. This is a strong simplification.
# A more detailed model would assign ratings based on conductor types or downstream load.
LINE_RATED_CURRENT_A = 100.0 # More conservative assumption for individual segments than 220A. Needs proper engineering values.
# For lines directly from substation, perhaps 220A is more appropriate.
# Let's use a dictionary for specific line ratings if known, else default.
DEFAULT_LINE_RATED_CURRENT_A = 100.0
# Tie Line Capacity
TIE_LINE_RATED_CURRENT_A = 150.0 # Assumption, should be based on tie switch/line capacity
# DG Locations (Node IDs from 1 to 62) - Based on Figure 1 interpretation
DG_LOCATIONS_KW = {
6: BASE_DG_CAPACITY_KW, 10: BASE_DG_CAPACITY_KW, 15: BASE_DG_CAPACITY_KW,
27: BASE_DG_CAPACITY_KW, 31: BASE_DG_CAPACITY_KW, 37: BASE_DG_CAPACITY_KW,
50: BASE_DG_CAPACITY_KW, 58: BASE_DG_CAPACITY_KW
}
# Tie Switches: (node1, node2, switch_id_text) - normally open
# Interpretation based on careful review of Figure 1:
# S13-1: (13, 22) - Intra-Feeder 1 (Connects two branches of Feeder 1)
# S29-2: (29, 42) - Intra-Feeder 2 (Connects two branches of Feeder 2)
# S62-3: (62, 19) - Inter-Feeder (Connects Feeder 3 (node 62) to Feeder 1 (node 19))
TIE_SWITCHES_INFO = [
{'nodes': (13, 22), 'id': 'S13-1', 'type': 'intra-F1', 'capacity_A': TIE_LINE_RATED_CURRENT_A},
{'nodes': (29, 42), 'id': 'S29-2', 'type': 'intra-F2', 'capacity_A': TIE_LINE_RATED_CURRENT_A},
{'nodes': (62, 19), 'id': 'S62-3', 'type': 'inter-F3_F1', 'capacity_A': TIE_LINE_RATED_CURRENT_A}
]
# This interpretation means Feeder 2 cannot directly receive support from F1 or F3.
# If problem implies all feeders can support each other, TIE_SWITCHES_INFO would need redefinition.
# Substation connection points (source nodes for feeders)
# CB1 -> Node 1, CB2 -> Node 23, CB3 -> Node 43
# Node 0 will represent the main grid / infinite source.
SOURCE_NODE = 0
SUBSTATION_CONNECTIONS = {
'CB1': (SOURCE_NODE, 1),
'CB2': (SOURCE_NODE, 23),
'CB3': (SOURCE_NODE, 43)
}
# Capacity of connection from source to substation nodes (effectively feeder capacity)
SUBSTATION_LINE_CAPACITY_A = FEEDER_RATED_CURRENT_A
# --- Data Loading Functions ---
def load_load_data(filename="C题附件：有源配电网62节点系统基本参数.xlsx - 表1 有源配电网62节点系统负荷参数.csv"):
df = pd.read_csv(filename)
df.columns = ['node_id', 'load_kw']
# Convert node_id to int if it's not already
df['node_id'] = df['node_id'].astype(int)
return df.set_index('node_id')['load_kw'].to_dict()
def load_topology_data(filename="C题附件：有源配电网62节点系统基本参数.xlsx - 表2 有源配电网62节点系统拓扑参数.csv"):
df = pd.read_csv(filename)
# Rename columns for easier access (assuming standard Chinese headers)
df.columns = ['line_num', 'from_node', 'to_node', 'length_km', 'resistance_ohm', 'reactance_ohm']
# Convert relevant columns to numeric
for col in ['from_node', 'to_node', 'length_km', 'resistance_ohm', 'reactance_ohm']:
df[col] = pd.to_numeric(df[col], errors='coerce')
return df
# --- Core Power Grid Model Class ---
class PowerGridModel:
def __init__(self, load_data, topology_data, dg_locations_kw, tie_switches_info, substation_connections):
self.loads_kw = load_data
self.topology_df = topology_data
self.dg_kw = dg_locations_kw.copy() # Allow modification for different scenarios
self.tie_switches_info = tie_switches_info
self.substation_connections = substation_connections
self.graph = self._build_graph()
self.feeder_info = self._identify_feeders()
def _build_graph(self):
G = nx.Graph() # Use Graph for undirected, or DiGraph if flow direction is fixed by sources
# Add nodes with load and DG info
all_nodes = set(self.topology_df['from_node']) | set(self.topology_df['to_node'])
for node_id in all_nodes:
node_id = int(node_id) # Ensure int
G.add_node(node_id,
load_kw=self.loads_kw.get(node_id, 0),
dg_kw=self.dg_kw.get(node_id, 0))
# Add lines from topology data
for _, row in self.topology_df.iterrows():
u, v = int(row['from_node']), int(row['to_node'])
G.add_edge(u, v,
id=row['line_num'],
length_km=row['length_km'],
resistance_ohm=row['resistance_ohm'],
# reactance_ohm=row['reactance_ohm'], # Ignoring reactance as per problem
rated_current_a=DEFAULT_LINE_RATED_CURRENT_A, # Default, can be refined
failed=False)
# Add substation connections (virtual lines from a common source)
# These represent the main feeder lines from CBs
G.add_node(SOURCE_NODE, type='source')
for cb_id, (src, dest_node) in self.substation_connections.items():
G.add_edge(src, dest_node, id=cb_id, length_km=0.01, resistance_ohm=0.001, # Minimal impedance
rated_current_a=SUBSTATION_LINE_CAPACITY_A, type='substation_link', failed=False)
return G
def _get_subgraph_with_operational_lines(self, graph_to_copy, faulty_line_edge=None):
"""Creates a subgraph considering only non-failed lines and open tie switches."""
g_op = graph_to_copy.copy()
# Remove failed lines
lines_to_remove = []
if faulty_line_edge: # faulty_line_edge is (u,v)
if g_op.has_edge(*faulty_line_edge):
lines_to_remove.append(faulty_line_edge)
for u, v, data in list(g_op.edges(data=True)):
if data.get('failed', False):
lines_to_remove.append((u,v))
g_op.remove_edges_from(lines_to_remove)
# Normally, tie switches are open. For restoration, specific ones might be closed.
# This base function assumes they are open unless explicitly handled by restoration logic.
return g_op
def _identify_feeders(self):
"""Identifies nodes belonging to each feeder under normal operation (tie switches open)."""
g_normal = self._get_subgraph_with_operational_lines(self.graph)
feeder_info = {} # {'CB1': {nodes}, 'CB2': {nodes}, ...}
for cb_id, (src_node, start_node) in self.substation_connections.items():
if g_normal.has_node(start_node) and g_normal.has_node(src_node) and nx.has_path(g_normal, src_node, start_node):
# Find all nodes reachable from start_node without passing through another substation's start_node
# or the main source node again, after removing other substation links.
temp_g = g_normal.copy()
other_cb_links = []
for other_cb, (s,d) in self.substation_connections.items():
if other_cb != cb_id and temp_g.has_edge(s,d):
other_cb_links.append((s,d))
temp_g.remove_edges_from(other_cb_links)
if nx.has_path(temp_g, src_node, start_node):
# All nodes in the component connected to start_node, excluding the source itself
component_nodes = nx.node_connected_component(temp_g.subgraph(
[n for n in temp_g.nodes if n != src_node or n == start_node] # Consider start_node part of feeder
), start_node)
feeder_info[cb_id] = component_nodes
else: # Should not happen if graph is built correctly
feeder_info[cb_id] = {start_node} if start_node in g_normal else set()
else:
feeder_info[cb_id] = {start_node} if start_node in g_normal else set()
return feeder_info
def _calculate_line_current_kw(self, power_kw):
"""Calculates current (A) given power (kW) at VOLTAGE_KV (line-to-line)."""
if VOLTAGE_KV <= 0: return float('inf')
return abs(power_kw) / (ROOT_3 * VOLTAGE_KV * 1.0) # Assumed PF=1 for current calculation from P
def _get_downstream_info(self, G, line_u, line_v, source_nodes_for_feeder):
"""
Calculates total load and DG power downstream of a directed line (u,v),
assuming v is further from the source_node for this path.
G: graph to operate on (can be a faulted graph)
source_nodes_for_feeder: list of possible source nodes for the current connected component.
"""
# Temporarily make graph directed from source to loads to find downstream nodes
# This is tricky if the graph is not purely radial or has loops after closing ties.
# A simpler approach for radial sections:
# Check connectivity from sources to line_v, if line_u is removed.
# If line_v is disconnected from all sources when (u,v) is cut, then everything
# in the component of line_v is downstream.
# Create a copy of G to modify
temp_g = G.copy()
if not temp_g.has_edge(line_u, line_v):
return [], 0, 0 # Line doesn't exist
temp_g.remove_edge(line_u, line_v)
downstream_nodes = set()
# Check which side (u or v) is disconnected from the source(s)
# Assume u is closer to source, v is further.
# If v is still connected to a source, then (u,v) might be part of a loop or fed from elsewhere.
# A robust way: find path from source to v. If (u,v) is on all paths, then v is downstream of u via this line.
# For radial feeders (normal operation):
# If we consider (u,v) where u is parent of v:
component_of_v = set()
q = [line_v]
visited = {line_u, line_v} # Start by marking u as visited (as if coming from u)
# If line_v is connected to any source_node without passing through line_u
v_connected_to_source_alt_path = False
for src in source_nodes_for_feeder:
if nx.has_path(temp_g, src, line_v):
v_connected_to_source_alt_path = True
break
if v_connected_to_source_alt_path: # (u,v) is part of a loop or v is fed from elsewhere
# This simple downstream logic is insufficient for meshed networks.
# For now, assume radial for this part of flow calculation.
# A more complex flow calculation (Newton-Raphson) would be needed for meshed.
# Given problem constraints, assume feeders are normally radial.
pass # This line might not have a clear "downstream" if looped.
# If v is disconnected from source when (u,v) is cut, then its component is downstream.
# Check connectivity for v in temp_g (where (u,v) is removed)
v_still_connected = any(nx.has_path(temp_g, src, line_v) for src in source_nodes_for_feeder if src in temp_g)
if not v_still_connected: # v is now isolated from source, so its component is downstream of (u,v)
component_of_v = nx.node_connected_component(temp_g, line_v)
else: # v is still connected, means (u,v) might be redundant or complex.
# Try to determine direction based on distance from source
dist_u = float('inf')
dist_v = float('inf')
for src in source_nodes_for_feeder:
if src not in G: continue
if nx.has_path(G, src, line_u): dist_u = min(dist_u, nx.shortest_path_length(G, src, line_u))
if nx.has_path(G, src, line_v): dist_v = min(dist_v, nx.shortest_path_length(G, src, line_v))
if dist_v > dist_u : # v is downstream of u
# Find component of v if (u,v) is removed and v is not connected to source
g_temp_removed_edge = G.copy()
g_temp_removed_edge.remove_edge(line_u,line_v)
is_v_conn_to_src = False
for src_node_feeder in source_nodes_for_feeder:
if src_node_feeder in g_temp_removed_edge and nx.has_path(g_temp_removed_edge, src_node_feeder, line_v):
is_v_conn_to_src = True
break
if not is_v_conn_to_src:
component_of_v = nx.node_connected_component(g_temp_removed_edge, line_v)
else: # v is still connected, (u,v) is likely a loop closing line. Flow is complex.
# For simplicity, this function will return 0 flow for loop lines if direction is ambiguous.
return [], 0, 0
elif dist_u > dist_v: # u is downstream of v (swap them)
# similar logic for u
g_temp_removed_edge = G.copy()
g_temp_removed_edge.remove_edge(line_u,line_v)
is_u_conn_to_src = False
for src_node_feeder in source_nodes_for_feeder:
if src_node_feeder in g_temp_removed_edge and nx.has_path(g_temp_removed_edge, src_node_feeder, line_u):
is_u_conn_to_src = True
break
if not is_u_conn_to_src:
component_of_v = nx.node_connected_component(g_temp_removed_edge, line_u) # component_of_v is actually comp of u
else:
return [], 0, 0
else: # Equidistant or complex, cannot determine simple downstream for this line
return [], 0, 0
downstream_nodes = component_of_v
total_downstream_load_kw = sum(G.nodes[n]['load_kw'] for n in downstream_nodes)
total_downstream_dg_kw = sum(G.nodes[n]['dg_kw'] for n in downstream_nodes if G.nodes[n]['dg_kw'] > 0)
return list(downstream_nodes), total_downstream_load_kw, total_downstream_dg_kw
def calculate_power_flows_and_currents(self, current_graph_state, active_dgs_kw):
"""
Simplified power flow for radial networks or parts of networks.
Returns dict of line flows and currents, and substation powers.
Flows: {(u,v): power_kw} where power_kw > 0 means u to v.
Currents: {(u,v): current_A}
Substation_powers: {'CB1': power_kw_drawn_from_substation}
"""
line_flows_kw = {}
line_currents_a = {}
substation_powers_kw = {}
# Update DG outputs in the graph state
for node_id, dg_val in active_dgs_kw.items():
if node_id in current_graph_state:
current_graph_state.nodes[node_id]['dg_kw'] = dg_val
for node_id in current_graph_state.nodes(): # Reset others if not in active_dgs_kw
if node_id not in active_dgs_kw and 'dg_kw' in current_graph_state.nodes[node_id]:
if current_graph_state.nodes[node_id].get('type') != 'source': # Don't zero out if it was never a DG
current_graph_state.nodes[node_id]['dg_kw'] = 0
# Determine connected components and their sources
# This is a very simplified load flow. It assumes power flows from sources (substations)
# down to loads. DG power reduces the load seen by upstream sections.
# It does not handle loops well without iterative methods (e.g. Hardy Cross or Newton-Raphson).
# For each feeder, calculate flows assuming radial structure
# This is an approximation. A full AC or DC power flow is more accurate.
# Initialize all line flows to 0
for u, v in current_graph_state.edges():
line_flows_kw[(u,v)] = 0
line_flows_kw[(v,u)] = 0 # For bi-directional calculation needs
line_currents_a[(u,v)] = 0
# Iterate multiple times for flow distribution in case of ties or complex paths
# This is a placeholder for a proper iterative flow solution.
# For now, a topological sort based flow for radial parts.
processed_nodes_for_flow_calc = set()
for cb_id, (src_node, start_node) in self.substation_connections.items():
if not current_graph_state.has_node(start_node) or not nx.is_connected(current_graph_state.subgraph([n for n in current_graph_state.nodes() if n != SOURCE_NODE])):
# This feeder might be entirely down if start_node is disconnected from actual nodes
substation_powers_kw[cb_id] = 0
continue
# Get nodes for this feeder (dynamic based on current_graph_state)
# Nodes connected to start_node, excluding SOURCE_NODE, if start_node is connected to SOURCE_NODE
feeder_nodes_component = set()
if current_graph_state.has_edge(src_node, start_node):
temp_g_for_feeder = current_graph_state.copy()
# Remove other substation links to isolate this feeder's component
other_links_to_remove = []
for other_cb, (s,d) in self.substation_connections.items():
if other_cb != cb_id and temp_g_for_feeder.has_edge(s,d):
other_links_to_remove.append((s,d))
temp_g_for_feeder.remove_edges_from(other_links_to_remove)
if nx.has_path(temp_g_for_feeder, src_node, start_node):
try:
# Consider only the part of the graph reachable from start_node, not crossing back to SOURCE_NODE
# except via the designated start_node path.
search_nodes = [n for n in temp_g_for_feeder.nodes if n != src_node]
sub_graph_feeder = temp_g_for_feeder.subgraph(search_nodes)
if start_node in sub_graph_feeder:
feeder_nodes_component = nx.node_connected_component(sub_graph_feeder, start_node)
except nx.NetworkXError: # if start_node not in subgraph
feeder_nodes_component = set()
if not feeder_nodes_component:
substation_powers_kw[cb_id] = 0
continue
# Order nodes from furthest to closest to substation for power accumulation
# This is for radial feeders. If loops exist, this is not sufficient.
# Using BFS layers from start_node
# Net load at each node (Load - DG)
node_net_power_kw = {}
for node in feeder_nodes_component:
node_net_power_kw[node] = current_graph_state.nodes[node]['load_kw'] - current_graph_state.nodes[node]['dg_kw']
# Accumulate power up towards the substation
# This requires a tree traversal (e.g., DFS post-order traversal) from leaves to root.
# For simplicity, if the feeder is a tree rooted at start_node:
if nx.is_tree(current_graph_state.subgraph(feeder_nodes_component | {start_node})): # Check if it's a tree
# Create a directed tree towards the source for easier traversal
# This part is complex if graph is not a tree.
# For now, sum all net loads on the feeder as the substation power (approximation)
total_feeder_net_load = sum(node_net_power_kw[n] for n in feeder_nodes_component)
substation_powers_kw[cb_id] = total_feeder_net_load
# Distribute this flow down the lines (highly simplified)
# A proper method: for each line, sum net_power of all nodes in subtree rooted by that line.
# This simplified flow calculation is a major placeholder.
# For overload risk, we need per-line flows.
# Simplified: assume current_graph_state is a tree rooted at start_node for this feeder
# Use BFS to assign flow from start_node downwards
# This is not a full power flow, but an estimation for line loading.
# Build a directed graph for this feeder based on BFS from start_node
# This is just for flow assignment direction.
# Actual flow needs to sum up demands from downstream.
# For each edge (u,v) in the feeder:
# Determine parent (closer to start_node) and child
# Power on (parent, child) = sum of net loads in subtree rooted at child.
# This is complex to implement robustly here without a full flow algorithm.
# Fallback: Use the _get_downstream_info logic if possible, iterate edges
# This needs to be called carefully to avoid double counting or misdirection in non-radial.
# For now, this function will primarily return substation_powers_kw and
# leave detailed line_flows_kw and line_currents_a for a more robust implementation
# or accept its high level of approximation.
# Let's try a slightly better approximation for line flows on a tree:
# For each edge (u,v) in the feeder tree (rooted at start_node)
# Assume u is parent of v. Power(u,v) = sum of net_power for all nodes in subtree of v.
if start_node in feeder_nodes_component: # Should be
try:
dfs_edges = list(nx.dfs_edges(nx.bfs_tree(current_graph_state.subgraph(feeder_nodes_component), start_node), source=start_node))
# Calculate power for each node including its children's power
node_total_subtree_power = node_net_power_kw.copy()
for u, v in reversed(dfs_edges): # From leaves up to root
if u in node_total_subtree_power and v in node_total_subtree_power:
node_total_subtree_power[u] += node_total_subtree_power[v]
for u,v in dfs_edges: # From root down
if v in node_total_subtree_power:
flow = node_total_subtree_power[v]
line_flows_kw[(u,v)] = flow
line_flows_kw[(v,u)] = -flow # Convention for direction
line_currents_a[(u,v)] = self._calculate_line_current_kw(flow)
except Exception as e:
# print(f"Warning: Could not perform tree-based flow for {cb_id} due to {e}")
pass # Keep substation power as sum, line flows might be inaccurate
else: # Not a tree, flow calculation is more complex.
# print(f"Warning: Feeder {cb_id} is not a tree. Simplified flow may be inaccurate.")
total_feeder_net_load = sum(node_net_power_kw[n] for n in feeder_nodes_component if n in node_net_power_kw)
substation_powers_kw[cb_id] = total_feeder_net_load
# Line flows in meshed networks require iterative solvers.
# For now, this part will be very approximate for meshed sections.
# Handle reverse power flow and inter-feeder DG adjustment
# "分布式能源不得向上级电网倒送功率"
# "可以在相邻馈线间进行调节"
for cb_id in list(substation_powers_kw.keys()):
if substation_powers_kw[cb_id] < 0: # Reverse power flow
excess_dg_on_feeder = -substation_powers_kw[cb_id]
# Try to transfer to other feeders via inter-feeder tie lines
# This logic is complex and needs careful state management.
# For Q1, a simpler approach might be DG curtailment on that feeder.
# Find inter-feeder tie switches connected to this feeder
# Example: S62-3 connects Feeder 3 (node 62) to Feeder 1 (node 19)
# If Feeder 1 has excess_dg, it might try to send to Feeder 3 via (19,62)
# This is an advanced feature. For now, assume DG curtailment if倒送.
# To implement curtailment: identify DGs on this feeder, reduce their output
# proportionally until substation_powers_kw[cb_id] >= 0.
# This would require re-calculating flows.
# For now, just flag it.
# print(f"Warning: Reverse power flow on {cb_id} of {substation_powers_kw[cb_id]} kW. DG curtailment or transfer needed.")
# A simple curtailment:
dG_on_feeder_nodes = [n for n in self.feeder_info.get(cb_id, set()) if n in active_dgs_kw and active_dgs_kw[n] > 0]
total_dg_cap_on_feeder = sum(active_dgs_kw[n] for n in dG_on_feeder_nodes)
if total_dg_cap_on_feeder > 0:
curtail_ratio = min(1.0, excess_dg_on_feeder / total_dg_cap_on_feeder) if total_dg_cap_on_feeder >0 else 0
for dg_node in dG_on_feeder_nodes:
active_dgs_kw[dg_node] *= (1 - curtail_ratio)
# Flows need to be recalculated after curtailment. This suggests an iterative solution.
# For this submission, we'll assume this check is done *before* final flow calc,
# or simply note the violation.
# To avoid recursion here, this function should ideally take DGs as fixed input.
# The adjustment logic should be outside or iterative.
pass
return line_flows_kw, line_currents_a, substation_powers_kw
def calculate_overload_risk(self):
"""
Calculates overload risk for the current DG setup.
Assumes DGs are at their BASE_DG_CAPACITY_KW.
"""
# Get current operational graph (no faults, ties normally open)
g_op = self._get_subgraph_with_operational_lines(self.graph)
# Calculate power flows and currents
# Need to handle DG outputs properly.
current_dg_outputs = self.dg_kw.copy() # Use the model's current DG settings
# Iterative step for DG curtailment if reverse power flow:
# This is a simplified loop. A more robust solution uses optimization or better heuristics.
for _iter in range(3): # Max 3 iterations for adjustment
line_flows, line_currents, substation_powers = self.calculate_power_flows_and_currents(g_op, current_dg_outputs)
reverse_power_detected = False
for cb_id, power_kw in substation_powers.items():
if power_kw < -1e-3: # Small threshold for倒送
reverse_power_detected = True
# print(f"Info: Reverse power on {cb_id} ({power_kw:.2f} kW). Attempting curtailment.")
excess_dg_on_feeder = abs(power_kw)
feeder_nodes_for_cb = self.feeder_info.get(cb_id, set())
dG_on_feeder_nodes = [n for n in feeder_nodes_for_cb if n in current_dg_outputs and current_dg_outputs[n] > 0]
total_dg_cap_on_feeder = sum(current_dg_outputs[n] for n in dG_on_feeder_nodes)
if total_dg_cap_on_feeder > 1e-3 : # Avoid division by zero
curtail_amount_total = excess_dg_on_feeder
for dg_node in dG_on_feeder_nodes:
# Proportional curtailment
proportion = current_dg_outputs[dg_node] / total_dg_cap_on_feeder
curtail_this_dg = proportion * curtail_amount_total
current_dg_outputs[dg_node] = max(0, current_dg_outputs[dg_node] - curtail_this_dg)
else: # No DG to curtail, reverse power might be from other sources or model issue
pass
if not reverse_power_detected:
break
# Final flows after potential curtailment
line_flows, line_currents, substation_powers = self.calculate_power_flows_and_currents(g_op, current_dg_outputs)
overloaded_lines = []
for u, v, data in g_op.edges(data=True):
if data.get('type') == 'substation_link': continue # Don't check substation links themselves for overload here
current = line_currents.get((u,v), 0)
# If flow is from v to u, current might be stored as current_uv = -current_vu
# Take absolute value of flow for current calculation, or ensure current is always positive.
# The _calculate_line_current_kw uses abs(power_kw) so current should be positive.
rated_current = data.get('rated_current_a', DEFAULT_LINE_RATED_CURRENT_A)
if current > 1.1 * rated_current:
overloaded_lines.append({'edge': (u,v), 'current': current, 'rated': rated_current, 'over_by_%': (current/(1.1*rated_current)-1)*100 if rated_current>0 else float('inf')})
if overloaded_lines:
# print(f"System Overload Detected. Overloaded lines: {overloaded_lines}")
# P_over = 1 (for this deterministic scenario)
# C_over = fixed penalty or sum of penalties
risk_overload = 1.0 * COST_PENALTY_FOR_ANY_OVERLOAD
# Or, sum of consequences for each overloaded line, if C_over is per line.
# risk_overload = sum(some_cost_function(ol['over_by_%']) for ol in overloaded_lines)
else:
# print("System is NOT overloaded in the base case.")
risk_overload = 0.0
return risk_overload, overloaded_lines, substation_powers, current_dg_outputs
def calculate_load_loss_risk(self):
"""
Calculates total load loss risk by considering single line faults.
R_loss = sum(P_fault_i * C_loss_i)
P_fault_i = annual probability of fault i
C_loss_i = consequence of fault i (e.g., unserved_load_kw * COST_VOLL_PER_KW)
"""
total_load_loss_risk = 0.0
detailed_fault_impacts = []
# Iterate through all operational lines (excluding substation virtual links for fault simulation)
original_edges = [ (u,v,data) for u,v,data in self.graph.edges(data=True)
if data.get('type') != 'substation_link' and not data.get('is_tie', False)]
for u_fault, v_fault, line_data_faulted in original_edges:
faulty_edge = (u_fault, v_fault)
line_length_km = line_data_faulted.get('length_km', 0)
# Probability of this specific line failing (annual)
# Assuming failure rates are independent and this is the probability of this line being the one to fail.
prob_line_fault = line_length_km * FAILURE_RATE_LINE_PER_KM
if prob_line_fault == 0: continue
# --- Simulate fault ---
g_faulted = self.graph.copy()
if not g_faulted.has_edge(*faulty_edge): continue
g_faulted.edges[faulty_edge]['failed'] = True # Mark as failed
g_after_fault_isolation = self._get_subgraph_with_operational_lines(g_faulted, faulty_edge=faulty_edge)
# --- Identify initial load loss ---
initial_unserved_load_kw = 0
disconnected_load_nodes = {} # {node: load_kw}
# Check connectivity for all load nodes
for node_id, load_kw in self.loads_kw.items():
if load_kw <= 0: continue
is_connected_to_source = False
for cb_id, (src_node, start_node) in self.substation_connections.items():
if nx.has_path(g_after_fault_isolation, src_node, node_id):
is_connected_to_source = True
break
if not is_connected_to_source:
initial_unserved_load_kw += load_kw
disconnected_load_nodes[node_id] = load_kw
if initial_unserved_load_kw == 0: # Fault does not cause load loss (e.g. redundant line)
detailed_fault_impacts.append({'fault': faulty_edge, 'unserved_kw_initial': 0, 'unserved_kw_final':0, 'restored_kw':0, 'risk_contrib':0})
continue
# --- Attempt restoration via tie lines ---
# This is a complex part. Needs to:
# 1. Identify disconnected areas and loads.
# 2. Identify available tie switches that can connect these areas to healthy feeders.
# 3. Check capacity of tie lines and the supporting feeder.
# 4. Prioritize restoration (e.g., maximize load restored).
# For this model, a simplified restoration:
# Iterate over available tie switches. If closing one helps, simulate it.
# This should be greedy or more optimized.
restored_load_kw_total_for_this_fault = 0
# Create a graph state for restoration attempts
g_for_restoration = g_after_fault_isolation.copy()
# Sort disconnected loads by size (optional, for prioritization)
sorted_disconnected_loads = sorted(disconnected_load_nodes.items(), key=lambda item: item[1], reverse=True)
# Try closing tie switches one by one (if they connect a live part to a dead part)
# This is a very simplified greedy approach.
# A proper approach would evaluate all combinations or use optimization.
# Identify current live sources/feeders
live_feeder_sources = [] # (cb_id, start_node_of_live_feeder)
for cb_id, (src,start) in self.substation_connections.items():
if nx.has_path(g_for_restoration, src, start): # Check if substation itself is connected
live_feeder_sources.append(start)
for tie in self.tie_switches_info:
tie_n1, tie_n2 = tie['nodes']
tie_capacity_a = tie['capacity_A']
if not g_for_restoration.has_node(tie_n1) or not g_for_restoration.has_node(tie_n2): continue
if g_for_restoration.has_edge(tie_n1, tie_n2): continue # Already closed or part of main graph (should not be for ties)
# Check if one end is live and other is dead (or part of the disconnected component)
tie_n1_is_live = any(nx.has_path(g_for_restoration, src, tie_n1) for src in live_feeder_sources)
tie_n2_is_live = any(nx.has_path(g_for_restoration, src, tie_n2) for src in live_feeder_sources)
if tie_n1_is_live == tie_n2_is_live: continue # Both live or both dead, closing doesn't bridge outage for now
live_tie_node, dead_tie_node = (tie_n1, tie_n2) if tie_n1_is_live else (tie_n2, tie_n1)
# Check if dead_tie_node is part of the current outage we are trying to fix
# This requires knowing which component dead_tie_node belongs to.
# For now, assume if it's not live, it's part of some outage.
# Simulate closing this tie switch
g_for_restoration.add_edge(live_tie_node, dead_tie_node, id=tie['id'], type='tie_closed',
rated_current_a=tie_capacity_a, resistance_ohm=0.001, length_km=0.01)
# Check how much load can be restored through this tie without overloading tie or new path
# This requires a flow calculation on g_for_restoration.
# Simplified: Check loads now connected.
newly_restored_load_kw_this_tie = 0
temp_restored_nodes_this_tie = []
for node_id, load_val in disconnected_load_nodes.items():
if node_id not in g_for_restoration: continue # Should not happen
# Check if this node is now connected to ANY source
is_now_connected = any(nx.has_path(g_for_restoration, src, node_id) for src in live_feeder_sources)
if is_now_connected and node_id not in temp_restored_nodes_this_tie: # And not already counted as restored by previous ties
# More checks needed:
# 1. Tie line capacity: Power through (live_tie_node, dead_tie_node) <= tie_capacity_a
# 2. Path capacity on the live feeder.
# This is where the simplified flow becomes a bottleneck.
# For now, assume if connected, it can be restored up to a certain limit.
# This is a MAJOR simplification.
newly_restored_load_kw_this_tie += load_val
temp_restored_nodes_this_tie.append(node_id)
# Here, we'd need to check if adding newly_restored_load_kw_this_tie overloads the tie or feeder.
# If current_through_tie > tie_capacity_a, then not all of this load can be restored.
# This part needs a proper constrained flow allocation.
# For now, let's assume a fraction can be restored if connected, or all if small.
# This is a placeholder for a more robust restoration algorithm.
# Let's assume, for now, if connected, it's restored. If this overloads things,
# the overload risk model should capture it (but that's for normal state).
# Here, the goal is to minimize unserved load.
# If this tie leads to overload, we shouldn't use it or only partially.
# For now, naively accept all newly connected load.
if newly_restored_load_kw_this_tie > 0:
restored_load_kw_total_for_this_fault += newly_restored_load_kw_this_tie
# Update disconnected_load_nodes:
for r_node in temp_restored_nodes_this_tie:
if r_node in disconnected_load_nodes:
del disconnected_load_nodes[r_node] # No longer disconnected
else: # Closing this tie didn't help, revert
if g_for_restoration.has_edge(live_tie_node, dead_tie_node):
g_for_restoration.remove_edge(live_tie_node, dead_tie_node)
# Final unserved load for this fault scenario
final_unserved_load_kw = initial_unserved_load_kw - restored_load_kw_total_for_this_fault
final_unserved_load_kw = max(0, final_unserved_load_kw) # Cannot be negative
consequence_c_loss = final_unserved_load_kw * COST_VOLL_PER_KW
risk_contribution = prob_line_fault * consequence_c_loss
total_load_loss_risk += risk_contribution
detailed_fault_impacts.append({
'fault_type': 'line', 'component_id': faulty_edge, 'prob_fault': prob_line_fault,
'unserved_kw_initial': initial_unserved_load_kw,
'restored_kw': restored_load_kw_total_for_this_fault,
'unserved_kw_final': final_unserved_load_kw,
'consequence_c_loss': consequence_c_loss,
'risk_contribution': risk_contribution
})
# TODO: Add DG faults, Switch faults, User faults if they cause wider outages.
# For DG faults: prob_dg_fault = FAILURE_RATE_DG_PERCENT
# A DG fault primarily impacts system's ability to meet load or avoid overload.
# It doesn't directly cause load loss unless it's islanded and the DG is the only source.
# The problem implies grid-connected DGs.
return total_load_loss_risk, detailed_fault_impacts
# --- Main Execution ---
if __name__ == '__main__':
print("--- 配电网风险评估模型 Q1 ---")
# 1. Load Data
print("\n1. 加载数据...")
loads = load_load_data()
topology = load_topology_data()
# print(f"负荷数据: {len(loads)} 点")
# print(f"拓扑数据: {len(topology)} 条线路")
# 2. Initialize Power Grid Model
print("\n2. 初始化电网模型...")
grid = PowerGridModel(loads, topology, DG_LOCATIONS_KW, TIE_SWITCHES_INFO, SUBSTATION_CONNECTIONS)
# print(f"电网图: {grid.graph.number_of_nodes()} 个节点, {grid.graph.number_of_edges()} 条边")
# print(f"馈线信息: {grid.feeder_info}")
# --- 问题1: 失负荷风险和过负荷风险计算模型 ---
print("\n--- 问题1: 风险计算 ---")
# A. 过负荷风险模型 (R_over = P_over * C_over)
# For Q1, DGs are at BASE_DG_CAPACITY_KW. This is a deterministic check for this state.
# P_over = 1 if overload occurs, 0 otherwise. C_over is the penalty.
print("\nA. 计算过负荷风险...")
# Note: The calculate_power_flows_and_currents is highly simplified.
# Results for overload depend heavily on its accuracy and line ratings.
try:
risk_overload, overloaded_lines_details, substation_p, final_dg_out = grid.calculate_overload_risk()
print(f" 计算得到的过负荷风险 (R_over): ${risk_overload:.2f}")
if overloaded_lines_details:
print(f" 检测到过负荷线路 ({len(overloaded_lines_details)} 条):")
# for ol in overloaded_lines_details[:3]: # Print first 3
# print(f" - 线路 {ol['edge']}, 电流: {ol['current']:.2f}A, 额定: {ol['rated']:.2f}A, 超出: {ol['over_by_%']:.2f}%")
else:
print(" 在当前DG配置下，未检测到线路过负荷。")
# print(f" 变电站出口功率 (kW): {substation_p}")
# print(f" 最终DG出力 (kW) (可能经过削减): {final_dg_out}")
except Exception as e:
print(f" 计算过负荷风险时发生错误: {e}")
risk_overload = -1 # Indicate error
# B. 失负荷风险模型 (R_loss = sum(P_fault_i * C_loss_i))
print("\nB. 计算失负荷风险...")
# Note: Restoration logic is simplified.
try:
total_r_loss, fault_details = grid.calculate_load_loss_risk()
print(f" 计算得到的总失负荷风险 (R_loss): ${total_r_loss:.2f} (基于所选成本)")
# print("\n 部分故障场景详情:")
# for fd in fault_details[:3]: # Print first 3
# print(f" - 故障线路: {fd.get('component_id')}, "
# f"初始失负荷: {fd.get('unserved_kw_initial'):.2f} kW, "
# f"最终失负荷: {fd.get('unserved_kw_final'):.2f} kW, "
# f"风险贡献: ${fd.get('risk_contribution'):.2f}")
except Exception as e:
print(f" 计算失负荷风险时发生错误: {e}")
total_r_loss = -1 # Indicate error
print("\n--- 模型执行完毕 ---")
print("注意: 此模型包含多项简化和假设 (如线路额定电流, 成本参数, 潮流计算简化, 恢复逻辑简化).")
print("结果的准确性取决于这些假设的合理性和参数的精确性。")
对此代码进行分析
扫一扫
7461
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
