opusonesolutions/carsons

View on GitHub
carsons/carsons.py

Summary

Maintainability
A
3 hrs
Test Coverage
from collections import defaultdict
from itertools import islice
from typing import Dict, Iterable, Iterator, Tuple

from numpy import arctan, array, cos, exp, log, ndarray
from numpy import pi as π
from numpy import sin, sqrt, zeros
from numpy.linalg import inv

alpha = exp(2j * π / 3)

# fmt: off
A = array([
    [1,            1,          1],  # noqa E202
    [1,     alpha**2,      alpha],  # noqa E202
    [1,        alpha,   alpha**2],  # noqa E202
])

Ainv = (1 / 3) * array([
    [1,            1,          1],  # noqa E202
    [1,        alpha,   alpha**2],  # noqa E202
    [1,     alpha**2,      alpha],  # noqa E202
])
# fmt: on


def convert_geometric_model(geometric_model) -> ndarray:
    carsons_model = CarsonsEquations(geometric_model)

    z_primitive = carsons_model.build_z_primitive()
    z_abc = perform_kron_reduction(z_primitive)
    return z_abc


def calculate_impedance(model) -> ndarray:
    z_primitive = model.build_z_primitive()
    z_abc = perform_kron_reduction(z_primitive, dimension=model.dimension)

    return z_abc


def perform_kron_reduction(z_primitive: ndarray, dimension=3) -> ndarray:
    """Reduces the primitive impedance matrix to an equivalent impedance
    matrix.

    We break z_primative up into four quadrants as follows:

              Ẑpp = [Ẑaa, Ẑab, Ẑac]   Ẑpn = [Ẑan]
                    [Ẑba, Ẑbb, Ẑbc]         [Ẑbn]
                    [Ẑca, Ẑcb, Ẑcc]         [Ẑcn]

              Ẑnp = [Ẑna, Ẑnb, Ẑnc]   Ẑnn = [Ẑnn]

    Ẑnn is of dimension mxm, where m is the number of neutrals. E.g. with m = 2:
                                          Ẑan = [Ẑan₁,  Ẑan₂]
                                                [Ẑbn₁,  Ẑbn₂]
                                                [Ẑcn₁,  Ẑcn₂]

              Ẑna = [Ẑn₁a, Ẑn₁b, Ẑn₁c]    Ẑnn = [Ẑn₁n₁, Ẑn₁n₂]
                    [Ẑn₂a, Ẑn₂b, Ẑn₂c]          [Ẑn₂n₁, Ẑn₂n₂]

    Definitions:
    Ẑ ----- "primative" impedance value, i.e. one that does not factor
            in the mutuals caused by neighboring neutral conductors.
    Z ----- a phase-phase impedance value that factors the mutual impedance
            of neighboring neutral conductors

    Returns:
    Z ----  a corrected impedance matrix in the form:

                     Zabc = [Zaa, Zab, Zac]
                            [Zba, Zbb, Zbc]
                            [Zca, Zcb, Zcc]
    """
    Ẑpp, Ẑpn = (
        z_primitive[0:dimension, 0:dimension],
        z_primitive[0:dimension, dimension:],
    )
    Ẑnp, Ẑnn = (
        z_primitive[dimension:, 0:dimension],
        z_primitive[dimension:, dimension:],
    )
    Z_abc = Ẑpp - Ẑpn @ inv(Ẑnn) @ Ẑnp
    return Z_abc


def calculate_sequence_impedance_matrix(Z):
    return Ainv @ Z @ A


def calculate_sequence_impedances(Z):
    Z012 = calculate_sequence_impedance_matrix(Z)
    return Z012[1, 1], Z012[0, 0]


