Thermal model

README.md

@@ -25,6 +25,7 @@ Disclaimer: This project is currently under development. Use at your own risk.
     <li><a href="#set-an-orbit-for-a-paseos-spacecraftactor">Set an orbit for a PASEOS SpacecraftActor</a></li>
     <li><a href="#how-to-add-a-communication-device">How to add a communication device</a></li>
     <li><a href="#how-to-add-a-power-device">How to add a power device</a></li>
+    <li><a href="#thermal-modeling">Thermal Modeling</a></li>
     </ul>
     <li><a href="#simulation-settings">Simulation Settings</a></li>
     <ul>
@@ -231,10 +232,55 @@ ActorBuilder.set_power_devices(actor=sat_actor,
                                # Charging rate in W
                                charging_rate_in_W=10)
 ```
+
+#### Thermal Modelling
+To model thermal constraints on spacecraft we utilize a model inspired by the one-node model described in [Martínez - Spacecraft Thermal Modelling and Test](http://imartinez.etsiae.upm.es/~isidoro/tc3/Spacecraft%20Thermal%20Modelling%20and%20Testing.pdf). Thus, we model the change in temperature as 
+
+$mc \, \frac{dT}{dt} = \dot{Q}_{solar} + \dot{Q}_{albedo} + \dot{Q}_{central_body_IR} - \dot{Q}_{dissipated} + \dot{Q}_{activity}.$
+
+This means your spacecraft will heat up due to being in sunlight, albedo reflections, infrared radiation emitted by the central body as well as due to power consumption of activities. It will cool down due to heat dissipation.
+
+The model is only available for a [SpacecraftActor](#spacecraftactor) and (like all the physical models) only evaluated for the [local actor](#local-actor).
+
+The following parameters have to be specified for this:
+* Spacecraft mass [kg], initial temperature [K], emissive area (for heat disspiation) and thermal capacity [J / (kg * K)]
+* Spacecraft absorptance of Sun light, infrared light. [0 to 1]
+* Spacecraft area [m^2] facing Sun and central body, respectively
+* Solar irradiance in this orbit [W] (defaults to 1360W)
+* Central body surface temperature [k] (defaults to 288K)
+* Central body emissivity and reflectance [0 to 1] (defaults to 0.6 and 0.3)
+* Ratio of power converted to heat (defaults to 0.5)
+
+To use it, simply equip your [SpacecraftActor](#spacecraftactor) with a thermal model with:
+
+```py
+from paseos import SpacecraftActor, ActorBuilder
+my_actor = ActorBuilder.get_actor_scaffold("my_actor", SpacecraftActor, pk.epoch(0))
+ActorBuilder.set_thermal_model(
+    actor=my_actor,
+    actor_mass=50.0, # Setting mass to 50kg
+    actor_initial_temperature_in_K=273.15, # Setting initialtemperature to 0°C
+    actor_sun_absorptance=1.0, # Depending on material, define absorptance
+    actor_infrared_absorptance=1.0, # Depending on material, define absorptance
+    actor_sun_facing_area=1.0, # Area in m2
+    actor_central_body_facing_area=1.0, # Area in m2
+    actor_emissive_area=1.0, # Area in m2
+    actor_thermal_capacity=1000, # Capacity in J / (kg * K)
+    # ... leaving out default valued parameters, see docs for details
+)
+```
+
+The model is evaluated automatically during [activities](#activity). You can check the spacecraft temperature with:
+
+```py
+print(my_actor.temperature_in_K)
+```
+
 ### Simulation Settings
 #### Initializing PASEOS
 We will now show how to create an instance of PASEOS. An instance of PASEOS shall be bounded to one PASEOS [actor](#actor) that we call [local actor](#local-actor). Please, notice that an orbit shall be placed for a [SpacecraftActor](#spacecraftactor) before being added to a PASEOS instance. <br>
 
+### How to instantiate PASEOS
 ```py 
 import pykep as pk
 import paseos
