The Ultimate Power Systems Cheatsheet

# Plant capacity factor comparison
plants = {
    "Nuclear":     0.92,
    "Coal":        0.72,
    "Natural Gas": 0.55,
    "Hydro":       0.44,
    "Wind":        0.35,
    "Solar PV":    0.25,
}
for name, cf in plants.items():
    # Energy per year for 1000 MW plant
    energy_mwh = 1000 * cf * 8760
    print(f"{name:14s}  CF={cf:.0%}  Annual={energy_mwh:,.0f} MWh")

import numpy as np

# Medium transmission line (π-model) per unit length
z = complex(0.1, 0.4)        # Ω/km  (R + jX)
y = complex(0, 2.8e-6)       # S/km  (jωC)
length = 200                  # km

Z = z * length               # total series impedance
Y = y * length               # total shunt admittance

# π-model ABCD
A = 1 + Y * Z / 2
B = Z
C = Y * (1 + Y * Z / 4)
D = A

print(f"A = {A:.6f}")
print(f"B = {B:.4f} Ω")
print(f"C = {C*1e6:.4f} µS")
print(f"D = {D:.6f}")
print(f"|AD - BC| = {abs(A*D - B*C):.6f}  (should be 1 for reciprocal)")

import math

# 3-phase radial feeder voltage drop
V_LL = 11e3          # 11 kV line-to-line
I = 200              # Amps
R_per_km = 0.5       # Ω/km
X_per_km = 0.35      # Ω/km
length = 10          # km
pf = 0.85            # lagging
phi = math.acos(pf)

R = R_per_km * length
X = X_per_km * length
VD = math.sqrt(3) * I * (R * math.cos(phi) + X * math.sin(phi))
V_phase = V_LL / math.sqrt(3)

pct_drop = VD / V_phase * 100
print(f"Voltage drop = {VD:.1f} V ({pct_drop:.2f}% of phase voltage)")
print(f"End voltage  = {V_LL - VD/math.sqrt(3)*math.sqrt(3):.0f} V line-to-line")

import numpy as np

# Symmetrical component transformation
a = complex(-0.5, np.sqrt(3)/2)  # 120° operator
A_mat = np.array([
    [1, 1, 1],
    [1, a**2, a],
    [1, a, a**2],
], dtype=complex)

# Single line-to-ground fault (phase A) on 11 kV system
E = 11000 / np.sqrt(3)  # pre-fault voltage (V)
Z1 = complex(2, 6)       # +ve seq impedance (Ω)
Z2 = complex(2, 6)       # -ve seq impedance (Ω)
Z0 = complex(2, 10)      # zero seq impedance (Ω)

Ia1 = E / (Z1 + Z2 + Z0)  # SLG fault: sequences in series
Ia  = 3 * Ia1
print(f"Fault current = {abs(Ia)/1000:.2f} kA")

# Unbalanced phase voltages after fault
V_seq = np.array([
    E - Z0 * Ia1,
    E - Z1 * Ia1,
    E - Z2 * Ia1,
], dtype=complex)
V_phases = A_mat @ V_seq
print(f"Va = {abs(V_phases[0]):.0f} V  (should ≈ 0 for bolted fault)")
print(f"Vb = {abs(V_phases[1]):.0f} V")
print(f"Vc = {abs(V_phases[2]):.0f} V")

# IEC 60255 inverse time overcurrent relay
# t = TMS * k / ((I/Is)^alpha - 1)

def relay_time(TMS, I, Is, curve="normal_inv"):
    """Calculate relay operating time per IEC 60255."""
    curves = {
        "normal_inv":   (0.14, 0.02),
        "very_inv":     (13.5,  1.0),
        "extremely_inv":(80.0,  2.0),
        "long_inv":     (120.0, 2.0),
    }
    K, alpha = curves[curve]
    M = I / Is  # multiple of plug setting
    if M <= 1.0:
        return float('inf')
    return TMS * K / (M**alpha - 1)

# Example: feeder protection
TMS = 0.3
Is = 400  # A plug setting
fault_currents = [800, 1200, 2000, 4000]

for If in fault_currents:
    t = relay_time(TMS, If, Is, "normal_inv")
    print(f"I_f = {If:5d} A  →  t = {t:.3f} s")

# Two-part tariff calculation with PF penalty
energy_rate = 8.0      # Rs/kWh
demand_rate = 300.0    # Rs/kVA per month
pf_actual = 0.82
pf_target = 0.90

