tidal_simulation

"""
Power Generation Sim

Refactored version of previous simulation.
"""
from mpl_toolkits.mplot3d import Axes3D
import matplotlib
from matplotlib import cm
import matplotlib.pyplot as plt
import pandas as pd
import math as math
import numpy as np

#function to find index of nearest value; used to find bay height
#for a given volume
def find_nearest(array,value):
    idx = (np.abs(array-value)).argmin()
    return idx

#function to calculate flow based on orifice and head
def flow_calc(tide_h, bay_h, orifice_area):
    if(tide_h > bay_h):
        return abs((0.65 * math.sqrt(19.62 * (orifice_area ** 2) * (tide_h - bay_h))))
    return abs((-0.65 * (math.sqrt(19.62 * (orifice_area ** 2) * (bay_h - tide_h)))))

#function to calculate instantaneous tide height
def tide_height(time):
    tidal_constants = [[0.131, 2.3621, 0.00438],[0.092, 6.1107, 0.00406],[1.390, 3.4478, 0.00843],[0.545, 4.0361, 0.00873]]
    tide = 2.893
    for element in tidal_constants:
        tide += (element[0] * math.cos((element[2] * time) + element[1]))
    return tide

#function to calculate efficiency of turbine based on given curve
def efficiency_calc(flow_percent):
    #values required for turbine calculations
    eff_values = np.array([0.5, 0.73, 0.82, 0.88, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.86, 0.75, 0.5, 0.4])
    flow_eff_values = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4])
    return eff_values[find_nearest(flow_eff_values, flow_percent)]

#function to calculate instantaneous power of turbine using equation P = 1000*Q*H*\eta, turbine_on value is due to way the control logic works
def inst_P(flow, head, efficiency, turbine_on, pumping = 0):
    if(pumping):
        return (1000.0 * flow * head * 9.81 * efficiency * turbine_on)
    return abs(1000.0 * flow * head * 9.81 * efficiency * turbine_on)

#function to control sleiuces, turbines, returns 2-element array, first element is whether
#turbine is on or off, second element is orifice size i.e. sluice or turbine inlet
def ebb_logic(tide_height, cut_in, change, turbine_area):
    #if tide is coming in, open sleuice gates
    if (change > 0):
        return [0, 10]
    #if tide is going out but still beneath cut in height, close sluice
    elif (change < 0 and (tide_height > cut_in)):
        return [0, 0]
    #otherwise open turbine gates, turbine on
    else:
        return [1, turbine_area]

def flood_logic(tide_height, cut_in, change, turbine_area):
    #if tide going out, open sluice gates
    if (change < 0):
        return [0, 10]
    #if tide coming in but less than cut in height, close sluice
    elif (change > 0 and (tide_height < cut_in)):
        return [0, 0]
    #else open turbine gates, turbine on
    else:
        return [1, turbine_area]

def two_way_logic(tide_height, cut_in, change, turbine_area):
    #if tide going our and head low, close gates
    if ((change < 0) and (tide_height > cut_in)):
        return [0, 0]
    #if tide going out an below cut in value, turbine on
    elif ((change < 0) and (tide_height < cut_in)):
        return [1, turbine_area]
    #if tide coming in but less that cut in, close gates
    elif ((change > 0) and (tide_height < cut_in)):
        return [0, 0]
    #if tide coming in and above cut in value
    else:
        return [1, turbine_area]

#function to calculate total power generated in kWh for a given
#rated flow and turbine size
def power_value_total(tide_list, rated_flow, turbine_area, control_function):
    #initial values
    previous_tide = 0
    previous_bay = 0
    cut_in = 3.4
    #local vars
    bay_h = 0
    bay_volume = 0
    inst_P_list = []
    for tide in tide_list:
        #calculate gradient
        grad = tide - previous_tide
        #calculate orifice/turbine values
        control_signal = control_function(tide, cut_in, grad, turbine_area)
        turbine_on = control_signal[0]
        orifice_area = control_signal[1]
        #calculate flow rate
        flow_rate = flow_calc(tide, previous_bay, orifice_area)
        #add fluid to bay volume and find new bay height
        bay_volume += (flow_rate * 60)
        bay_h = (find_nearest(bay_vol_ref, bay_volume) / 100.0)


        head = abs(tide - bay_h)

        flow_percent = abs(flow_rate/rated_flow)
        efficiency = efficiency_calc(flow_percent)

        instantaneous_P = inst_P(flow_rate, head, efficiency, turbine_on)
        inst_P_list.append(instantaneous_P)
        previous_tide = tide
        previous_bay = bay_h
        total_energy = [power*60 for power in inst_P_list]

    return (sum(total_energy) / 3600000)

