line_design_v_vs_p

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
from functools import partial
import pandapower
import pandapower.networks as pn
import matplotlib.pyplot as plt

"""
Code to explore the question raised at
https://electronics.stackexchange.com/questions/654161/what-is-the-relationship-between-a-power-lines-voltage-and-the-amount-of-power/

Environment requirements:
matplotlib, numba, numpy, scipy, pandapower
"""


matpower_loader = pandapower.converter.from_mpc

"""
References for test cases:
https://pandapower.readthedocs.io/en/latest/networks/power_system_test_cases.html

https://egriddata.org/group/electric-grid-test-case-repository

"""

# https://egriddata.org/sites/default/files/case_ACTIVSg2000.m
ACTIVSg2000 = partial(matpower_loader, 'case_ACTIVSg2000.m')


# https://egriddata.org/sites/default/files/case_ACTIVSg70k.m
ACTIVSg70k = partial(matpower_loader, 'case_ACTIVSg70k.m')


case_loaders = [pn.case2869pegase,
         pn.case3120sp,
         # pn.GBnetwork,  # No rating information
         # pn.iceland,    # No rating information
         ACTIVSg2000,
         ACTIVSg70k]
#case_loaders = [ACTIVSg70k]


def extract_line_info(net):
    lines_w_ratings = net.line[net.line['max_i_ka'] > 0]
    df = pd.merge(lines_w_ratings, net.bus, left_on='from_bus', right_index=True)
    df = df[['max_i_ka', 'vn_kv']]
    # Some datasets use kA value of 99999 presumably for no info.
    df = df[(df['vn_kv'] > 0) & (df['max_i_ka'] < 99998)]
    df['max_mva'] = df['max_i_ka']*df['vn_kv']*1.732

    # Some datasets use MVA value of 9900 presumably for no info.
    df = df[df['max_mva'] < 9000]

    print(f"Min KA rating: {df['max_i_ka'].min()}")
    print(f"Max KA rating: {df['max_i_ka'].max()}")
    print(f"Min MVA rating: {df['max_mva'].min()}")
    print(f"Max MVA rating: {df['max_mva'].max()}")
    print(f"Number of lines: {len(df)}")
    return df


nets = []
for i, loader in enumerate(case_loaders):
    print(f'Loading case {i}')
    nets.append(extract_line_info(loader()))

data = pd.concat(nets)

plt.plot(data['vn_kv'], data['max_mva'], '.', label='Individual Lines')

print('Voltage levels:')
print(np.sort(data['vn_kv'].unique()))

data_grouped = data.groupby('vn_kv')

print(data_grouped.describe())

mean_by_voltage = data_grouped['max_mva'].mean()
plt.plot(mean_by_voltage.index, mean_by_voltage, '*', label='Mean by voltage')


def power_curve(x, m):
    return x**m


[m_fit], cov = curve_fit(power_curve, mean_by_voltage.index, mean_by_voltage)

print(f'Best fit for power curve: m={m_fit}')

plot_xs = np.linspace(mean_by_voltage.index.min(), mean_by_voltage.index.max())

plt.plot(plot_xs, power_curve(plot_xs, m_fit), label=f'Power curve, m={m_fit:.2f}')


"""
def exp_curve(x, a, b):
    return a*np.exp(b*x)

y1, y2 = mean_by_voltage.iloc[[0, -1]]
x1, x2 = mean_by_voltage.index[[0, -1]]
b = (np.log(y1) - np.log(y2)) / (x1 - x2)
a = (np.log(y1) + np.log(y2) - b*(x1 + x2)) / 2
p0 = (a, b)
params, cov = curve_fit(exp_curve, mean_by_voltage.index, mean_by_voltage, p0=p0)
plt.plot(plot_xs, exp_curve(plot_xs, *params))
"""

"""
def quadratic_curve(x, a, b, c):
    return a*x**2 + b*x + c


params, cov = curve_fit(quadratic_curve, mean_by_voltage.index, mean_by_voltage)

print(f'Best fit for quadratic curve: a={params[0]}, b={params[1]}, c={params[2]}')
plt.plot(plot_xs, quadratic_curve(plot_xs, *params))
"""

def quadratic_curve2(x, a, b):
    c = 0
    return a*x**2 + b*x + c


params, cov = curve_fit(quadratic_curve2, mean_by_voltage.index, mean_by_voltage)
(a, b), c = params, 0

print(f'Best fit for quadratic curve: a={a}, b={b}, c={c}')
plt.plot(plot_xs, quadratic_curve2(plot_xs, *params), label=f'Quadratic curve, a={a:.3f}, b={b:.2f}, c={c:.2f}')

plt.title('Overhead Transmission Lines')
plt.xlabel('Line Voltage (kV)')
plt.ylabel('Line Capacity (MVA)')
plt.grid(True)
plt.legend()
plt.show()