From 545622a486355ca8f8fdae8d3580ee5d2a7fbd47 Mon Sep 17 00:00:00 2001 From: marcostfermin <90121492+marcostfermin@users.noreply.github.com> Date: Sun, 11 Jan 2026 01:31:15 +0000 Subject: [PATCH] Refactor simulation, computation, fault analysis, and machine modules --- electricpy/compute.py | 37 +- electricpy/fault.py | 1225 ++-------------------------------------- electricpy/machines.py | 1079 +++++++---------------------------- electricpy/sim.py | 323 +++-------- 4 files changed, 351 insertions(+), 2313 deletions(-) diff --git a/electricpy/compute.py b/electricpy/compute.py index b2549535..37747cb6 100644 --- a/electricpy/compute.py +++ b/electricpy/compute.py @@ -54,6 +54,11 @@ def largest_integer(numBits, signed=True): >>> cmp.largest_integer(32, signed=True) 2147483647 """ + if numBits is None: + raise ValueError("numBits must be a positive integer.") + numBits = int(numBits) + if numBits <= 0: + raise ValueError("numBits must be a positive integer.") # Use Signed or Unsigned Formula if signed: return int(2 ** (numBits - 1) - 1) @@ -123,7 +128,7 @@ def mod2div(divident, divisor): # part used in each step) is 0, the step cannot # use the regular divisor; we need to use an # all-0s divisor. - tmp = xor('0' * pick, tmp) + divident[pick] + tmp = xor('0' * len(divisor), tmp) + divident[pick] # increment pick to move further pick += 1 @@ -134,7 +139,7 @@ def mod2div(divident, divisor): if tmp[0] == '1': tmp = xor(divisor, tmp) else: - tmp = xor('0' * pick, tmp) + tmp = xor('0' * len(divisor), tmp) checkword = tmp return checkword @@ -143,6 +148,17 @@ def mod2div(divident, divisor): data = str(data) # Condition Key key = str(key) + + # Basic validation + if len(key) < 2: + raise ValueError("CRC key must be at least 2 bits long.") + if set(data) - set('01'): + raise ValueError("CRC data must be a string of bits containing only '0' and '1'.") + if set(key) - set('01'): + raise ValueError("CRC key must be a string of bits containing only '0' and '1'.") + if key[0] != '1': + raise ValueError("CRC key must start with '1' (highest-order term present).") + l_key = len(key) # Appends n-1 zeroes at end of data @@ -216,7 +232,7 @@ def mod2div(divident, divisor): # part used in each step) is 0, the step cannot # use the regular divisor; we need to use an # all-0s divisor. - tmp = xor('0' * pick, tmp) + divident[pick] + tmp = xor('0' * len(divisor), tmp) + divident[pick] # increment pick to move further pick += 1 @@ -227,7 +243,7 @@ def mod2div(divident, divisor): if tmp[0] == '1': tmp = xor(divisor, tmp) else: - tmp = xor('0' * pick, tmp) + tmp = xor('0' * len(divisor), tmp) checkword = tmp return checkword @@ -236,6 +252,17 @@ def mod2div(divident, divisor): data = str(data) # Condition Key key = str(key) + + # Basic validation + if len(key) < 2: + raise ValueError("CRC key must be at least 2 bits long.") + if set(data) - set('01'): + raise ValueError("CRC data must be a string of bits containing only '0' and '1'.") + if set(key) - set('01'): + raise ValueError("CRC key must be a string of bits containing only '0' and '1'.") + if key[0] != '1': + raise ValueError("CRC key must start with '1' (highest-order term present).") + l_key = len(key) # Appends n-1 zeroes at end of data @@ -264,7 +291,7 @@ def string_to_bits(str): The binary representation of the input string. """ - data = (''.join(format(ord(x), 'b') for x in str)) + data = (''.join(format(ord(x), '08b') for x in str)) return data # END diff --git a/electricpy/fault.py b/electricpy/fault.py index db20b0f1..afba4f85 100644 --- a/electricpy/fault.py +++ b/electricpy/fault.py @@ -32,26 +32,6 @@ def single_phase_to_ground_fault(Vth, Zseq, Rf=0, sequence=True, reference='A'): .. math:: I_2 = I_1 .. math:: I_0 = I_1 - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - Rf: complex, optional - The fault resistance, default=0 - sequence: bool, optional - Control argument to force return into symmetrical- - or phase-domain values. - reference: {'A', 'B', 'C'} - Single character denoting the reference, - default='A' - - Returns - ------- - Ifault: list of complex, - The Array of Fault Currents as (If0, If1, If2) """ # Decompose Reactance Tuple X0, X1, X2 = Zseq @@ -82,32 +62,6 @@ def double_phase_to_ground_fault(Vth, Zseq, Rf=0, sequence=True, reference='A'): This function will evaluate the Zero, Positive, and Negative sequence currents for a double-line-to-ground fault. - - .. math:: I_1 = \frac{V_{th}}{Z_1+\frac{Z_2*(Z_0+3*R_f)}{Z_0+Z_2+3*R_f}} - - .. math:: I_2 = -\frac{V_{th}-Z_1*I_1}{X_2} - - .. math:: I_0 = -\frac{V_{th}-Z_1*I_1}{X_0+3*R_f} - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - Rf: complex, optional - The fault resistance, default=0 - sequence: bool, optional - Control argument to force return into symmetrical- - or phase-domain values. - reference: {'A', 'B', 'C'} - Single character denoting the reference, - default='A' - - Returns - ------- - Ifault: list of complex, - The Array of Fault Currents as (If0, If1, If2) """ # Decompose Reactance Tuple X0, X1, X2 = Zseq @@ -139,32 +93,6 @@ def phase_to_phase_fault(Vth, Zseq, Rf=0, sequence=True, reference='A'): This function will evaluate the Zero, Positive, and Negative sequence currents for a phase-to-phase fault. - - .. math:: I_1 = \frac{V_{th}}{Z_1+Z_2+R_f} - - .. math:: I_2 = -I_1 - - .. math:: I_0 = 0 - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - Rf: complex, optional - The fault resistance, default=0 - sequence: bool, optional - Control argument to force return into symmetrical- - or phase-domain values. - reference: {'A', 'B', 'C'} - Single character denoting the reference, - default='A' - - Returns - ------- - Ifault: list of complex, - The Array of Fault Currents as (If0, If1, If2) """ # Decompose Reactance Tuple X0, X1, X2 = Zseq @@ -196,32 +124,6 @@ def three_phase_fault(Vth, Zseq, Rf=0, sequence=True, reference='A'): This function will evaluate the Zero, Positive, and Negative sequence currents for a three-phase fault. - - .. math:: I_1 = \frac{V_{th}}{Z_1+R_1} - - .. math:: I_2 = 0 - - .. math:: I_0 = 0 - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - Rf: complex, optional - The fault resistance, default=0 - sequence: bool, optional - Control argument to force return into symmetrical- - or phase-domain values. - reference: {'A', 'B', 'C'} - Single character denoting the reference, - default='A' - - Returns - ------- - Ifault: list of complex - The Fault Current, equal for 0, pos., and neg. seq. """ # Decompose Reactance Tuple X0, X1, X2 = Zseq @@ -247,30 +149,6 @@ def poleopen1(Vth, Zseq, sequence=True, reference='A'): This function will evaluate the Zero, Positive, and Negative sequence currents for a single pole open fault. - - .. math:: I_1 = \frac{V_{th}}{Z_1+(\frac{1}{Z_2}+\frac{1}{Z_0})^-1} - - .. math:: I_2 = -I_1 * \frac{Z_0}{Z_2+Z_0} - - .. math:: I_0 = -I_1 * \frac{Z_2}{Z_2+Z_0} - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - sequence: bool, optional - Control argument to force return into symmetrical- - or phase-domain values. - reference: {'A', 'B', 'C'} - Single character denoting the reference, or the - faulted phase indicator; default='A' - - Returns - ------- - Ifault: list of complex, - The Array of Fault Currents as (If0, If1, If2) """ # Decompose Reactance Tuple X0, X1, X2 = Zseq @@ -295,34 +173,10 @@ def poleopen1(Vth, Zseq, sequence=True, reference='A'): # Define Double Pole Open Calculator def poleopen2(Vth, Zseq, sequence=True, reference='A'): r""" - Single Pole Open Fault Calculator. + Double Pole Open Fault Calculator. This function will evaluate the Zero, Positive, and Negative - sequence currents for a single pole open fault. - - .. math:: I_1 = \frac{V_{th}}{Z_1+Z_2+Z_0} - - .. math:: I_2 = I_1 - - .. math:: I_0 = I_1 - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - sequence: bool, optional - Control argument to force return into symmetrical- - or phase-domain values. - reference: {'A', 'B', 'C'} - Single character denoting the reference, or the - faulted phase indicator; default='A' - - Returns - ------- - Ifault: list of complex, - The Array of Fault Currents as (If0, If1, If2) + sequence currents for a double pole open fault. """ # Decompose Reactance Tuple X0, X1, X2 = Zseq @@ -348,40 +202,9 @@ def poleopen2(Vth, Zseq, sequence=True, reference='A'): def short_circuit_mva(Zth=None, Isc=None, Vth=1): r""" Short-Circuit MVA Calculator. - - Function defines a method of interpretively - calculating the short-circuit MVA value - given two of the three arguments. The formulas - are all based around the following: - - .. math:: MVA_{sc} = V_{th}*I_{sc} - - .. math:: V_{th} = I_{sc}*Z_{th} - - Parameters - ---------- - Zth: float - The Thevenin-Equivalent-Impedance - Isc: float, optional - Short-Circuit-Current, if left as - None, will force function to use - default setting for Vth. - default=None - Vth: float, optional - The Thevenin-Equivalent-Voltage, - defaults to a 1-per-unit value. - default=1 - - Returns - ------- - MVA: float - Short-Circuit MVA, not described - as three-phase or otherwise, such - determination is dependent upon - inputs. """ # Test for too few inputs - if not any((Zth, Isc)): + if Zth is None and Isc is None: raise ValueError("Either Zth or Isc must be specified.") # Condition Inputs if Zth is not None: @@ -391,9 +214,9 @@ def short_circuit_mva(Zth=None, Isc=None, Vth=1): if Vth != 1: Vth = abs(Vth) # Calculate MVA from one of the available methods - if all((Zth, Isc)): + if Zth is not None and Isc is not None: MVA = Isc ** 2 * Zth - elif all((Zth, Vth)): + elif Zth is not None: MVA = Vth ** 2 / Zth else: MVA = Vth * Isc @@ -408,27 +231,6 @@ def short_circuit_mva(Zth=None, Isc=None, Vth=1): def phs3mvasc(Vth, Zseq, Rf=0, Sbase=1): r""" Three-Phase MVA Short-Circuit Calculator. - - Calculator to evaluate the Short-Circuit MVA of a three-phase fault given the system - parameters of Vth, Zseq, and an optional Rf. Uses the formula as follows: - - .. math:: MVA_{sc} = \frac{\left|V_{th}^2\right|}{|Z_1|} * Sbase - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - Rf: complex, optional - The fault resistance, default=0 - Sbase: real, optional - The per-unit base for power. default=1 - - Returns - ------- - MVA: real - Three-Phase Short-Circuit MVA. """ # Calculate Three-Phase MVA MVA = abs(Vth) ** 2 / abs(Zseq[1]) * Sbase @@ -443,31 +245,6 @@ def phs3mvasc(Vth, Zseq, Rf=0, Sbase=1): def phs1mvasc(Vth, Zseq, Rf=0, Sbase=1): r""" Single-Phase MVA Short-Circuit Calculator. - - Calculator to evaluate the Short-Circuit MVA of a single-phase fault given the system - parameters of Vth, Zseq, and an optional Rf. Uses the formula as follows: - - .. math:: MVA_{sc} = \left|I_1^2\right|*|Z_1| * Sbase - - where: - - .. math:: I_1 = \frac{V_{th}}{Z_0+Z_1+Z_2+3*R_f} - - Parameters - ---------- - Vth: complex - The Thevenin-Equivalent-Voltage - Zseq: list of complex - Tupple of sequence reactances as (Z0, Z1, Z2) - Rf: complex, optional - The fault resistance, default=0 - Sbase: real, optional - The per-unit base for power. default=1 - - Returns - ------- - MVA: real - Single-Phase Short-Circuit MVA. """ # Decompose Reactance Tuple X0, X1, X2 = Zseq @@ -493,38 +270,6 @@ def phs1mvasc(Vth, Zseq, Rf=0, Sbase=1): def busvolt(k, n, Vpf, Z0, Z1, Z2, If, sequence=True, reference='A'): """ Faulted Bus Voltage Calculator. - - This function is designed to calculate the bus voltage(s) - given a specific set of fault characteristics. - - Parameters - ---------- - k: float - Bus index at which to calculate faulted voltage - n: float - Bus index at which fault occurred - Vpf: complex - Voltage Pre-Fault, Singular Number - Z0: ndarray - Zero-Sequence Impedance Matrix - Z1: ndarray - Positive-Sequence Impedance Matrix - Z2: ndarray - Negative-Sequence Impedance Matrix - If: complex - Sequence Fault Current Evaluated at Bus *n* - sequence: bool, optional - Control argument to force return into symmetrical- - or phase-domain values. - reference: {'A', 'B', 'C'} - Single character denoting the reference, - default='A' - - Returns - ------- - Vf: complex - The Fault Voltage, set of sequence or phase voltages as - specified by *sequence* """ # Condition Inputs k = k - 1 @@ -550,49 +295,6 @@ def ct_saturation(XoR, Imag, Vrated, Irated, CTR, Rb, Xb, remnance=0, freq=60, ALF=20): r""" Electrical Current Transformer Saturation Calculator. - - A function to determine the saturation value and a boolean indicator - showing whether or not CT is -in fact- saturated. - - To perform this evaluation, we must satisfy the equation: - - .. math:: - 20\geq(1+\frac{X}{R})*\frac{|I_{mag}|}{I_{rated}*CTR} - *\frac{\left|R_{burden}+j*\omega*\frac{X_{burden}} - {\omega}\right|*100}{V_{rated}*(1-remnanc)} - - Parameters - ---------- - XoR: float - The X-over-R ratio of the system. - Imag: float - The (maximum) current magnitude to use for calculation, - typically the fault current. - Vrated: float - The rated voltage (accompanying the C-Class value) of - the CT. - Irated: float - The rated secondary current for the CT. - CTR: float - The CT Ratio (primary/secondary, N) to be used. - Rb: float - The total burden resistance in ohms. - Xb: float - The total burden reactance in ohms. - remnance: float, optional - The system flux remnance, default=0. - freq: float, optional - The system frequency in Hz, default=60. - ALF: float, optional - The Saturation Constant which must be satisfied, - default=20. - - Returns - ------- - result: float - The calculated Saturation value. - saturation: bool - Boolean indicator to mark presence of saturation. """ # Define omega w = 2 * _np.pi * freq @@ -616,55 +318,6 @@ def ct_saturation(XoR, Imag, Vrated, Irated, CTR, Rb, Xb, remnance=0, freq=60, def ct_cclass(XoR, Imag, Irated, CTR, Rb, Xb, remnance=0, sat_crit=20): r""" Electrical Current Transformer (CT) C-Class Function. - - A function to determine the C-Class rated voltage for a CT. - The formula shown below demonstrates the standard formula - which is normally used to evaluate the saturation criteria. - Worth noting here, is the fact that :math:`V_{rated}` is the - CT C-Class. - - .. math:: - \text{Saturation Criteria}=\frac{(1+\frac{X}{R})\cdot - \frac{|I_{mag}|}{I_{rated}\cdot CTR}\cdot\frac{\left| - R_{burden}+j\cdot X_{burden}\right|\cdot100}{V_{rated}}} - {1-remnance} - - For the purposes of this function, the above formula is applied - as follows to evaluate the CT C-Class such as to satisfy the - saturation criteria defined. - - .. math:: - \text{CT C-Class}=\frac{(1+\frac{X}{R})\cdot - \frac{|I_{mag}|}{I_{rated}\cdot CTR}\cdot\frac{ - \left|R_{burden}+j\cdot X_{burden}\right|\cdot100} - {\text{Saturation Criteria (i.e., 20)}}}{1-remnance} - - Parameters - ---------- - XoR: float - The X-over-R ratio of the system. - Imag: float - The (maximum) current magnitude to use for calculation, - typically the fault current. - Irated: float - The rated secondary current for the CT. - CTR: float - The CT Ratio (primary/secondary, N) to be used. - Rb: float - The total burden resistance in ohms. - Xb: float - The total burden reactance in ohms. - remnance: float, optional - The system flux remnance, default=0. - sat_crit: float, optional - The saturation criteria which must be satisfied, - typically such that CT saturation will not occur, - default=20. - - Returns - ------- - c_class: float - The calculated C-Class rated voltage. """ # Calculate each "term" (multiple) t1 = (1 + XoR) @@ -678,33 +331,9 @@ def ct_cclass(XoR, Imag, Irated, CTR, Rb, Xb, remnance=0, sat_crit=20): # Define Saturation Voltage at Rated Burden -def ct_satratburden(Inom, VArat=None, ANSIv=None, ALF=20, ): +def ct_satratburden(Inom, VArat=None, ANSIv=None, ALF=20): r""" Electrical Current Transformer (CT) Saturation at Rated Burden Calculator. - - A function to determine the Saturation at rated burden. - - .. math:: V_{saturated}=ALF*\frac{VA_{rated}}{I_{nominal}} - - where: - - .. math:: VA_{rated}=I_{nominal}*\frac{ANSI_{voltage}}{20} - - Parameters - ---------- - Inom: float - Nominal Current - VArat: float, optional, exclusive - The apparent power (VA) rating of the CT. - ANSIv: float, optional, exclusive - The ANSI voltage requirement to meet. - ALF: float, optional - Accuracy Limit Factor, default=20. - - Returns - ------- - Vsat: float - The saturated voltage. """ # Validate Inputs if VArat is None and ANSIv is None: @@ -721,24 +350,6 @@ def ct_satratburden(Inom, VArat=None, ANSIv=None, ALF=20, ): def ct_vpeak(Zb, Ip, CTR): r""" Electrical Current Transformer (CT) Peak Voltage Calculator. - - Simple formula to calculate the Peak Voltage of a CT. - - .. math:: \sqrt{3.5*|Z_burden|*I_{peak}*CTR} - - Parameters - ---------- - Zb: float - The burden impedance magnitude (in ohms). - Ip: float - The peak current for the CT. - CTR: float - The CTR turns ratio of the CT. - - Returns - ------- - Vpeak: float - The peak voltage. """ return _np.sqrt(3.5 * abs(Zb) * Ip * CTR) @@ -748,34 +359,6 @@ def ct_timetosat(Vknee, XoR, Rb, CTR, Imax, ts=None, npts=100, freq=60, plot=False): r""" Electrical Current Transformer (CT) Time to Saturation Function. - - Function to determine the "time to saturate" for an underrated C-Class - CT using three standard curves described by Juergen Holbach. - - Parameters - ---------- - Vknee: float - The knee-voltage for the CT. - XoR: float - The X-over-R ratio of the system. - Rb: float - The total burden resistance in ohms. - CTR: float - The CT Ratio (primary/secondary, N) to be used. - Imax: float - The (maximum) current magnitude to use for calculation, - typically the fault current. - ts: numpy.ndarray or float, optional - The time-array or particular (floatint point) time at which - to calculate the values. default=_np.linspace(0,0.1,freq*npts) - npts: float, optional - The number of points (per cycle) to calculate if ts is not - specified, default=100. - freq: float, optional - The system frequency in Hz, default=60. - plot: bool, optional - Control argument to enable plotting of calculated curves, - default=False. """ # Calculate omega w = 2 * _np.pi * freq @@ -822,41 +405,10 @@ def ct_timetosat(Vknee, XoR, Rb, CTR, Imax, ts=None, npts=100, freq=60, def pktransrecvolt(C, L, R=0, VLL=None, VLN=None, freq=60): """ Peak Transient Recovery Function. - - Peak Transient Recovery Voltage calculation function, evaluates the peak - transient recovery voltage (restriking voltage) and the - Rate-of-Rise-Recovery Voltage. - - Parameters - ---------- - C: float - Capacitance Value in Farads. - L: float - Inductance in Henries. - R: float, optional - The resistance of the system used for - calculation, default=0. - VLL: float, exclusive - Line-to-Line voltage, exclusive - optional argument. - VLN: float, exclusive - Line-to-Neutral voltage, exclusive - optional argument. - freq: float, optional - System frequency in Hz. - - Returns - ------- - Vcpk: float - Peak Transient Recovery Voltage in volts. - RRRV: float - The RRRV (Rate-of-Rise-Recovery Voltage) - calculated given the parameters in volts - per second. """ # Evaluate alpha, omega-n, and fn alpha = R / (2 * L) - wn = 1 / _np.sqrt(L * C) - alpha + wn = _np.sqrt(max(0.0, 1 / (L * C) - alpha ** 2)) fn = wn / (2 * _np.pi) # Evaluate Vm if VLL is not None: @@ -866,7 +418,11 @@ def pktransrecvolt(C, L, R=0, VLL=None, VLN=None, freq=60): else: raise ValueError("One voltage must be specified.") # Evaluate Vcpk (worst case) - Vcpk = wn ** 2 / (wn ** 2 - 2 * _np.pi * freq) * Vm * 2 + we = 2 * _np.pi * freq + denom = (wn ** 2 - we ** 2) + if denom == 0: + raise ZeroDivisionError("Invalid parameters: wn equals system frequency (resonance).") + Vcpk = (wn ** 2 / denom) * Vm * 2 # Evaluate RRRV RRRV = 2 * Vm * fn / 0.5 return Vcpk, RRRV @@ -876,33 +432,6 @@ def pktransrecvolt(C, L, R=0, VLL=None, VLN=None, freq=60): def trvresistor(C, L, reduction, Rd0=500, wd0=260e3, tpk0=10e-6): """ Transient Recovery Voltage (TRV) Reduction Resistor Function. - - Function to find the resistor value that will reduce the TRV by a specified - percentage. - - Parameters - ---------- - C: float - Capacitance Value in Farads. - L: float - Inductance in Henries. - reduction: float - The percentage that the TRV should be reduced by. - Rd0: float, optional - Damping Resistor Evaluation Starting Point, default=500 - wd0: float, optional - Omega-d evaluation starting point, default=260*k - tpk0: float, optional - Time of peak voltage evaluation starting point, default=10*u - - Returns - ------- - Rd: float - Damping resistor value, in ohms. - wd: float - Omega-d - tpk: float - Time of peak voltage. """ # Evaluate omega-n wn = 1 / _np.sqrt(L * C) @@ -925,29 +454,6 @@ def equations(data): def toctriptime(I, Ipickup, TD, curve="U1", CTR=1): """ Time OverCurrent Trip Time Function. - - Time-OverCurrent Trip Time Calculator, evaluates the time - to trip for a specific TOC (51) element given the curve - type, current characteristics and time-dial setting. - - Parameters - ---------- - I: float - Measured Current in Amps - Ipickup: float - Fault Current Pickup Setting (in Amps) - TD: float - Time Dial Setting - curve: string, optional - Name of specified TOC curve, may be entry from set: - {U1,U2,U3,U4,U5,C1,C2,C3,C4,C5}, default=U1 - CTR: float, optional - Current Transformer Ratio, default=1 - - Returns - ------- - tt: float - Time-to-Trip for characterized element. """ # Condition Inputs curve = curve.upper() @@ -977,28 +483,6 @@ def toctriptime(I, Ipickup, TD, curve="U1", CTR=1): def tocreset(I, Ipickup, TD, curve="U1", CTR=1): """ Time OverCurrent Reset Time Function. - - Function to calculate the time to reset for a TOC - (Time-OverCurrent, 51) element. - - Parameters - ---------- - I: float - Measured Current in Amps - Ipickup: float - Fault Current Pickup Setting (in Amps) - TD: float - Time Dial Setting - curve: string, optional - Name of specified TOC curve, may be entry from set: - {U1,U2,U3,U4,U5,C1,C2,C3,C4,C5}, default=U1 - CTR: float, optional - Current Transformer Ratio, default=1 - - Returns - ------- - tr: float - Time-to-Reset for characterized element. """ # Condition Inputs curve = curve.upper() @@ -1018,32 +502,6 @@ def tocreset(I, Ipickup, TD, curve="U1", CTR=1): def pickup(Iloadmax, Ifaultmin, scale=0, printout=False, units="A"): """ Electrical Current Pickup Selection Assistant. - - Used to assist in evaluating an optimal phase-over-current pickup - setting. Uses maximum load and minimum fault current to provide - user assistance. - - Parameters - ---------- - Iloadmax: float - The maximum load current in amps. - Ifaultmin: float - The minimum fault current in amps. - scale: int, optional - Control scaling to set number of significant figures. - default=0 - printout: boolean, optional - Control argument to enable printing of intermediate - stages, default=False. - units: string, optional - String to be appended to any printed output denoting - the units of which are being printed, default="A" - - Returns - ------- - setpoint: float - The evaluated setpoint at which the function suggests - the phase-over-current pickup setting be placed. """ IL2 = 2 * Iloadmax IF2 = Ifaultmin / 2 @@ -1065,54 +523,6 @@ def tdradial(I, CTI, Ipu_up, Ipu_dn=0, TDdn=0, curve="U1", scale=2, freq=60, CTR_up=1, CTR_dn=1, tfixed=None): """ Radial Time Dial Coordination Function. - - Function to evaluate the Time-Dial (TD) setting in radial schemes - where the Coordinating Time Interval (CTI) and the up/downstream - pickup settings are known along with the TD setting for the - downstream protection. - - Parameters - ---------- - I: float - Measured fault current in Amps, typically set using the - maximum fault current available. - CTI: float - Coordinating Time Interval in cycles. - Ipu_up: float - Pickup setting for upstream protection, - specified in amps - Ipu_dn: float, optional - Pickup setting for downstream protection, - specified in amps, default=0 - TDdn: float, optional - Time-Dial setting for downstream protection, - specified in seconds, default=0 - curve: string, optional - Name of specified TOC curve, may be entry from set: - {U1,U2,U3,U4,U5,C1,C2,C3,C4,C5}, default=U1 - scale: int, optional - Scaling value used to evaluate a practical TD - setting, default=2 - freq: float, optional - System operating frequency, default=60 - CTR_up: float, optional - Current Transformer Ratio for upstream relay. - default=1 - CTR_dn: float, optional - Current Transformer Ratio for downstream relay. - default=1 - tfixed: float, optional - Used to specify a fixed time delay for coordinated - protection elements, primarily used for coordinating - TOC elements (51) with OC elements (50) with a fixed - tripping time. Overrides downstream TOC arguments - including *Ipu_dn* and *TDdn*. - - Returns - ------- - TD: float - Calculated Time-Dial setting according to radial - scheme logical analysis. """ # Condition Inputs curve = curve.upper() @@ -1142,8 +552,8 @@ def tdradial(I, CTI, Ipu_up, Ipu_dn=0, TDdn=0, curve="U1", scale=2, freq=60, tpu_desired = tfixed + CTI # Re-Evaluate M M = I / (CTR_up * Ipu_up) - # Calculate TD setting - TD = tpu_desired / (A / (M ** 2 - 1) + B) + # Calculate TD setting (use curve exponent P, not hard-coded 2) + TD = tpu_desired / (A / (M ** P - 1) + B) # Scale and Round TD = _np.floor(TD * 10 ** scale) / 10 ** scale return TD @@ -1153,28 +563,6 @@ def tdradial(I, CTI, Ipu_up, Ipu_dn=0, TDdn=0, curve="U1", scale=2, freq=60, def protectiontap(S, CTR=1, VLN=None, VLL=None): """ Protection TAP Setting Calculator. - - Evaluates the required TAP setting based on the rated power of - a transformer (the object being protected) and the voltage - (either primary or secondary) in conjunction with the CTR - (current transformer ratio) for the side in question (primary/ - secondary). - - Parameters - ---------- - CTR: float - The Current Transformer Ratio. - S: float - Rated apparent power magnitude (VA/VAR/W). - VLN: float, exclusive - Line-to-Neutral voltage in volts. - VLL: float, exclusive - Line-to-Line voltage in volts. - - Returns - ------- - TAP: float - The TAP setting required to meet the specifications. """ # Condition Voltage(s) if VLL is not None: @@ -1192,30 +580,6 @@ def protectiontap(S, CTR=1, VLN=None, VLL=None): def correctedcurrents(Ipri, TAP, correction="Y", CTR=1): """ Electrical Transformer Current Correction Function. - - Function to evaluate the currents as corrected for microprocessor- - based relay protection schemes. - - Parameters - ---------- - Ipri: list of complex - Three-phase set (IA, IB, IC) of primary currents. - TAP: float - Relay's TAP setting. - correction: string, optional - String defining correction factor, may be one of: - (Y, D+, D-, Z); Y denotes Y (Y0) connection, D+ - denotes Dab (D1) connection, D- denotes Dac (D11) - connection, and Z (Z12) denotes zero-sequence - removal. default="Y" - CTR: float - Current Transformer Ratio, default=1 - - Returns - ------- - Isec_corr: list of complex - The corrected currents to perform operate/restraint - calculations with. """ # Define Matrix Lookup MAT = {"Y": XFMY0, @@ -1244,53 +608,6 @@ def iopirt(IpriHV, IpriLV, TAPHV, TAPLV, corrHV="Y", corrLV="Y", CTRHV=1, CTRLV=1): """ Operate/Restraint Current Calculator. - - Calculates the operating current (Iop) and the restraint - current (Irt) as well as the slope. - - Parameters - ---------- - IpriHV: list of complex - Three-phase set (IA, IB, IC) of primary currents - on the high-voltage side of power transformer. - IpriLV list of complex - Three-phase set (IA, IB, IC) of primary currents - on the low-voltage side of power transformer. - TAPHV float - Relay's TAP setting for high-voltage side of - power transformer. - TAPLV float - Relay's TAP setting for low-voltage side of - power transformer. - corrHV string, optional - String defining correction factor on high-voltage - side of power transformer, may be one of: - (Y, D+, D-, Z); Y denotes Y (Y0) connection, D+ - denotes Dab (D1) connection, D- denotes Dac (D11) - connection, and Z (Z12) denotes zero-sequence - removal. default="Y" - corrLV string, optional - String defining correction factor on low-voltage - side of power transformer, may be one of: - (Y, D+, D-, Z); Y denotes Y (Y0) connection, D+ - denotes Dab (D1) connection, D- denotes Dac (D11) - connection, and Z (Z12) denotes zero-sequence - removal. default="Y" - CTRHV float - Current Transformer Ratio for high-voltage side - of power transformer, default=1 - CTRLV float - Current Transformer Ratio for low-voltage side - of power transformer, default=1 - - Returns - ------- - Iop: list of float - The operating currents for phases A, B, and C. - Irt: list of float - The restraint currents for phases A, B, and C. - slope: list of float - The calculated slopes for phases A, B, and C. """ # Calculate Corrected Currents IcorHV = correctedcurrents(IpriHV, TAPHV, corrHV, CTRHV) @@ -1307,33 +624,6 @@ def iopirt(IpriHV, IpriLV, TAPHV, TAPLV, corrHV="Y", corrLV="Y", CTRHV=1, def symrmsfaultcur(V, R, X, t=1 / 60, freq=60): """ Symmetrical/RMS Current Calculator. - - Function to evaluate the time-constant tau, the symmetrical fault current, - and the RMS current for a faulted circuit. - - Parameters - ---------- - V: float - Voltage magnitude at fault point, - not described as line-to-line or - line-to-neutral. - R: float - The fault resistance in ohms. - X: float - The fault impedance in ohms. - t: float, optional - The time in seconds. - freq: float, optional - The system frequency in Hz. - - Returns - ------- - tau: float - The time-constant tau in seconds. - Isym: float - Symmetrical fault current in amps. - Irms: float - RMS fault current in amps. """ # Calculate Z and tau Z = _np.sqrt(R ** 2 + X ** 2) @@ -1349,25 +639,6 @@ def symrmsfaultcur(V, R, X, t=1 / 60, freq=60): def faultratio(I, Ipickup, CTR=1): """ Fault Multiple of Pickup (Ratio) Calculator. - - Evaluates the CTR-scaled pickup measured to pickup current ratio. - - M = meas / pickup - - Parameters - ---------- - I: float - Measured Current in Amps - Ipickup: float - Fault Current Pickup Setting (in Amps) - CTR: float, optional - Current Transformer Ratio for relay, - default=1 - - Returns - ------- - M: float - The measured-to-pickup ratio """ M = I / (CTR * Ipickup) return M @@ -1377,31 +648,6 @@ def faultratio(I, Ipickup, CTR=1): def residcomp(z1, z0, linelength=1): """ Residual Compensation Factor Function. - - Evaluates the residual compensation factor based on the line's positive and - zero sequence impedance characteristics. - - Parameters - ---------- - z1: complex - The positive-sequence impedance - characteristic of the line, specified in - ohms-per-unit where the total line length - (of same unit) is specified in - *linelength* argument. - z0: complex - The zero-sequence impedance characteristic - of the line, specified in ohms-per-unit - where the total line length (of same unit) - is specified in *linelength* argument. - linelength: float, optional - The length (in same base unit as impedance - characteristics) of the line. default=1 - - Returns - ------- - k0: complex - The residual compensation factor. """ # Evaluate Z1L and Z0L Z1L = z1 * linelength @@ -1416,52 +662,9 @@ def distmeasz(VLNmeas, If, Ip, Ipp, CTR=1, VTR=1, k0=None, z1=None, z0=None, linelength=1): """ Distance Element Measured Impedance Function. - - Function to evaluate the Relay-Measured-Impedance as calculated from - the measured voltage, current, and line parameters. - - Parameters - ---------- - VLNmeas: complex - Measured Line-to-Neutral voltage for the - faulted phase in primary volts. - If: complex - Faulted phase current measured in primary - amps. - Ip: complex - Secondary phase current measured in primary - amps. - Ipp: complex - Terchiary phase current measured in primary - amps. - CTR: float, optional - Current transformer ratio, default=1 - VTR: float, optional - Voltage transformer ratio, default=1 - k0: complex, optional - Residual Compensation Factor - z1: complex, optional - The positive-sequence impedance - characteristic of the line, specified in - ohms-per-unit where the total line length - (of same unit) is specified in - *linelength* argument. - z0: complex, optional - The zero-sequence impedance characteristic - of the line, specified in ohms-per-unit - where the total line length (of same unit) - is specified in *linelength* argument. - linelength: float, optional - The length (in same base unit as impedance - characteristics) of the line. default=1 - - Returns - ------- - Zmeas: complex - The "measured" impedance as calculated by the relay. """ # Validate Residual Compensation Inputs - if k0 == z1 == z0 is None: + if k0 is None and z1 is None and z0 is None: raise ValueError("Residual compensation arguments must be set.") if k0 is None and (z1 is None or z0 is None): raise ValueError("Both *z1* and *z0* must be specified.") @@ -1481,25 +684,6 @@ def distmeasz(VLNmeas, If, Ip, Ipp, CTR=1, VTR=1, k0=None, z1=None, z0=None, def transmismatch(I1, I2, tap1, tap2): """ Electrical Transformer TAP Mismatch Function. - - Function to evaluate the transformer ratio mismatch for protection. - - Parameters - ---------- - I1: complex - Current (in amps) on transformer primary side. - I2: complex - Current (in amps) on transformer secondary. - tap1: float - Relay TAP setting on the primary side. - tap2: float - Relay TAP setting on the secondary side. - - Returns - ------- - mismatch: float - Transformer CT mismatch value associated with - relay. """ # Evaluate MR MR = min(abs(I1 / I2), abs(tap1 / tap2)) @@ -1513,40 +697,6 @@ def highzvpickup(I, RL, Rct, CTR=1, threephase=False, Ks=1.5, Vstd=400, Kd=0.5): """ High Impedance Pickup Setting Function. - - Evaluates the voltage pickup setting for a high - impedance bus protection system. - - Parameters - ---------- - I: float - Fault current on primary side (in amps) - RL: float - One-way line resistance in ohms - Rct: float - Current Transformer Resistance in ohms - CTR: float, optional - Current Transformer Ratio, default=1 - threephase: bool, optional - Control argument to set the function to - evaluate the result for a three-phase - fault or unbalanced fault. default=False - Ks: float, optional - Security Factor for secure voltage pickup - setting, default=1.5 - Vstd: float, optional - C-Class Voltage rating (i.e. C-400), - default=400 - Kd: float, optional - The dependability factor for dependable - voltage pickup setting, default=0.5 - - Returns - ------- - Vsens: float - The calculated sensetive voltage-pickup. - Vdep: float - The calculated dependable voltage-pickup. """ # Condition Based on threephase Argument n = 2 @@ -1562,36 +712,9 @@ def highzvpickup(I, RL, Rct, CTR=1, threephase=False, Ks=1.5, def highzmini(N, Ie, Irly=None, Vset=None, Rrly=2000, Imov=0, CTR=1): """ Minimum Current for High Impedance Protection Calculator. - - Evaluates the minimum pickup current required to cause - high-impedance bus protection element pickup. - - Parameters - ---------- - N: int - Number of Current Transformers included in scheme - Ie: float - The excitation current at the voltage setting - Irly: float, optional - The relay current at voltage setting - Vset: float, optional - The relay's voltage pickup setting in volts. - Rrly: float, optional - The relay's internal resistance in ohms, default=2000 - Imov: float, optional - The overvoltage protection current at the - voltage setting. default=0.0 - CTR: float, optional - Current Transformer Ratio, default=1 - - Returns - ------- - Imin: float - Minimum current required to cause high-impedance - bus protection element pickup. """ # Validate Inputs - if Irly == Vset is None: + if Irly is None and Vset is None: raise ValueError("Relay Current Required.") # Condition Inputs Ie = abs(Ie) @@ -1610,24 +733,6 @@ def highzmini(N, Ie, Irly=None, Vset=None, Rrly=2000, Imov=0, CTR=1): def instoc(Imin, CTR=1, Ki=0.5): """ Instantaneous OverCurrent Pickup Calculator. - - Using a sensetivity factor and the CTR, evaluates the secondary-level pickup - setting for an instantaneous overcurrent element. - - Parameters - ---------- - Imin: float - The minimum fault current in primary amps. - CTR: float, optional - Current Transformer Ratio, default=1 - Ki: Sensetivity factor, default=0.5 - - Returns - ------- - Ipu: float - The pickup setting for the instantaneous - overcurrent element as referred to the - secondary side. """ # Evaluate Overcurrent Pickup Setting Ipu = Ki * abs(Imin) / CTR @@ -1638,36 +743,6 @@ def instoc(Imin, CTR=1, Ki=0.5): def genlossfield(Xd, Xpd, Zbase=1, CTR=1, VTR=1): """ Electric Generator Loss of Field Function. - - Generates the Loss-of-Field Element settings for a generator using the Xd - value and per-unit base information. - - Parameters - ---------- - Xd: float - The Transient Reactance (Xd) term. May be - specified in ohms or Per-Unit ohms if - *Zbase* is set. - Xpd: float - The Sub-Transient Reactance (X'd) term. May - be specified in ohms or Per-Unit ohms if - *Zbase* is set. - Zbase: float, optional - Base impedance, used to convert per-unit - Xd and Xpd to secondary values. default=1 - CTR: float, optional - Current Transformer Ratio, default=1 - VTR: float, optional - Voltage Transformer Ratio, default=1 - - Returns - ------- - ZoneOff: float - Zone Offset in ohms. - Z1dia: float - Zone 1 diameter in ohms. - Z2dia: float - Zone 2 diameter in ohms. """ # Condition Inputs Xd = abs(Xd) @@ -1675,7 +750,7 @@ def genlossfield(Xd, Xpd, Zbase=1, CTR=1, VTR=1): Zbase = abs(Zbase) # Evaluate Xd_sec and Xpd_sec Xd_sec = Xd * Zbase * (CTR / VTR) - Xpd_sec = Xd * Zbase * (CTR / VTR) + Xpd_sec = Xpd * Zbase * (CTR / VTR) # Determine Zone Offset ZoneOff = Xpd_sec / 2 # Evaluate Z1 Diameter and Z2 Diameter @@ -1689,31 +764,6 @@ def genlossfield(Xd, Xpd, Zbase=1, CTR=1, VTR=1): def thermaltime(In, Ibase, tbase): r""" Thermal Time Limit Calculator. - - Computes the maximum allowable time for a specified current `In` given - parameters for a maximum current and time at some other level, (`Ibase`, - `tbase`). - - Uses the following formula: - - .. math:: t_n=\frac{I_{base}^2*t_{base}}{I_n^2} - - Parameters - ---------- - In: float - Current at which to calculate max time. - Ibase: float - Base current, at which maximum time - `tbase` is allowable. - tbase: float - Base time for which a maximum allowable - current `Ibase` is specified. Unitless. - - Returns - ------- - tn: float - Time allowable for specified current, - `In`. """ # Perform Calculation tn = (Ibase ** 2 * tbase) / (In ** 2) @@ -1724,42 +774,6 @@ def thermaltime(In, Ibase, tbase): def synmach_Isym(t, Eq, Xd, Xdp, Xdpp, Tdp, Tdpp): r""" Synch. Machine Symmetrical Fault Current Calc. - - Determines the Symmetrical Fault Current of a synchronous - machine given the machine parameters, the internal voltage, - and the time for which to calculate. - - .. math:: I_a(t)=\sqrt{2}\left|E_q\right|\left[ - \frac{1}{X_d}+\left(\frac{1}{X'_d}-\frac{1}{X_d} - \right)\cdot e^{\frac{-t}{T'_d}}+\left(\frac{1} - {X"_d}-\frac{1}{X'_d}\right)\cdot e^{\frac{-t}{T"_d}} - \right] - - Parameters - ---------- - t: float - Time at which to calculate the fault current - Eq: float - The internal machine voltage in per-unit-volts - Xd: float - The Xd (d-axis) reactance in per-unit-ohms - Xdp: float - The X"d (d-axis transient) reactance in - per-unit-ohms - Xdpp: float - The X"d (d-axis subtransient) reactance in - per-unit-ohms - Tdp: float - The T'd (d-axis transient) time constant of the - machine in seconds - Tdpp: float - The T"d (d-axis subtransient) time constant of - the machine in seconds - - Returns - ------- - Ia: float - Peak symmetrical fault current in per-unit-amps """ # Calculate Time-Constant Term t_c = ( @@ -1776,34 +790,6 @@ def synmach_Isym(t, Eq, Xd, Xdp, Xdpp, Tdp, Tdpp): def synmach_Iasym(t, Eq, Xdpp, Xqpp, Ta): r""" Synch. Machine Asymmetrical Fault Current Calc. - - Determines the asymmetrical fault current of a synchronous - machine given the machine parameters, the internal voltage, - and the time for which to calculate. - - .. math:: I_{asym}=\sqrt{2}\left|E_q\right|\frac{1}{2} - \left(\frac{1}{X"_d}+\frac{1}{X"_q}\right)e^{\frac{-t} - {T_a}} - - Parameters - ---------- - t: float - Time at which to calculate the fault current - Eq: float - The internal machine voltage in per-unit-volts - Xdpp: float - The X"d (d-axis subtransient) reactance in - per-unit-ohms - Xqpp: float - The X"q (q-axis subtransient) reactance in - per-unit-ohms - Ta: float - Armature short-circuit (DC) time constant in seconds - - Returns - ------- - Iasym: float - Peak asymmetrical fault current in per-unit-amps """ # Calculate Time Constant Term t_c = 1 / Xdpp + 1 / Xqpp @@ -1816,37 +802,6 @@ def synmach_Iasym(t, Eq, Xdpp, Xqpp, Ta): def indmacheigenvalues(Lr, Ls, Lm, Rr, Rs, wrf=0, freq=60): """ Induction Machine Eigenvalue Calculator. - - Calculates the pertinent eigenvalues for an unloaded - induction machine given a specific set of machine - parameters. - - Parameters - ---------- - Lr: float - Inductance of the Rotor (in Henrys). - Ls: float - Inductance of the Stator (in Henrys). - Lm: float - Inductance of the Magnetizing branch - (in Henrys). - Rr: float - Resistance of the Rotor (in Ohms). - Rs: float - Resistance of the Stator (in Ohms). - wrf: float, optional - Frequency (in radians/sec) of the rotor slip. - default=0 - freq: float, optional - Base frequency of the system (in Hertz). - default=60 - - Returns - ------- - lam1: complex - The First Eigenvalue - lam2: complex - The Second Eigenvalue """ # Calculate Required Values omega_e_base = 2 * _np.pi * freq @@ -1875,44 +830,6 @@ def indmacheigenvalues(Lr, Ls, Lm, Rr, Rs, wrf=0, freq=60): def indmachphs3sc(t, Is0, Lr, Ls, Lm, Rr, Rs, wrf=0, freq=60, real=True): """ Induction Machine 3-Phase SC Calculator. - - Determines the short-circuit current at a specified time for a three-phase - fault on an unloaded induction machine. - - Parameters - ---------- - t: array_like - The time at which to find the - current, may be int, float, or - numpy array. - Is0: complex - The initial (t=0) current on - the stator. - Lr: float - Inductance of the Rotor (in Henrys). - Ls: float - Inductance of the Stator (in Henrys). - Lm: float - Inductance of the Magnetizing branch - (in Henrys). - Rr: float - Resistance of the Rotor (in Ohms). - Rs: float - Resistance of the Stator (in Ohms). - wrf: float, optional - Frequency (in radians/sec) of the rotor slip. - default=0 - freq: float, optional - Base frequency of the system (in Hertz). - default=60 - real: bool, optional - Control argument to force returned value - to be real part only. default=True - - Returns - ------- - ias: array_like - Fault Current """ # Calculate Required Values omega_r = 2 * _np.pi * freq @@ -1936,43 +853,6 @@ def indmachphs3sc(t, Is0, Lr, Ls, Lm, Rr, Rs, wrf=0, freq=60, real=True): def indmachphs3torq(t, Is0, Lr, Ls, Lm, Rr, Rs, wrf=0, freq=60): """ Induction Machine 3-Phase Torque Calculator. - - Determines the torque exerted during a three-phase fault on an induction - machine. - - Parameters - ---------- - t: array_like - The time at which to find the - current, may be int, float, or - numpy array. - Is0: complex - The initial (t=0) current on - the stator. - Lr: float - Inductance of the Rotor (in Henrys). - Ls: float - Inductance of the Stator (in Henrys). - Lm: float - Inductance of the Magnetizing branch - (in Henrys). - Rr: float - Resistance of the Rotor (in Ohms). - Rs: float - Resistance of the Stator (in Ohms). - p: int - Number of electrical poles. - wrf: float, optional - Frequency (in radians/sec) of the rotor slip. - default=0 - freq: float, optional - Base frequency of the system (in Hertz). - default=60 - - Returns - ------- - Tem: array_like - Induction machine torque in N*m """ # Calculate Required Values omega_r = 2 * _np.pi * freq @@ -2000,51 +880,6 @@ def synmach_ifault(t, Ea, alpha, Xd, Xdp, Xdpp, Xqpp, Tdp, Tdpp, Ta, freq=60): # noqa: D401 "Synchronous" is intentional descriptor """ Synchronous Machine Fault Current Calculator. - - Given machine parameters, fault inception angle, and time at - which to calculate fault current, this function will identify - the complete (symmetrical, asymmetrical, and double frequency) - fault current. - - .. image:: /static/synmach_ifault_formula.png - - Parameters - ---------- - t: float - Time at which to calculate the fault current - Eq: float - The internal machine voltage in per-unit-volts - alpha: float - Fault inception angle (in degrees) - Xd: float - The Xd (d-axis) reactance in per-unit-ohms - Xdp: float - The X"d (d-axis transient) reactance in - per-unit-ohms - Xdpp: float - The X"d (d-axis sub-transient) reactance in - per-unit-ohms - Xqpp: float - The X"q (q-axis sub-transient) reactance in - per-unit-ohms - Tdp: float - The T'd (d-axis transient) time constant of the - machine in seconds - Tdpp: float - The T"d (d-axis sub-transient) time constant of - the machine in seconds - Ta: float - Armature short-circuit (DC) time constant in seconds - freq: float, optional - System (electrical) frequency (in degrees), - default=60 - - Returns - ------- - ias: float - Synchronous machine fault current (symmetrical, - asymmetrical, and double frequency component) in - amps """ # Calculate we Component we = 2 * _np.pi * freq @@ -2053,20 +888,30 @@ def synmach_ifault(t, Ea, alpha, Xd, Xdp, Xdpp, Xqpp, Tdp, Tdpp, Ta, freq=60): alpha = _np.radians(alpha) # Define Constant Term const = _np.sqrt(2) * Ea + + # Symmetrical AC component magnitude term + isym_mag = ( + (1 / Xd) + + (1 / Xdp - 1 / Xd) * _np.exp(-t / Tdp) + + (1 / Xdpp - 1 / Xdp) * _np.exp(-t / Tdpp) + ) + + # DC offset / asymmetry magnitude term if Xqpp != 0: val = 1 / Xqpp else: val = 0 - asym = 1 / 2 * (1 / Xdpp + val) * _np.exp(t / Ta) + asym_mag = 1 / 2 * (1 / Xdpp + val) * _np.exp(-t / Ta) + # Define Symmetrical Portion - isym = ( - const * (1 / Xd + (1 / Xdp - 1 / Xd) * _np.exp(-t / Tdp) + - (1 / Xdpp - 1 / Xdp) * _np.exp(-t / Tdpp)) * _np.sin(we * t + alpha) - ) - # Define Asymmetrical Portion - iasym = const * asym * _np.sin(alpha) + isym = const * isym_mag * _np.sin(we * t + alpha) + + # Define Asymmetrical (DC offset) Portion + iasym = const * asym_mag * _np.sin(alpha) + # Define Double Frequency Term - idbl = const * 1 / 2 * asym * _np.sin(2 * we * t + alpha) + idbl = const * 1 / 2 * asym_mag * _np.sin(2 * we * t + alpha) + # Compose Complete Current Value ias = isym - iasym - idbl return ias diff --git a/electricpy/machines.py b/electricpy/machines.py index efadc749..08acc2b7 100644 --- a/electricpy/machines.py +++ b/electricpy/machines.py @@ -13,9 +13,46 @@ from electricpy.phasors import compose +# Internal Helpers +def _as_array(x): + return _np.asarray(x) if isinstance(x, (list, tuple)) else x + + +def _validate_pf(PF): + # Clamp for numerical safety + if PF is None: + raise ValueError("Power factor (PF) must be specified.") + PF = float(PF) + if PF > 1.0: + PF = 1.0 + if PF < -1.0: + PF = -1.0 + return PF + + +def _derive_leakage(Lm, Lls=0, Llr=0, Ls=None, Lr=None): + # Preserve original semantics: if Ls/Lr provided, they override leakage via Ls-Lm / Lr-Lm + if Ls is not None: + Lls = Ls - Lm + if Lr is not None: + Llr = Lr - Lm + return Lls, Llr + + +def _to_reactances(Lm, Lls=0, Llr=0, freq=60, calcX=True): + # If calcX True: inputs are inductances (H) and convert to reactances (ohms) + # If calcX False: inputs are already reactances (ohms) and should be used directly + if calcX: + w = 2 * _np.pi * freq + return (Lm * w), (Lls * w), (Llr * w), w + else: + w = 2 * _np.pi * freq + return (Lm), (Lls), (Llr), w + + # Define Transformer Short-Circuit/Open-Circuit Function -def transformertest(Poc=False, Voc=False, Ioc=False, Psc=False, Vsc=False, - Isc=False): +def transformertest(Poc=None, Voc=None, Ioc=None, Psc=None, Vsc=None, + Isc=None): """ Electrical Transformer Rated Test Evaluator. @@ -26,48 +63,38 @@ def transformertest(Poc=False, Voc=False, Ioc=False, Psc=False, Vsc=False, { Psc, Vsc, Isc }. All values given must be given as absolute value, not complex. All values returned are given with respect to primary. - - Parameters - ---------- - Poc: float, optional - The open-circuit measured power (real power), default=None - Voc: float, optional - The open-circuit measured voltage (measured on X), - default=None - Ioc: float, optional - The open-circuit measured current (measured on primary), - default=None - Psc: float, optional - The short-circuit measured power (real power), default=None - Vsc: float, optional - The short-circuit measured voltage (measured on X), - default=None - Isc: float, optional - The short-circuit measured current (measured on X), - default=None - - Returns - ------- - {Req,Xeq,Rc,Xm}: Given all optional args - {Rc, Xm}: Given open-circuit parameters - {Req, Xeq}: Given short-circuit parameters """ SC = False OC = False + # Given Open-Circuit Values if (Poc is not None) and (Voc is not None) and (Ioc is not None): + if Voc == 0 or Ioc == 0: + raise ValueError("Voc and Ioc must be non-zero for open-circuit test.") PF = Poc / (Voc * Ioc) + PF = _validate_pf(PF) + # Admittance magnitude = I/V, angle = -acos(PF) (lagging magnetizing current) Y = _c.rect(Ioc / Voc, -_np.arccos(PF)) + if Y.real == 0: + raise ValueError("Invalid open-circuit data produced zero conductance.") Rc = 1 / Y.real + if Y.imag == 0: + raise ValueError("Invalid open-circuit data produced zero susceptance.") Xm = -1 / Y.imag OC = True + # Given Short-Circuit Values if (Psc is not None) and (Vsc is not None) and (Isc is not None): + if Vsc == 0 or Isc == 0: + raise ValueError("Vsc and Isc must be non-zero for short-circuit test.") PF = Psc / (Vsc * Isc) - Zeq = _c.rect(Vsc / Isc, _np.arccos(PF)) + PF = _validate_pf(PF) + # Impedance magnitude = V/I, angle = acos(PF) (assume inductive) + Zeq = _c.rect(Vsc / Isc, _np.arccos(abs(PF))) Req = Zeq.real Xeq = Zeq.imag SC = True + # Return All if Found if OC and SC: return (Req, Xeq, Rc, Xm) @@ -89,32 +116,6 @@ def phase_shift_transformer(style="DY", shift=30): Use with transformer orientation to evaluate the phase-shift across a transformer. For example, find the phase shift for a Delta-Wye transformer as seen from the delta side. - - Parameters - ---------- - style: {'DY','YD','DD','YY'}, optional - String that denotes the transformer - orientation. default='DY' - shift: float, optional - Transformer angle shift, default=30 - - Returns - ------- - phase: complex - Phasor including the phase shift and - positive or negative characteristic. - - Examples - -------- - >>> import electricpy as ep - >>> import electricpy.machines as machines - >>> # Find shift of Delta-Wye Transformer w/ 30° shift - >>> shift = machines.xfmphs(style="DY", shift=30) - >>> ep.cprint(shift) - 1.0 ∠ 30.0° - >>> shift = machines.phase_shift_transformer(style="DY", shift=30) - >>> ep.cprint(shift) - 1.0 ∠ 30.0° """ # Define Direction Dictionary orientation = { @@ -123,10 +124,13 @@ def phase_shift_transformer(style="DY", shift=30): "DD": 0, "YY": 0, } + key = str(style).upper() + if key not in orientation: + raise ValueError("Invalid transformer style. Use one of: 'DY','YD','DD','YY'.") # Find Direction - v = orientation[style.upper()] + v = orientation[key] # Calculate Shift - phase = _np.exp(1j * _np.radians(v * abs(shift))) + phase = _np.exp(1j * _np.radians(v * shift)) # Return return (phase) @@ -139,60 +143,16 @@ def indmachvth(Vas, Rs, Lm, Lls=0, Ls=None, freq=60, calcX=True): r""" Induction Machine Thevenin Voltage Calculator. - Function to calculate the Thevenin equivalent voltage of an - induction machine given a specific set of parameters. - .. math:: V_{th}=\frac{j\omega L_m}{R_s+j\omega(L_{ls}+L_m)}V_{as} - - where: - - .. math:: \omega = \omega_{es} = 2\pi\cdot f_{\text{electric}} - - Parameters - ---------- - Vas: complex - Terminal Stator Voltage in Volts - Rs: float - Stator resistance in ohms - Lm: float - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - Vth: complex - Thevenin-Equivalent voltage (in volts) of induction - machine described. - - See Also - -------- - indmachzth: Induction Machine Thevenin Impedance Calculator - indmachpem: Induction Machine Electro-Mechanical Power Calculator - indmachtem: Induction Machine Electro-Mechanical Torque Calculator - indmachpkslip: Induction Machine Peak Slip Calculator - indmachpktorq: Induction Machine Peak Torque Calculator - indmachiar: Induction Machine Phase-A Rotor Current Calculator - indmachstarttorq: Induction Machine Starting Torque Calculator """ - # Condition Inputs + # Condition Inputs (do not mutate user-provided variables) if Ls is not None: # Use Ls instead of Lls Lls = Ls - Lm - if calcX: # Convert Inductances to Reactances - w = 2 * _np.pi * freq - Lm *= w - Lls *= w + + Xlm, Xls, _, _ = _to_reactances(Lm, Lls, 0, freq, calcX) + # Calculate Thevenin Voltage, Return - Vth = 1j * Lm / (Rs + 1j * (Lls + Lm)) * Vas + Vth = 1j * Xlm / (Rs + 1j * (Xls + Xlm)) * Vas return (Vth) @@ -201,67 +161,17 @@ def indmachzth(Rs, Lm, Lls=0, Llr=0, Ls=None, Lr=None, freq=60, calcX=True): r""" Induction Machine Thevenin Impedance Calculator. - Function to calculate the Thevenin equivalent impedance of an - induction machine given a specific set of parameters. - .. math:: Z_{th} = \frac{(R_s+j\omega L_{ls})j\omega L_m} {R_s+j\omega(L_{ls}+L_m)}+j\omega L_{lr} - - where: - - .. math:: \omega = \omega_{es} = 2\pi\cdot f_{\text{electric}} - - Parameters - ---------- - Rs: float - Stator resistance in ohms - Lm: float - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Llr: float, optional - Rotor leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - Lr: float, optional - Rotor inductance in Henrys - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - Zth: complex - Thevenin-Equivalent impedance (in ohms) of induction - machine described. - - See Also - -------- - indmachvth: Induction Machine Thevenin Voltage Calculator - indmachpem: Induction Machine Electro-Mechanical Power Calculator - indmachtem: Induction Machine Electro-Mechanical Torque Calculator - indmachpkslip: Induction Machine Peak Slip Calculator - indmachpktorq: Induction Machine Peak Torque Calculator - indmachiar: Induction Machine Phase-A Rotor Current Calculator - indmachstarttorq: Induction Machine Starting Torque Calculator """ # Condition Inputs - if Ls is not None: # Use Ls instead of Lls - Lls = Ls - Lm - if Lr is not None: # Use Lr instead of Llr - Llr = Lr - Lm - if calcX: # Convert Inductances to Reactances - w = 2 * _np.pi * freq - Lm *= w - Lls *= w - Llr *= w + Lls, Llr = _derive_leakage(Lm, Lls, Llr, Ls, Lr) + + Xlm, Xls, Xlr, _ = _to_reactances(Lm, Lls, Llr, freq, calcX) + # Calculate Thevenin Impedance - Zth = (Rs + 1j * Lls) * (1j * Lm) / (Rs + 1j * (Lls + Lm)) + 1j * Llr + Zth = (Rs + 1j * Xls) * (1j * Xlm) / (Rs + 1j * (Xls + Xlm)) + 1j * Xlr return (Zth) @@ -271,92 +181,31 @@ def indmachpem(slip, Rr, Vth=None, Zth=None, Vas=0, Rs=0, Lm=0, Lls=0, r""" Mechanical Power Calculator for Induction Machines. - Function to calculate the mechanical power using the thevenin - equivalent circuit terms. - - .. math:: - P_{em}=\frac{|V_{th_{\text{stator}}}|^2\cdot\frac{R_r}{slip}} - {\left[\left(\frac{R_r}{slip}+R_{th_{\text{stator}}}\right)^2 - +X_{th_{\text{stator}}}^2\right]\cdot\omega_{es}}\cdot(1-slip) - - Parameters - ---------- - slip: float - The mechanical/electrical slip factor of the - induction machine. - Rr: float - Rotor resistance in ohms - Vth: complex, optional - Thevenin-equivalent stator voltage of the - induction machine, may be calculated internally - if given stator voltage and machine parameters. - Zth: complex, optional - Thevenin-equivalent inductance (in ohms) of the - induction machine, may be calculated internally - if given machine parameters. - Vas: complex, optional - Terminal Stator Voltage in Volts - Rs: float, optional - Stator resistance in ohms - Lm: float, optional - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Llr: float, optional - Rotor leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - Lr: float, optional - Rotor inductance in Henrys - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - Pem: float - Power (in watts) that is produced or consumed - by the mechanical portion of the induction machine. - - See Also - -------- - indmachvth: Induction Machine Thevenin Voltage Calculator - indmachzth: Induction Machine Thevenin Impedance Calculator - indmachtem: Induction Machine Electro-Mechanical Torque Calculator - indmachpkslip: Induction Machine Peak Slip Calculator - indmachpktorq: Induction Machine Peak Torque Calculator - indmachiar: Induction Machine Phase-A Rotor Current Calculator - indmachstarttorq: Induction Machine Starting Torque Calculator + Uses thevenin-equivalent terms. """ - # Condition Inputs + if slip == 0: + raise ValueError("slip must be non-zero.") + if slip < 0: + raise ValueError("slip must be positive for motoring conventions used here.") + + # Electrical rad/sec (used in power conversion term) w = 2 * _np.pi * freq - if Ls is not None: # Use Ls instead of Lls - Lls = Ls - Lm - if Lr is not None: # Use Lr instead of Llr - Llr = Lr - Lm - if calcX: # Convert Inductances to Reactances - Lm *= w - Lls *= w - Llr *= w - # Test for Valid Input Set + + # Test for Valid Input Set and Calculate Vth/Zth if needed if Vth is None: - if not all((Vas, Rs, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Vth + if not all((Vas, Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Vth.") Vth = indmachvth(Vas, Rs, Lm, Lls, Ls, freq, calcX) + if Zth is None: - if not all((Rs, Llr, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Zth + if not all((Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Zth.") Zth = indmachzth(Rs, Lm, Lls, Llr, Ls, Lr, freq, calcX) + # Use Terms to Calculate Pem Rth = Zth.real Xth = Zth.imag - Pem = (abs(Vth) ** 2 * Rr / slip) / (((Rr / slip + Rth) ** 2 + Xth ** 2) * w) * (1 - slip) + Pem = (abs(Vth) ** 2 * (Rr / slip)) / (((Rr / slip + Rth) ** 2 + (Xth ** 2)) * w) * (1 - slip) return (Pem) @@ -365,108 +214,36 @@ def indmachtem(slip, Rr, p=0, Vth=None, Zth=None, Vas=0, Rs=0, Lm=0, Lls=0, Llr=0, Ls=None, Lr=None, wsyn=None, freq=60, calcX=True): r""" Induction Machine Torque Calculator. - - Calculate the torque generated or consumed by an induction - machine given the machine parameters of Vth and Zth by use - of the equation below. - - .. math:: - T_{em}=\frac{3|V_{th_{\text{stator}}}|^2} - {\left[\left(\frac{R_r}{slip}+R_{th_{\text{stator}}}\right)^2 - +X_{th_{\text{stator}}}\right]}\frac{R_r}{slip*\omega_{sync}} - - where: - - .. math:: - \omega_{sync}=\frac{\omega_{es}}{\left(\frac{poles}{2}\right)} - - Parameters - ---------- - slip: float - The mechanical/electrical slip factor of the - induction machine. - Rr: float - Rotor resistance in ohms - p: int, optional - Number of poles in the induction machine - Vth: complex, optional - Thevenin-equivalent stator voltage of the - induction machine, may be calculated internally - if given stator voltage and machine parameters. - Zth: complex, optional - Thevenin-equivalent inductance (in ohms) of the - induction machine, may be calculated internally - if given machine parameters. - Vas: complex, optional - Terminal Stator Voltage in Volts - Rs: float, optional - Stator resistance in ohms - Lm: float, optional - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Llr: float, optional - Rotor leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - Lr: float, optional - Rotor inductance in Henrys - wsync: float, optional - Synchronous speed in rad/sec, may be specified - directly as a replacement of p (number of poles). - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - Tem: float - Torque (in Newton-meters) that is produced or consumed - by the mechanical portion of the induction machine. - - See Also - -------- - indmachvth: Induction Machine Thevenin Voltage Calculator - indmachzth: Induction Machine Thevenin Impedance Calculator - indmachpem: Induction Machine Electro-Mechanical Power Calculator - indmachpkslip: Induction Machine Peak Slip Calculator - indmachpktorq: Induction Machine Peak Torque Calculator - indmachiar: Induction Machine Phase-A Rotor Current Calculator - indmachstarttorq: Induction Machine Starting Torque Calculator """ - # Condition Inputs - w = 2 * _np.pi * freq - if Ls is not None: # Use Ls instead of Lls - Lls = Ls - Lm - if Lr is not None: # Use Lr instead of Llr - Llr = Lr - Lm - if p != 0: # Calculate Sync. Speed from Num. Poles - wsyn = w / (p / 2) - if calcX: # Convert Inductances to Reactances - Lm *= w - Lls *= w - Llr *= w - # Test for Valid Input Set + if slip == 0: + raise ValueError("slip must be non-zero.") if not any((p, wsyn)): raise ValueError("Poles or Synchronous Speed must be specified.") + + # Electrical rad/sec + w = 2 * _np.pi * freq + + # Determine synchronous mechanical speed + if wsyn is None: + if p == 0: + raise ValueError("If wsyn is not specified, p must be non-zero.") + wsyn = w / (p / 2) + + # Calculate Vth/Zth if needed if Vth is None: - if not all((Vas, Rs, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Vth + if not all((Vas, Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Vth.") Vth = indmachvth(Vas, Rs, Lm, Lls, Ls, freq, calcX) + if Zth is None: - if not all((Rs, Llr, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Zth + if not all((Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Zth.") Zth = indmachzth(Rs, Lm, Lls, Llr, Ls, Lr, freq, calcX) - # Use Terms to Calculate Pem + + # Use Terms to Calculate Tem Rth = Zth.real Xth = Zth.imag - Tem = 3 * abs(Vth) ** 2 / ((Rr / slip + Rth) ** 2 + Xth) * Rr / (slip * wsyn) + Tem = 3 * (abs(Vth) ** 2) * (Rr / slip) / (((Rr / slip + Rth) ** 2 + (Xth ** 2)) * wsyn) return Tem @@ -476,78 +253,18 @@ def indmachpkslip(Rr, Zth=None, Rs=0, Lm=0, Lls=0, Llr=0, Ls=None, r""" Induction Machine Slip at Peak Torque Calculator. - Function to calculate the slip encountered by an induction machine - with the parameters specified when the machine is generating peak - torque. Uses formula as shown below. - .. math:: \text{slip} = \frac{R_r}{|Z_{th}|} - - where: - - .. math:: - Z_{th} = \frac{(R_s+j\omega L_{ls})j\omega L_m} - {R_s+j\omega(L_{ls}+L_m)}+j\omega L_{lr} - - Parameters - ---------- - Rr: float - Rotor resistance in ohms - Zth: complex, optional - Thevenin-equivalent inductance (in ohms) of the - induction machine, may be calculated internally - if given machine parameters. - Rs: float, optional - Stator resistance in ohms - Lm: float, optional - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Llr: float, optional - Rotor leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - Lr: float, optional - Rotor inductance in Henrys - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - s_peak: float - The peak slip for the induction machine described. - - See Also - -------- - indmachvth: Induction Machine Thevenin Voltage Calculator - indmachzth: Induction Machine Thevenin Impedance Calculator - indmachpem: Induction Machine Electro-Mechanical Power Calculator - indmachtem: Induction Machine Electro-Mechanical Torque Calculator - indmachpktorq: Induction Machine Peak Torque Calculator - indmachiar: Induction Machine Phase-A Rotor Current Calculator - indmachstarttorq: Induction Machine Starting Torque Calculator """ - # Condition Inputs - w = 2 * _np.pi * freq - if Ls is not None: # Use Ls instead of Lls - Lls = Ls - Lm - if Lr is not None: # Use Lr instead of Llr - Llr = Lr - Lm - if calcX: # Convert Inductances to Reactances - Lm *= w - Lls *= w - Llr *= w - # Test for Valid Input Set + if Rr is None: + raise ValueError("Rr must be specified.") + if Zth is None: - if not all((Rs, Llr, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Zth + if not all((Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Zth.") Zth = indmachzth(Rs, Lm, Lls, Llr, Ls, Lr, freq, calcX) - # Calculate Peak Slip + + if abs(Zth) == 0: + raise ValueError("Invalid Zth: magnitude is zero.") s_peak = Rr / abs(Zth) return s_peak @@ -558,97 +275,22 @@ def indmachiar(poles=0, Vth=None, Zth=None, Vas=0, Rs=0, Lm=0, Lls=0, Llr=0, r""" Induction Machine Rotor Current Calculator. - Calculation function to find the phase-A, rotor current for an - induction machine given the thevenin voltage and impedance. - - This current is calculated using the following formulas: - .. math:: I_{a_{\text{rotor}}} = \frac{V_{th}}{|Z_{th}|+Z_{th}} - - where: - - .. math:: V_{th}=\frac{j\omega L_m}{R_s+j\omega(L_{ls}+L_m)}V_{as} - - .. math:: - Z_{th} = \frac{(R_s+j\omega L_{ls})j\omega L_m} - {R_s+j\omega(L_{ls}+L_m)}+j\omega L_{lr} - - .. math:: \omega = \omega_{es} = 2\pi\cdot f_{\text{electric}} - - Parameters - ---------- - poles: int, optional - Number of poles for the induction machine. - Vth: complex, optional - Thevenin-equivalent stator voltage of the - induction machine, may be calculated internally - if given stator voltage and machine parameters. - Zth: complex, optional - Thevenin-equivalent inductance (in ohms) of the - induction machine, may be calculated internally - if given machine parameters. - Vas: complex, optional - Terminal Stator Voltage in Volts - Rs: float, optional - Stator resistance in ohms - Lm: float, optional - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Llr: float, optional - Rotor leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - Lr: float, optional - Rotor inductance in Henrys - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - Iar: complex - The rotor, phase-A current in amps. - - See Also - -------- - indmachvth: Induction Machine Thevenin Voltage Calculator - indmachzth: Induction Machine Thevenin Impedance Calculator - indmachpem: Induction Machine Electro-Mechanical Power Calculator - indmachtem: Induction Machine Electro-Mechanical Torque Calculator - indmachpkslip: Induction Machine Peak Slip Calculator - indmachpktorq: Induction Machine Peak Torque Calculator - indmachstarttorq: Induction Machine Starting Torque Calculator """ - # Condition Inputs - w = 2 * _np.pi * freq - if Ls is not None: # Use Ls instead of Lls - Lls = Ls - Lm - if Lr is not None: # Use Lr instead of Llr - Llr = Lr - Lm - if poles != 0: # Calculate Sync. Speed from Num. Poles - wsyn = w / (poles / 2) - if calcX: # Convert Inductances to Reactances - Lm *= w - Lls *= w - Llr *= w - # Test for Valid Input Set if Vth is None: - if not all((Vas, Rs, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Vth + if not all((Vas, Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Vth.") Vth = indmachvth(Vas, Rs, Lm, Lls, Ls, freq, calcX) + if Zth is None: - if not all((Rs, Llr, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Zth + if not all((Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Zth.") Zth = indmachzth(Rs, Lm, Lls, Llr, Ls, Lr, freq, calcX) - # Calculate Rotor Current - Iar = Vth / (Zth.real + Zth) + + den = abs(Zth) + Zth + if den == 0: + raise ValueError("Invalid operating point: |Zth| + Zth equals zero.") + Iar = Vth / den return Iar @@ -658,122 +300,25 @@ def indmachpktorq(Rr, poles=0, s_pk=None, Iar=None, Vth=None, Zth=None, Vas=0, calcX=True): r""" Induction Machine Peak Torque Calculator. - - Calculation function to find the peak torque for an - induction machine given the thevenin voltage and impedance. - - This current is calculated using the following formulas: - - .. math:: - T_{em}=(|I_{a_{\text{rotor}}}|)^2\cdot\frac{R_r} - {\text{slip}_{\text{peak}}} - - where: - - .. math:: I_{a_{\text{rotor}}} = \frac{V_{th}}{|Z_{th}|+Z_{th}} - - .. math:: V_{th}=\frac{j\omega L_m}{R_s+j\omega(L_{ls}+L_m)}V_{as} - - .. math:: - Z_{th} = \frac{(R_s+j\omega L_{ls})j\omega L_m} - {R_s+j\omega(L_{ls}+L_m)}+j\omega L_{lr} - - .. math:: \omega = \omega_{es} = 2\pi\cdot f_{\text{electric}} - - Parameters - ---------- - Rr: float - Rotor resistance in Ohms - poles: int, optional - Number of poles for the induction machine. - s_pk: float, optional - Peak induction machine slip, may be calculated - internally if remaining machine characteristics are - provided. - Iar: complex, optional - Phase-A, Rotor Current in Amps, may be calculated - internally if remaining machine characteristics are - provided. - Vth: complex, optional - Thevenin-equivalent stator voltage of the - induction machine, may be calculated internally - if given stator voltage and machine parameters. - Zth: complex, optional - Thevenin-equivalent inductance (in ohms) of the - induction machine, may be calculated internally - if given machine parameters. - Vas: complex, optional - Terminal Stator Voltage in Volts - Rs: float, optional - Stator resistance in ohms - Lm: float, optional - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Llr: float, optional - Rotor leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - Lr: float, optional - Rotor inductance in Henrys - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - Tpk: float - Peak torque of specified induction machine in - newton-meters. - - See Also - -------- - indmachvth: Induction Machine Thevenin Voltage Calculator - indmachzth: Induction Machine Thevenin Impedance Calculator - indmachpem: Induction Machine Electro-Mechanical Power Calculator - indmachtem: Induction Machine Electro-Mechanical Torque Calculator - indmachpkslip: Induction Machine Peak Slip Calculator - indmachiar: Induction Machine Phase-A Rotor Current Calculator - indmachstarttorq: Induction Machine Starting Torque Calculator """ - # Condition Inputs - w = 2 * _np.pi * freq - if Ls is not None: # Use Ls instead of Lls - Lls = Ls - Lm - if Lr is not None: # Use Lr instead of Llr - Llr = Lr - Lm - if poles != 0: # Calculate Sync. Speed from Num. Poles - wsyn = w / (poles / 2) - if calcX: # Convert Inductances to Reactances - Lm *= w - Lls *= w - Llr *= w - # Test for Valid Input Set if Vth is None: - if not all((Vas, Rs, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Vth + if not all((Vas, Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Vth.") Vth = indmachvth(Vas, Rs, Lm, Lls, Ls, freq, calcX) + if Zth is None: - if not all((Rs, Llr, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Zth + if not all((Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Zth.") Zth = indmachzth(Rs, Lm, Lls, Llr, Ls, Lr, freq, calcX) + if Iar is None: - if not all((Vth, Zth)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Ias Iar = indmachiar(Vth=Vth, Zth=Zth) + if s_pk is None: - if not all((Rr, Zth)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Peak Slip s_pk = indmachpkslip(Rr=Rr, Zth=Zth) - # Use Terms to Calculate Peak Torque + + if s_pk == 0: + raise ValueError("s_pk must be non-zero.") Tpk = abs(Iar) ** 2 * Rr / s_pk return Tpk @@ -783,118 +328,22 @@ def indmachstarttorq(Rr, poles=0, Iar=None, Vth=None, Zth=None, Vas=0, Rs=0, Lm=0, Lls=0, Llr=0, Ls=None, Lr=None, freq=60, calcX=True): r""" Induction Machine Starting Torque Calculator. - - Calculation function to find the starting torque for an - induction machine given the thevenin voltage and impedance. - - This current is calculated using the following formulas: - - .. math:: - T_{em}=(|I_{a_{\text{rotor}}}|)^2\cdot\frac{R_r} - {\text{slip}_{\text{peak}}} - - where: - - .. math:: \text{slip} = 1 - - .. math:: - I_{a_{\text{rotor}}} = \frac{V_{th}}{\frac{R_r}{\text{slip}}+Z_{th}} - - .. math:: V_{th}=\frac{j\omega L_m}{R_s+j\omega(L_{ls}+L_m)}V_{as} - - .. math:: - Z_{th} = \frac{(R_s+j\omega L_{ls})j\omega L_m} - {R_s+j\omega(L_{ls}+L_m)}+j\omega L_{lr} - - .. math:: \omega = \omega_{es} = 2\pi\cdot f_{\text{electric}} - - Parameters - ---------- - Rr: float - Rotor resistance in Ohms - poles: int, optional - Number of poles for the induction machine. - Iar: complex, optional - Phase-A, Rotor Current in Amps, may be calculated - internally if remaining machine characteristics are - provided. - Vth: complex, optional - Thevenin-equivalent stator voltage of the - induction machine, may be calculated internally - if given stator voltage and machine parameters. - Zth: complex, optional - Thevenin-equivalent inductance (in ohms) of the - induction machine, may be calculated internally - if given machine parameters. - Vas: complex, optional - Terminal Stator Voltage in Volts - Rs: float, optional - Stator resistance in ohms - Lm: float, optional - Magnetizing inductance in Henrys - Lls: float, optional - Stator leakage inductance in Henrys, default=0 - Llr: float, optional - Rotor leakage inductance in Henrys, default=0 - Ls: float, optional - Stator inductance in Henrys - Lr: float, optional - Rotor inductance in Henrys - freq: float, optional - System (electrical) frequency in Hz, default=60 - calcX: bool, optional - Control argument to force system to calculate - system reactances with system frequency, or to - treat them as previously-calculated reactances. - default=True - - Returns - ------- - Tstart: float - Peak torque of specified induction machine in - newton-meters. - - See Also - -------- - indmachvth: Induction Machine Thevenin Voltage Calculator - indmachzth: Induction Machine Thevenin Impedance Calculator - indmachpem: Induction Machine Electro-Mechanical Power Calculator - indmachtem: Induction Machine Electro-Mechanical Torque Calculator - indmachpkslip: Induction Machine Peak Slip Calculator - indmachpktorq: Induction Machine Peak Torque Calculator - indmachiar: Induction Machine Phase-A Rotor Current Calculator """ - # Condition Inputs - w = 2 * _np.pi * freq - if Ls is not None: # Use Ls instead of Lls - Lls = Ls - Lm - if Lr is not None: # Use Lr instead of Llr - Llr = Lr - Lm - if poles != 0: # Calculate Sync. Speed from Num. Poles - wsyn = w / (poles / 2) - if calcX: # Convert Inductances to Reactances - Lm *= w - Lls *= w - Llr *= w - # Slip is 1 (one) for starting slip = 1 - # Test for Valid Input Set + if Vth is None: - if not all((Vas, Rs, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Vth + if not all((Vas, Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Vth.") Vth = indmachvth(Vas, Rs, Lm, Lls, Ls, freq, calcX) + if Zth is None: - if not all((Rs, Llr, Lm, Lls)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Zth + if not all((Rs, Lm)): + raise ValueError("Invalid Argument Set, too few provided to compute Zth.") Zth = indmachzth(Rs, Lm, Lls, Llr, Ls, Lr, freq, calcX) + if Iar is None: - if not all((Vth, Zth)): - raise ValueError("Invalid Argument Set, too few provided.") - # Valid Argument Set, Calculate Ias Iar = Vth / (Rr / slip + Zth) - # Use Terms to Calculate Peak Torque + Tstart = abs(Iar) ** 2 * Rr / slip return Tstart @@ -904,29 +353,10 @@ def pstator(Pem, slip): r""" Stator Power Calculator for Induction Machine. - Given the electromechanical power and the slip, - this function will calculate the power related to the - stator (provided or consumed). - .. math:: P_s=\frac{P_{em}}{1-\text{slip}} - - Parameters - ---------- - Pem: float - Electromechanical power in watts. - slip: float - Slip factor in rad/sec. - - Returns - ------- - Ps: float - Power related to the stator in watts. - - See Also - -------- - protor: Rotor Power Calculator for Induction Machines """ - # Calculate and Return + if (1 - slip) == 0: + raise ValueError("Invalid slip: 1 - slip equals zero.") Ps = Pem / (1 - slip) return Ps @@ -936,29 +366,10 @@ def protor(Pem, slip): r""" Rotor Power Calculator for Induction Machine. - Given the electromechanical power and the slip, - this function will calculate the power related to the - rotor (provided or consumed). - .. math:: P_r=-\text{slip}\cdot\frac{P_{em}}{1-\text{slip}} - - Parameters - ---------- - Pem: float - Electromechanical power in watts. - slip: float - Slip factor in rad/sec. - - Returns - ------- - Pr: float - Power related to the rotor in watts. - - See Also - -------- - pstator: Stator Power Calculator for Induction Machines """ - # Calculate and Return + if (1 - slip) == 0: + raise ValueError("Invalid slip: 1 - slip equals zero.") Pr = -slip * (Pem / (1 - slip)) return Pr @@ -968,54 +379,6 @@ def indmachfocratings(Rr, Rs, Lm, Llr=0, Lls=0, Lr=None, Ls=None, Vdqs=1, Tem=1, wes=1): r""" FOC Ind. Machine Rated Operation Calculator. - - Determines the parameters and characteristics of a Field- - Oriented-Controlled Induction Machine operating at its - rated limits. - - Parameters - ---------- - Rr: float - Rotor resistance in per-unit-ohms - Rs: float - Stator resistance in per-unit-ohms - Lm: float - Magnetizing inductance in per-unit-Henrys - Llr: float, optional - Rotor leakage inductance in per-unit-Henrys, - default=0 - Lls: float, optional - Stator leakage inductance in per-unit-Henrys, - default=0 - Lr: float, optional - Rotor inductance in per-unit-Henrys - Ls: float, optional - Stator inductance in per-unit-Henrys - Vdqs: complex, optional - The combined DQ-axis voltage required for rated - operation, in per-unit-volts, default=1+j0 - Tem: float, optional - The mechanical torque required for rated operation, - in per-unit-newton-meters, default=1 - wes: float, optional - The per-unit electrical system frequency, default=1 - - Returns - ------- - Idqr: complex - Combined DQ-axis Rotor Current in per-unit-amps - Idqs: complex - Combined DQ-axis Stator Current in per-unit-amps - LAMdqr: complex - Combined DQ-axis Rotor Flux in per-unit - LAMdqs: complex - Combined DQ-axis Stator Flux in per-unit - slip_rat: float - Rated Slip as percent of rotational and system frequencies - w_rat: float - Rated System frequency in per-unit-rad/sec - lamdr_rat: float - Rated D-axis rotor flux in per-unit """ # Condition Inputs: if Ls is None: # Use Lls instead of Ls @@ -1026,10 +389,14 @@ def indmachfocratings(Rr, Rs, Lm, Llr=0, Lls=0, Lr=None, # Define Equations Function as Solver def equations(val): Idr, Iqr, Ids, Iqs, LAMdr, LAMqr, LAMds, LAMqs, wr = val - A = (Rs * Ids - wes * LAMqs) - Vdqs - B = Rs * Iqs - wes * LAMds - C = Rr * Idr - (wes - wr) * LAMqr - D = Rr * Iqr + (wes - wr) * LAMdr + # Force Vdqs to be treated as complex: Vd (real) + j Vq (imag) + Vds_cmd = _np.real(Vdqs) + Vqs_cmd = _np.imag(Vdqs) + + A = (Rs * Ids - wes * LAMqs) - Vds_cmd + B = (Rs * Iqs + wes * LAMds) - Vqs_cmd + C = (Rr * Idr - (wes - wr) * LAMqr) + D = (Rr * Iqr + (wes - wr) * LAMdr) E = (Ls * Ids + Lm * Idr) - LAMds F = (Ls * Iqs + Lm * Iqr) - LAMqs G = (Lm * Ids + Lr * Idr) - LAMdr @@ -1047,10 +414,14 @@ def equations(val): LAMds0 = Ls * Ids0 + Lm * Idr0 LAMqs0 = Ls * Iqs0 + Lm * Iqr0 wr = 1 + # Use Iterative Solver to Find Results Idr, Iqr, Ids, Iqs, LAMdr, LAMqr, LAMds, LAMqs, wr = _fsolve(equations, ( Idr0, Iqr0, Ids0, Iqs0, LAMdr0, LAMqr0, LAMds0, LAMqs0, wr)) + # Calculate Remaining Rating Terms + if wes == 0: + raise ValueError("wes must be non-zero.") slip_rated = (wes - wr) / wes w_rated = wr lamdr_rated = abs(LAMdr + 1j * LAMqr) @@ -1070,72 +441,38 @@ def imfoc_control(Tem_cmd, LAMdr_cmd, wr_cmd, Rr, Rs, Lm, Llr=0, Lls=0, Lr=None, Ls=None, s_err=0): """ FOC Ind. Machine Rated Operation Calculator. - - Determines the parameters and characteristics of a Field- - Oriented-Controlled Induction Machine operating at its - rated limits. - - Parameters - ---------- - Tem_cmd: float - Mechanical torque setpoint in per-unit-newton-meters - LAMdr_cmd: float - D-axis flux setpoint in per-unit - wr_cmd: float - Mechanical (rotor) speed in per-unit-rad/sec - Rr: float - Rotor resistance in per-unit-ohms - Rs: float - Stator resistance in per-unit-ohms - Lm: float - Magnetizing inductance in per-unit-Henrys - Llr: float, optional - Rotor leakage inductance in per-unit-Henrys, - default=0 - Lls: float, optional - Stator leakage inductance in per-unit-Henrys, - default=0 - Lr: float, optional - Rotor inductance in per-unit-Henrys - Ls: float, optional - Stator inductance in per-unit-Henrys - s_err: float, optional - Error in slip calculation as a percent (e.g. 0.25), - default=0 - - Returns - ------- - Vdqs: complex - Combined DQ-axis Stator Voltage in per-unit volts - Idqr: complex - Combined DQ-axis Rotor Current in per-unit-amps - Idqs: complex - Combined DQ-axis Stator Current in per-unit-amps - LAMdqr: complex - Combined DQ-axis Rotor Flux in per-unit - LAMdqs: complex - Combined DQ-axis Stator Flux in per-unit - wslip: float - Machine Slip frequency in per-unit-rad/sec - wes: float - The electrical system frequency in per-unit-rad/sec """ # Condition Inputs: if Ls is None: # Use Lls instead of Ls Ls = Lls + Lm if Lr is None: # Use Llr instead of Lr Lr = Llr + Lm + + if Ls == 0 or Lr == 0: + raise ValueError("Ls and Lr must be non-zero.") + # Calculate Additional Constraints sigma = (1 - Lm ** 2 / (Ls * Lr)) accuracy = 1 + s_err + if accuracy == 0: + raise ValueError("Invalid s_err: 1 + s_err equals zero.") + # Command Values (Transient and Steady State) + if Lm == 0: + raise ValueError("Lm must be non-zero.") Ids = LAMdr_cmd / Lm + + if LAMdr_cmd == 0: + raise ValueError("LAMdr_cmd must be non-zero to compute Iqs and slip.") Iqs = Tem_cmd / ((Lm / Lr) * LAMdr_cmd) - wslip = Rr / (Lr * accuracy) * (Lm * Iqs) / LAMdr_cmd + + wslip = (Rr / (Lr * accuracy)) * (Lm * Iqs) / LAMdr_cmd wes = wslip + wr_cmd + # Stator dq Voltages (Steady State) Vds = Rs * Ids - wes * sigma * Ls * Iqs - Vqs = Rs * Iqs - wes * Ls * Ids + Vqs = Rs * Iqs + wes * Ls * Ids + # Remaining Steady State Iqr = -Lm / Lr * Iqs Idr = 0 @@ -1158,55 +495,27 @@ def synmach_Eq(Vt_pu, Itmag, PF, Ra, Xd, Xq): # noqa: D401 "Synchronous" is an intentional descriptor r""" Synchronous Machine Eq Calculator. - - Given specified parameter set, will calculate - the internal voltage on the q-axis (Eq). - - .. math:: E_q=V_{t_{pu}}-\left[R_a\cdot I_{t_{pu}}+ - j\cdot X_q\cdot I_{t_{pu}}+j(X_d-X_q)\cdot I_{ad}\right] - - where: - - .. math:: I_{t_{pu}}=I_{t_{mag}}\cdot e^{-j( - \angle{V_{t_{pu}}}-\cos^{-1}(PF))} - - .. math:: \theta_q=\angle{V_{t_{pu}}-\left(R_a - I_{t_{pu}}+j\cdot X_qI_{t_{pu}}\right) - - .. math:: I_{ad}=\left|I_{t_{pu}}\cdot\sin( - -\cos^{-1}(PF)+\theta_q)\right|e^{j(\theta_q - -90°)} - - Parameters - ---------- - Vt_pu: complex - Terminal voltage in per-unit-volts - Itmag: float - Terminal current magnitude in per- - unit-amps - PF: float - Machine Power Factor, (+)ive values denote - leading power factor, (-)ive values denote - lagging power factor - Ra: float - AC resistance in per-unit-ohms - Xd: float - D-axis reactance in per-unit-ohms - Xq: float - Q-axis reactance in per-unit-ohms - - Returns - ------- - Eq: complex - Internal Synchronous Machine Voltage - in per-unit-volts """ - # Calculate Required Terms - phi = _np.arccos(PF) + PF = _validate_pf(PF) Itmag = abs(Itmag) - It_pu = Itmag * _np.exp(-1j * (_np.angle(Vt_pu) + phi)) + + # PF sign convention (as stated in docstring): + # (+) leading, (-) lagging + phi = _np.arccos(abs(PF)) + phi_signed = -phi if PF >= 0 else phi + + # Current angle relative to terminal voltage + angV = _np.angle(Vt_pu) + angI = angV - phi_signed # lagging: angI = angV - phi; leading: angI = angV + phi + It_pu = Itmag * _np.exp(1j * angI) + + # Quadrature-axis angle th_q = _np.angle(Vt_pu - (Ra * It_pu + 1j * Xq * It_pu)) - Iad = (abs(It_pu) * _np.sin(phi + th_q)) * _np.exp(1j * (th_q - _np.pi / 2)) + + # Approximate direct-axis component current phasor (per provided formulation intent) + Iad_mag = abs(It_pu) * abs(_np.sin(-phi_signed + th_q)) + Iad = Iad_mag * _np.exp(1j * (th_q - _np.pi / 2)) + # Calculate Eq Eq = Vt_pu - (Ra * It_pu + 1j * Xq * It_pu + 1j * (Xd - Xq) * Iad) return Eq diff --git a/electricpy/sim.py b/electricpy/sim.py index f5b44618..d51119b3 100644 --- a/electricpy/sim.py +++ b/electricpy/sim.py @@ -93,6 +93,9 @@ def digifiltersim(fin, filter, freqs, NN=1000, dt=0.01, title="", for k in range(NN): x[k] = fin(k * dt, freq) + # Ensure filter is a numpy array for .size and .shape usage + filter = _np.asarray(filter) + # Identify how many rows were provided sz = filter.size if (sz < 5): @@ -100,8 +103,12 @@ def digifiltersim(fin, filter, freqs, NN=1000, dt=0.01, title="", "Refer to documentation for proper format.") elif (sz == 5): rows = 1 + filter = _np.reshape(filter, (1, 5)) else: rows, cols = filter.shape + if cols != 5: + raise ValueError("ERROR: Invalid filter shape, expected Nx5.") + # Operate with each individual filter set x_tmp = _np.copy(x) nsteps = NN - 4 @@ -120,6 +127,10 @@ def digifiltersim(fin, filter, freqs, NN=1000, dt=0.01, title="", B0 * x_tmp[T] + B1 * x_tmp[T - 1] + B2 * x_tmp[T - 2]) # Copy New output into temporary input x_tmp = _np.copy(y) + + # Reset y for next cascade stage to avoid carryover from prior stage + y = _np.zeros(NN) + # Copy finalized output into *ytime* for plotting ytime = _np.copy(x_tmp) # Plot Filtered Output @@ -193,9 +204,8 @@ def step_response(system, npts=1000, dt=0.01, combine=True, xlim=False, for i in range(npts): step[i] = 1.0 - # Simulate Response for each input (step, ramp, parabola) - # All 'x' values are variables that are considered don't-care - x, y1, x = _sig.lsim((system), step, TT) + # Simulate Response + t_out, y1, x = _sig.lsim((system), step, TT) # Calculate error over all points for k in range(npts): @@ -278,9 +288,8 @@ def ramp_response(system, npts=1000, dt=0.01, combine=True, xlim=False, for i in range(npts): ramp[i] = (dt * i) - # Simulate Response for each input (step, ramp, parabola) - # All 'x' values are variables that are considered don't-care - x, y2, x = _sig.lsim((system), ramp, TT) + # Simulate Response + t_out, y2, x = _sig.lsim((system), ramp, TT) # Calculate error over all points for k in range(npts): @@ -362,9 +371,8 @@ def parabolic_response(system, npts=1000, dt=0.01, combine=True, xlim=False, for i in range(npts): parabola[i] = (dt * i) ** (2) - # Simulate Response for each input (step, ramp, parabola) - # All 'x' values are variables that are considered don't-care - x, y3, x = _sig.lsim((system), parabola, TT) + # Simulate Response + t_out, y3, x = _sig.lsim((system), parabola, TT) # Calculate error over all points for k in range(npts): @@ -499,6 +507,12 @@ def func_c(self, x): # Concatenated Function "simpts=" + str(simpts), " Autocorrecting simpts to be NN-1.") simpts = NN - 1 + # Normalize plotting flags: default to False if not explicitly True + if plotforcing is None: + plotforcing = False + if plotresult is None: + plotresult = False + # Test for C and D matricies if isinstance(C, typetest) and isinstance(D, typetest): solution = 3 # Set to solve and plot complete output @@ -520,7 +534,7 @@ def func_c(self, x): # Concatenated Function if callable(func): # if f is a function, test as one mF = func(1) # f should return: int, float, tuple, _np.arr, _np.matrix elif isinstance(func, (tuple, list)): # if f is tupple of arguments - if callable(func[0]): # if first argument is a function + if len(func) > 0 and callable(func[0]): # if first argument is a function c_funcs = c_func_concat(func) # concatinate functions into one mF = "MultiFunctions" # label as multiple concatenated functions else: @@ -564,6 +578,23 @@ def func_c(self, x): # Concatenated Function "\n or numpy.matrixlib.defmatrix.matrix. Nor does function " + "\ncontain tuple of function handles. Please review function.") + # Normalize forcing function to always return a column matrix for consistent math + def _fn_col_matrix(T): + y = fn(T) + if isinstance(y, _np.matrixlib.defmatrix.matrix): + m = _np.asmatrix(y) + elif isinstance(y, _np.ndarray): + m = _np.asmatrix(y) + elif isinstance(y, tuple): + m = _np.asmatrix(y) + elif isinstance(y, (int, float, _np.float64)): + m = _np.asmatrix([y]) + else: + m = _np.asmatrix(y) + if m.shape[1] != 1: + m = _np.matrix.reshape(m, (m.size, 1)) + return m + # Test for size correlation between matricies if (cA != rA): # A isn't nxn matrix raise ValueError("Matrix 'A' is not NxN matrix.") @@ -571,22 +602,25 @@ def func_c(self, x): # Concatenated Function if (B.size % rA) == 0: # Elements in B divisible by rows in A _warn("WARNING: Reshaping 'B' matrix to match 'A' matrix.") B = _np.matrix.reshape(B, (rA, int(B.size / rA))) # Reshape Matrix + rB, cB = B.shape else: raise ValueError("'A' matrix dimensions don't match 'B' matrix dimensions.") elif (rA != rx): # A and x matricies don't have same number of rows if (x.size % rA) == 0: # Elements in x divisible by rows in A _warn("WARNING: Reshaping 'x' matrix to match 'A' matrix.") x = _np.matrix.reshape(x, (rA, 1)) # Reshape Matrix + rx, cx = x.shape else: raise ValueError("'A' matrix dimensions don't match 'B' matrix dimensions.") elif (cB != rF) or (cF != 1): # Forcing Function matrix doesn't match B matrix raise ValueError("'B' matrix dimensions don't match forcing function dimensions.") - elif (solution == 3) and (cC != cA) or (rC != 1): # Number of elements in C don't meet requirements + elif ((solution == 3) and ((cC != cA) or (rC != 1))): # Number of elements in C don't meet requirements raise ValueError("'C' matrix dimensions don't match state-space variable dimensions.") - elif (solution == 3) and ((cD != rF) or (rD != 1)): # Number of elements in D don't meet requirements + elif ((solution == 3) and ((cD != rF) or (rD != 1))): # Number of elements in D don't meet requirements if (cD == rD) and (cD == 1) and (D[0] == 0): # D matrix is set to [0] D = _np.asmatrix(_np.zeros(rF)) # Re-create D to meet requirements _warn("WARNING: Autogenerating 'D' matrix of zeros to match forcing functions.") + rD, cD = D.shape else: raise ValueError("'D' matrix dimensions don't match forcing function dimensions.") @@ -597,7 +631,7 @@ def func_c(self, x): # Concatenated Function # Start by defining Constants T = 0 TT = _np.arange(0, (dt * (NN)), dt) - yout = 0 + yout = _np.zeros(NN) # Define list of strings for plot output soltype = ["(Zero-Input)", "(Zero-State)", "(Complete Simulation)", "(Complete Sim., Combined Output)"] @@ -631,20 +665,28 @@ def func_c(self, x): # Concatenated Function for i in range(0, simpts): for n in range(xtim_len): xtim[n][i] = x[n] # xtim[state-variable][domain] = x[state-variable] + # Create Forcing Function output + if solution == 0: + fcol = None + else: + fcol = _fn_col_matrix(T) if fnc > 1: # More than one forcing function for n in range(fnc): - fn_arr[n][i] = _np.asarray(fn(T))[n][0] + fn_arr[n][i] = _np.asarray(fcol)[n][0] else: # only one forcing function - fn_arr[i] = fn(T) + if solution != 0: + fn_arr[i] = _np.asarray(fcol)[0][0] + else: + fn_arr[i] = 0 if solution == 0: # Zero-input, no added function input x = x + dt * A * x else: # Zero-state or Total, add function input - x = x + dt * A * x + dt * B * fn(T) + x = x + dt * A * x + dt * B * fcol if solution == 3: - yout = yout + dt * D * fn(T) + yout[i] = _np.asarray((D * fcol))[0][0] T = T + dt # Add discrete increment to T @@ -690,10 +732,11 @@ def func_c(self, x): # Concatenated Function if (plotresult and solution == 3): cofig = _plt.figure("Combined Output") C = _np.asarray(C) # convert back to array for operation + ysum = _np.zeros(NN) for i in range(cC): - yout = yout + xtim[i] * C[0][i] # Sum all st-space var mult. by their coeff - yout = _np.asarray(yout) # convert output to array for plotting purposes - _plt.plot(TT, yout[0]) + ysum = ysum + xtim[i] * C[0][i] # Sum all st-space var mult. by their coeff + ysum = ysum + yout + _plt.plot(TT, ysum) if xlim != False: _plt.xlim(xlim) if ylim != False: @@ -750,43 +793,26 @@ def NewtonRaphson(F, J, X0, eps=1e-4, mxiter=100, lsq_eps=0.25): The number of iterations completed before returning either due to solution being found, or max iterations being surpassed. - - Examples - -------- - >>> # doctest: +SKIP - >>> import numpy as np - >>> from electricpy import sim # Import Simulation Submodule - >>> def F(x): - ... matr = np.array([[x[1]*10*np.sin(x[0])+2], - ... [x[1]*(-10)*np.cos(x[0])+x[1]**2*10+1]]) - ... return(matr) - >>> def J(x): - ... matr = np.array([[10*x[1]*np.cos(x[0]), 10*np.sin(x[0])], - ... [10*x[1]*np.sin(x[0]), -10*np.cos(x[0])+20*x[1]]]) - ... return(matr) - >>> # Now compute Newton-Raphson - >>> X0 = [0, 1] - >>> results, iter = sim.NewtonRaphson(F,J,X0) - >>> print(results) - [-0.236,0.8554] - >>> print(iter) # Iteration Counter - 4 - - See Also - -------- - nr_pq: Newton-Raphson System Generator - mbuspowerflow: Multi-Bus Power Flow Calculator - """ + # Ensure X0 is numpy array for vector math + X0 = _np.asarray(X0, dtype=_np.float64) + # Test for one-variable inputs if isinstance(F(X0), (int, float, _np.float64)): # System is size-1 if not isinstance(J(X0), (int, float, _np.float64)): # Jacobian isn't size-1 raise ValueError("ERROR: The Jacobian isn't size-1.") return newton(F, X0, J) + # Condition outputs to arrays + def _asvec(v): + v = _np.asarray(v, dtype=_np.float64) + return v.reshape(-1) + # Test for valid argument sizes - f0sz = len(F(X0)) - j0sz = len(J(X0)) + F_value = _asvec(F(X0)) + J_value = _np.asarray(J(X0), dtype=_np.float64) + f0sz = F_value.size + j0sz = J_value.shape[0] if f0sz != j0sz: # Size mismatch raise ValueError("ERROR: The arguments return arrays or lists" + " of different sizes: f0=" + str(f0sz) + "; j0=" + str(j0sz)) @@ -799,27 +825,27 @@ def inv(m): i = _np.eye(a, a) return _np.linalg.lstsq(m, i, rcond=None)[0] - F_value = F(X0) F_norm = _np.linalg.norm(F_value, ord=2) # L2 norm of vector iteration_counter = 0 teps = eps + userhasbeenwarned = False + while abs(F_norm) > teps and iteration_counter < mxiter: try: # Try Solve Operation - delta = _np.linalg.solve(J(X0), -F_value) + delta = _np.linalg.solve(_np.asarray(J(X0), dtype=_np.float64), -F_value) teps = eps # Change Test Epsilon if Needed except _np.linalg.LinAlgError: # Use Least Square if Error # Warn User, but only Once - try: - tst = userhasbeenwarned - except NameError: + if not userhasbeenwarned: userhasbeenwarned = True _warn("WARNING: Singular matrix, attempting LSQ method.") # Calculate Delta Using Least-Squares Inverse - delta = - inv(J(X0)).dot(F_value) + delta = - inv(_np.asarray(J(X0), dtype=_np.float64)).dot(F_value) # Change Epsilon Test teps = lsq_eps + X0 = X0 + delta - F_value = F(X0) + F_value = _asvec(F(X0)) F_norm = _np.linalg.norm(F_value, ord=2) iteration_counter += 1 @@ -836,85 +862,6 @@ def inv(m): def nr_pq(Ybus, V_set, P_set, Q_set, extend=True, argshape=False, verbose=False): """ Newton Raphson Real/Reactive Power Function Generator. - - Given specified parameters, will generate the necessary real and reactive - power functions necessary to compute the system's power flow. - - Parameters - ---------- - Ybus: array_like - Postitive Sequence Y-Bus Matrix for Network. - V_set: list of list of float - List of known and unknown voltages. - Known voltages should be provided as - a list of floating values in the form of: - [mag, ang], unknown voltages should - be provided as None. - P_set: list of float - List of known and unknown real powers. - Known powers should be provided as - a floating point value, unknown powers - should be provided as None. Generated power - should be noted as positive, consumed power - should be noted as negative. - Q_set: list of float - List of known and unknown reactive powers. - Known powers should be provided as - a floating point value, unknown powers - should be provided as None. Generated power - should be noted as positive, consumed power - should be noted as negative. - extend: bool, optional - Control argument to format returned value as - singular function handle or lists of function - handles. default=True; singular function - argshape: bool, optional - Control argument to force return of the voltage - argument array as a tuple of: (Θ-len, V-len). - default=False - verbose: bool, optional - Control argument to print verbose information - about function generation, useful for debugging. - default=False - - Returns - ------- - retset: array_like - An array of function handles corresponding to - the appropriate voltage magnitude and angle - calculation functions based on the real and - reactive power values. Function(s) will accept - an argument of the form: - [Θ1, Θ2,..., Θn, V1, V2,..., Vm] - where n is the number of busses with unknown - voltage angle, and m is the number of busses - with unknown voltage magnitude. - - Examples - -------- - >>> # doctest: +SKIP - >>> import numpy as np - >>> from electricpy import sim # Import Simulation Submodule - >>> ybustest = [[-10j,10j], - ... [10j,-10j]] - >>> Vlist = [[1,0],[None,None]] # We don't know the voltage or angle at bus 2 - >>> Plist = [None,-2.0] # 2pu watts consumed - >>> Qlist = [None,-1.0] # 1pu vars consumed - >>> F = nr_pq(ybustest,Vlist,Plist,Qlist) - >>> X0 = [0,1] # Define Initial Conditions - >>> J = sim.jacobian(F) # Find Jacobian - >>> # Now use Newton-Raphson to Solve - >>> results, iter = sim.NewtonRaphson(F,J,X0) - >>> print(results) - [-0.236,0.8554] - >>> print(iter) # Iteration Counter - 4 - - See Also - -------- - NewtonRaphson: Newton-Raphson System Solver - mbuspowerflow: Multi-Bus Power Flow Calculator - """ # Condition Inputs Ybus = _np.asarray(Ybus) @@ -1081,103 +1028,6 @@ def mbuspowerflow(Ybus, Vknown, Pknown, Qknown, X0='flatstart', eps=1e-4, slackbus=0, lsq_eps=0.25): """ Multi-Bus Power Flow Calculator. - - Function wrapper to simplify the evaluation of a power flow calculation. - Determines the function array (F) and the Jacobian array (J) and uses the - Newton-Raphson method to iteratively evaluate the system to converge to a - solution. - - Parameters - ---------- - Ybus: array_like - Postitive Sequence Y-Bus Matrix for Network. - Vknown: list of list of float - List of known and unknown voltages. - Known voltages should be provided as - a list of floating values in the form of: - [mag, ang], unknown voltages should - be provided as None. - Pknown: list of float - List of known and unknown real powers. - Known powers should be provided as - a floating point value, unknown powers - should be provided as None. Generated power - should be noted as positive, consumed power - should be noted as negative. - Qknown: list of float - List of known and unknown reactive powers. - Known powers should be provided as - a floating point value, unknown powers - should be provided as None. Generated power - should be noted as positive, consumed power - should be noted as negative. - X0: {'flatstart', list of float}, optional - Initial conditions/Initial guess. May be set - to 'flatstart' to force function to generate - flat voltages and angles of 1∠0°. Must be - specified in the form: - [Θ1, Θ2,..., Θn, V1, V2,..., Vm] - where n is the number of busses with unknown - voltage angle, and m is the number of busses - with unknown voltage magnitude. - eps: float, optional - Epsilon - The error value, default=0.0001 - mxiter: int, optional - Maximum Iterations - The highest number of - iterations allowed, default=100 - returnct: bool, optional - Control argument to force function to return - the iteration counter from the Newton-Raphson - solution. default=False - degrees: bool, optional - Control argument to force returned angles to - degrees. default=True - split: bool, optional - Control argument to force returned array to - split into lists of magnitudes and angles. - default=False - slackbus: int, optional - Control argument to specify the bus index for - the slack bus. If the slack bus is not positioned - at bus index 1 (default), this control can be - used to reformat the data sets to a format - necessary for proper generation and Newton - Raphson computation. Must be zero-based. - default=0 - lsq_eps: float, optional - Least Squares Method (Failover) Epsilon - the error value. - default=0.25 - - - .. image:: /static/mbuspowerflow_example.png - - Examples - -------- - >>> # doctest: +SKIP - >>> # Perform Power-Flow Analysis for Figure - >>> import numpy as np - >>> from electricpy import sim # Import Simulation Submodule - >>> ybustest = [[-10j,10j], - ... [10j,-10j]] - >>> Vlist = [[1,0],[None,None]] # We don't know the voltage or angle at bus 2 - >>> Plist = [None,-2.0] # 2pu watts consumed - >>> Qlist = [None,-1.0] # 1pu vars consumed - >>> sim.mbuspowerflow( - ... ybustest, - ... Vlist, - ... Plist, - ... Qlist, - ... degrees=True, - ... split=True, - ... returnct=True - ... ) - ([array([-13.52185223]), array([ 0.85537271])], 4) - - See Also - -------- - NewtonRaphson: Newton-Raphson System Solver - nr_pq: Newton-Raphson System Generator - electricpy.powerflow: Simple (2-bus) Power Flow Calculator """ # Identify Lack of Support if not __NUMDIFFTOOL_SUPPORT__: @@ -1198,21 +1048,28 @@ def mbuspowerflow(Ybus, Vknown, Pknown, Qknown, X0='flatstart', eps=1e-4, Qknown = _np.roll(Qknown, (len(Qknown) - slackbus), 0).tolist() # Generate F Function Array F, shp = nr_pq(Ybus, Vknown, Pknown, Qknown, True, True, False) + + ang_len, mag_len = shp + # Handle Flat-Start Condition if X0 == 'flatstart': - ang_len, mag_len = shp X0 = _np.append(_np.zeros(ang_len), _np.ones(mag_len)) + # Evaluate Jacobian J = jacobian(F) + # Compute Newton-Raphson nr_result, iter_count = NewtonRaphson(F, J, X0, eps, mxiter, lsq_eps) + # Convert to Degrees if Necessary if degrees: for i in range(ang_len): nr_result[i] = _np.degrees(nr_result[i]) + # Split into Mag/Ang Arrays if Necessary if split: nr_result = [nr_result[:ang_len], nr_result[-mag_len:]] + # Return with Iteration Counter if returnct: return nr_result, iter_count