#function takes two arrays and brute forces the best configuration of flow/turbine sizes by calculating the total power output for all combinations, then returning the indices of the maximal optimal values - Note that this is a very inefficient method, with run times for a full month of data approaching 12 hours
def optimisation(turbine_array, flow_array, control_logic):
    #initialise C style array
    power_array = np.empty([len(turbine_array),len(flow_array)])
    #calculate power output for all values of turbine size and flow rate
    for t_idx, turbine_area in enumerate(turbine_array):
        for r_idx, rated_flow in enumerate(flow_array):
            power_array[t_idx][r_idx] = power_value_total(tide_list, rated_flow, turbine_area, control_logic)
    #find maximum power value stored in array
    max_power = np.amax(power_array)
    #find index of maximum value
    optimum_tup = np.unravel_index(power_array.argmax(), power_array.shape)
    optimum_turbine = turbine_area_list[optimum_tup[0]]
    optimum_flow = rated_flow_list[optimum_tup[1]]


    return power_array, max_power, optimum_turbine, optimum_flow

#function to save a plot of power by turbine size and flow rate
def surface_map(rated_flow_list, turbine_area_list, power_array, save_name, title_string):
    rated_flow_list, turbine_area_list = np.meshgrid(rated_flow_list, turbine_area_list)
    plot_title = 'Total Energy Generation By Turbine Size and Rated Flow: ' + title_string
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(rated_flow_list, turbine_area_list, power_array, rstride=1, cstride=1, cmap=cm.YlOrRd_r)
    ax.set_xlabel(r'Rated Flow $(m^{3}s^{-1})$')
    ax.set_ylabel(r'Turbine Swept Area $(m^2)$')
    ax.set_zlabel(r'Total Energy Generated $(kWh)$')
    plt.suptitle(plot_title)
    plt.savefig(save_name)

#set parameters for simulation
sim_length = 44000
sim_time_division = 2
#time index for all calculations, calculate all tides
time_values = np.arange(0, sim_length, sim_time_division)
tide_list = [tide_height(time) for time in time_values]

#values for turbine swept area and turbine rating for optimisation
#purposes
rated_flow_list = [x for x in np.arange(1, 100, 1)]
turbine_area_list = [x for x in np.arange(0.2, 20, 0.2)]

#import bay volume values and read trapeziodal as bay model
bay_data_df = pd.read_csv('bay_volume_data')
bay_vol_ref = bay_data_df['trapezoidal_bay']

#simulate total power generated for turbine sizes, flow rate and control logic
ebb_power_gen, ebb_power_max, ebb_turbine, ebb_flow = optimisation(turbine_area_list, rated_flow_list, ebb_logic)
flood_power_gen, flood_power_max, flood_turbine, flood_flow = optimisation(turbine_area_list, rated_flow_list, flood_logic)
two_way_power_gen, two_way_power_max, two_way_turbine, two_way_flow = optimisation(turbine_area_list, rated_flow_list, two_way_logic)

#save surface map plots of previously calculated values
surface_map(rated_flow_list, turbine_area_list, ebb_power_gen, 'ebb_power_optimisation.png', 'Ebb Generation')
surface_map(rated_flow_list, turbine_area_list, flood_power_gen, 'flood_power_optimisation.png', 'Flood Generation')
surface_map(rated_flow_list, turbine_area_list, two_way_power_gen, 'two_way_power_optimisation', 'Two Way Generation')

#formatted string for writing to file
desc_string = 'power_total, optimal_turbine, optimal_flow\n'
formatted_string = '%.2f, %.2f, %.2f\n%.2f, %.2f, %.2f\n%.2f, %.2f, %.2f' % (ebb_power_max, ebb_turbine, ebb_flow, flood_power_max, flood_turbine, flood_flow, two_way_power_max, two_way_turbine, two_way_flow)

#write to and save file
with open('optimisation_file.dat', 'w') as data_file:
    data_file.write(desc_string)
    data_file.write(formatted_string)