Probabilistic Graph Neural Inference for satellite anomaly response operations for low-power autonomous deployments

Probabilistic Graph Neural Inference for satellite anomaly response operations for low-power autonomous deployments

Introduction: A Noisy Signal from the Edge of Space

It was 3 AM, and I was staring at a stream of telemetry data from a cubesat prototype in my lab. The satellite simulation was running on a Raspberry Pi 4, deliberately power-constrained to mimic orbital conditions. Suddenly, the temperature sensor readings from the propulsion module began oscillating wildly—not enough to trigger traditional threshold alarms, but enough to suggest something was wrong. The onboard rule-based system remained silent. As I manually traced the dependencies between systems, I realized what I was seeing: a cascading anomaly propagating through the satellite's interconnected subsystems. This moment crystallized a fundamental insight from my research: in complex autonomous systems, anomalies aren't isolated events—they're network phenomena.

Through studying distributed satellite constellations and their failure modes, I learned that traditional anomaly detection approaches fail precisely because they treat sensors as independent. In reality, a thermal anomaly in the power system might manifest as communication latency, which then affects attitude control. My exploration of graph neural networks revealed their potential to capture these relationships, but standard GNNs lacked the uncertainty quantification crucial for autonomous decision-making in low-power environments where false positives carry severe consequences.

Technical Background: From Deterministic Graphs to Probabilistic Inference

The Graph Representation Problem

During my investigation of satellite telemetry systems, I found that the most challenging aspect wasn't detecting anomalies—it was representing the system in a way that captured both structural and temporal dependencies. A satellite isn't just a collection of sensors; it's a dynamic network where:

1. Physical connections (power buses, data buses, thermal paths) create hard dependencies
2. Functional dependencies (attitude control needing power, communication needing thermal stability) create soft constraints
3. Temporal patterns (orbital periods, thermal cycles, communication windows) create time-varying relationships

While exploring probabilistic graphical models, I discovered that traditional Bayesian networks struggled with the high-dimensional, time-series nature of telemetry data. Markov logic networks offered more flexibility but became computationally intractable for real-time inference on edge hardware.

Enter Probabilistic Graph Neural Networks

My research into modern graph learning revealed that Probabilistic Graph Neural Networks (PGNNs) combine the representational power of GNNs with the uncertainty quantification of probabilistic models. The key insight I gained from studying papers like "Uncertainty in Graph Neural Networks" (Zhu et al., 2021) was that we could model two types of uncertainty:

1. Aleatoric uncertainty: Inherent noise in the observations (sensor noise, cosmic radiation effects)
2. Epistemic uncertainty: Model uncertainty due to limited training data (especially important for rare anomalies)

One interesting finding from my experimentation with different uncertainty quantification methods was that Monte Carlo Dropout, while computationally efficient, often underestimated uncertainty in out-of-distribution scenarios common in space operations. This led me to explore Bayesian neural network approaches adapted for graph structures.

Implementation Details: Building a PGNN for Satellite Operations

Graph Construction from Telemetry Data

In my implementation, I approached graph construction as a multi-modal problem. Each subsystem becomes a node with features including:

• Current sensor readings (normalized)
• Historical statistics (mean, variance over sliding windows)
• Operational mode flags
• Time since last maintenance event

Edges are constructed based on:

• Physical connectivity matrices (from system diagrams)
• Correlation patterns learned from historical data
• Known functional dependencies

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch_geometric.nn import GCNConv, MessagePassing
from torch_geometric.data import Data
import numpy as np