class CarsonsEquations:

    ρ = 100  # resistivity, ohms/meter^3
    μ = 4 * π * 1e-7  # permeability, Henry / meter

    def __init__(self, model):
        self.phases: Iterable[str] = model.phases
        self.phase_positions: Dict[str, Tuple[float, float]] = model.wire_positions
        self.gmr: Dict[str, float] = model.geometric_mean_radius
        self.r: Dict[str, float] = model.resistance

        self.ƒ = getattr(model, "frequency", 60)
        self.ω = 2.0 * π * self.ƒ  # angular frequency radians / second

    def build_z_primitive(self) -> ndarray:
        dimension = len(self.conductors)
        z_primitive = zeros(shape=(dimension, dimension), dtype=complex)

        for index_i, phase_i in enumerate(self.conductors):
            for index_j, phase_j in enumerate(self.conductors):
                if phase_i not in self.phases or phase_j not in self.phases:
                    continue
                R = self.compute_R(phase_i, phase_j)
                X = self.compute_X(phase_i, phase_j)
                z_primitive[index_i, index_j] = complex(R, X)

        return z_primitive

    def compute_R(self, i, j) -> float:
        rᵢ = self.r[i]
        ΔR = self.μ * self.ω / π * self.compute_P(i, j)

        if i == j:
            return rᵢ + ΔR
        else:
            return ΔR

    def compute_X(self, i, j) -> float:
        Qᵢⱼ = self.compute_Q(i, j)
        ΔX = self.μ * self.ω / π * Qᵢⱼ

        # calculate geometry ratio 𝛥G
        if i != j:
            Dᵢⱼ = self.compute_D(i, j)
            dᵢⱼ = self.compute_d(i, j)
            𝛥G = Dᵢⱼ / dᵢⱼ
        else:
            hᵢ = self.get_h(i)
            gmrⱼ = self.gmr[j]
            𝛥G = 2.0 * hᵢ / gmrⱼ

        X_o = self.ω * self.μ / (2 * π) * log(𝛥G)

        return X_o + ΔX

    def compute_P(self, i, j, number_of_terms=1) -> float:
        terms = islice(self.compute_P_terms(i, j), number_of_terms)
        return sum(terms)

    def compute_P_terms(self, i, j) -> Iterator[float]:
        yield π / 8.0

        kᵢⱼ = self.compute_k(i, j)
        θᵢⱼ = self.compute_θ(i, j)

        yield -kᵢⱼ / (3 * sqrt(2)) * cos(θᵢⱼ)
        yield kᵢⱼ**2 / 16 * (0.6728 + log(2 / kᵢⱼ)) * cos(2 * θᵢⱼ)
        yield kᵢⱼ**2 / 16 * θᵢⱼ * sin(2 * θᵢⱼ)
        yield kᵢⱼ**3 / (45 * sqrt(2)) * cos(3 * θᵢⱼ)
        yield -π * kᵢⱼ**4 * cos(4 * θᵢⱼ) / 1536

    def compute_Q(self, i, j, number_of_terms=2) -> float:
        terms = islice(self.compute_Q_terms(i, j), number_of_terms)
        return sum(terms)

    def compute_Q_terms(self, i, j) -> Iterator[float]:
        yield -0.0386

        kᵢⱼ = self.compute_k(i, j)
        yield 0.5 * log(2 / kᵢⱼ)

        θᵢⱼ = self.compute_θ(i, j)
        yield kᵢⱼ / (3 * sqrt(2)) * cos(θᵢⱼ)
        yield -π * kᵢⱼ**2 / 64 * cos(2 * θᵢⱼ)
        yield kᵢⱼ**3 / (45 * sqrt(2)) * cos(3 * θᵢⱼ)
        yield -(kᵢⱼ**4) / 384 * θᵢⱼ * sin(4 * θᵢⱼ)
        yield -(kᵢⱼ**4) / 384 * cos(4 * θᵢⱼ) * (log(2 / kᵢⱼ) + 1.0895)

    def compute_k(self, i, j) -> float:
        Dᵢⱼ = self.compute_D(i, j)
        return Dᵢⱼ * sqrt(self.ω * self.μ / self.ρ)

    def compute_θ(self, i, j) -> float:
        xᵢ, _ = self.phase_positions[i]
        xⱼ, _ = self.phase_positions[j]
        xᵢⱼ = abs(xⱼ - xᵢ)
        hᵢ, hⱼ = self.get_h(i), self.get_h(j)

        return arctan(xᵢⱼ / (hᵢ + hⱼ))

    def compute_d(self, i, j) -> float:
        return self.calculate_distance(
            self.phase_positions[i],
            self.phase_positions[j],
        )

    def compute_D(self, i, j) -> float:
        xⱼ, yⱼ = self.phase_positions[j]

        return self.calculate_distance(self.phase_positions[i], (xⱼ, -yⱼ))

    @staticmethod
    def calculate_distance(positionᵢ, positionⱼ) -> float:
        xᵢ, yᵢ = positionᵢ
        xⱼ, yⱼ = positionⱼ
        return sqrt((xᵢ - xⱼ) ** 2 + (yᵢ - yⱼ) ** 2)

    def get_h(self, i):
        _, yᵢ = self.phase_positions[i]
        return yᵢ

    @property
    def dimension(self):
        return 2 if getattr(self, "is_secondary", False) else 3

    @property
    def conductors(self):
        neutral_conductors = sorted([ph for ph in self.phases if ph.startswith("N")])

        return ["A", "B", "C"] + neutral_conductors