@@ -680,4 +726,4 @@ Created by $\Phi$[-lab@Sweden](https://www.ai.se/en/data-factory/f-lab-sweden).
 * Pablo Gómez - pablo.gomez at esa.int, pablo.gomez at ai.se
 * Gabriele Meoni - gabriele.meoni at esa.int, gabriele.meoni at ai.se
 * Johan Östman - johan.ostman at ai.se
-* Vinutha Magal Shreenath - vinutha at ai.se
+* Vinutha Magal Shreenath - vinutha at ai.se

paseos/activities/activity_processor.py

@@ -107,6 +107,16 @@ async def _update(self, elapsed_time: float):
         if self._advance_paseos_clock:
             self._paseos_instance.advance_time(elapsed_time)
 
+        # Update actor temperature
+        if (
+            hasattr(self._paseos_instance.local_actor, "_thermal_model")
+            and self._paseos_instance.local_actor._thermal_model is not None
+        ):
+            self._paseos_instance.local_actor._thermal_model.update_temperature(
+                elapsed_time, self._power_consumption_in_watt
+            )
+
+        # Update state of charge
         self._paseos_instance.local_actor.discharge(
             self._power_consumption_in_watt, elapsed_time
         )

paseos/actors/actor_builder.py

@@ -6,6 +6,7 @@
 from .base_actor import BaseActor
 from .spacecraft_actor import SpacecraftActor
 from .ground_station_actor import GroundstationActor
+from ..thermal.thermal_model import ThermalModel
 
 
 class ActorBuilder:
@@ -165,6 +166,102 @@ def set_power_devices(
             + f"ChargingRate={charging_rate_in_W}W to actor {actor}"
         )
 
+    def set_thermal_model(
+        actor: SpacecraftActor,
+        actor_mass: float,
+        actor_initial_temperature_in_K: float,
+        actor_sun_absorptance: float,
+        actor_infrared_absorptance: float,
+        actor_sun_facing_area: float,
+        actor_central_body_facing_area: float,
+        actor_emissive_area: float,
+        actor_thermal_capacity: float,
+        body_solar_irradiance: float = 1360,
+        body_surface_temperature_in_K: float = 288,
+        body_emissivity: float = 0.6,
+        body_reflectance: float = 0.3,
+        power_consumption_to_heat_ratio: float = 0.5,
+    ):
+        """Add a thermal model to the actor to model temperature based on
+        heat flux from sun, central body albedo, central body IR, actor IR
+        emission and due to actor activities.
+        For the moment, it is a slightly simplified version
+        of the single node model from "Spacecraft Thermal Control" by Prof. Isidoro Martínez
+        available at http://imartinez.etsiae.upm.es/~isidoro/tc3/Spacecraft%20Thermal%20Modelling%20and%20Testing.pdf
+
+        Args:
+            actor (SpacecraftActor): Actor to model.
+            actor_mass (float): Actor's mass in kg.
+            actor_initial_temperature_in_K (float): Actor's initial temperature in K.
+            actor_sun_absorptance (float): Actor's absorptance ([0,1]) of solar light
+            actor_infrared_absorptance (float): Actor's absportance ([0,1]) of IR.
+            actor_sun_facing_area (float): Actor area facing the sun in m^2.
+            actor_central_body_facing_area (float): Actor area facing central body in m^2.
+            actor_emissive_area (float): Actor area emitting (radiating) heat.
+            actor_thermal_capacity (float): Actor's thermal capacity in J / (kg * K).
+            body_solar_irradiance (float, optional): Irradiance from the sun in W. Defaults to 1360.
+            body_surface_temperature_in_K (float, optional): Central body surface temperature. Defaults to 288.
+            body_emissivity (float, optional): Centrla body emissivity [0,1] in IR. Defaults to 0.6.
+            body_reflectance (float, optional): Central body reflectance of sun light. Defaults to 0.3.
+            power_consumption_to_heat_ratio (float, optional): Conversion ratio for activities.
+            0 leads to know heat-up due to activity. Defaults to 0.5.
+        """
+        # check for spacecraft actor
+        assert isinstance(
+            actor, SpacecraftActor
+        ), "Thermal models are only supported for SpacecraftActors"
+
+        assert actor_mass > 0, "Actor mass has to be positive."
+
+        assert (
+            0 <= power_consumption_to_heat_ratio
+            and power_consumption_to_heat_ratio <= 1.0
+        ), "Heat ratio has to be 0 to 1."
+
+        logger.trace("Checking actor thermal values for sensibility.")
+        assert (
+            0 <= actor_initial_temperature_in_K
+        ), "Actor initial temperature cannot be below 0K."
+        assert (
+            0 <= actor_sun_absorptance and actor_sun_absorptance <= 1.0
+        ), "Absorptance has to be 0 to 1."
+        assert (
+            0 <= actor_infrared_absorptance and actor_infrared_absorptance <= 1.0
+        ), "Absorptance has to be 0 to 1."
+        assert 0 < actor_sun_facing_area, "Sun-facing area has to be > 0."
+        assert 0 < actor_central_body_facing_area, "Body-facing area has to be > 0."
+        assert 0 < actor_emissive_area, "Actor emissive area has to be > 0."
+        assert 0 < actor_thermal_capacity, "Thermal capacity has to be > 0"
+
+        logger.trace("Checking body thermal values for sensibility.")
+        assert 0 < body_solar_irradiance, "Solar irradiance has to be > 0."
+        assert (
+            0 <= body_surface_temperature_in_K
+        ), "Body surface temperature cannot be below 0K."
+        assert (
+            0 <= body_emissivity and body_emissivity <= 1.0
+        ), "Body emissivity has to be 0 to 1"
+        assert (
+            0 <= body_reflectance and body_reflectance <= 1.0
+        ), "Body reflectance has to be 0 to 1"
+
+        actor._mass = actor_mass
+        actor._thermal_model = ThermalModel(
+            local_actor=actor,
+            actor_initial_temperature_in_K=actor_initial_temperature_in_K,
+            actor_sun_absorptance=actor_sun_absorptance,
+            actor_infrared_absorptance=actor_infrared_absorptance,
+            actor_sun_facing_area=actor_sun_facing_area,
+            actor_central_body_facing_area=actor_central_body_facing_area,
+            actor_emissive_area=actor_emissive_area,
+            actor_thermal_capacity=actor_thermal_capacity,
+            body_solar_irradiance=body_solar_irradiance,
+            body_surface_temperature_in_K=body_surface_temperature_in_K,
+            body_emissivity=body_emissivity,
+            body_reflectance=body_reflectance,
+            power_consumption_to_heat_ratio=power_consumption_to_heat_ratio,
+        )
+
     def add_comm_device(actor: BaseActor, device_name: str, bandwidth_in_kbps: float):
         """Creates a communication device.
 

paseos/actors/base_actor.py

@@ -1,11 +1,11 @@
+from abc import ABC
+
 from loguru import logger
 import pykep as pk
-
+import numpy as np
 from skspatial.objects import Sphere
-
-from abc import ABC
-
 from dotmap import DotMap
+
 from ..communication.get_communication_window import get_communication_window
 from ..communication.is_in_line_of_sight import is_in_line_of_sight
 from ..power.is_in_eclipse import is_in_eclipse
@@ -44,9 +44,12 @@ class BaseActor(ABC):
     _current_activity = None
 
     # The following variables are used to track last evaluated state vectors to avoid recomputation.
-    _last_position = None
-    _last_velocity = None
-    _last_eclipse_status = None
+    _previous_position = None
+    _time_of_previous_position = None
+    _previous_velocity = None
+    _previous_eclipse_status = None
+    _time_of_previous_eclipse_status = None
+    _previous_altitude = None
 
     def __init__(self, name: str, epoch: pk.epoch) -> None:
         """Constructor for a base actor
@@ -134,6 +137,32 @@ def discharge(self, consumption_rate_in_W: float, duration_in_s: float):
         """
         pass
 
+    @property
+    def altitude(
+        self,
+        t0: pk.epoch = None,
+    ) -> float:
+        """Returns altitude above [0,0,0]. Will only compute if not computed for this timestep.
+
+        Args:
+            t0 (pk.epoch): Epoch to get altitude at. Defaults to local time.
+
+        Returns:
+            float: Altitude in meters.
+        """
+        if t0 is None:
+            t0 = self._local_time
+        if (
+            t0.mjd2000 == self._time_of_previous_position
+            and self._previous_altitude is not None
+        ):
+            return self._previous_altitude
+        else:
+            self._previous_altitude = np.sqrt(
+                np.sum(np.power(self.get_position(t0), 2))
+            )
+            return self._previous_altitude
+
     def get_position(self, epoch: pk.epoch):
         """Compute the position of this actor at a specific time. Requires orbital parameters or position set.
 
@@ -159,7 +188,8 @@ def get_position(self, epoch: pk.epoch):
         # If the actor has no orbit, return position
         if self._orbital_parameters is None:
             if self._position is not None:
-                self._last_position = self._position
+                self._previous_position = self._position
+                self._time_of_previous_position = epoch.mjd2000
                 return self._position
         else:
             return self._orbital_parameters.eph(epoch)[0]
@@ -190,8 +220,9 @@ def get_position_velocity(self, epoch: pk.epoch):
             + " (mjd2000)."
         )
         pos, vel = self._orbital_parameters.eph(epoch)