# Monthly consumption
monthly_kwh = 500_000
peak_kw = 1200
peak_kva = peak_kw / pf_actual

# Base charges
demand_charge = peak_kva * demand_rate
energy_charge = monthly_kwh * energy_rate

# PF penalty (reactive charge)
pf_penalty = peak_kw * (pf_target / pf_actual - 1) * demand_rate

total = demand_charge + energy_charge + pf_penalty
print(f"Demand charge (kVA billing): Rs {demand_charge:,.0f}")
print(f"Energy charge:               Rs {energy_charge:,.0f}")
print(f"PF penalty:                  Rs {pf_penalty:,.0f}")
print(f"Total monthly bill:          Rs {total:,.0f}")
print(f"Cost per kWh:                Rs {total/monthly_kwh:.2f}")

# If PF improved to 0.95 via capacitor bank:
new_kva = peak_kw / 0.95
new_total = new_kva * demand_rate + energy_charge
savings = total - new_total
print(f"\nWith PF = 0.95: savings = Rs {savings:,.0f}/month")

# HVDC power flow — LCC converter
Vd_ideal = 500e3  # rated DC voltage (V)
Id = 2000          # DC current (A)
alpha = 15          # firing angle (degrees)
gamma = 18          # extinction angle (degrees)
Xc = 15             # commutating reactance per bridge (Ω)
n_bridges = 2       # 12-pulse = 2 × 6-pulse

import math
alpha_r = math.radians(alpha)
gamma_r = math.radians(gamma)

# Rectifier
Vd_r = n_bridges * (3*math.sqrt(2)/math.pi * 33e3 * math.cos(alpha_r)
        - 3 * Xc / math.pi * Id)
# Inverter
Vd_i = n_bridges * (3*math.sqrt(2)/math.pi * 33e3 * math.cos(gamma_r)
        - 3 * Xc / math.pi * Id)

Pdc = Vd_r * Id
losses = (Vd_r - Vd_i) * Id
efficiency = Vd_i / Vd_r * 100

print(f"Rectifier V_dc = {Vd_r/1000:.1f} kV")
print(f"Inverter  V_dc = {Vd_i/1000:.1f} kV")
print(f"DC Power      = {Pdc/1e6:.0f} MW")
print(f"DC Losses     = {losses/1e6:.1f} MW ({losses/Pdc*100:.2f}%)")
print(f"Efficiency    = {efficiency:.2f}%")

import numpy as np

# Solar PV annual energy yield estimation
# Using NREL PVWatts simplified model
dc_rating_kw = 1000           # 1 MW DC array
dc_ac_ratio = 1.25            # panel-to-inverter ratio
losses = 0.14                 # 14% system losses (soiling, wiring, mismatch)
temp_coeff = -0.004           # -0.4%/°C

# Monthly solar irradiance (kWh/m²/day) — typical mid-latitude
monthly_ghi = [3.2, 4.1, 5.2, 6.3, 6.8, 6.5, 6.7, 6.1, 5.3, 4.0, 3.0, 2.8]

# Monthly average ambient temperature (°C)
monthly_temp = [5, 7, 12, 17, 22, 27, 30, 29, 24, 18, 11, 6]

print(f"{'Month':>8s}  {'GHI':>5s}  {'Tamb':>5s}  {'Temp Corr':>10s}  {'kWh':>8s}")
print("-" * 50)
total_kwh = 0
for i, (ghi, temp) in enumerate(zip(monthly_ghi, monthly_temp)):
    days = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31][i]
    # Temperature derating of STC output
    t_corr = 1 + temp_coeff * (temp - 25)
    t_corr = max(t_corr, 0.7)  # floor at 70%
    # Simplified energy: P_dc × hours × DC/AC × losses × temp
    energy = dc_rating_kw * ghi * days * (1/dc_ac_ratio) * (1 - losses) * t_corr
    total_kwh += energy
    print(f"{i+1:>8d}  {ghi:>5.1f}  {temp:>5.1f}  {t_corr:>10.3f}  {energy:>8.0f}")

cf = total_kwh / (dc_rating_kw * 8760)
print(f"
Annual yield:  {total_kwh:,.0f} kWh ({total_kwh/1000:,.1f} MWh)")
print(f"Capacity factor: {cf:.1%}")