class ModifiedCarsonsEquations(CarsonsEquations):
    """
    Modified Carson's Equation. Two approximations are made:
    only the first term of P and the first two terms of Q are considered.
    """

    number_of_P_terms = 1

    def compute_P(self, i, j, number_of_terms=1) -> float:
        return super().compute_P(i, j, self.number_of_P_terms)

    def compute_X(self, i, j) -> float:
        Q_first_term = super().compute_Q(i, j, 1)

        # Simplify equations and don't compute Dᵢⱼ explicitly
        kᵢⱼ_Dᵢⱼ_ratio = sqrt(self.ω * self.μ / self.ρ)
        ΔX = Q_first_term * 2 + log(2)

        if i == j:
            X_o = -log(self.gmr[i]) - log(kᵢⱼ_Dᵢⱼ_ratio)
        else:
            X_o = -log(self.compute_d(i, j)) - log(kᵢⱼ_Dᵢⱼ_ratio)

        return (X_o + ΔX) * self.ω * self.μ / (2 * π)


class ConcentricNeutralCarsonsEquations(ModifiedCarsonsEquations):
    def __init__(self, model, *args, **kwargs):
        super().__init__(model)
        self.neutral_strand_gmr: Dict[str, float] = model.neutral_strand_gmr
        self.neutral_strand_count: Dict[str, float] = defaultdict(
            lambda: None, model.neutral_strand_count
        )
        self.neutral_strand_resistance: Dict[
            str, float
        ] = model.neutral_strand_resistance
        # fmt: off
        self.radius: Dict[str, float] = defaultdict(
            lambda: None,
            {
                phase: (diameter_over_neutral - model.neutral_strand_diameter[phase]) / 2
                for phase, diameter_over_neutral in model.diameter_over_neutral.items()
            },
        )
        self.phase_positions.update({
            f"N{phase}": self.phase_positions[phase]
            for phase in self.phase_positions.keys()
        })
        self.gmr.update({
            phase: self.GMR_cn(phase)
            for phase in model.diameter_over_neutral.keys()
        })
        self.r.update({
            phase: resistance / model.neutral_strand_count[phase]
            for phase, resistance in model.neutral_strand_resistance.items()
        })
        # fmt: on
        return

    def compute_d(self, i, j) -> float:
        I, J = set(i), set(j)
        r = self.radius[i] or self.radius[j]

        one_neutral_same_phase = I ^ J == set("N")
        different_phase = not I & J
        one_neutral = "N" in I ^ J

        if one_neutral_same_phase:
            # Distance between a neutral/phase conductor of same phase
            return r

        distance_ij = self.calculate_distance(
            self.phase_positions[i], self.phase_positions[j]
        )
        if different_phase and one_neutral:
            # Distance between a neutral/phase conductor of different phase
            # approximate by modelling the concentric neutral cables as one
            # equivalent conductor directly above the phase conductor
            return (distance_ij**2 + r**2) ** 0.5
        else:
            # Distance between two neutral/phase conductors
            return distance_ij

    def GMR_cn(self, phase) -> float:
        GMR_s = self.neutral_strand_gmr[phase]
        k = self.neutral_strand_count[phase]
        R = self.radius[phase]
        return (GMR_s * k * R ** (k - 1)) ** (1 / k)


