diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index ad6317c13b..5580a94340 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp @@ -16,6 +16,7 @@ #include "../common/exception.hpp" #include +#include #include #include #include @@ -23,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -256,10 +258,377 @@ inline void complete_bidirectional_neighbourhood_info(std::vector(n)), rank_(static_cast(n), Idx{0}) { + std::iota(parent_.begin(), parent_.end(), Idx{0}); + } + + // Find the representative of x, applying path compression. + Idx find(Idx x) { + while (parent_[x] != x) { + parent_[x] = parent_[parent_[x]]; + x = parent_[x]; + } + return x; + } + + // Unite the sets containing a and b. Returns true if they were distinct. + bool unite(Idx node_a, Idx node_b) { + Idx const root_a = find(node_a); + Idx const root_b = find(node_b); + if (root_a == root_b) { + return false; + } + // form union by rank + if (rank_[root_a] < rank_[root_b]) { + parent_[root_a] = root_b; + } else if (rank_[root_a] > rank_[root_b]) { + parent_[root_b] = root_a; + } else { + parent_[root_b] = root_a; + ++rank_[root_a]; + } + return true; + } + + private: + std::vector parent_; + std::vector rank_; +}; + +// Contract every branch that carries its own native flow measurement, merging +// its two end buses. The resulting partition is the set of observable +// components and is independent of bus or neighbour ordering. +inline DisjointSet contract_branch_measured_edges(std::vector const& neighbour_list) { + auto const n_bus = static_cast(neighbour_list.size()); + DisjointSet components{n_bus}; + for (Idx bus = 0; bus < n_bus; ++bus) { + for (auto const& neighbour : neighbour_list[bus].direct_neighbours) { + if (branch_has_native_measurement(neighbour.status)) { + components.unite(bus, neighbour.bus); + } + } + } + return components; +} + +// Count the number of distinct components in the disjoint set over [0, n_bus). +inline Idx count_components(DisjointSet& components, Idx n_bus) { + Idx count = 0; + for (Idx bus = 0; bus < n_bus; ++bus) { + if (components.find(bus) == bus) { + ++count; + } + } + return count; +} + +// An unmeasured branch that still connects two distinct observable components. +// It can be added to the spanning tree only by consuming a nodal (injection) +// measurement at one of its two end buses. Branches that carry their own flow +// measurement do not appear here: they were already contracted into components. +struct CandidateEdge { + Idx from_bus; // original from-bus + Idx to_bus; // original to-bus + Idx from_component; // compact component index of from_bus + Idx to_component; // compact component index of to_bus +}; + +// The network reduced to its observable components. Branch-measured edges have +// been contracted away, leaving the super-nodes (components), the unmeasured +// branches that still join distinct components, and the buses that hold an +// unused nodal measurement available for assignment. +struct ContractedNetwork { + Idx n_components{}; + std::vector component_of_bus; // bus -> compact component index in [0, n_components) + std::vector candidate_edges; + std::vector bus_has_injection; // 1 if the bus has an unused nodal measurement +}; + +// Reduce the network to its observable components: compact the contracted +// component representatives to dense indices, record per-bus injection +// availability, and collect the unmeasured branches that join distinct +// components. The result depends only on the graph and its measurements, not on +// bus or neighbour ordering. +inline ContractedNetwork build_contracted_network(std::vector const& neighbour_list, + DisjointSet& components) { + auto const n_bus = static_cast(neighbour_list.size()); + ContractedNetwork result; + result.component_of_bus.assign(static_cast(n_bus), -1); + result.bus_has_injection.assign(static_cast(n_bus), 0); + + // Compact component representatives to dense indices [0, n_components). + std::vector representative_to_compact(static_cast(n_bus), -1); + for (Idx bus = 0; bus < n_bus; ++bus) { + Idx const representative = components.find(bus); + if (representative_to_compact[representative] == -1) { + representative_to_compact[representative] = result.n_components++; + } + result.component_of_bus[bus] = representative_to_compact[representative]; + if (neighbour_list[bus].status == ConnectivityStatus::node_measured) { + result.bus_has_injection[bus] = 1; + } + } + + // Collect candidate edges. The neighbour list is bidirectional, so each + // undirected branch is taken once (from the lower-indexed bus). Branches + // whose endpoints already share a component would be self-loops in the + // contracted graph and are dropped. + for (Idx bus = 0; bus < n_bus; ++bus) { + for (auto const& neighbour : neighbour_list[bus].direct_neighbours) { + if (bus >= neighbour.bus || branch_has_native_measurement(neighbour.status)) { + continue; + } + Idx const from_component = result.component_of_bus[bus]; + Idx const to_component = result.component_of_bus[neighbour.bus]; + if (from_component == to_component) { + continue; + } + result.candidate_edges.push_back({.from_bus = bus, + .to_bus = neighbour.bus, + .from_component = from_component, + .to_component = to_component}); + } + } + + return result; +} + +// Graphic-matroid independence: do the given candidate edges form a forest over +// the components? +inline bool contracted_edges_form_forest(ContractedNetwork const& net, std::vector const& edge_indices) { + DisjointSet components{net.n_components}; + for (Idx const idx : edge_indices) { + auto const& edge = net.candidate_edges[idx]; + if (!components.unite(edge.from_component, edge.to_component)) { + return false; // closing a cycle + } + } + return true; +} + +// Transversal-matroid augmenting step: try to assign the given candidate edge to +// an injection bus at one of its two endpoints, reassigning earlier edges along +// an alternating path if necessary. +inline bool assign_edge_to_injection(ContractedNetwork const& net, Idx edge_idx, std::vector& bus_to_edge, + std::vector& visited_bus) { + auto const& edge = net.candidate_edges[edge_idx]; + std::array const endpoint_buses{edge.from_bus, edge.to_bus}; + for (Idx const bus : endpoint_buses) { + if (net.bus_has_injection[bus] == 0 || visited_bus[bus] != 0) { + continue; + } + visited_bus[bus] = 1; + if (bus_to_edge[bus] == -1 || assign_edge_to_injection(net, bus_to_edge[bus], bus_to_edge, visited_bus)) { + bus_to_edge[bus] = edge_idx; + return true; + } + } + return false; +} + +// Transversal-matroid independence: can every candidate edge be matched to a +// distinct injection bus at one of its endpoints? +inline bool contracted_edges_transversal_independent(ContractedNetwork const& net, + std::vector const& edge_indices) { + auto const n_bus = static_cast(net.component_of_bus.size()); + std::vector bus_to_edge(static_cast(n_bus), -1); + for (Idx const idx : edge_indices) { + std::vector visited_bus(static_cast(n_bus), 0); + if (!assign_edge_to_injection(net, idx, bus_to_edge, visited_bus)) { + return false; + } + } + return true; +} + +// Build the edge-index set I, optionally adding one edge and removing one. +inline std::vector edge_set_with(std::vector const& in_set, Idx n_edges, Idx add_edge, + Idx remove_edge) { + std::vector result; + for (Idx e = 0; e < n_edges; ++e) { + if (e != remove_edge && (in_set[e] != 0 || e == add_edge)) { + result.push_back(e); + } + } + return result; +} + +// Sources X1: edges that can extend the forest; sinks X2: edges that can extend +// the injection assignment. +struct ExchangeGraphSources { + std::vector can_extend_forest; + std::vector can_extend_assignment; +}; + +// Classify every edge not in I by whether adding it keeps a forest and/or a +// valid injection assignment. +inline ExchangeGraphSources classify_exchange_edges(ContractedNetwork const& net, + std::vector const& in_set, Idx n_edges) { + ExchangeGraphSources sources{.can_extend_forest = std::vector(static_cast(n_edges), 0), + .can_extend_assignment = + std::vector(static_cast(n_edges), 0)}; + for (Idx x = 0; x < n_edges; ++x) { + if (in_set[x] != 0) { + continue; + } + sources.can_extend_forest[x] = contracted_edges_form_forest(net, edge_set_with(in_set, n_edges, x, -1)) ? 1 : 0; + sources.can_extend_assignment[x] = + contracted_edges_transversal_independent(net, edge_set_with(in_set, n_edges, x, -1)) ? 1 : 0; + } + return sources; +} + +// Length-zero augmentation: directly add an edge independent in both matroids. +// Returns true and updates in_set if such an edge exists. +inline bool try_length_zero_augmentation(std::vector& in_set, Idx n_edges, + ExchangeGraphSources const& sources) { + for (Idx x = 0; x < n_edges; ++x) { + if (in_set[x] == 0 && sources.can_extend_forest[x] != 0 && sources.can_extend_assignment[x] != 0) { + in_set[x] = 1; + return true; + } + } + return false; +} + +// Breadth-first search for a shortest augmenting path from X1 to X2 in the +// exchange graph. Arcs leaving an edge x not in I go to edges y in I such that +// I - y + x stays assignable; arcs leaving an edge y in I go to edges x not in I +// such that I - y + x stays a forest. Fills predecessor and returns the sink +// edge index, or -1 if no augmenting path exists. +inline Idx find_augmenting_sink(ContractedNetwork const& net, std::vector const& in_set, Idx n_edges, + ExchangeGraphSources const& sources, std::vector& predecessor) { + std::vector visited(static_cast(n_edges), 0); + std::queue bfs_queue; + for (Idx x = 0; x < n_edges; ++x) { + if (in_set[x] == 0 && sources.can_extend_forest[x] != 0) { + visited[x] = 1; + bfs_queue.push(x); + } + } + + // Whether an arc current -> other exists in the exchange graph. + auto is_reachable = [&](Idx current, Idx other) { + if (in_set[current] == 0) { + return in_set[other] != 0 && + contracted_edges_transversal_independent(net, edge_set_with(in_set, n_edges, current, other)); + } + return in_set[other] == 0 && contracted_edges_form_forest(net, edge_set_with(in_set, n_edges, other, current)); + }; + + // An out-of-set edge that can extend the assignment and was actually reached + // (pure sources are already known not to be sinks) terminates the path. + auto is_sink = [&](Idx current) { + return in_set[current] == 0 && sources.can_extend_assignment[current] != 0 && predecessor[current] != -1; + }; + + while (!bfs_queue.empty()) { + Idx const current = bfs_queue.front(); + bfs_queue.pop(); + if (is_sink(current)) { + return current; + } + for (Idx other = 0; other < n_edges; ++other) { + if (visited[other] == 0 && is_reachable(current, other)) { + visited[other] = 1; + predecessor[other] = current; + bfs_queue.push(other); + } + } + } + return -1; +} + +// Flip membership of every edge along the augmenting path ending at sink. +inline void flip_augmenting_path(std::vector& in_set, std::vector const& predecessor, Idx sink) { + for (Idx node = sink; node != -1; node = predecessor[node]) { + in_set[node] = (in_set[node] == 0) ? std::uint8_t{1} : std::uint8_t{0}; + } +} + +// Try to enlarge the current common independent set I (candidate edges that +// simultaneously form a forest over the components and can be matched to +// distinct injection buses) by one element, using a shortest augmenting path in +// the matroid-intersection exchange graph. Returns true and updates in_set on +// success. +inline bool grow_assignable_forest(ContractedNetwork const& net, std::vector& in_set) { + auto const n_edges = static_cast(net.candidate_edges.size()); + ExchangeGraphSources const sources = classify_exchange_edges(net, in_set, n_edges); + + // First try a direct (length-zero) augmentation. + if (try_length_zero_augmentation(in_set, n_edges, sources)) { + return true; + } + + // Otherwise search for a shortest augmenting path and flip it. + std::vector predecessor(static_cast(n_edges), -1); + Idx const sink = find_augmenting_sink(net, in_set, n_edges, sources, predecessor); + if (sink == -1) { + return false; + } + flip_augmenting_path(in_set, predecessor, sink); + return true; +} + +// Decide observability of the contracted network: does a spanning tree over the +// components exist whose every edge can be assigned a distinct injection bus? +inline bool is_contracted_network_observable(ContractedNetwork const& net) { + if (net.n_components <= 1) { + return true; + } + std::vector in_set(net.candidate_edges.size(), 0); + Idx const target = net.n_components - 1; + Idx selected = 0; + while (selected < target && grow_assignable_forest(net, in_set)) { + ++selected; + } + return selected == target; +} + +// Order-independent meshed observability check based on matroid intersection of +// the network's graphic matroid and the measurement-assignment transversal +// matroid. +inline bool meshed_observable_matroid_intersection(std::vector const& neighbour_list) { + auto components = contract_branch_measured_edges(neighbour_list); + auto const net = build_contracted_network(neighbour_list, components); + return is_contracted_network_observable(net); +} + +// To be deprecated post matroid intersection implementation. inline void prepare_starting_nodes(std::vector const& neighbour_list, Idx n_bus, std::vector& starting_candidates) { - // First find a list of starting points. These are nodes without measurements and all edges connecting to it has no - // edge measurements. + // First collect nodes without nodal measurements and no measured incident edges. + // Then append nodes without nodal measurements that do have measured incident edges. + std::vector secondary_candidates; + secondary_candidates.reserve(static_cast(n_bus)); + for (Idx bus = 0; bus < n_bus; ++bus) { if (neighbour_list[bus].status == ConnectivityStatus::has_no_measurement) { bool all_neighbours_no_edge_measurement = true; @@ -271,18 +640,13 @@ inline void prepare_starting_nodes(std::vector const& neig } if (all_neighbours_no_edge_measurement) { starting_candidates.push_back(bus); + } else { + secondary_candidates.push_back(bus); } } } - // If no such starting point, find nodes without measurements - if (starting_candidates.empty()) { - for (Idx bus = 0; bus < n_bus; ++bus) { - if (neighbour_list[bus].status == ConnectivityStatus::has_no_measurement) { - starting_candidates.push_back(bus); - } - } - } + starting_candidates.insert(starting_candidates.end(), secondary_candidates.begin(), secondary_candidates.end()); // If no nodes without measurements, start from first node // (but network should be observable, so this is just a fallback) @@ -300,6 +664,7 @@ struct StatusModification { ConnectivityStatus old_value; // original value before modification }; +// To be deprecated post matroid intersection implementation. // Context struct to hold shared state during spanning tree search struct SpanningTreeContext { std::vector* neighbour_list; @@ -319,6 +684,7 @@ struct SpanningTreeContext { } }; +// To be deprecated post matroid intersection implementation. // Helper function: Try to traverse edges with native measurements inline bool try_native_edge_measurements(SpanningTreeContext& ctx, bool& step_success) { if ((*ctx.visited)[ctx.current_bus] == std::to_underlying(BusVisited::NotVisited)) { @@ -352,6 +718,7 @@ inline bool try_native_edge_measurements(SpanningTreeContext& ctx, bool& step_su return false; } +// To be deprecated post matroid intersection implementation. // Helper function: Try to use downwind measurement inline bool try_downwind_measurement(SpanningTreeContext& ctx, bool& step_success, bool current_bus_no_measurement) { if (!current_bus_no_measurement && ctx.downwind) { @@ -389,6 +756,7 @@ inline bool try_downwind_measurement(SpanningTreeContext& ctx, bool& step_succes return false; } +// To be deprecated post matroid intersection implementation. // Helper function: Process a single edge during general connection rules inline bool process_edge(SpanningTreeContext& ctx, BusNeighbourhoodInfo::neighbour& neighbour, ConnectivityStatus neighbour_status, ConnectivityStatus reverse_status, bool use_current_node, @@ -425,6 +793,7 @@ inline bool process_edge(SpanningTreeContext& ctx, BusNeighbourhoodInfo::neighbo return true; } +// To be deprecated post matroid intersection implementation. // Helper function: Try general connection rules inline bool try_general_connection_rules(SpanningTreeContext& ctx, bool& step_success, bool current_bus_no_measurement) { @@ -450,6 +819,7 @@ inline bool try_general_connection_rules(SpanningTreeContext& ctx, bool& step_su return false; } +// To be deprecated post matroid intersection implementation. // Helper function: Reassign nodal measurement between two connected nodes inline void reassign_nodal_measurement(SpanningTreeContext& ctx, Idx from_node, Idx to_node) { // no reassignment possible if reached via edge measurement @@ -487,6 +857,7 @@ inline void reassign_nodal_measurement(SpanningTreeContext& ctx, Idx from_node, } } +// To be deprecated post matroid intersection implementation. // Helper function: Try backtracking inline bool try_backtrack(SpanningTreeContext& ctx, bool& step_success) { if (!ctx.edge_track->empty()) { @@ -516,6 +887,7 @@ inline bool try_backtrack(SpanningTreeContext& ctx, bool& step_success) { return false; } +// To be deprecated post matroid intersection implementation. inline bool find_spanning_tree_from_node_impl(Idx start_bus, Idx n_bus, std::vector& neighbour_list, std::vector& modifications, @@ -581,6 +953,7 @@ inline bool find_spanning_tree_from_node_impl(Idx start_bus, Idx n_bus, return ctx.visited_count == n_bus; } +// To be deprecated post matroid intersection implementation. // Backward-compatible overload for tests - makes a copy and creates modification tracker inline bool find_spanning_tree_from_node(Idx start_bus, Idx n_bus, std::vector const& neighbour_list) { @@ -594,6 +967,7 @@ inline bool find_spanning_tree_from_node(Idx start_bus, Idx n_bus, edge_track_buffer); } +// To be deprecated post matroid intersection implementation. inline bool sufficient_condition_meshed_without_voltage_phasor(std::vector& neighbour_list) { auto const n_bus = static_cast(neighbour_list.size()); std::vector starting_candidates; @@ -626,6 +1000,7 @@ inline bool sufficient_condition_meshed_without_voltage_phasor(std::vector const& neighbour_list) { @@ -685,8 +1060,11 @@ inline ObservabilityResult observability_check(MeasuredValues const& measur is_sufficient_condition_met = detail::sufficient_condition_radial_with_voltage_phasor( y_bus_structure, observability_sensors, n_voltage_phasor_sensors); } else { - is_sufficient_condition_met = - detail::sufficient_condition_meshed_without_voltage_phasor(bus_neighbourhood_info); + // Order-independent meshed sufficient condition based on matroid intersection. + if (!detail::meshed_observable_matroid_intersection(bus_neighbourhood_info)) { + throw NotObservableError{"Meshed observability check fail. Network unobservable.\n"}; + } + is_sufficient_condition_met = true; } return ObservabilityResult{.is_observable = is_necessary_condition_met && is_sufficient_condition_met, diff --git a/tests/cpp_unit_tests/math_solver/test_observability.cpp b/tests/cpp_unit_tests/math_solver/test_observability.cpp index 0a3e5d3827..fce4950485 100644 --- a/tests/cpp_unit_tests/math_solver/test_observability.cpp +++ b/tests/cpp_unit_tests/math_solver/test_observability.cpp @@ -740,6 +740,291 @@ TEST_CASE("Test Observability - complete_bidirectional_neighbourhood_info") { } } +// Robust (order-independent) meshed observability check: +// contracting branch-measured edges into observable components via union-find. +TEST_CASE("Test Observability - contract_branch_measured_edges") { + using power_grid_model::math_solver::detail::BusNeighbourhoodInfo; + using power_grid_model::math_solver::detail::contract_branch_measured_edges; + using power_grid_model::math_solver::detail::count_components; + using enum power_grid_model::math_solver::detail::ConnectivityStatus; + + SUBCASE("No branch measurements - every bus is its own component") { + std::vector neighbour_list(3); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = has_no_measurement}}; + neighbour_list[1].direct_neighbours = {{.bus = 0, .status = has_no_measurement}, + {.bus = 2, .status = has_no_measurement}}; + neighbour_list[2].direct_neighbours = {{.bus = 1, .status = has_no_measurement}}; + + auto components = contract_branch_measured_edges(neighbour_list); + CHECK(count_components(components, 3) == 3); + CHECK(components.find(0) != components.find(1)); + CHECK(components.find(1) != components.find(2)); + } + + SUBCASE("Measured edges merge their end buses") { + // 0 =meas= 1 --- 2 =meas= 3 + std::vector neighbour_list(4); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}}; + neighbour_list[1].direct_neighbours = {{.bus = 0, .status = branch_native_measurement_unused}, + {.bus = 2, .status = has_no_measurement}}; + neighbour_list[2].direct_neighbours = {{.bus = 1, .status = has_no_measurement}, + {.bus = 3, .status = branch_native_measurement_unused}}; + neighbour_list[3].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + + auto components = contract_branch_measured_edges(neighbour_list); + CHECK(count_components(components, 4) == 2); + CHECK(components.find(0) == components.find(1)); + CHECK(components.find(2) == components.find(3)); + CHECK(components.find(1) != components.find(2)); + } + + SUBCASE("Chain of measured edges collapses to a single component") { + std::vector neighbour_list(4); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}}; + neighbour_list[1].direct_neighbours = {{.bus = 0, .status = branch_native_measurement_unused}, + {.bus = 2, .status = branch_native_measurement_unused}}; + neighbour_list[2].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}, + {.bus = 3, .status = branch_native_measurement_unused}}; + neighbour_list[3].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + + auto components = contract_branch_measured_edges(neighbour_list); + CHECK(count_components(components, 4) == 1); + } + + SUBCASE("Partition is independent of neighbour ordering") { + // Same graph (0 =meas= 1, 2 =meas= 3, plain edge 0--2), but neighbours + // listed in ascending vs. shuffled order. The contraction must agree. + std::vector ascending(4); + ascending[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}, + {.bus = 2, .status = has_no_measurement}}; + ascending[1].direct_neighbours = {{.bus = 0, .status = branch_native_measurement_unused}}; + ascending[2].direct_neighbours = {{.bus = 0, .status = has_no_measurement}, + {.bus = 3, .status = branch_native_measurement_unused}}; + ascending[3].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + + std::vector shuffled(4); + shuffled[0].direct_neighbours = {{.bus = 2, .status = has_no_measurement}, + {.bus = 1, .status = branch_native_measurement_unused}}; + shuffled[1].direct_neighbours = {{.bus = 0, .status = branch_native_measurement_unused}}; + shuffled[2].direct_neighbours = {{.bus = 3, .status = branch_native_measurement_unused}, + {.bus = 0, .status = has_no_measurement}}; + shuffled[3].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + + auto comp_a = contract_branch_measured_edges(ascending); + auto comp_b = contract_branch_measured_edges(shuffled); + + CHECK(count_components(comp_a, 4) == count_components(comp_b, 4)); + CHECK((comp_a.find(0) == comp_a.find(1)) == (comp_b.find(0) == comp_b.find(1))); + CHECK((comp_a.find(2) == comp_a.find(3)) == (comp_b.find(2) == comp_b.find(3))); + CHECK((comp_a.find(0) == comp_a.find(2)) == (comp_b.find(0) == comp_b.find(2))); + } + + SUBCASE("Empty network") { + std::vector const neighbour_list; + auto components = contract_branch_measured_edges(neighbour_list); + CHECK(count_components(components, 0) == 0); + } +} + +TEST_CASE("Test Observability - build_contracted_network") { + using power_grid_model::math_solver::detail::build_contracted_network; + using power_grid_model::math_solver::detail::BusNeighbourhoodInfo; + using power_grid_model::math_solver::detail::contract_branch_measured_edges; + using enum power_grid_model::math_solver::detail::ConnectivityStatus; + + SUBCASE("Square with no measured branches keeps every bus as its own component") { + // 0--1--2--3--0 ring, injections at bus 0 and bus 2 + std::vector neighbour_list(4); + neighbour_list[0].status = node_measured; + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = has_no_measurement}, + {.bus = 3, .status = has_no_measurement}}; + neighbour_list[1].direct_neighbours = {{.bus = 0, .status = has_no_measurement}, + {.bus = 2, .status = has_no_measurement}}; + neighbour_list[2].status = node_measured; + neighbour_list[2].direct_neighbours = {{.bus = 1, .status = has_no_measurement}, + {.bus = 3, .status = has_no_measurement}}; + neighbour_list[3].direct_neighbours = {{.bus = 2, .status = has_no_measurement}, + {.bus = 0, .status = has_no_measurement}}; + + auto components = contract_branch_measured_edges(neighbour_list); + auto const net = build_contracted_network(neighbour_list, components); + + CHECK(net.n_components == 4); + CHECK(net.candidate_edges.size() == 4); // the four ring branches + CHECK(net.bus_has_injection[0] == 1); + CHECK(net.bus_has_injection[1] == 0); + CHECK(net.bus_has_injection[2] == 1); + CHECK(net.bus_has_injection[3] == 0); + // every candidate edge joins two distinct components + for (auto const& edge : net.candidate_edges) { + CHECK(edge.from_component != edge.to_component); + CHECK(net.component_of_bus[edge.from_bus] == edge.from_component); + CHECK(net.component_of_bus[edge.to_bus] == edge.to_component); + } + } + + SUBCASE("Measured branches contract endpoints and reduce candidate edges") { + // 0 =meas= 1 --- 2 =meas= 3, injection at bus 1 + std::vector neighbour_list(4); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}}; + neighbour_list[1].status = node_measured; + neighbour_list[1].direct_neighbours = {{.bus = 0, .status = branch_native_measurement_unused}, + {.bus = 2, .status = has_no_measurement}}; + neighbour_list[2].direct_neighbours = {{.bus = 1, .status = has_no_measurement}, + {.bus = 3, .status = branch_native_measurement_unused}}; + neighbour_list[3].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + + auto components = contract_branch_measured_edges(neighbour_list); + auto const net = build_contracted_network(neighbour_list, components); + + CHECK(net.n_components == 2); + CHECK(net.component_of_bus[0] == net.component_of_bus[1]); + CHECK(net.component_of_bus[2] == net.component_of_bus[3]); + CHECK(net.component_of_bus[1] != net.component_of_bus[2]); + // only the unmeasured 1--2 branch remains as a candidate edge + REQUIRE(net.candidate_edges.size() == 1); + CHECK(net.candidate_edges[0].from_bus == 1); + CHECK(net.candidate_edges[0].to_bus == 2); + CHECK(net.bus_has_injection[1] == 1); + } + + SUBCASE("Unmeasured branch inside a single component is dropped") { + // triangle: 0 =meas= 1 =meas= 2, plus unmeasured 0--2 (a self-loop after contraction) + std::vector neighbour_list(3); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}, + {.bus = 2, .status = has_no_measurement}}; + neighbour_list[1].direct_neighbours = {{.bus = 0, .status = branch_native_measurement_unused}, + {.bus = 2, .status = branch_native_measurement_unused}}; + neighbour_list[2].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}, + {.bus = 0, .status = has_no_measurement}}; + + auto components = contract_branch_measured_edges(neighbour_list); + auto const net = build_contracted_network(neighbour_list, components); + + CHECK(net.n_components == 1); + CHECK(net.candidate_edges.empty()); + } + + SUBCASE("Empty network") { + std::vector const neighbour_list; + auto components = contract_branch_measured_edges(neighbour_list); + auto const net = build_contracted_network(neighbour_list, components); + CHECK(net.n_components == 0); + CHECK(net.candidate_edges.empty()); + CHECK(net.component_of_bus.empty()); + } +} + +TEST_CASE("Test Observability - meshed_observable_matroid_intersection") { + using power_grid_model::math_solver::detail::BusNeighbourhoodInfo; + using power_grid_model::math_solver::detail::complete_bidirectional_neighbourhood_info; + using power_grid_model::math_solver::detail::meshed_observable_matroid_intersection; + using enum power_grid_model::math_solver::detail::ConnectivityStatus; + + SUBCASE("All branches measured and connected is observable") { + std::vector neighbour_list(3); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}}; + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + complete_bidirectional_neighbourhood_info(neighbour_list); + CHECK(meshed_observable_matroid_intersection(neighbour_list) == true); + } + + SUBCASE("Line of unmeasured branches with enough injections is observable") { + // 0 -- 1 -- 2, injections at buses 0 and 1 + std::vector neighbour_list(3); + neighbour_list[0].status = node_measured; + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = has_no_measurement}}; + neighbour_list[1].status = node_measured; + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = has_no_measurement}}; + complete_bidirectional_neighbourhood_info(neighbour_list); + CHECK(meshed_observable_matroid_intersection(neighbour_list) == true); + } + + SUBCASE("Injection reassignment is found") { + // 0 -- 1 -- 2, injections at buses 1 and 2: bus 1 must hand its injection + // to edge 0--1 while edge 1--2 takes the injection at bus 2 + std::vector neighbour_list(3); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = has_no_measurement}}; + neighbour_list[1].status = node_measured; + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = has_no_measurement}}; + neighbour_list[2].status = node_measured; + complete_bidirectional_neighbourhood_info(neighbour_list); + CHECK(meshed_observable_matroid_intersection(neighbour_list) == true); + } + + SUBCASE("Too few injections is not observable") { + // 0 -- 1 -- 2, single injection at bus 1 but two merges are needed + std::vector neighbour_list(3); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = has_no_measurement}}; + neighbour_list[1].status = node_measured; + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = has_no_measurement}}; + complete_bidirectional_neighbourhood_info(neighbour_list); + CHECK(meshed_observable_matroid_intersection(neighbour_list) == false); + } + + SUBCASE("Measured triangle with an uncovered pendant is not observable") { + // measured triangle 0-1-2; pendant bus 3 hangs off bus 0 by an unmeasured + // branch and no injection can cover it + std::vector neighbour_list(4); + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}, + {.bus = 2, .status = branch_native_measurement_unused}, + {.bus = 3, .status = has_no_measurement}}; + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + complete_bidirectional_neighbourhood_info(neighbour_list); + CHECK(meshed_observable_matroid_intersection(neighbour_list) == false); + } + + SUBCASE("Contracted component keeps both internal injections for two ports") { + // A=B=C blob (measured branches 0-1 and 1-2) with injections at A(0) and + // C(2). Bus 3 hangs off A and bus 4 hangs off C via unmeasured branches. + // Both injections must survive the merge to cover the two external ports. + std::vector neighbour_list(5); + neighbour_list[0].status = node_measured; + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}, + {.bus = 3, .status = has_no_measurement}}; + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + neighbour_list[2].status = node_measured; + neighbour_list[2].direct_neighbours = {{.bus = 4, .status = has_no_measurement}}; + complete_bidirectional_neighbourhood_info(neighbour_list); + CHECK(meshed_observable_matroid_intersection(neighbour_list) == true); + } + + SUBCASE("A single internal injection cannot serve two ports") { + // Same blob as above but only A(0) carries an injection: the port at C + // can no longer be covered. + std::vector neighbour_list(5); + neighbour_list[0].status = node_measured; + neighbour_list[0].direct_neighbours = {{.bus = 1, .status = branch_native_measurement_unused}, + {.bus = 3, .status = has_no_measurement}}; + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = branch_native_measurement_unused}}; + neighbour_list[2].direct_neighbours = {{.bus = 4, .status = has_no_measurement}}; + complete_bidirectional_neighbourhood_info(neighbour_list); + CHECK(meshed_observable_matroid_intersection(neighbour_list) == false); + } + + SUBCASE("Verdict is independent of neighbour ordering") { + + auto build = [](bool reversed) { + std::vector neighbour_list(3); + neighbour_list[1].status = node_measured; + neighbour_list[2].status = node_measured; + if (reversed) { + neighbour_list[1].direct_neighbours = {{.bus = 2, .status = has_no_measurement}, + {.bus = 0, .status = has_no_measurement}}; + } else { + neighbour_list[1].direct_neighbours = {{.bus = 0, .status = has_no_measurement}, + {.bus = 2, .status = has_no_measurement}}; + } + complete_bidirectional_neighbourhood_info(neighbour_list); + return neighbour_list; + }; + auto ascending = build(false); + auto reversed = build(true); + CHECK(meshed_observable_matroid_intersection(ascending) == true); + CHECK(meshed_observable_matroid_intersection(ascending) == meshed_observable_matroid_intersection(reversed)); + } +} + // TODO: properly clean up after y-bus access refactoring TEST_CASE("Test Observability - assign_independent_sensors_radial") { using power_grid_model::math_solver::YBusStructure; diff --git a/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/input.json b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/input.json new file mode 100644 index 0000000000..97a8128780 --- /dev/null +++ b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/input.json @@ -0,0 +1,567 @@ +{ + "version": "1.0", + "type": "input", + "is_batch": false, + "attributes": {}, + "data": { + "source": [ + { + "id": 0, + "node": 1, + "status": 1, + "u_ref": 1.0 + } + ], + "node": [ + { + "id": 1, + "u_rated": 10500 + }, + { + "id": 2, + "u_rated": 10500 + }, + { + "id": 3, + "u_rated": 10500 + }, + { + "id": 4, + "u_rated": 10500 + }, + { + "id": 5, + "u_rated": 10500 + }, + { + "id": 6, + "u_rated": 10500 + }, + { + "id": 7, + "u_rated": 10500 + }, + { + "id": 8, + "u_rated": 10500 + }, + { + "id": 9, + "u_rated": 10500 + }, + { + "id": 10, + "u_rated": 10500 + }, + { + "id": 11, + "u_rated": 10500 + }, + { + "id": 12, + "u_rated": 10500 + }, + { + "id": 13, + "u_rated": 10500 + }, + { + "id": 14, + "u_rated": 10500 + } + ], + "line": [ + { + "id": 101, + "from_node": 1, + "to_node": 2, + "from_status": 1, + "to_status": 1, + "r1": 0.1118, + "x1": 0.0112, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 102, + "from_node": 1, + "to_node": 5, + "from_status": 1, + "to_status": 1, + "r1": 0.1015, + "x1": 0.0102, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 103, + "from_node": 2, + "to_node": 5, + "from_status": 1, + "to_status": 1, + "r1": 0.1508, + "x1": 0.0151, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 104, + "from_node": 2, + "to_node": 4, + "from_status": 1, + "to_status": 1, + "r1": 0.3980, + "x1": 0.0398, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 105, + "from_node": 2, + "to_node": 3, + "from_status": 1, + "to_status": 1, + "r1": 0.1309, + "x1": 0.0131, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 106, + "from_node": 3, + "to_node": 4, + "from_status": 1, + "to_status": 1, + "r1": 0.2927, + "x1": 0.0293, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 107, + "from_node": 4, + "to_node": 5, + "from_status": 1, + "to_status": 1, + "r1": 0.5853, + "x1": 0.0585, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 108, + "from_node": 4, + "to_node": 7, + "from_status": 1, + "to_status": 1, + "r1": 0.2117, + "x1": 0.0212, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 109, + "from_node": 4, + "to_node": 9, + "from_status": 1, + "to_status": 1, + "r1": 0.4975, + "x1": 0.0498, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 110, + "from_node": 5, + "to_node": 6, + "from_status": 1, + "to_status": 1, + "r1": 0.1951, + "x1": 0.0195, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 111, + "from_node": 6, + "to_node": 12, + "from_status": 1, + "to_status": 1, + "r1": 0.3292, + "x1": 0.4829, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 112, + "from_node": 6, + "to_node": 13, + "from_status": 1, + "to_status": 1, + "r1": 0.1485, + "x1": 0.0149, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 113, + "from_node": 6, + "to_node": 10, + "from_status": 1, + "to_status": 1, + "r1": 0.3015, + "x1": 0.0302, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 114, + "from_node": 7, + "to_node": 8, + "from_status": 1, + "to_status": 1, + "r1": 0.1716, + "x1": 0.0172, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 115, + "from_node": 7, + "to_node": 9, + "from_status": 1, + "to_status": 1, + "r1": 0.1327, + "x1": 0.0133, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 116, + "from_node": 9, + "to_node": 11, + "from_status": 1, + "to_status": 1, + "r1": 0.1026, + "x1": 0.0103, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 117, + "from_node": 10, + "to_node": 11, + "from_status": 1, + "to_status": 1, + "r1": 0.1401, + "x1": 0.0140, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 118, + "from_node": 12, + "to_node": 13, + "from_status": 1, + "to_status": 1, + "r1": 0.1070, + "x1": 0.0107, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 119, + "from_node": 13, + "to_node": 14, + "from_status": 1, + "to_status": 1, + "r1": 0.1843, + "x1": 0.0184, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + }, + { + "id": 120, + "from_node": 9, + "to_node": 14, + "from_status": 1, + "to_status": 1, + "r1": 0.2369, + "x1": 0.0237, + "c1": 0.0, + "tan1": 0.0, + "i_n": 500 + } + ], + "sym_power_sensor": [ + { + "id": 401, + "measured_object": 4, + "measured_terminal_type": 9, + "p_sigma": 1e6, + "q_sigma": 1e6, + "p_measured": 0, + "q_measured": 0 + }, + { + "id": 402, + "measured_object": 5, + "measured_terminal_type": 9, + "p_sigma": 1e6, + "q_sigma": 1e6, + "p_measured": 0, + "q_measured": 0 + }, + { + "id": 403, + "measured_object": 6, + "measured_terminal_type": 9, + "p_sigma": 1e6, + "q_sigma": 1e6, + "p_measured": 0, + "q_measured": 0 + }, + { + "id": 404, + "measured_object": 7, + "measured_terminal_type": 9, + "p_sigma": 1e6, + "q_sigma": 1e6, + "p_measured": 0, + "q_measured": 0 + }, + { + "id": 405, + "measured_object": 10, + "measured_terminal_type": 9, + "p_sigma": 1e6, + "q_sigma": 1e6, + "p_measured": 0, + "q_measured": 0 + }, + { + "id": 406, + "measured_object": 13, + "measured_terminal_type": 9, + "p_sigma": 1e6, + "q_sigma": 1e6, + "p_measured": 0, + "q_measured": 0 + } + ], + "sym_current_sensor": [ + { + "id": 301, + "measured_object": 101, + "measured_terminal_type": 0, + "angle_measurement_type": 0, + "i_sigma": 1.0, + "i_measured": 1.0, + "i_angle_sigma": 1.0, + "i_angle_measured": 1.0 + }, + { + "id": 302, + "measured_object": 105, + "measured_terminal_type": 0, + "angle_measurement_type": 0, + "i_sigma": 1.0, + "i_measured": 1.0, + "i_angle_sigma": 1.0, + "i_angle_measured": 1.0 + }, + { + "id": 303, + "measured_object": 104, + "measured_terminal_type": 0, + "angle_measurement_type": 0, + "i_sigma": 1.0, + "i_measured": 1.0, + "i_angle_sigma": 1.0, + "i_angle_measured": 1.0 + }, + { + "id": 304, + "measured_object": 103, + "measured_terminal_type": 0, + "angle_measurement_type": 0, + "i_sigma": 1.0, + "i_measured": 1.0, + "i_angle_sigma": 1.0, + "i_angle_measured": 1.0 + }, + { + "id": 305, + "measured_object": 111, + "measured_terminal_type": 0, + "angle_measurement_type": 0, + "i_sigma": 1.0, + "i_measured": 1.0, + "i_angle_sigma": 1.0, + "i_angle_measured": 1.0 + }, + { + "id": 306, + "measured_object": 120, + "measured_terminal_type": 0, + "angle_measurement_type": 0, + "i_sigma": 1.0, + "i_measured": 1.0, + "i_angle_sigma": 1.0, + "i_angle_measured": 1.0 + }, + { + "id": 307, + "measured_object": 117, + "measured_terminal_type": 0, + "angle_measurement_type": 0, + "i_sigma": 1.0, + "i_measured": 1.0, + "i_angle_sigma": 1.0, + "i_angle_measured": 1.0 + } + ], + "sym_voltage_sensor": [ + { + "id": 601, + "measured_object": 6, + "u_sigma": 1, + "u_measured": 10000, + "u_angle_measured": 0 + } + ], + "sym_load": [ + { + "id": 501, + "node": 3, + "status": 1, + "type": 0, + "p_specified": 15e6, + "q_specified": 3e6 + }, + { + "id": 502, + "node": 4, + "status": 1, + "type": 0, + "p_specified": 12e6, + "q_specified": 2e6 + }, + { + "id": 503, + "node": 10, + "status": 1, + "type": 0, + "p_specified": 8e6, + "q_specified": 1.5e6 + }, + { + "id": 504, + "node": 1, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 505, + "node": 2, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 506, + "node": 5, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 507, + "node": 6, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 508, + "node": 7, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 509, + "node": 8, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 510, + "node": 9, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 511, + "node": 11, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 512, + "node": 12, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 513, + "node": 13, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + }, + { + "id": 514, + "node": 14, + "status": 1, + "type": 0, + "p_specified": 5e6, + "q_specified": 1e6 + } + ] + } +} diff --git a/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/input.json.license b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/input.json.license new file mode 100644 index 0000000000..7601059167 --- /dev/null +++ b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/input.json.license @@ -0,0 +1,3 @@ +SPDX-FileCopyrightText: Contributors to the Power Grid Model project + +SPDX-License-Identifier: MPL-2.0 diff --git a/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/params.json b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/params.json new file mode 100644 index 0000000000..1b473f2467 --- /dev/null +++ b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/params.json @@ -0,0 +1,8 @@ +{ + "calculation_method": ["newton_raphson", "iterative_linear"], + "rtol": 1e-8, + "atol": { + "default": 1e-8, + ".+_residual": 5e-4 + } +} diff --git a/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/params.json.license b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/params.json.license new file mode 100644 index 0000000000..7601059167 --- /dev/null +++ b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/params.json.license @@ -0,0 +1,3 @@ +SPDX-FileCopyrightText: Contributors to the Power Grid Model project + +SPDX-License-Identifier: MPL-2.0 diff --git a/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/sym_output.json b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/sym_output.json new file mode 100644 index 0000000000..83bbecf194 --- /dev/null +++ b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/sym_output.json @@ -0,0 +1,8 @@ +{ + "version": "1.0", + "type": "sym_output", + "is_batch": false, + "attributes": {}, + "data": { + } +} \ No newline at end of file diff --git a/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/sym_output.json.license b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/sym_output.json.license new file mode 100644 index 0000000000..7601059167 --- /dev/null +++ b/tests/data/state_estimation/meshed-network-observability/08-starting-candidate-dominance/sym_output.json.license @@ -0,0 +1,3 @@ +SPDX-FileCopyrightText: Contributors to the Power Grid Model project + +SPDX-License-Identifier: MPL-2.0 diff --git a/tests/data/state_estimation/meshed-network-observability/test_network_diagrams.svg b/tests/data/state_estimation/meshed-network-observability/test_network_diagrams.svg index 96cbcaf5dd..47b303141f 100644 --- a/tests/data/state_estimation/meshed-network-observability/test_network_diagrams.svg +++ b/tests/data/state_estimation/meshed-network-observability/test_network_diagrams.svg @@ -1,6 +1,6 @@ - - + + \ No newline at end of file