class SatelliteGraphBuilder:
    """Constructs dynamic graph from satellite telemetry"""

    def __init__(self, n_subsystems, history_window=100):
        self.n_subsystems = n_subsystems
        self.history_window = history_window
        self.physical_adjacency = self.load_physical_connectivity()

    def telemetry_to_graph(self, telemetry_batch):
        """Convert time-series telemetry to graph representation"""
        # Node features: [current_value, mean_10s, std_10s, trend]
        node_features = []

        for subsystem_id in range(self.n_subsystems):
            current = telemetry_batch[subsystem_id, -1]
            history = telemetry_batch[subsystem_id, -self.history_window:]

            features = [
                current,
                np.mean(history),
                np.std(history),
                self.calculate_trend(history)
            ]
            node_features.append(features)

        # Dynamic edge weights based on correlation
        edge_index, edge_attr = self.compute_dynamic_edges(telemetry_batch)

        return Data(
            x=torch.tensor(node_features, dtype=torch.float32),
            edge_index=edge_index,
            edge_attr=edge_attr
        )

    def compute_dynamic_edges(self, telemetry_batch):
        """Compute edges based on recent correlations"""
        # Implementation of dynamic correlation-based edges
        correlations = np.corrcoef(telemetry_batch)

        # Threshold and create edge list
        edge_index = []
        edge_weights = []

        for i in range(self.n_subsystems):
            for j in range(i + 1, self.n_subsystems):
                if abs(correlations[i, j]) > 0.7:  # Strong correlation threshold
                    edge_index.append([i, j])
                    edge_index.append([j, i])  # Undirected graph

                    # Combine physical and correlation information
                    weight = (self.physical_adjacency[i, j] * 0.3 +
                             abs(correlations[i, j]) * 0.7)
                    edge_weights.extend([weight, weight])

        return (torch.tensor(edge_index, dtype=torch.long).t().contiguous(),
                torch.tensor(edge_weights, dtype=torch.float32))

Probabilistic Graph Neural Network Architecture

Through experimenting with different architectures, I developed a PGNN that balances expressiveness with computational efficiency—critical for low-power deployments. The key innovation was separating the deterministic feature extraction from probabilistic uncertainty estimation.

class ProbabilisticGNNLayer(MessagePassing):
    """Single layer of probabilistic graph neural network"""

    def __init__(self, in_channels, out_channels, dropout_rate=0.1):
        super().__init__(aggr='mean')

        # Deterministic transformation
        self.deterministic_lin = nn.Linear(in_channels, out_channels)

        # Probabilistic components
        self.mu_lin = nn.Linear(in_channels, out_channels)
        self.log_var_lin = nn.Linear(in_channels, out_channels)

        self.dropout = nn.Dropout(dropout_rate)
        self.activation = nn.ReLU()

    def forward(self, x, edge_index, edge_weight=None):
        # Message passing
        aggregated = self.propagate(edge_index, x=x, edge_weight=edge_weight)

        # Deterministic features
        deterministic = self.activation(self.deterministic_lin(aggregated))
        deterministic = self.dropout(deterministic)

        # Probabilistic parameters
        mu = self.mu_lin(aggregated)
        log_var = self.log_var_lin(aggregated)

        # Reparameterization trick for training
        if self.training:
            std = torch.exp(0.5 * log_var)
            eps = torch.randn_like(std)
            probabilistic = mu + eps * std
        else:
            probabilistic = mu  # Use mean during inference

        return deterministic, mu, log_var, probabilistic