class TapeShieldedCableCarsonsEquations(ModifiedCarsonsEquations):

    ρ_tape_shield = 1.7721e-8  # copper resistivity at 20 degrees, ohm-meter

    def __init__(self, model, *args, **kwargs):
        super().__init__(model)
        self.ds = model.tape_shield_outer_diameter
        self.thickness = model.tape_shield_thickness

        self.phases.extend(f"{ph}t" for ph in model.tape_shield_outer_diameter.keys())
        # fmt: off
        self.r.update({
            f"{ph}t": self.compute_shield_r(self.ds[ph], self.thickness[ph])
            for ph in model.tape_shield_outer_diameter.keys()
        })
        self.gmr.update({
            f"{ph}t": self.compute_shield_gmr(self.ds[ph], self.thickness[ph])
            for ph in model.tape_shield_outer_diameter.keys()
        })
        # fmt: on

        # set shield position to be the same as its phase conductor position
        tape_shield_positions = {
            f"{ph}t": self.phase_positions[ph]
            for ph in model.tape_shield_outer_diameter.keys()
        }
        self.phase_positions.update(tape_shield_positions)

    @staticmethod
    def compute_shield_gmr(ds, thickness) -> float:
        return (ds - thickness) / 2

    @classmethod
    def compute_shield_r(cls, ds, thickness) -> float:
        area = (π * (ds / 2) ** 2) - (π * (ds / 2 - thickness) ** 2)
        return cls.ρ_tape_shield / area

    def compute_d(self, i, j) -> float:
        I, J = set(i), set(j)

        one_tape_shield_same_phase = I ^ J == set("t")

        if one_tape_shield_same_phase:
            return self.gmr[i] if "t" in i else self.gmr[j]
        else:
            return super().compute_d(i, j)

    @property
    def conductors(self):
        neutral_conductors = sorted([ph for ph in self.phases if ph.startswith("N")])
        shield_conductors = sorted([ph for ph in self.phases if ph.endswith("t")])

        return ["A", "B", "C"] + shield_conductors + neutral_conductors


class MultiConductorCarsonsEquations(ModifiedCarsonsEquations):
    def __init__(self, model):
        super().__init__(model)
        self.outside_radius: Dict[str, float] = model.outside_radius

    def compute_d(self, i, j) -> float:
        # Assumptions:
        # 1. All conductors in the cable are touching each other and
        #    therefore equidistant.
        # 2. In case of quadruplex cables, the space between conductors
        #    which are diagonally positioned is neglected.
        return self.outside_radius[i] + self.outside_radius[j]

    @property
    def conductors(self):
        neutral_conductors = sorted([ph for ph in self.phases if ph.startswith("N")])
        if self.is_secondary:
            conductors = ["S1", "S2"] + neutral_conductors
        else:
            conductors = ["A", "B", "C"] + neutral_conductors

        return conductors

    @property
    def is_secondary(self):
        phase_conductors = [ph for ph in self.phases if not ph.startswith("N")]
        if phase_conductors == ["S1", "S2"]:
            return True
        else:
            return False