deepspace/sim.py

"""Spacecraft simulator for oracle7-deepspace. Deterministic (seeded), stdlib only, one tick = 60 s.

Four mission PROFILES modelled on real deep-space craft (Voyager 1, New Horizons, Juno, Europa Clipper). A profile sets
the power source (RTG with a small capacitor buffer and no battery, or solar array + battery), the attitude actuators
(reaction wheels, or thrusters on A/B branches), the thermal environment (colder with distance), the power budget and
the one-way light time. The physics are deliberately simple and documented; the profiles set which subsystems exist
and which faults are possible, not a claim of engineering-level fidelity to the real vehicles.

Faults are injected into the TRUE state; the onboard software only ever sees TELEMETRY. Two faults are TRAPS where
acting would be wrong: a biased sensor (the zone is fine) and a single-tick radiation upset. The same class, built from
telemetry with perturbed parameters, is the onboard DIGITAL TWIN.
"""
from __future__ import annotations

import copy
import math
import random

DT = 60.0
ZONES = ("battery", "avionics", "tank")
BANDS = {"battery": (5.0, 30.0), "avionics": (0.0, 40.0), "tank": (10.0, 35.0)}        # operating limits (deg C)
SURVIVAL = {"battery": (-10.0, 45.0), "avionics": (-20.0, 60.0), "tank": (2.0, 50.0)}  # beyond these = damage
THERMOSTAT = {"battery": (12.0, 22.0), "avionics": (8.0, 25.0), "tank": (16.0, 26.0)}   # heater on below lo, off above hi
HEATER_ON_EQ = {"battery": 49.0, "avionics": 20.0, "tank": 40.0}                        # where a stuck-on heater settles

# one-way light time in minutes is approximate and for narrative only (distances change; see the site's ephemeris)
PROFILES = {
    "voyager1": {"label": "Voyager 1", "power": "rtg", "actuator": "thrusters", "distance_au": 167.0, "t_env_shift": -30.0,
                 "source_amp": 7.2, "base_load": 3.9, "science_amp": 1.2, "heater_amp": 0.8, "cap_ah": 2.0,
                 "faults": ("nominal", "heater_stuck_on", "heater_stuck_off", "sensor_bias", "seu_spike",
                            "thruster_degradation", "rtg_decline", "star_tracker_dropout")},
    "new_horizons": {"label": "New Horizons", "power": "rtg", "actuator": "thrusters", "distance_au": 61.0, "t_env_shift": -25.0,
                     "source_amp": 6.6, "base_load": 3.5, "science_amp": 1.1, "heater_amp": 0.8, "cap_ah": 2.0,
                     "faults": ("nominal", "heater_stuck_on", "heater_stuck_off", "sensor_bias", "seu_spike",
                                "thruster_degradation", "rtg_decline", "star_tracker_dropout")},
    "juno": {"label": "Juno", "power": "solar", "actuator": "thrusters", "distance_au": 5.2, "t_env_shift": -10.0,
             "source_amp": 14.0, "base_load": 6.0, "science_amp": 2.5, "heater_amp": 1.5, "cap_ah": 110.0,
             "faults": ("nominal", "heater_stuck_on", "heater_stuck_off", "sensor_bias", "seu_spike",
                        "thruster_degradation", "solar_degradation", "battery_short", "star_tracker_dropout")},
    "europa_clipper": {"label": "Europa Clipper", "power": "solar", "actuator": "wheels", "distance_au": 3.0, "t_env_shift": -4.0,
                       "source_amp": 12.0, "base_load": 5.0, "science_amp": 2.0, "heater_amp": 1.5, "cap_ah": 100.0,
                       "faults": ("nominal", "heater_stuck_on", "heater_stuck_off", "wheel_friction", "sensor_bias",
                                  "star_tracker_dropout", "solar_degradation", "battery_short", "seu_spike")},
}
FAULTS = ("nominal", "heater_stuck_on", "heater_stuck_off", "wheel_friction", "sensor_bias", "star_tracker_dropout",
          "solar_degradation", "battery_short", "seu_spike", "thruster_degradation", "rtg_decline")