-        self._last_position = pos
-        self._last_velocity = vel
+        self._previous_position = pos
+        self._previous_velocity = vel
+        self._time_of_previous_position = epoch.mjd2000
         return pos, vel
 
     def is_in_line_of_sight(
@@ -226,8 +257,12 @@ def is_in_eclipse(self, t: pk.epoch = None):
         """
         if t is None:
             t = self._local_time
-        self._last_eclipse_status = is_in_eclipse(self, self._central_body, t)
-        return self._last_eclipse_status
+        if t.mjd2000 == self._time_of_previous_eclipse_status:
+            return self._previous_eclipse_status
+        else:
+            self._previous_eclipse_status = is_in_eclipse(self, self._central_body, t)
+            self._time_of_previous_eclipse_status = t.mjd2000
+        return self._previous_eclipse_status
 
     def get_communication_window(
         self,

paseos/actors/spacecraft_actor.py

@@ -16,6 +16,11 @@ class SpacecraftActor(BaseActor):
     _max_battery_level_in_Ws = None
     _charging_rate_in_W = None
 
+    # Actor's mass in kg
+    _mass = None
+
+    _thermal_model = None
+
     def __init__(
         self,
         name: str,
@@ -48,6 +53,24 @@ def battery_level_in_Ws(self):
         """
         return self._battery_level_in_Ws
 
+    @property
+    def mass(self) -> float:
+        """Returns actor's mass in kg.
+
+        Returns:
+            float: Mass
+        """
+        return self._mass
+
+    @property
+    def temperature_in_K(self) -> float:
+        """Returns the current temperature of the actor in K.
+
+        Returns:
+            float: Temperature in Kelvin.
+        """
+        return self._thermal_model.temperature_in_K
+
     @property
     def state_of_charge(self):
         """Get the current battery level as ratio of maximum.

paseos/paseos.py

@@ -36,7 +36,7 @@ class PASEOS:
 
     # Used to monitor the local actor over execution and write performance stats
     _operations_monitor = None
-    _time_since_last_log = sys.float_info.max
+    _time_since_previous_log = sys.float_info.max
 
     # Use automatic clock (default on for now)
     use_automatic_clock = True
@@ -113,11 +113,11 @@ def advance_time(self, time_to_advance: float):
             self._local_actor.set_time(pk.epoch(self._state.time * pk.SEC2DAY))
 
             # Check if we should update the status log
-            if self._time_since_last_log > self._cfg.io.logging_interval:
+            if self._time_since_previous_log > self._cfg.io.logging_interval:
                 self.log_status()
-                self._time_since_last_log = 0
+                self._time_since_previous_log = 0
             else:
-                self._time_since_last_log += dt
+                self._time_since_previous_log += dt
 
         logger.debug("New time is: " + str(self._state.time) + " s.")
 

paseos/power/is_in_eclipse.py

@@ -21,7 +21,9 @@ def is_in_eclipse(actor, central_body: pk.planet, t: pk.epoch, plot=False) -> bo
     try:
         logger.debug(f"Checking whether {actor} is in eclipse at {t}.")
     except RuntimeError:
-        logger.debug(f"Checking whether {actor} is in eclipse at {t.mjd2000} (mjd2000).")
+        logger.debug(
+            f"Checking whether {actor} is in eclipse at {t.mjd2000} (mjd2000)."
+        )
 
     # Compute central body position in solar reference frame
     r_central_body_heliocentric, _ = np.array(central_body.eph(t))
@@ -31,7 +33,7 @@ def is_in_eclipse(actor, central_body: pk.planet, t: pk.epoch, plot=False) -> bo
     )
 
     # Compute satellite / actor position in solar reference frame
-    r_sat_central_body_frame = np.array(actor.get_position_velocity(t)[0])
+    r_sat_central_body_frame = np.array(actor.get_position(t))
     logger.trace("r_sat_central_body_frame is" + str(r_sat_central_body_frame))
     r_sat_heliocentric = r_central_body_heliocentric + r_sat_central_body_frame
     logger.trace("r_sat_heliocentric is" + str(r_sat_heliocentric))