class SatellitePGNN(nn.Module):
    """Complete PGNN for satellite anomaly detection"""

    def __init__(self, input_dim, hidden_dim, output_dim, n_layers=3):
        super().__init__()

        self.input_proj = nn.Linear(input_dim, hidden_dim)

        # Stack of PGNN layers
        self.layers = nn.ModuleList([
            ProbabilisticGNNLayer(
                hidden_dim if i == 0 else hidden_dim * 2,
                hidden_dim
            ) for i in range(n_layers)
        ])

        # Output heads
        self.anomaly_head = nn.Linear(hidden_dim * 2, output_dim)
        self.uncertainty_head = nn.Linear(hidden_dim * 2, output_dim)

        # For uncertainty calibration
        self.temperature = nn.Parameter(torch.ones(1))

    def forward(self, data, n_samples=5):
        x, edge_index, edge_weight = data.x, data.edge_index, data.edge_attr

        x = F.relu(self.input_proj(x))

        deterministic_features = []
        all_mus = []
        all_log_vars = []

        # Forward through layers
        for layer in self.layers:
            det, mu, log_var, prob = layer(x, edge_index, edge_weight)
            deterministic_features.append(det)
            all_mus.append(mu)
            all_log_vars.append(log_var)
            x = torch.cat([det, prob], dim=-1)

        # Monte Carlo sampling for uncertainty estimation
        anomaly_scores = []
        uncertainties = []

        for _ in range(n_samples):
            # Sample from final layer distribution
            final_mu = all_mus[-1]
            final_log_var = all_log_vars[-1]

            if self.training or n_samples > 1:
                std = torch.exp(0.5 * final_log_var)
                eps = torch.randn_like(std)
                sampled = final_mu + eps * std
            else:
                sampled = final_mu

            # Anomaly score
            anomaly_score = torch.sigmoid(self.anomaly_head(sampled))
            anomaly_scores.append(anomaly_score)

            # Uncertainty estimation
            uncertainty = torch.sigmoid(self.uncertainty_head(sampled))
            uncertainties.append(uncertainty)

        # Aggregate samples
        anomaly_scores = torch.stack(anomaly_scores)
        uncertainties = torch.stack(uncertainties)

        mean_anomaly = anomaly_scores.mean(dim=0)
        mean_uncertainty = uncertainties.mean(dim=0)
        score_variance = anomaly_scores.var(dim=0)

        return {
            'anomaly_score': mean_anomaly,
            'uncertainty': mean_uncertainty,
            'score_variance': score_variance,
            'all_scores': anomaly_scores
        }

Low-Power Optimization Techniques

During my experimentation with edge deployment, I discovered several critical optimizations for low-power environments:

class LowPowerPGNN(SatellitePGNN):
    """Optimized version for low-power deployment"""

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

        # Quantization-aware training setup
        self.quant = torch.quantization.QuantStub()
        self.dequant = torch.quantization.DeQuantStub()

        # Adaptive computation - skip layers when confident
        self.confidence_threshold = 0.9

    def forward(self, data, max_layers=None):
        # Quantize input
        x = self.quant(data.x)

        # Adaptive depth computation
        if max_layers is None:
            max_layers = len(self.layers)

        early_exit_scores = []
        layer_uncertainties = []

        x = F.relu(self.input_proj(x))

        for i, layer in enumerate(self.layers[:max_layers]):
            det, mu, log_var, prob = layer(x, data.edge_index, data.edge_attr)

            # Early exit check
            if i > 0:  # Check after at least one layer
                current_score = torch.sigmoid(
                    self.anomaly_head(torch.cat([det, prob], dim=-1))
                )
                current_uncertainty = torch.sigmoid(
                    self.uncertainty_head(torch.cat([det, prob], dim=-1))
                )

                # Check if we can exit early
                high_confidence = (current_uncertainty < 0.1).all()
                if high_confidence and i >= 1:
                    # Use this layer's output
                    x = torch.cat([det, prob], dim=-1)
                    break

            x = torch.cat([det, prob], dim=-1)

        # Dequantize for output
        x = self.dequant(x)

        # Final computation with reduced Monte Carlo samples
        result = super().forward(
            Data(x=x, edge_index=data.edge_index, edge_attr=data.edge_attr),
            n_samples=3  # Reduced for power savings
        )

        return result