# the response that resolves each fault (ground truth for SCORING only; the onboard software never sees it)
CORRECT = {"nominal": {"no_action"}, "heater_stuck_on": {"switch_heater_backup"}, "heater_stuck_off": {"switch_heater_backup"},
           "wheel_friction": {"switch_to_spare_wheel"}, "sensor_bias": {"mark_sensor_invalid"},
           "star_tracker_dropout": {"switch_attitude_sensor"}, "solar_degradation": {"shed_loads"},
           "battery_short": {"safe_mode"}, "seu_spike": {"no_action"}, "thruster_degradation": {"switch_thruster_branch"},
           "rtg_decline": {"shed_loads"}}
HARMFUL = {"sensor_bias": {"switch_heater_backup", "safe_mode", "shed_loads"},
           "seu_spike": {"switch_heater_backup", "safe_mode", "shed_loads", "mark_sensor_invalid", "switch_to_spare_wheel",
                         "switch_thruster_branch"},
           "nominal": {"switch_heater_backup", "safe_mode", "shed_loads", "switch_to_spare_wheel", "switch_attitude_sensor",
                       "switch_thruster_branch"}}


def light_time_min(profile):
    return PROFILES[profile]["distance_au"] * 8.317


class Spacecraft:
    def __init__(self, seed=0, params=None, profile="europa_clipper"):
        self.profile = profile
        pr = PROFILES[profile]
        self.rng = random.Random(seed)
        k_env = {"battery": 0.0022, "avionics": 0.0030, "tank": 0.0018}
        t_env = {z: v + pr["t_env_shift"] for z, v in {"battery": -5.0, "avionics": 5.0, "tank": -15.0}.items()}
        avionics_self = 0.020 + 0.0030 * max(0.0, -pr["t_env_shift"])   # colder craft keep avionics warm with internal heat
        heater_w = {z: max(0.0, k_env[z] * (HEATER_ON_EQ[z] - t_env[z]) - (avionics_self if z == "avionics" else 0.0))
                    for z in ZONES}
        p = {"cap_ah": pr["cap_ah"], "source_amp": pr["source_amp"], "base_load": pr["base_load"],
             "heater_amp": pr["heater_amp"], "science_amp": pr["science_amp"], "wheel_amp": 0.3, "k_env": k_env,
             "t_env": t_env, "heater_w": heater_w, "avionics_self_heat": avionics_self, "noise": 0.25}
        if params:
            p.update(params)
        self.p = p
        self.power, self.actuator = pr["power"], pr["actuator"]
        self.t = 0
        self.soc = 0.85 if self.power == "solar" else 0.9
        self.temp = {"battery": 18.0, "avionics": 18.0, "tank": 21.0}
        self.heater_primary_ok = {z: "ok" for z in ZONES}
        self.heater_circuit = {z: "primary" for z in ZONES}
        self.heater_cmd = {z: False for z in ZONES}
        self.sensor_bias = {(z, s): 0.0 for z in ZONES for s in "AB"}
        self.sensor_invalid = set()
        self.wheel_friction = [1.0, 1.0, 1.0, 1.0]
        self.wheels_active = [0, 1, 2] if self.actuator == "wheels" else []
        self.thruster_branch = "A"
        self.thruster_eff = {"A": 1.0, "B": 1.0}
        self.star_tracker_ok = True
        self.attitude_sensor = "star_tracker"
        self.source_health = 1.0
        self.short_amp = 0.0
        self.science_on = True
        self.safe = False
        self.pointing_err = 0.2
        self.seu = None
        self.fault = "nominal"
        self.damage = []
        self.i_source = self.i_load = 0.0
        self.duty = 0.05

    # ---- fault injection (TRUE state) ----
    def inject(self, fault, zone="battery", wheel=1, sensor="A", magnitude=1.0):
        self.fault = fault
        if fault == "heater_stuck_on":
            self.heater_primary_ok[zone] = "stuck_on"
        elif fault == "heater_stuck_off":
            self.heater_primary_ok[zone] = "stuck_off"
        elif fault == "wheel_friction":
            self.wheel_friction[wheel] = 6.0 + magnitude
        elif fault == "sensor_bias":
            self.sensor_bias[(zone, sensor)] = 18.0 + 10.0 * magnitude
        elif fault == "star_tracker_dropout":
            self.star_tracker_ok = False
        elif fault == "solar_degradation":
            self.source_health = 0.78 - 0.08 * magnitude
        elif fault == "rtg_decline":
            self.source_health = 0.86 - 0.04 * magnitude            # a thermocouple string lost: ~0.9-1.3 A on a tight budget
        elif fault == "battery_short":
            self.short_amp = 5.0 + magnitude
        elif fault == "thruster_degradation":
            self.thruster_eff["A"] = 0.18 - 0.08 * magnitude         # clogged branch-A attitude thrusters
        elif fault == "seu_spike":
            self.seu = (self.t + 1, zone, sensor)

    # ---- commands ----
    def command(self, action, **kw):
        if action == "switch_heater_backup":
            self.heater_circuit[kw["zone"]] = "backup"
        elif action == "mark_sensor_invalid":
            self.sensor_invalid.add((kw["zone"], kw["sensor"]))
        elif action == "switch_to_spare_wheel":
            w = kw["wheel"]
            if w in self.wheels_active and 3 not in self.wheels_active:
                self.wheels_active = [x if x != w else 3 for x in self.wheels_active]
        elif action == "switch_thruster_branch":
            self.thruster_branch = "B" if self.thruster_branch == "A" else "A"
        elif action == "switch_attitude_sensor":
            self.attitude_sensor = "sun_sensor"
        elif action == "shed_loads":
            self.science_on = False
        elif action == "safe_mode":
            self.safe = True
            self.science_on = False

    # ---- physics ----
    def _heater_on(self, z):
        if self.heater_circuit[z] == "primary":
            st = self.heater_primary_ok[z]
            if st == "stuck_on":
                return True
            if st == "stuck_off":
                return False
        return self.heater_cmd[z]

    def step(self):
        p = self.p
        self.t += 1
        for z in ZONES:
            est = self.zone_estimate(z)
            lo, hi = THERMOSTAT[z]
            if self.safe:
                lo, hi = lo - 4, hi - 4
            if est < lo:
                self.heater_cmd[z] = True
            elif est > hi:
                self.heater_cmd[z] = False
        for z in ZONES:
            heat = p["heater_w"][z] if self._heater_on(z) else 0.0
            if z == "avionics":
                heat += p["avionics_self_heat"]
            self.temp[z] += DT * (heat - p["k_env"][z] * (self.temp[z] - p["t_env"][z])) * 0.05
        # attitude
        if self.actuator == "wheels":
            authority = sum(1.0 / self.wheel_friction[w] for w in self.wheels_active) / 3.0
            self.duty = 0.05
        else:
            authority = self.thruster_eff[self.thruster_branch]
            self.duty = min(1.0, 0.05 / max(0.05, authority))       # weak thrusters fire longer for the same control
        knowledge = 1.0 if (self.attitude_sensor == "star_tracker" and self.star_tracker_ok) else (
            0.6 if self.attitude_sensor == "sun_sensor" else 0.0)
        target = 0.2 if knowledge >= 1.0 else (1.2 if knowledge > 0 else 6.0)
        if authority < 0.8:
            target += (0.8 - authority) * 25.0
        if self.safe:
            target = min(target, 2.0)
        self.pointing_err += (target - self.pointing_err) * 0.3 + self.rng.gauss(0, 0.05)
        self.pointing_err = max(0.05, self.pointing_err)
        # power
        if self.power == "solar":
            sun = math.cos(math.radians(min(self.pointing_err, 80.0)))
            src = p["source_amp"] * self.source_health * sun * (1.05 if self.safe else 1.0)
        else:
            src = p["source_amp"] * self.source_health * (1.0 - 1e-6 * self.t)   # RTG: pointing-independent, slow decay
        load = p["base_load"] + sum(p["heater_amp"] for z in ZONES if self._heater_on(z)) \
            + (p["science_amp"] if self.science_on else 0.0) + self.short_amp
        if self.actuator == "wheels":
            load += p["wheel_amp"] * sum(self.wheel_friction[w] for w in self.wheels_active)
        else:
            load += 0.6 * self.duty
        if self.safe:
            load -= 0.5 + 0.6 * self.short_amp + (0.3 if self.power == "rtg" else 0.5)
        self.i_source, self.i_load = src, load
        self.soc = min(1.0, max(0.0, self.soc + (src - load) * DT / 3600.0 / p["cap_ah"]))
        for z in ZONES:
            lo, hi = SURVIVAL[z]
            if not lo <= self.temp[z] <= hi:
                self.damage.append((self.t, "%s_temp" % z))
        if self.soc < 0.15:
            self.damage.append((self.t, "bus_undervoltage" if self.power == "rtg" else "battery_depleted"))

    def zone_estimate(self, z):
        vals = [self._read(z, s) for s in "AB" if (z, s) not in self.sensor_invalid]
        return sum(vals) / len(vals) if vals else self.temp[z]

    def _read(self, z, s):
        v = self.temp[z] + self.sensor_bias[(z, s)] + self.rng.gauss(0, self.p["noise"])
        if self.seu and self.seu[0] == self.t and self.seu[1] == z and self.seu[2] == s:
            v += 90.0
        return v

    def nominal_source(self):
        """what the ONBOARD model expects the power source to deliver now (known design value, no fault knowledge)"""
        if self.power == "solar":
            return self.p["source_amp"] * math.cos(math.radians(min(self.pointing_err, 80.0)))
        return self.p["source_amp"] * (1.0 - 1e-6 * self.t)

    def telemetry(self):
        tm = {"t": self.t, "profile": self.profile, "power_source": self.power, "actuator": self.actuator,
              "soc": round(self.soc, 4),
              "bus_v": round(26.0 + 6.0 * self.soc - 0.25 * self.short_amp + self.rng.gauss(0, 0.05), 3),
              "i_source": round(self.i_source, 3), "i_source_nominal": round(self.nominal_source(), 3),
              "i_load": round(self.i_load, 3), "pointing_err": round(self.pointing_err, 3),
              "star_tracker_valid": self.star_tracker_ok, "attitude_sensor": self.attitude_sensor,
              "wheels_active": list(self.wheels_active), "thruster_branch": self.thruster_branch,
              "thruster_duty": round(self.duty + self.rng.gauss(0, 0.004), 4),
              "science_on": self.science_on, "safe_mode": self.safe,
              "sensor_invalid": sorted("%s.%s" % zs for zs in self.sensor_invalid)}
        tm["i_solar"] = tm["i_source"]                                  # compatibility name used by v0 code
        for z in ZONES:
            tm["%s_A" % z] = round(self._read(z, "A"), 2)
            tm["%s_B" % z] = round(self._read(z, "B"), 2)
            tm["%s_heater_cmd" % z] = self.heater_cmd[z]
            tm["%s_heater_circuit" % z] = self.heater_circuit[z]
        for w in range(4):
            tm["wheel%d_current" % w] = round(0.3 * self.wheel_friction[w] * (1 if w in self.wheels_active else 0.1)
                                              + self.rng.gauss(0, 0.01), 3)
        return tm

    def clone(self, rng_seed=None, perturb=0.0):
        """digital twin: a copy with physics parameters perturbed by +-perturb (thermal and source constants)"""
        c = copy.deepcopy(self)
        c.rng = random.Random(rng_seed if rng_seed is not None else 12345)
        if perturb:
            r = random.Random(777)
            for key in ("k_env", "heater_w"):
                c.p[key] = {z: v * (1 + r.uniform(-perturb, perturb)) for z, v in c.p[key].items()}
            c.p["source_amp"] *= 1 + r.uniform(-perturb, perturb)
        return c