From c2f564d0a51b585bed03e6d00eab255880101028 Mon Sep 17 00:00:00 2001 From: erkkoi <49675834+ErkkoI@users.noreply.github.com> Date: Sat, 7 Mar 2026 14:12:59 +0200 Subject: [PATCH 1/8] Fixed heat gain to I*I*R instead of I * R --- linerate/models/cigre601.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/linerate/models/cigre601.py b/linerate/models/cigre601.py index 7ab9c6e..72edf69 100644 --- a/linerate/models/cigre601.py +++ b/linerate/models/cigre601.py @@ -138,7 +138,7 @@ def compute_temperature_gradient( I = current / n # noqa R = self.compute_resistance(conductor_temperature=T_c, current=I) return cigre601.convective_cooling.compute_temperature_gradient( - total_heat_gain=I * R, + total_heat_gain=I * I * R, conductor_thermal_conductivity=self.span.conductor.thermal_conductivity, # type: ignore # noqa core_diameter=self.span.conductor.core_diameter, conductor_diameter=self.span.conductor.conductor_diameter, From 84ad94a5be2d67cebc514eef9b75ceddd924cf1b Mon Sep 17 00:00:00 2001 From: erkkoi <49675834+ErkkoI@users.noreply.github.com> Date: Tue, 10 Mar 2026 21:11:30 +0200 Subject: [PATCH 2/8] Implemented temperature tracking and temporary ampacity solver based on it. --- linerate/equations/heat_capacity.py | 20 +++ linerate/models/thermal_model.py | 78 +++++++++++- linerate/solver.py | 60 ++++++++- linerate/types.py | 56 ++++++++- linerate/units.py | 4 + .../test_cigre_transient_cases.py | 114 ++++++++++++++++++ tests/equations/test_heat_capacity.py | 55 +++++++++ tests/integration_tests/test_vectorization.py | 100 ++++++++------- 8 files changed, 434 insertions(+), 53 deletions(-) create mode 100644 linerate/equations/heat_capacity.py create mode 100644 tests/acceptance_tests/test_cigre_transient_cases.py create mode 100644 tests/equations/test_heat_capacity.py diff --git a/linerate/equations/heat_capacity.py b/linerate/equations/heat_capacity.py new file mode 100644 index 0000000..f93581f --- /dev/null +++ b/linerate/equations/heat_capacity.py @@ -0,0 +1,20 @@ +from linerate.units import ( + Celsius, + JoulePerKilogramPerKelvin, + KilogramPerMeter, + PerKelvin, + JoulePerMeter, +) + + +def calculate_heat_capacity_per_unit_length( + mass_per_unit_length: KilogramPerMeter, + specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin, + specific_heat_capacity_coefficient: PerKelvin, + conductor_temperature: Celsius, +) -> JoulePerMeter: + m = mass_per_unit_length + c = specific_heat_capacity_at_20_celsius + T = conductor_temperature + beta = specific_heat_capacity_coefficient + return m * c * (1 + (T - 20) * beta) diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py index 61bab94..a9e86bf 100644 --- a/linerate/models/thermal_model.py +++ b/linerate/models/thermal_model.py @@ -1,10 +1,12 @@ from abc import ABC, abstractmethod from typing import Dict +import numpy as np +from functools import partial from linerate import solver -from linerate.equations import joule_heating, radiative_cooling -from linerate.types import BaseWeather, Span -from linerate.units import Ampere, Celsius, OhmPerMeter, WattPerMeter +from linerate.equations import joule_heating, radiative_cooling, heat_capacity +from linerate.types import BaseWeather, Span, conductor_heat_capacity_defined +from linerate.units import Ampere, Celsius, Duration, OhmPerMeter, WattPerMeter def _copy_method_docstring(parent_class): @@ -269,3 +271,73 @@ def compute_conductor_temperature( tolerance=tolerance, ) return T + + def compute_heat_capacity_per_unit_length(self, conductor_temperature: Celsius): + conductor = self.span.conductor + if conductor_heat_capacity_defined(conductor): + mc_s = heat_capacity.calculate_heat_capacity_per_unit_length( + mass_per_unit_length=conductor.steel_mass_per_unit_length, + specific_heat_capacity_at_20_celsius=conductor.steel_specific_heat_capacity_at_20_celsius, + specific_heat_capacity_coefficient=conductor.steel_specific_heat_capacity_temperature_coefficient, + conductor_temperature=conductor_temperature, + ) + mc_a = heat_capacity.calculate_heat_capacity_per_unit_length( + mass_per_unit_length=conductor.aluminum_mass_per_unit_length, + specific_heat_capacity_at_20_celsius=conductor.aluminum_specific_heat_capacity_at_20_celsius, + specific_heat_capacity_coefficient=conductor.aluminum_specific_heat_capacity_temperature_coefficient, + conductor_temperature=conductor_temperature, + ) + return mc_s + mc_a + else: + raise RuntimeError("Heat capacity data must be defined to calculate heat capacity.") + + def compute_final_temperature( + self, + initial_conductor_temperature: Celsius, + heating_time: Duration, + current: Ampere, + time_step: Duration = np.timedelta64(60, "s"), + ) -> Ampere: + time_step = np.broadcast_to(time_step, heating_time.shape) + step_count = heating_time // time_step + remainder = heating_time % time_step + dt = time_step / np.timedelta64(1, "s") + modification_mask = step_count > 0 + temperature = initial_conductor_temperature + while np.any(modification_mask): + heat_capacity_ = self.compute_heat_capacity_per_unit_length(temperature) + heat_balance = self.compute_heat_balance(temperature, current=current) + dT = (heat_balance / heat_capacity_) * dt + temperature = temperature + modification_mask * dT + step_count -= 1 + modification_mask = step_count > 0 + if np.any(remainder > 0): + heat_capacity_ = self.compute_heat_capacity_per_unit_length(temperature) + heat_balance = self.compute_heat_balance(temperature, current=current) + dT = (heat_balance / heat_capacity_) * remainder + temperature = temperature + dT + return temperature + + def compute_temporary_ampacity( + self, + max_conductor_temperature: Celsius, + heating_time: Duration, + initial_conductor_temperature: Celsius, + time_step: Duration = np.timedelta64(60, "s"), + min_ampacity: Ampere = 0, + max_ampacity: Ampere = 5000, + tolerance: float = 1.0, + accept_invalid_values: bool = False, + ): + n = self.span.num_conductors + I = solver.compute_conductor_transient_ampacity( # noqa + partial(self.compute_final_temperature, time_step=time_step), + max_conductor_temperature, + initial_conductor_temperature, + heating_time, + min_ampacity, + max_ampacity, + tolerance, + accept_invalid_values, + ) + return I * n diff --git a/linerate/solver.py b/linerate/solver.py index e6fc755..92dfd39 100644 --- a/linerate/solver.py +++ b/linerate/solver.py @@ -3,7 +3,7 @@ import numpy as np -from .units import Ampere, Celsius, FloatOrFloatArray, WattPerMeter +from .units import Ampere, Celsius, Duration, FloatOrFloatArray, WattPerMeter __all__ = ["bisect", "compute_conductor_temperature", "compute_conductor_ampacity"] @@ -164,3 +164,61 @@ def compute_conductor_ampacity( return bisect( f, min_ampacity, max_ampacity, tolerance, accept_invalid_values=accept_invalid_values ) + + +def compute_conductor_transient_ampacity( + final_temperature: Callable[[Celsius, Duration, Ampere], WattPerMeter], + max_conductor_temperature: Celsius, + initial_conductor_temperature: Celsius, + heating_time: Duration, + min_ampacity: Ampere = 0, + max_ampacity: Ampere = 5_000, + tolerance: float = 1, # Ampere + accept_invalid_values: bool = False, +) -> Ampere: + r"""Use the bisection method to compute the temporary thermal rating (ampacity). + + Parameters + ---------- + final_temperature: + :math:`f(T, A) = P_J + P_s - P_c - P_r~\left[\text{W}~\text{m}^{-1}\right]`. A function of + temperature, heating time and current that returns the conductor temperature at the end of heating time. + initial_conductor_temperature: + :math:`T_\text{inital}~\left[^\circ\text{C}\right]`. Initial conductor temperature + heating_time: + :math:`t_\text{heating}~\left[\text{s}\right]`. Time the conductor is heated. + max_conductor_temperature: + :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature + min_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Lower bound for the numerical scheme for + computing the ampacity + max_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Upper bound for the numerical scheme for + computing the ampacity + tolerance: + :math:`\Delta I~\left[\text{A}\right]`. The numerical accuracy of the ampacity. The + bisection iterations will stop once the numerical ampacity uncertainty is below + :math:`\Delta I`. The bisection method will run for + :math:`\left\lceil\frac{I_\text{max} - I_\text{min}}{\Delta I}\right\rceil` iterations. + accept_invalid_values: + If True, np.nan is returned whenever the current cannot be found within the provided + search interval. If False, a ValueError will be raised instead. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I~\left[\text{A}\right]`. The thermal rating. + """ + + def temperature_difference(current: Ampere): + return max_conductor_temperature - final_temperature( + initial_conductor_temperature, heating_time, current + ) + + return bisect( + temperature_difference, + min_ampacity, + max_ampacity, + tolerance, + accept_invalid_values=accept_invalid_values, + ) diff --git a/linerate/types.py b/linerate/types.py index 599e158..253e297 100644 --- a/linerate/types.py +++ b/linerate/types.py @@ -1,6 +1,6 @@ from dataclasses import dataclass from functools import cached_property -from typing import Optional +from typing import Optional, TypeGuard import numpy as np import pygeodesy @@ -8,18 +8,29 @@ from .units import ( Celsius, Degrees, + JoulePerKilogramPerKelvin, + KilogramPerMeter, Meter, MeterPerSecond, OhmPerMeter, + PerKelvin, Radian, SquareMeter, SquareMeterPerAmpere, Unitless, - WattPerMeterPerKelvin, WattPerSquareMeter, + WattPerMeterPerKelvin, ) -__all__ = ["Conductor", "Weather", "Tower", "Span", "WeatherWithSolarRadiation"] +__all__ = [ + "Conductor", + "BaseWeather", + "Weather", + "Tower", + "Span", + "WeatherWithSolarRadiation", + "conductor_heat_capacity_defined", +] @dataclass(frozen=True) @@ -74,6 +85,45 @@ class Conductor: #: with aluminium strands under a tension of at least 40 N :cite:p:`cigre601`. thermal_conductivity: Optional[WattPerMeterPerKelvin] = None + steel_mass_per_unit_length: Optional[KilogramPerMeter] = None + + steel_specific_heat_capacity_at_20_celsius: Optional[JoulePerKilogramPerKelvin] = None + + steel_specific_heat_capacity_temperature_coefficient: Optional[PerKelvin] = None + + aluminum_mass_per_unit_length: Optional[KilogramPerMeter] = None + + aluminum_specific_heat_capacity_at_20_celsius: Optional[JoulePerKilogramPerKelvin] = None + + aluminum_specific_heat_capacity_temperature_coefficient: Optional[PerKelvin] = None + + +class ConductorWithTransientData(Conductor): + steel_mass_per_unit_length: KilogramPerMeter + + steel_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin + + steel_specific_heat_capacity_temperature_coefficient: PerKelvin + + aluminum_mass_per_unit_length: KilogramPerMeter + + aluminum_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin + + aluminum_specific_heat_capacity_temperature_coefficient: PerKelvin + + +def conductor_heat_capacity_defined(conductor: Conductor) -> TypeGuard[ConductorWithTransientData]: + return all( + [ + conductor.steel_mass_per_unit_length is not None, + conductor.steel_specific_heat_capacity_at_20_celsius is not None, + conductor.steel_specific_heat_capacity_temperature_coefficient is not None, + conductor.aluminum_mass_per_unit_length is not None, + conductor.aluminum_specific_heat_capacity_at_20_celsius is not None, + conductor.aluminum_specific_heat_capacity_temperature_coefficient is not None, + ] + ) + @dataclass(frozen=True) class Tower: diff --git a/linerate/units.py b/linerate/units.py index af44c46..a833aae 100644 --- a/linerate/units.py +++ b/linerate/units.py @@ -32,5 +32,9 @@ SquareMeterPerSecond = Annotated[FloatOrFloatArray, "m²/s"] KilogramPerMeterPerSecond = Annotated[FloatOrFloatArray, "kg/(m s)"] KilogramPerCubeMeter = Annotated[FloatOrFloatArray, "kg/m³"] +KilogramPerMeter = Annotated[FloatOrFloatArray, "kg/m"] +PerKelvin = Annotated[FloatOrFloatArray, "1/K"] +JoulePerMeter = Annotated[FloatOrFloatArray, "J/m"] Date = Union[np.datetime64, npt.NDArray[np.datetime64]] +Duration = Union[np.timedelta64, npt.NDArray[np.timedelta64]] diff --git a/tests/acceptance_tests/test_cigre_transient_cases.py b/tests/acceptance_tests/test_cigre_transient_cases.py new file mode 100644 index 0000000..4bf866b --- /dev/null +++ b/tests/acceptance_tests/test_cigre_transient_cases.py @@ -0,0 +1,114 @@ +import numpy as np +import pytest +from pytest import approx +from linerate import Conductor, Weather, Tower, Span, Cigre601 + + +@pytest.fixture +def drake_conductor_transient() -> Conductor: + return Conductor( + conductor_diameter=28.143e-3, + core_diameter=10.4e-3, + outer_layer_strand_diameter=4.4e-3, + emissivity=0.8, + solar_absorptivity=0.8, + temperature1=25, + temperature2=75, + resistance_at_temperature1=0.0727e-3, + resistance_at_temperature2=0.0872e-3, + aluminium_cross_section_area=float("nan"), + constant_magnetic_effect=1, + current_density_proportional_magnetic_effect=0, + max_magnetic_core_relative_resistance_increase=1, + steel_mass_per_unit_length=0.5119, + steel_specific_heat_capacity_at_20_celsius=481, + steel_specific_heat_capacity_temperature_coefficient=1e-4, + aluminum_mass_per_unit_length=1.116, + aluminum_specific_heat_capacity_at_20_celsius=897, + aluminum_specific_heat_capacity_temperature_coefficient=3.8e-4, + ) + + +@pytest.fixture +def initial_weather() -> Weather: + return Weather( + air_temperature=24, + wind_speed=1.9, + wind_direction=np.radians(35), + ground_albedo=1.0, + ) + + +@pytest.fixture +def weather_1() -> Weather: + return Weather( + air_temperature=23.7, + wind_speed=1.7, + wind_direction=np.radians(28), + ground_albedo=1.0, + ) + + +@pytest.fixture +def weather_2() -> Weather: + return Weather( + air_temperature=23.5, + wind_speed=0.8, + wind_direction=np.radians(53), + ground_albedo=1.0, + ) + + +@pytest.fixture +def example_span_transient(drake_conductor_transient: Conductor) -> Span: + start_tower = Tower(latitude=50, longitude=0.0001, altitude=0) + end_tower = Tower(latitude=50, longitude=-0.0001, altitude=0) + return Span( + conductor=drake_conductor_transient, + start_tower=start_tower, + end_tower=end_tower, + num_conductors=1, + ) + + +@pytest.fixture +def transient_initial_model(example_span_transient, initial_weather) -> Cigre601: + return Cigre601( + span=example_span_transient, weather=initial_weather, time=np.datetime64("2016-06-10 00:00") + ) + + +@pytest.fixture +def transient_model_1(example_span_transient, weather_1) -> Cigre601: + return Cigre601( + span=example_span_transient, weather=weather_1, time=np.datetime64("2016-06-10 00:10") + ) + + +@pytest.fixture +def transient_model_2(example_span_transient, weather_2) -> Cigre601: + return Cigre601( + span=example_span_transient, weather=weather_2, time=np.datetime64("2016-06-10 00:20") + ) + + +def test_initial_condition(transient_initial_model): + assert transient_initial_model.compute_conductor_temperature(802, tolerance=0.01) == approx( + 42.01, abs=0.01 + ) + + +def test_tracking(transient_initial_model, transient_model_1, transient_model_2): + start_temperature = transient_initial_model.compute_conductor_temperature(802, tolerance=0.01) + results_1 = [ + transient_model_1.compute_final_temperature(start_temperature, np.timedelta64(t, "m"), 819) + for t in range(1, 11) + ] + results_2 = [ + transient_model_2.compute_final_temperature(results_1[-1], np.timedelta64(t, "m"), 856) + for t in range(1, 11) + ] + correct_1 = [42.175, 42.321, 42.449, 42.562, 42.662, 42.750, 42.828, 42.897, 42.958, 43.011] + correct_2 = [44.147, 45.199, 46.174, 47.075, 47.910, 48.682, 49.396, 50.057, 50.668, 51.233] + assert all(result == approx(correct, abs=0.01) for result, correct in zip(results_1, correct_1)) + assert all(result == approx(correct, abs=0.01) for result, correct in zip(results_2, correct_2)) diff --git a/tests/equations/test_heat_capacity.py b/tests/equations/test_heat_capacity.py new file mode 100644 index 0000000..fe1c2ea --- /dev/null +++ b/tests/equations/test_heat_capacity.py @@ -0,0 +1,55 @@ +from pytest import approx + +import hypothesis +import hypothesis.strategies as st + +from linerate.units import Celsius, KilogramPerMeter +from linerate.equations import heat_capacity + + +@hypothesis.given(mass_per_unit_length=st.floats(min_value=0, max_value=1e10, allow_nan=False)) +def test_heat_capacity_is_linear_in_mass_per_unit_length(mass_per_unit_length: KilogramPerMeter): + m = mass_per_unit_length + c = 1 + beta = 1 + T = 20 + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(m) + + +@hypothesis.given( + specific_heat_capacity_at_20_celsius=st.floats(min_value=0, max_value=1e10, allow_nan=False) +) +def test_heat_capacity_is_linear_in_specific_heat_capacity( + specific_heat_capacity_at_20_celsius: KilogramPerMeter, +): + m = 1 + c = specific_heat_capacity_at_20_celsius + beta = 1 + T = 20 + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(c) + + +@hypothesis.given( + specific_heat_capacity_coefficient=st.floats(min_value=0, max_value=1e10, allow_nan=False) +) +def test_heat_capacity_is_linear_in_heat_capacity_coefficient( + specific_heat_capacity_coefficient: KilogramPerMeter, +): + m = 1 + c = 1 + beta = specific_heat_capacity_coefficient + T = 21 + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(beta + 1) + + +@hypothesis.given(temperature=st.floats(min_value=-1e10, max_value=1e10, allow_nan=False)) +def test_heat_capacity_is_linear_in_temperature(temperature: Celsius): + m = 1 + c = 1 + beta = 1 + T = temperature + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(T - 19) diff --git a/tests/integration_tests/test_vectorization.py b/tests/integration_tests/test_vectorization.py index 4adb809..65dc400 100644 --- a/tests/integration_tests/test_vectorization.py +++ b/tests/integration_tests/test_vectorization.py @@ -1,67 +1,75 @@ import hypothesis -import hypothesis.strategies as st +from hypothesis.extra.numpy import mutually_broadcastable_shapes import numpy as np from linerate.models.cigre601 import Cigre601 from linerate.types import Conductor, Span, Tower, Weather +NUMBER_OF_VARIABLES = 33 -@hypothesis.given( - vectorization_indices=st.sets(st.integers(min_value=1, max_value=24), max_size=10) -) -def test_vectorization(vectorization_indices: set[int]): - def get_zeros(index, length): - out = np.ones(length, dtype=int) - if index in vectorization_indices: - out[index] = 2 - return np.zeros(out) + +class Shapes: + def __init__(self, shapes): + self.shapes = shapes + self.index = 0 + + def get_zeros(self, dtype=None): + out = np.zeros(self.shapes[self.index], dtype=dtype) + self.index += 1 + return out + + +@hypothesis.given(array_shapes=mutually_broadcastable_shapes(num_shapes=NUMBER_OF_VARIABLES)) +def test_vectorization(array_shapes): + shapes = Shapes(array_shapes.input_shapes) conductor = Conductor( - core_diameter=get_zeros(0, 25) + 1, - conductor_diameter=get_zeros(1, 25) + 2, - outer_layer_strand_diameter=get_zeros(2, 25), - emissivity=get_zeros(3, 25), - solar_absorptivity=get_zeros(4, 25), - temperature1=get_zeros(5, 25), - temperature2=get_zeros(6, 25) + 1, - resistance_at_temperature1=get_zeros(7, 25), - resistance_at_temperature2=get_zeros(8, 25), - aluminium_cross_section_area=1 + get_zeros(9, 25), - constant_magnetic_effect=get_zeros(10, 25), - current_density_proportional_magnetic_effect=get_zeros(11, 25) + 1, - max_magnetic_core_relative_resistance_increase=get_zeros(12, 25), + core_diameter=1, # Core diameter is not used in any calculations so its shape will not be broadcast + conductor_diameter=shapes.get_zeros() + 2, + outer_layer_strand_diameter=shapes.get_zeros(), + emissivity=shapes.get_zeros(), + solar_absorptivity=shapes.get_zeros(), + temperature1=shapes.get_zeros(), + temperature2=shapes.get_zeros() + 1, + resistance_at_temperature1=shapes.get_zeros(), + resistance_at_temperature2=shapes.get_zeros(), + aluminium_cross_section_area=1 + shapes.get_zeros(), + constant_magnetic_effect=shapes.get_zeros(), + current_density_proportional_magnetic_effect=shapes.get_zeros() + 1, + max_magnetic_core_relative_resistance_increase=shapes.get_zeros(), + steel_mass_per_unit_length=shapes.get_zeros(), + steel_specific_heat_capacity_at_20_celsius=shapes.get_zeros(), + steel_specific_heat_capacity_temperature_coefficient=shapes.get_zeros(), + aluminum_mass_per_unit_length=shapes.get_zeros(), + aluminum_specific_heat_capacity_at_20_celsius=shapes.get_zeros(), + aluminum_specific_heat_capacity_temperature_coefficient=shapes.get_zeros(), ) weather = Weather( - air_temperature=get_zeros(13, 25), - wind_direction=get_zeros(14, 25), - wind_speed=get_zeros(15, 25), - ground_albedo=get_zeros(23, 25), - clearness_ratio=get_zeros(16, 25), + air_temperature=shapes.get_zeros(), + wind_direction=shapes.get_zeros(), + wind_speed=shapes.get_zeros(), + ground_albedo=shapes.get_zeros(), + clearness_ratio=shapes.get_zeros(), ) span = Span( conductor=conductor, start_tower=Tower( - get_zeros(17, 25), - get_zeros(18, 25), - get_zeros(19, 25), + shapes.get_zeros(), + shapes.get_zeros(), + shapes.get_zeros(), ), end_tower=Tower( - 0.1 + get_zeros(20, 25), - 0.1 + get_zeros(21, 25), - 10 + get_zeros(22, 25), + 0.1 + shapes.get_zeros(), + 0.1 + shapes.get_zeros(), + 10 + shapes.get_zeros(), ), num_conductors=1, ) - time = np.array([np.datetime64("2022-01-01")]).reshape([1] * 25) - if 24 in vectorization_indices: - time = np.concatenate([time, time], axis=-1) - + time = shapes.get_zeros(dtype="datetime64[D]") model = Cigre601(span=span, weather=weather, time=time) - temperature = model.compute_conductor_temperature(100) - shape = [1] * 25 - for vectorisation_index in vectorization_indices: - shape[vectorisation_index] = 2 - shape[0] = 1 - - assert isinstance(temperature, np.ndarray) - assert temperature.shape == tuple(shape) + final_temperature = model.compute_final_temperature( + shapes.get_zeros() + 20, + shapes.get_zeros(dtype=" Date: Sun, 10 May 2026 13:18:35 +0300 Subject: [PATCH 3/8] Addressed comments - Added docstrings - Added some type information - Made ConductorWithTransientData to a proper type --- linerate/models/cigre601.py | 3 +- linerate/models/thermal_model.py | 94 ++++++++++++++----- linerate/solver.py | 10 +- linerate/types.py | 36 ++----- .../test_cigre_transient_cases.py | 16 ++-- tests/integration_tests/test_vectorization.py | 89 +++++++++++++++--- 6 files changed, 173 insertions(+), 75 deletions(-) diff --git a/linerate/models/cigre601.py b/linerate/models/cigre601.py index 72edf69..7fbcc1b 100644 --- a/linerate/models/cigre601.py +++ b/linerate/models/cigre601.py @@ -136,9 +136,8 @@ def compute_temperature_gradient( n = self.span.num_conductors T_c = conductor_temperature I = current / n # noqa - R = self.compute_resistance(conductor_temperature=T_c, current=I) return cigre601.convective_cooling.compute_temperature_gradient( - total_heat_gain=I * I * R, + total_heat_gain=self.compute_joule_heating(conductor_temperature=T_c, current=I), conductor_thermal_conductivity=self.span.conductor.thermal_conductivity, # type: ignore # noqa core_diameter=self.span.conductor.core_diameter, conductor_diameter=self.span.conductor.conductor_diameter, diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py index a9e86bf..df67d3f 100644 --- a/linerate/models/thermal_model.py +++ b/linerate/models/thermal_model.py @@ -5,7 +5,7 @@ from functools import partial from linerate import solver from linerate.equations import joule_heating, radiative_cooling, heat_capacity -from linerate.types import BaseWeather, Span, conductor_heat_capacity_defined +from linerate.types import BaseWeather, ConductorWithTransientData, Span from linerate.units import Ampere, Celsius, Duration, OhmPerMeter, WattPerMeter @@ -274,30 +274,47 @@ def compute_conductor_temperature( def compute_heat_capacity_per_unit_length(self, conductor_temperature: Celsius): conductor = self.span.conductor - if conductor_heat_capacity_defined(conductor): - mc_s = heat_capacity.calculate_heat_capacity_per_unit_length( - mass_per_unit_length=conductor.steel_mass_per_unit_length, - specific_heat_capacity_at_20_celsius=conductor.steel_specific_heat_capacity_at_20_celsius, - specific_heat_capacity_coefficient=conductor.steel_specific_heat_capacity_temperature_coefficient, - conductor_temperature=conductor_temperature, - ) - mc_a = heat_capacity.calculate_heat_capacity_per_unit_length( - mass_per_unit_length=conductor.aluminum_mass_per_unit_length, - specific_heat_capacity_at_20_celsius=conductor.aluminum_specific_heat_capacity_at_20_celsius, - specific_heat_capacity_coefficient=conductor.aluminum_specific_heat_capacity_temperature_coefficient, - conductor_temperature=conductor_temperature, - ) - return mc_s + mc_a - else: - raise RuntimeError("Heat capacity data must be defined to calculate heat capacity.") - - def compute_final_temperature( + if not isinstance(conductor, ConductorWithTransientData): + raise ValueError("Heat capacity data must be defined to calculate heat capacity.") + mc_s = heat_capacity.calculate_heat_capacity_per_unit_length( + mass_per_unit_length=conductor.steel_mass_per_unit_length, + specific_heat_capacity_at_20_celsius=conductor.steel_specific_heat_capacity_at_20_celsius, + specific_heat_capacity_coefficient=conductor.steel_specific_heat_capacity_temperature_coefficient, + conductor_temperature=conductor_temperature, + ) + mc_a = heat_capacity.calculate_heat_capacity_per_unit_length( + mass_per_unit_length=conductor.aluminum_mass_per_unit_length, + specific_heat_capacity_at_20_celsius=conductor.aluminum_specific_heat_capacity_at_20_celsius, + specific_heat_capacity_coefficient=conductor.aluminum_specific_heat_capacity_temperature_coefficient, + conductor_temperature=conductor_temperature, + ) + return mc_s + mc_a + + def compute_temperature_after_heating_time( self, initial_conductor_temperature: Celsius, heating_time: Duration, current: Ampere, time_step: Duration = np.timedelta64(60, "s"), - ) -> Ampere: + ) -> Celsius: + r"""Compute conductor temperature after it is heated for a time using an iterative method. + + Parameters + ---------- + initial_conductor_temperature: + :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating. + heating_time: + :math:`t`. Duration of heating. Unit is changed to seconds in the function. + current: + :math:`I~\left[\text{A}\right]`. Current heating the conductor + time_step: + :math:`\Delta t`. Time step of heating. Unit is changed to seconds in the function. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`T~\left[^\circ\text{C}\right]`. Conductor temperature. + """ time_step = np.broadcast_to(time_step, heating_time.shape) step_count = heating_time // time_step remainder = heating_time % time_step @@ -318,7 +335,7 @@ def compute_final_temperature( temperature = temperature + dT return temperature - def compute_temporary_ampacity( + def compute_transient_ampacity( self, max_conductor_temperature: Celsius, heating_time: Duration, @@ -329,9 +346,42 @@ def compute_temporary_ampacity( tolerance: float = 1.0, accept_invalid_values: bool = False, ): + r"""Use the bisection method to compute the transient thermal rating (ampacity). + + Parameters + ---------- + max_conductor_temperature: + :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature + heating_time: + :math:`t`. Duration of heating. Unit is changed to seconds in the function. + initial_conductor_temperature: + :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating. + time_step: + :math:`\Delta t`. Time step of heating. Unit is changed to seconds in the function. + min_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Lower bound for the numerical scheme for + computing the ampacity + max_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Upper bound for the numerical scheme for + computing the ampacity + tolerance: + :math:`\Delta I~\left[\text{A}\right]`. The numerical accuracy of the ampacity. The + bisection iterations will stop once the numerical ampacity uncertainty is below + :math:`\Delta I`. The bisection method will run for + :math:`\left\lceil\frac{I_\text{min} - I_\text{min}}{\Delta I}\right\rceil` iterations. + accept_invalid_values: + If True, np.nan is returned whenever the current cannot be found within the provided + search interval. If False, a ValueError will be raised instead. + + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I~\left[\text{A}\right]`. The transient thermal rating. + """ n = self.span.num_conductors I = solver.compute_conductor_transient_ampacity( # noqa - partial(self.compute_final_temperature, time_step=time_step), + partial(self.compute_temperature_after_heating_time, time_step=time_step), max_conductor_temperature, initial_conductor_temperature, heating_time, diff --git a/linerate/solver.py b/linerate/solver.py index 92dfd39..f890995 100644 --- a/linerate/solver.py +++ b/linerate/solver.py @@ -170,7 +170,7 @@ def compute_conductor_transient_ampacity( final_temperature: Callable[[Celsius, Duration, Ampere], WattPerMeter], max_conductor_temperature: Celsius, initial_conductor_temperature: Celsius, - heating_time: Duration, + heating_duration: Duration, min_ampacity: Ampere = 0, max_ampacity: Ampere = 5_000, tolerance: float = 1, # Ampere @@ -185,8 +185,8 @@ def compute_conductor_transient_ampacity( temperature, heating time and current that returns the conductor temperature at the end of heating time. initial_conductor_temperature: :math:`T_\text{inital}~\left[^\circ\text{C}\right]`. Initial conductor temperature - heating_time: - :math:`t_\text{heating}~\left[\text{s}\right]`. Time the conductor is heated. + heating_duration: + :math:`\delta t_\text{heating}~\left[\text{s}\right]`. Duration the conductor is heated. max_conductor_temperature: :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature min_ampacity: @@ -210,9 +210,9 @@ def compute_conductor_transient_ampacity( :math:`I~\left[\text{A}\right]`. The thermal rating. """ - def temperature_difference(current: Ampere): + def temperature_difference(current: Ampere) -> Celsius: return max_conductor_temperature - final_temperature( - initial_conductor_temperature, heating_time, current + initial_conductor_temperature, heating_duration, current ) return bisect( diff --git a/linerate/types.py b/linerate/types.py index 253e297..bbfc1ed 100644 --- a/linerate/types.py +++ b/linerate/types.py @@ -1,6 +1,6 @@ from dataclasses import dataclass from functools import cached_property -from typing import Optional, TypeGuard +from typing import Optional import numpy as np import pygeodesy @@ -29,7 +29,7 @@ "Tower", "Span", "WeatherWithSolarRadiation", - "conductor_heat_capacity_defined", + "ConductorWithTransientData", ] @@ -85,46 +85,28 @@ class Conductor: #: with aluminium strands under a tension of at least 40 N :cite:p:`cigre601`. thermal_conductivity: Optional[WattPerMeterPerKelvin] = None - steel_mass_per_unit_length: Optional[KilogramPerMeter] = None - - steel_specific_heat_capacity_at_20_celsius: Optional[JoulePerKilogramPerKelvin] = None - - steel_specific_heat_capacity_temperature_coefficient: Optional[PerKelvin] = None - - aluminum_mass_per_unit_length: Optional[KilogramPerMeter] = None - - aluminum_specific_heat_capacity_at_20_celsius: Optional[JoulePerKilogramPerKelvin] = None - - aluminum_specific_heat_capacity_temperature_coefficient: Optional[PerKelvin] = None - +@dataclass(frozen=True, kw_only=True) class ConductorWithTransientData(Conductor): + #: :math:`m_s`. Mass of steel (core material) per unit length of conductor :math:`\left[\text{kg}\text{m}^{-1}\right]` steel_mass_per_unit_length: KilogramPerMeter + #: :math:`c_s^20`. Heat capacity of steel (core material) at 20 ^\circ C :math:`\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\right]` steel_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin + #: :math:`\beta_20^s`. Heat capacity temperature coefficient of the steel (core material) at 20 ^\circ C :math:`\left[text{K}^{-1}\right]` steel_specific_heat_capacity_temperature_coefficient: PerKelvin + #: :math:`m_a`. Mass of aluminium (outer material) per unit length of conductor :math:`\left[\text{kg}\text{m}^{-1}\right]` aluminum_mass_per_unit_length: KilogramPerMeter + #: :math:`c_a^20`. Heat capacity of aluminium (outer material) at 20 ^\circ C :math:`\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\right]` aluminum_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin + #: :math:`\beta_a^20`. Heat capacity temperature coefficient of the aluminium (outer material) at 20 ^\circ C :math:`\left[text{K}^{-1}\right]` aluminum_specific_heat_capacity_temperature_coefficient: PerKelvin -def conductor_heat_capacity_defined(conductor: Conductor) -> TypeGuard[ConductorWithTransientData]: - return all( - [ - conductor.steel_mass_per_unit_length is not None, - conductor.steel_specific_heat_capacity_at_20_celsius is not None, - conductor.steel_specific_heat_capacity_temperature_coefficient is not None, - conductor.aluminum_mass_per_unit_length is not None, - conductor.aluminum_specific_heat_capacity_at_20_celsius is not None, - conductor.aluminum_specific_heat_capacity_temperature_coefficient is not None, - ] - ) - - @dataclass(frozen=True) class Tower: """Container for a tower (span end-point).""" diff --git a/tests/acceptance_tests/test_cigre_transient_cases.py b/tests/acceptance_tests/test_cigre_transient_cases.py index 4bf866b..237f415 100644 --- a/tests/acceptance_tests/test_cigre_transient_cases.py +++ b/tests/acceptance_tests/test_cigre_transient_cases.py @@ -1,12 +1,12 @@ import numpy as np import pytest from pytest import approx -from linerate import Conductor, Weather, Tower, Span, Cigre601 +from linerate import ConductorWithTransientData, Weather, Tower, Span, Cigre601 @pytest.fixture -def drake_conductor_transient() -> Conductor: - return Conductor( +def drake_conductor_transient() -> ConductorWithTransientData: + return ConductorWithTransientData( conductor_diameter=28.143e-3, core_diameter=10.4e-3, outer_layer_strand_diameter=4.4e-3, @@ -60,7 +60,7 @@ def weather_2() -> Weather: @pytest.fixture -def example_span_transient(drake_conductor_transient: Conductor) -> Span: +def example_span_transient(drake_conductor_transient: ConductorWithTransientData) -> Span: start_tower = Tower(latitude=50, longitude=0.0001, altitude=0) end_tower = Tower(latitude=50, longitude=-0.0001, altitude=0) return Span( @@ -101,11 +101,15 @@ def test_initial_condition(transient_initial_model): def test_tracking(transient_initial_model, transient_model_1, transient_model_2): start_temperature = transient_initial_model.compute_conductor_temperature(802, tolerance=0.01) results_1 = [ - transient_model_1.compute_final_temperature(start_temperature, np.timedelta64(t, "m"), 819) + transient_model_1.compute_temperature_after_heating_time( + start_temperature, np.timedelta64(t, "m"), 819 + ) for t in range(1, 11) ] results_2 = [ - transient_model_2.compute_final_temperature(results_1[-1], np.timedelta64(t, "m"), 856) + transient_model_2.compute_temperature_after_heating_time( + results_1[-1], np.timedelta64(t, "m"), 856 + ) for t in range(1, 11) ] correct_1 = [42.175, 42.321, 42.449, 42.562, 42.662, 42.750, 42.828, 42.897, 42.958, 43.011] diff --git a/tests/integration_tests/test_vectorization.py b/tests/integration_tests/test_vectorization.py index 65dc400..6bb4daf 100644 --- a/tests/integration_tests/test_vectorization.py +++ b/tests/integration_tests/test_vectorization.py @@ -1,29 +1,92 @@ import hypothesis -from hypothesis.extra.numpy import mutually_broadcastable_shapes +import hypothesis.strategies as st +from hypothesis.extra.numpy import mutually_broadcastable_shapes, BroadcastableShapes import numpy as np +import numpy.typing as npt from linerate.models.cigre601 import Cigre601 -from linerate.types import Conductor, Span, Tower, Weather +from linerate.types import Conductor, ConductorWithTransientData, Span, Tower, Weather NUMBER_OF_VARIABLES = 33 class Shapes: - def __init__(self, shapes): + def __init__(self, shapes: tuple[tuple[int, ...], ...]): self.shapes = shapes self.index = 0 - def get_zeros(self, dtype=None): + def get_zeros(self, dtype: np.dtype = np.dtype("f8")) -> npt.NDArray: out = np.zeros(self.shapes[self.index], dtype=dtype) self.index += 1 return out +@hypothesis.given( + vectorization_indices=st.sets(st.integers(min_value=1, max_value=24), max_size=10) +) +def test_vectorization(vectorization_indices: set[int]): + def get_zeros(index, length) -> npt.NDArray[np.floating]: + out = np.ones(length, dtype=int) + if index in vectorization_indices: + out[index] = 2 + return np.zeros(out, dtype=np.dtype("f8")) + + conductor = Conductor( + core_diameter=get_zeros(0, 25) + 1, + conductor_diameter=get_zeros(1, 25) + 2, + outer_layer_strand_diameter=get_zeros(2, 25), + emissivity=get_zeros(3, 25), + solar_absorptivity=get_zeros(4, 25), + temperature1=get_zeros(5, 25), + temperature2=get_zeros(6, 25) + 1, + resistance_at_temperature1=get_zeros(7, 25), + resistance_at_temperature2=get_zeros(8, 25), + aluminium_cross_section_area=1 + get_zeros(9, 25), + constant_magnetic_effect=get_zeros(10, 25), + current_density_proportional_magnetic_effect=get_zeros(11, 25) + 1, + max_magnetic_core_relative_resistance_increase=get_zeros(12, 25), + ) + weather = Weather( + air_temperature=get_zeros(13, 25), + wind_direction=get_zeros(14, 25), + wind_speed=get_zeros(15, 25), + ground_albedo=get_zeros(23, 25), + clearness_ratio=get_zeros(16, 25), + ) + span = Span( + conductor=conductor, + start_tower=Tower( + get_zeros(17, 25), + get_zeros(18, 25), + get_zeros(19, 25), + ), + end_tower=Tower( + 0.1 + get_zeros(20, 25), + 0.1 + get_zeros(21, 25), + 10 + get_zeros(22, 25), + ), + num_conductors=1, + ) + time = np.array([np.datetime64("2022-01-01")]).reshape([1] * 25) + if 24 in vectorization_indices: + time = np.concatenate([time, time], axis=-1) + + model = Cigre601(span=span, weather=weather, time=time) + temperature = model.compute_conductor_temperature(100) + shape = [1] * 25 + for vectorisation_index in vectorization_indices: + shape[vectorisation_index] = 2 + shape[0] = 1 + + assert isinstance(temperature, np.ndarray) + assert temperature.shape == tuple(shape) + + @hypothesis.given(array_shapes=mutually_broadcastable_shapes(num_shapes=NUMBER_OF_VARIABLES)) -def test_vectorization(array_shapes): +def test_transient_vectorization(array_shapes: BroadcastableShapes): shapes = Shapes(array_shapes.input_shapes) - conductor = Conductor( + conductor = ConductorWithTransientData( core_diameter=1, # Core diameter is not used in any calculations so its shape will not be broadcast conductor_diameter=shapes.get_zeros() + 2, outer_layer_strand_diameter=shapes.get_zeros(), @@ -37,11 +100,11 @@ def test_vectorization(array_shapes): constant_magnetic_effect=shapes.get_zeros(), current_density_proportional_magnetic_effect=shapes.get_zeros() + 1, max_magnetic_core_relative_resistance_increase=shapes.get_zeros(), - steel_mass_per_unit_length=shapes.get_zeros(), - steel_specific_heat_capacity_at_20_celsius=shapes.get_zeros(), + steel_mass_per_unit_length=shapes.get_zeros() + 0.5, + steel_specific_heat_capacity_at_20_celsius=shapes.get_zeros() + 500, steel_specific_heat_capacity_temperature_coefficient=shapes.get_zeros(), - aluminum_mass_per_unit_length=shapes.get_zeros(), - aluminum_specific_heat_capacity_at_20_celsius=shapes.get_zeros(), + aluminum_mass_per_unit_length=shapes.get_zeros() + 1, + aluminum_specific_heat_capacity_at_20_celsius=shapes.get_zeros() + 1000, aluminum_specific_heat_capacity_temperature_coefficient=shapes.get_zeros(), ) weather = Weather( @@ -65,11 +128,11 @@ def test_vectorization(array_shapes): ), num_conductors=1, ) - time = shapes.get_zeros(dtype="datetime64[D]") + time = shapes.get_zeros(dtype=np.dtype("datetime64[D]")) model = Cigre601(span=span, weather=weather, time=time) - final_temperature = model.compute_final_temperature( + final_temperature = model.compute_temperature_after_heating_time( shapes.get_zeros() + 20, - shapes.get_zeros(dtype=" Date: Sun, 10 May 2026 13:19:13 +0300 Subject: [PATCH 4/8] Added .* to pytest configuration parameter norecursedirs to fix hypothesis warning --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 4e0368c..5218f7b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -84,7 +84,7 @@ source = ["linerate"] log_file = "pytest.log" log_level = "DEBUG" log_file_format = "%(asctime)s - %(name)s - %(levelname)s: %(message)s" -norecursedirs = "_build tmp* __pycache__ src prof wheel_files" +norecursedirs = ".* _build tmp* __pycache__ src prof wheel_files" markers = [ "integration: Mark test as an integration test", ] From f96c918f40672acaeb49f31212b7317f470f4629 Mon Sep 17 00:00:00 2001 From: erkkoi <49675834+ErkkoI@users.noreply.github.com> Date: Tue, 19 May 2026 22:16:43 +0300 Subject: [PATCH 5/8] Changed .* to .hypothesis to only fix hypothesis warning and not skip all hidden folders. --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 5218f7b..aa8035a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -84,7 +84,7 @@ source = ["linerate"] log_file = "pytest.log" log_level = "DEBUG" log_file_format = "%(asctime)s - %(name)s - %(levelname)s: %(message)s" -norecursedirs = ".* _build tmp* __pycache__ src prof wheel_files" +norecursedirs = ".hypothesis _build tmp* __pycache__ src prof wheel_files" markers = [ "integration: Mark test as an integration test", ] From ff9bb48e8dc85a18a6f25a72f09ca0ad11cd30b7 Mon Sep 17 00:00:00 2001 From: Erkko Ihalainen Date: Wed, 20 May 2026 09:15:50 +0300 Subject: [PATCH 6/8] Changed heating_time -> heating_duration and compute_temperature_after_heating_time -> compute_temperature_after_heating --- linerate/models/thermal_model.py | 20 +++++++++---------- .../test_cigre_transient_cases.py | 4 ++-- tests/integration_tests/test_vectorization.py | 2 +- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py index df67d3f..884d047 100644 --- a/linerate/models/thermal_model.py +++ b/linerate/models/thermal_model.py @@ -290,10 +290,10 @@ def compute_heat_capacity_per_unit_length(self, conductor_temperature: Celsius): ) return mc_s + mc_a - def compute_temperature_after_heating_time( + def compute_temperature_after_heating( self, initial_conductor_temperature: Celsius, - heating_time: Duration, + heating_duration: Duration, current: Ampere, time_step: Duration = np.timedelta64(60, "s"), ) -> Celsius: @@ -303,7 +303,7 @@ def compute_temperature_after_heating_time( ---------- initial_conductor_temperature: :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating. - heating_time: + heating_duration: :math:`t`. Duration of heating. Unit is changed to seconds in the function. current: :math:`I~\left[\text{A}\right]`. Current heating the conductor @@ -315,9 +315,9 @@ def compute_temperature_after_heating_time( Union[float, float64, ndarray[Any, dtype[float64]]] :math:`T~\left[^\circ\text{C}\right]`. Conductor temperature. """ - time_step = np.broadcast_to(time_step, heating_time.shape) - step_count = heating_time // time_step - remainder = heating_time % time_step + time_step = np.broadcast_to(time_step, heating_duration.shape) + step_count = heating_duration // time_step + remainder = heating_duration % time_step dt = time_step / np.timedelta64(1, "s") modification_mask = step_count > 0 temperature = initial_conductor_temperature @@ -338,7 +338,7 @@ def compute_temperature_after_heating_time( def compute_transient_ampacity( self, max_conductor_temperature: Celsius, - heating_time: Duration, + heating_duration: Duration, initial_conductor_temperature: Celsius, time_step: Duration = np.timedelta64(60, "s"), min_ampacity: Ampere = 0, @@ -352,7 +352,7 @@ def compute_transient_ampacity( ---------- max_conductor_temperature: :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature - heating_time: + heating_duration: :math:`t`. Duration of heating. Unit is changed to seconds in the function. initial_conductor_temperature: :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating. @@ -381,10 +381,10 @@ def compute_transient_ampacity( """ n = self.span.num_conductors I = solver.compute_conductor_transient_ampacity( # noqa - partial(self.compute_temperature_after_heating_time, time_step=time_step), + partial(self.compute_temperature_after_heating, time_step=time_step), max_conductor_temperature, initial_conductor_temperature, - heating_time, + heating_duration, min_ampacity, max_ampacity, tolerance, diff --git a/tests/acceptance_tests/test_cigre_transient_cases.py b/tests/acceptance_tests/test_cigre_transient_cases.py index 237f415..38310f1 100644 --- a/tests/acceptance_tests/test_cigre_transient_cases.py +++ b/tests/acceptance_tests/test_cigre_transient_cases.py @@ -101,13 +101,13 @@ def test_initial_condition(transient_initial_model): def test_tracking(transient_initial_model, transient_model_1, transient_model_2): start_temperature = transient_initial_model.compute_conductor_temperature(802, tolerance=0.01) results_1 = [ - transient_model_1.compute_temperature_after_heating_time( + transient_model_1.compute_temperature_after_heating( start_temperature, np.timedelta64(t, "m"), 819 ) for t in range(1, 11) ] results_2 = [ - transient_model_2.compute_temperature_after_heating_time( + transient_model_2.compute_temperature_after_heating( results_1[-1], np.timedelta64(t, "m"), 856 ) for t in range(1, 11) diff --git a/tests/integration_tests/test_vectorization.py b/tests/integration_tests/test_vectorization.py index 6bb4daf..ae25aa6 100644 --- a/tests/integration_tests/test_vectorization.py +++ b/tests/integration_tests/test_vectorization.py @@ -130,7 +130,7 @@ def test_transient_vectorization(array_shapes: BroadcastableShapes): ) time = shapes.get_zeros(dtype=np.dtype("datetime64[D]")) model = Cigre601(span=span, weather=weather, time=time) - final_temperature = model.compute_temperature_after_heating_time( + final_temperature = model.compute_temperature_after_heating( shapes.get_zeros() + 20, shapes.get_zeros(dtype=np.dtype(" Date: Wed, 20 May 2026 12:56:46 +0300 Subject: [PATCH 7/8] Split IVP solver to a more generalized (still basic) function in solver.py. --- linerate/models/thermal_model.py | 44 ++++++++++--------- linerate/solver.py | 34 ++++++++++++++ .../test_cigre_transient_cases.py | 5 ++- 3 files changed, 61 insertions(+), 22 deletions(-) diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py index 884d047..62297b0 100644 --- a/linerate/models/thermal_model.py +++ b/linerate/models/thermal_model.py @@ -5,8 +5,9 @@ from functools import partial from linerate import solver from linerate.equations import joule_heating, radiative_cooling, heat_capacity +from linerate.solver import solve_ivp_forward_euler from linerate.types import BaseWeather, ConductorWithTransientData, Span -from linerate.units import Ampere, Celsius, Duration, OhmPerMeter, WattPerMeter +from linerate.units import Ampere, Celsius, Duration, OhmPerMeter, WattPerMeter, JoulePerKilogramPerKelvin def _copy_method_docstring(parent_class): @@ -272,7 +273,20 @@ def compute_conductor_temperature( ) return T - def compute_heat_capacity_per_unit_length(self, conductor_temperature: Celsius): + def compute_heat_capacity_per_unit_length(self, conductor_temperature: Celsius) -> JoulePerKilogramPerKelvin: + r""" + Compute the heat capacity of the conductor per unit length. + Parameters + ---------- + conductor_temperature: + :math:`T~\left[^\circ\text{C}\right]` Conductor temperature. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`c~\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\text{m}^{-1}\right]`. Conductor heat capacity per unit length. + + """ conductor = self.span.conductor if not isinstance(conductor, ConductorWithTransientData): raise ValueError("Heat capacity data must be defined to calculate heat capacity.") @@ -295,7 +309,7 @@ def compute_temperature_after_heating( initial_conductor_temperature: Celsius, heating_duration: Duration, current: Ampere, - time_step: Duration = np.timedelta64(60, "s"), + time_step: Duration = np.timedelta64(1, "s"), ) -> Celsius: r"""Compute conductor temperature after it is heated for a time using an iterative method. @@ -315,25 +329,15 @@ def compute_temperature_after_heating( Union[float, float64, ndarray[Any, dtype[float64]]] :math:`T~\left[^\circ\text{C}\right]`. Conductor temperature. """ - time_step = np.broadcast_to(time_step, heating_duration.shape) - step_count = heating_duration // time_step - remainder = heating_duration % time_step - dt = time_step / np.timedelta64(1, "s") - modification_mask = step_count > 0 - temperature = initial_conductor_temperature - while np.any(modification_mask): + def conductor_heating(temperature): heat_capacity_ = self.compute_heat_capacity_per_unit_length(temperature) heat_balance = self.compute_heat_balance(temperature, current=current) - dT = (heat_balance / heat_capacity_) * dt - temperature = temperature + modification_mask * dT - step_count -= 1 - modification_mask = step_count > 0 - if np.any(remainder > 0): - heat_capacity_ = self.compute_heat_capacity_per_unit_length(temperature) - heat_balance = self.compute_heat_balance(temperature, current=current) - dT = (heat_balance / heat_capacity_) * remainder - temperature = temperature + dT - return temperature + return heat_balance / heat_capacity_ + # Convert durations to floating point numbers + dt = time_step / np.timedelta64(1, "s") + tf = heating_duration / np.timedelta64(1, "s") + final_temperature = solve_ivp_forward_euler(conductor_heating, initial_conductor_temperature, tf, dt) + return final_temperature def compute_transient_ampacity( self, diff --git a/linerate/solver.py b/linerate/solver.py index f890995..164850e 100644 --- a/linerate/solver.py +++ b/linerate/solver.py @@ -222,3 +222,37 @@ def temperature_difference(current: Ampere) -> Celsius: tolerance, accept_invalid_values=accept_invalid_values, ) + + +def solve_ivp_forward_euler( + f: Callable[[FloatOrFloatArray], FloatOrFloatArray], + y0: FloatOrFloatArray, + duration: FloatOrFloatArray, + step: FloatOrFloatArray, + ) -> FloatOrFloatArray: + """ + A very basic vectorized forward Euler solver time-invariant systems. + Parameters + ---------- + f: Time derivative of the system, assumed to only depend on state. dy / dt = f(y) + y0: Initial state of the system + duration: Time for which to apply the function + step: Time step + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + Approximation of function f at the end of the given duration. + """ + y = y0 + step = np.broadcast_to(step, np.shape(duration)) + step_count = duration // step + remainder = duration % step + modification_mask = step_count > 0 + while np.any(modification_mask): + y = y + f(y) * step + step_count -= 1 + modification_mask = step_count > 0 + if np.any(remainder > 0): + y = y + f(y) * remainder + return y diff --git a/tests/acceptance_tests/test_cigre_transient_cases.py b/tests/acceptance_tests/test_cigre_transient_cases.py index 38310f1..6ec6e7c 100644 --- a/tests/acceptance_tests/test_cigre_transient_cases.py +++ b/tests/acceptance_tests/test_cigre_transient_cases.py @@ -100,15 +100,16 @@ def test_initial_condition(transient_initial_model): def test_tracking(transient_initial_model, transient_model_1, transient_model_2): start_temperature = transient_initial_model.compute_conductor_temperature(802, tolerance=0.01) + # 60s time step is used since the calculations in the brochure use the same step size. results_1 = [ transient_model_1.compute_temperature_after_heating( - start_temperature, np.timedelta64(t, "m"), 819 + start_temperature, np.timedelta64(t, "m"), 819, time_step=np.timedelta64(60, "s") ) for t in range(1, 11) ] results_2 = [ transient_model_2.compute_temperature_after_heating( - results_1[-1], np.timedelta64(t, "m"), 856 + results_1[-1], np.timedelta64(t, "m"), 856, time_step=np.timedelta64(60, "s") ) for t in range(1, 11) ] From 152aaa20282c356ee3edb74b8b6fb6640c253951 Mon Sep 17 00:00:00 2001 From: Erkko Ihalainen Date: Wed, 20 May 2026 12:57:36 +0300 Subject: [PATCH 8/8] Made vectorization tests more similar and increased maximum number of dimensions to 5. --- tests/integration_tests/test_vectorization.py | 75 ++++++++----------- 1 file changed, 32 insertions(+), 43 deletions(-) diff --git a/tests/integration_tests/test_vectorization.py b/tests/integration_tests/test_vectorization.py index ae25aa6..2ac2ee1 100644 --- a/tests/integration_tests/test_vectorization.py +++ b/tests/integration_tests/test_vectorization.py @@ -1,5 +1,4 @@ import hypothesis -import hypothesis.strategies as st from hypothesis.extra.numpy import mutually_broadcastable_shapes, BroadcastableShapes import numpy as np import numpy.typing as npt @@ -7,7 +6,8 @@ from linerate.models.cigre601 import Cigre601 from linerate.types import Conductor, ConductorWithTransientData, Span, Tower, Weather -NUMBER_OF_VARIABLES = 33 +NUMBER_OF_VARIABLES_TRANSIENT = 33 +NUMBER_OF_VARIABLES_STEADY_STATE = 24 class Shapes: @@ -22,67 +22,56 @@ def get_zeros(self, dtype: np.dtype = np.dtype("f8")) -> npt.NDArray: @hypothesis.given( - vectorization_indices=st.sets(st.integers(min_value=1, max_value=24), max_size=10) + array_shapes=mutually_broadcastable_shapes(num_shapes=NUMBER_OF_VARIABLES_STEADY_STATE, max_dims=5) ) -def test_vectorization(vectorization_indices: set[int]): - def get_zeros(index, length) -> npt.NDArray[np.floating]: - out = np.ones(length, dtype=int) - if index in vectorization_indices: - out[index] = 2 - return np.zeros(out, dtype=np.dtype("f8")) +def test_vectorization(array_shapes: BroadcastableShapes): + shapes = Shapes(array_shapes.input_shapes) conductor = Conductor( - core_diameter=get_zeros(0, 25) + 1, - conductor_diameter=get_zeros(1, 25) + 2, - outer_layer_strand_diameter=get_zeros(2, 25), - emissivity=get_zeros(3, 25), - solar_absorptivity=get_zeros(4, 25), - temperature1=get_zeros(5, 25), - temperature2=get_zeros(6, 25) + 1, - resistance_at_temperature1=get_zeros(7, 25), - resistance_at_temperature2=get_zeros(8, 25), - aluminium_cross_section_area=1 + get_zeros(9, 25), - constant_magnetic_effect=get_zeros(10, 25), - current_density_proportional_magnetic_effect=get_zeros(11, 25) + 1, - max_magnetic_core_relative_resistance_increase=get_zeros(12, 25), + core_diameter=1, + conductor_diameter=shapes.get_zeros() + 2, + outer_layer_strand_diameter=shapes.get_zeros(), + emissivity=shapes.get_zeros(), + solar_absorptivity=shapes.get_zeros(), + temperature1=shapes.get_zeros(), + temperature2=shapes.get_zeros() + 1, + resistance_at_temperature1=shapes.get_zeros(), + resistance_at_temperature2=shapes.get_zeros(), + aluminium_cross_section_area=1 + shapes.get_zeros(), + constant_magnetic_effect=shapes.get_zeros(), + current_density_proportional_magnetic_effect=shapes.get_zeros() + 1, + max_magnetic_core_relative_resistance_increase=shapes.get_zeros(), ) weather = Weather( - air_temperature=get_zeros(13, 25), - wind_direction=get_zeros(14, 25), - wind_speed=get_zeros(15, 25), - ground_albedo=get_zeros(23, 25), - clearness_ratio=get_zeros(16, 25), + air_temperature=shapes.get_zeros(), + wind_direction=shapes.get_zeros(), + wind_speed=shapes.get_zeros(), + ground_albedo=shapes.get_zeros(), + clearness_ratio=shapes.get_zeros(), ) span = Span( conductor=conductor, start_tower=Tower( - get_zeros(17, 25), - get_zeros(18, 25), - get_zeros(19, 25), + shapes.get_zeros(), + shapes.get_zeros(), + shapes.get_zeros(), ), end_tower=Tower( - 0.1 + get_zeros(20, 25), - 0.1 + get_zeros(21, 25), - 10 + get_zeros(22, 25), + 0.1 + shapes.get_zeros(), + 0.1 + shapes.get_zeros(), + 10 + shapes.get_zeros(), ), num_conductors=1, ) - time = np.array([np.datetime64("2022-01-01")]).reshape([1] * 25) - if 24 in vectorization_indices: - time = np.concatenate([time, time], axis=-1) + time = shapes.get_zeros(dtype=np.dtype("datetime64[D]")) model = Cigre601(span=span, weather=weather, time=time) temperature = model.compute_conductor_temperature(100) - shape = [1] * 25 - for vectorisation_index in vectorization_indices: - shape[vectorisation_index] = 2 - shape[0] = 1 - assert isinstance(temperature, np.ndarray) - assert temperature.shape == tuple(shape) + assert temperature.shape == array_shapes.result_shape -@hypothesis.given(array_shapes=mutually_broadcastable_shapes(num_shapes=NUMBER_OF_VARIABLES)) +@hypothesis.given(array_shapes=mutually_broadcastable_shapes(num_shapes=NUMBER_OF_VARIABLES_TRANSIENT, max_dims=5)) def test_transient_vectorization(array_shapes: BroadcastableShapes): shapes = Shapes(array_shapes.input_shapes)