I'm a novice programmer, but know my way around excel. However, I'm trying to teach myself python to enable myself to work with much larger datasets and, primarily, because I'm finding it really interesting and enjoyable.
I'm trying to figure out how to recreate the Excel goal seek function (I believe SciPy newton should be equivalent) within the script I have written. However, instead of being able to define a simple function f(x) for which to find the root of, I'm trying to find the root of the sum of a dataframe column. And I have no idea how to approach this.
My code up until the goal seek part is as follows:
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import weibull_min
# need to use a gamma function later on, so import math
import math
%matplotlib inline
# create dataframe using lidar experimental data
df = pd.read_csv(r'C:\Users\Latitude\Documents\Coursera\Wind Resource\Proj' \
'ect\Wind_Lidar_40and140.txt',
sep=' ',
header=None,
names=['Year','Month','Day','Hour','v_40','v_140'])
# add in columns for velocity cubed
df['v_40_cubed'] = df['v_40']**3
df['v_140_cubed'] = df['v_140']**3
# calculate mean wind speed, mean cubed wind speed, mean wind speed cubed
# use these to calculate energy patter factor, c and k
v_40_bar = df['v_40'].mean()
v_40_cubed_bar = df['v_40_cubed'].mean()
v_40_bar_cubed = v_40_bar ** 3
# energy pattern factor = epf
epf = v_40_cubed_bar / v_40_bar_cubed
# shape parameter = k
k_40 = 1 + 3.69/epf**2
# scale factor = c
# use imported math library to use gamma function math.gamma
c_40 = v_40_bar / math.gamma(1+1/k_40)
# create new dataframe from current, using bins of 0.25 and generate frequency for these
#' bins'
bins_1 = np.linspace(0,16,65,endpoint=True)
freq_df = df.apply(pd.Series.value_counts, bins=bins_1)
# tidy up the dataframe by dropping superfluous columns and adding in a % time column for
# frequency
freq_df_tidy = freq_df.drop(['Year','Month','Day','Hour','v_40_cubed','v_140_cubed'], axis=1)
freq_df_tidy['v_40_%time'] = freq_df_tidy['v_40']/freq_df_tidy['v_40'].sum()
# add in usable bin value for potential calculation of weibull
freq_df_tidy['windspeed_bin'] = np.linspace(0,16,64,endpoint=False)
# calculate weibull column and wind power density from the weibull fit
freq_df_tidy['Weibull_40'] = weibull_min.pdf(freq_df_tidy['windspeed_bin'], k_40, loc=0, scale=c_40)/4
freq_df_tidy['Wind_Power_Density_40'] = 0.5 * 1.225 * freq_df_tidy['Weibull_40'] * freq_df_tidy['windspeed_bin']**3
# calculate wind power density from experimental data
df['Wind_Power_Density_40'] = 0.5 * 1.225 * df['v_40']**3
At this stage, the result from the Weibull data, round(freq_df_tidy['Wind_Power_Density_40'].sum(),2) gives 98.12.
The result from the experimental data, round(df['Wind_Power_Density_40'].mean(),2) gives 101.14.
My aim now is to optimise the parameter c_40, which is used in the weibull calculated weibull power density (98.12), so that the result of the function round(freq_df_tidy['Wind_Power_Density_40'].sum(),2) is close to equal the experimental wind power density (101.14).
Any help on this would be hugely appreciated. Apologies if I've entered too much code into the request - I wanted to provide as much detail as possible. From my research, I think the SciPy newton method should do the trick, but I can't figure out how to apply it here.