# Power-aware inference scheduler
class PowerAwareInferenceScheduler:
    """Dynamically adjusts model complexity based on available power"""

    def __init__(self, battery_capacity, power_budget):
        self.battery_capacity = battery_capacity
        self.power_budget = power_budget
        self.complexity_levels = {
            'critical': {'n_samples': 10, 'n_layers': 3},
            'high': {'n_samples': 5, 'n_layers': 3},
            'medium': {'n_samples': 3, 'n_layers': 2},
            'low': {'n_samples': 1, 'n_layers': 1}
        }

    def select_complexity(self, battery_level, anomaly_confidence):
        """Select appropriate model complexity"""
        power_ratio = battery_level / self.battery_capacity

        if power_ratio < 0.2:
            return 'low'
        elif anomaly_confidence < 0.7:
            return 'critical'  # Need more certainty
        elif power_ratio < 0.5:
            return 'medium'
        else:
            return 'high'

Real-World Applications: Autonomous Anomaly Response

Decision Making Under Uncertainty

One of the most valuable insights from my research was that anomaly detection is only half the battle. The real challenge is deciding what to do when you detect something unusual. Through studying autonomous satellite operations, I developed a probabilistic decision framework:

class AnomalyResponsePlanner:
    """Plans responses based on anomaly probabilities and uncertainties"""

    def __init__(self, action_costs, action_effectiveness):
        self.action_costs = action_costs  # Power, time, risk
        self.action_effectiveness = action_effectiveness

    def plan_response(self, anomaly_scores, uncertainties, system_state):
        """Select optimal response using expected utility"""

        possible_actions = self.generate_actions(anomaly_scores, system_state)

        best_action = None
        best_expected_utility = -float('inf')

        for action in possible_actions:
            # Calculate expected utility considering uncertainty
            expected_utility = self.calculate_expected_utility(
                action, anomaly_scores, uncertainties, system_state
            )

            if expected_utility > best_expected_utility:
                best_expected_utility = expected_utility
                best_action = action

        return best_action, best_expected_utility

    def calculate_expected_utility(self, action, anomaly_scores, uncertainties, state):
        """Monte Carlo estimation of expected utility"""

        n_samples = 100
        total_utility = 0

        for _ in range(n_samples):
            # Sample from anomaly distribution
            sampled_anomalies = []
            for score, uncertainty in zip(anomaly_scores, uncertainties):
                # Sample from beta distribution (bounded [0,1])
                alpha = score * (1 - uncertainty) * 10 + 1
                beta = (1 - score) * (1 - uncertainty) * 10 + 1
                sampled = np.random.beta(alpha.item(), beta.item())
                sampled_anomalies.append(sampled)

            # Simulate action outcome
            outcome = self.simulate_action(action, sampled_anomalies, state)
            utility = self.utility_function(outcome, action)
            total_utility += utility

        return total_utility / n_samples

# Example action templates
RESPONSE_ACTIONS = {
    'monitor': {
        'cost': {'power': 0.1, 'time': 60},
        'effect': 'increases_observation'
    },
    'reset_subsystem': {
        'cost': {'power': 0.5, 'time': 5, 'risk': 0.1},
        'effect': 'clears_software_state'
    },
    'safe_mode': {
        'cost': {'power': 0.3, 'time': 30, 'risk': 0.05},
        'effect': 'minimal_operations'
    },
    'reconfigure': {
        'cost': {'power': 0.8, 'time': 120, 'risk': 0.2},
        'effect': 'hardware_reconfiguration'
    }
}

Distributed Constellation Coordination

While exploring multi-satellite systems, I realized that anomalies often propagate across constellations. My research into distributed graph learning led to a federated approach:

python
class FederatedPGNNTrainer:
    """Federated learning for constellation-wide anomaly detection"""

    def __init__(self, constellation_size, communication_budget):
        self.constellation_size = constellation_size
        self.communication_budget = communication_budget
        self.global_model = None
        self.client_models = []

    def federated_round(self, client_data, rounds=10):
        """Execute one round of federated learning"""

        # 1. Send global model to clients (within communication budget)
        compressed_model = self.compress_model(self.global_model)

        # 2. Clients train locally on their anomaly data
        client_updates = []
        for client_id, data