AEP Calculate Problem

[donghyeoncheon]

The code below is part of the code I am currently working on to calculate AEP (Annual Energy Production). However, I am facing issues with setting wind direction and wind speed, specifically regarding FlorisInterface.reinitialize for wind directions. Can you suggest a solution?

Is there anything specific to check when saving yaw optimization values to a CSV file and importing them later? The reason is that when calculating AEP, Geometrics, Scipy, and SR all produce the same values.

D = 40.0
F = D * 3  

fi.reinitialize(
    layout_x=[882.51*F, 883.71*F, 884.82*F, 886.34*F, 888.48*F, 888.07*F, 886.76*F, 885.32*F, 887.92*F, 889.53*F, 890.96*F, 892.56*F, 890.73*F], 
    layout_y=[1758.3*F, 1755.87*F, 1753.26*F, 1751.23*F, 1752.14*F, 1755.57*F, 1758.35*F, 1760.55*F, 1761.95*F, 1759.4*F, 1757.2*F, 1754.7*F, 1753.35*F],
    wind_directions=np.arange(270.0, 280.0, 1.0),
    wind_speeds=[7.0],
)

...
#### AEP (Annual Energy Production) Calculation ####

# Path to FLORIS configuration file
floris_config_path = "examples\inputs\gch.yaml"

# Path to yaw angle CSV file
yaw_angles_csv_path = "optimized_yaw_angles_geo_v8.csv"

# Initialize the FLORIS interface
fi = wfct.floris_interface.FlorisInterface(floris_config_path)

# Load yaw angle data
df_yaw = pd.read_csv(yaw_angles_csv_path)

# Select yaw angles for a specific wind direction (270 degrees)
wind_direction = 270.0
yaw_angles = df_yaw[df_yaw['Wind Direction'] == wind_direction].iloc[0, 1:].values

# Set wind speed (e.g., 8.0 m/s)
wind_speed = 7.0

# Run the FLORIS model

# Set wind direction and wind speed
fi.reinitialize(wind_direction=[wind_direction], wind_speed=[wind_speed])
fi.calculate_wake(yaw_angles=yaw_angles)

# Calculate wake effects
fi.calculate_wake(yaw_angles=np.array([yaw_angles]))

# Calculate power output of each turbine
turbine_powers = fi.get_turbine_powers()

# Calculate Annual Energy Production (AEP)
# Assumption: Annual operating hours are 8760 hours
annual_hours = 8760
turbine_aep = turbine_powers * annual_hours / 1e9  # Converted to GWh
total_aep = np.sum(turbine_aep)

# Output the results
print("Turbines AEP (MWh):", turbine_aep)
print("Total AEP (MWh):", total_aep)

Below is the code that produces identical values. It may be a somewhat frustrating question, but I would appreciate your advice

# For Geometric Optimization
print("\nOptimized Yaw Angles - Geometric Method:")
for wd, yaws in zip(df_opt_geo.wind_direction, yaw_angles_opt_geo):
    print(f"Wind Direction {wd}: {yaws}")

# For Serial Refine Optimization
print("\nOptimized Yaw Angles - Serial Refine Method:")
for wd, yaws in zip(df_opt_sr.wind_direction, yaw_angles_opt_sr):
    print(f"Wind Direction {wd}: {yaws}")

# For SciPy Optimization
print("\nOptimized Yaw Angles - SciPy Method:")
for wd, yaws in zip(df_opt_scipy.wind_direction, yaw_angles_opt_scipy):
    print(f"Wind Direction {wd}: {yaws}")

# Create a DataFrame for each optimization method
df_yaw_angles_geo = pd.DataFrame(data=yaw_angles_opt_geo, columns=[f"Turbine_{t+1}" for t in range(13)])
df_yaw_angles_geo.insert(0, 'Wind Direction', df_opt_geo.wind_direction)

df_yaw_angles_sr = pd.DataFrame(data=yaw_angles_opt_sr, columns=[f"Turbine_{t+1}" for t in range(13)])
df_yaw_angles_sr.insert(0, 'Wind Direction', df_opt_sr.wind_direction)

df_yaw_angles_scipy = pd.DataFrame(data=yaw_angles_opt_scipy, columns=[f"Turbine_{t+1}" for t in range(13)])
df_yaw_angles_scipy.insert(0, 'Wind Direction', df_opt_scipy.wind_direction)

# Save to CSV
df_yaw_angles_geo.to_csv('optimized_yaw_angles_geo_v10.csv', index=False)
df_yaw_angles_sr.to_csv('optimized_yaw_angles_sr_v10.csv', index=False)
df_yaw_angles_scipy.to_csv('optimized_yaw_angles_scipy_v10.csv', index=False)

print("Yaw angles saved to CSV files.")

# Create a frequency array for specific condition
n_wind_directions = 10  # Adjust according to the number of wind directions
n_wind_speeds = 1       # Assuming one wind speed in the model

# Create a frequency array with equal probability for each wind direction
freq = np.full((len(wind_directions), 1), 1.0 / len(wind_directions))

# Compute AEP for each optimization method
# Geometric Optimization
fi.calculate_wake(yaw_angles=np.expand_dims(yaw_angles_opt_geo, axis=1))
aep_geo = fi.get_farm_AEP(freq=freq)

# Serial Refine Optimization
fi.calculate_wake(yaw_angles=np.expand_dims(yaw_angles_opt_sr, axis=1))
aep_sr = fi.get_farm_AEP(freq=freq)

# SciPy Optimization
fi.calculate_wake(yaw_angles=np.expand_dims(yaw_angles_opt_scipy, axis=1))
aep_scipy = fi.get_farm_AEP(freq=freq)

# Output AEP results
print("AEP - Geometric Method: {:.3f} GWh".format(aep_geo / 1.0e9))
print("AEP - Serial Refine Method: {:.3f} GWh".format(aep_sr / 1.0e9))
print("AEP - SciPy Method: {:.3f} GWh".format(aep_scipy / 1.0e9))

Also, how can I overcome the "dimension 0, 1, 2" error in the third case? Could I get some advice?

fi.calculate_wake(yaw_angles=yaw_angles)
axis 2 is out of bounds for array of dimension 0 axis 2 is out of bounds for array of dimension 1 axis 2 is out of bounds for array of dimension 2

I've tried using "reshape" and working with arrays, but the "dimension 0, 1, 2" error keeps occurring. I'm sorry for asking simple questions.

[rafmudaf]

@donghyeoncheon can you share the optimized yaw angle file or describe how you created it?

[donghyeoncheon]

optimized_yaw_angles_scipy_v10.csv
This is one of the Optimized file. I have referred to the Floris documentation to implement the AEP calculation, but I'm still finding it a bit challenging to apply AEP values for specific yaw angles from the CSV file. Thank youu

[rafmudaf]

@donghyeoncheon I see you closed this discussion. Were you able to resolve your problem? If so, could you please share the answer for the benefit of other users?

[donghyeoncheon]

@rafmudaf Yes, sir and It's necessary to verify exactly how AEP calculations are conducted. Could you please check it when you have time?

[donghyeoncheon]

N_TURBINES = 13 
# Calculate baseline AEP without yaw adjustments
fi.calculate_wake()
turbine_powers_baseline = fi.get_turbine_powers() / 1e6  # Convert from W to MW
average_power_output_per_hour_baseline = np.mean(turbine_powers_baseline, axis=(0, 2))
scaled_total_AEP_baseline = np.sum(average_power_output_per_hour_baseline) * 8760 / len(np.arange(150.0, 180.0, 5.0)) *0.1#270*0.1# * 0.1 #  # Scale down for simplification #######################

yaw_angles_geo = df_yaw_angles_geo.iloc[0, 1:].values
yaw_angles_sr = df_yaw_angles_sr.iloc[0, 1:].values
yaw_angles_scipy = df_yaw_angles_scipy.iloc[0, 1:].values
yaw_angles_bayesian = df_yaw_angles_bay.iloc[0, 1:].values
# Calculate baseline AEP without yaw adjustments
fi.calculate_wake()
turbine_powers_baseline = fi.get_turbine_powers() / 1e6  # Convert from W to MW
average_power_output_per_hour_baseline = np.mean(turbine_powers_baseline, axis=(0, 2))
total_hours_in_year = 8760  # Total hours in a year
scaled_total_AEP_baseline = np.sum(average_power_output_per_hour_baseline) * total_hours_in_year / len(wind_directions) * 0.1  # Scale down for simplification


# Function to calculate AEP with given yaw angles and return power difference
def calculate_aep_and_power_difference(fi, wind_directions, yaw_angles, total_hours_in_year, baseline_power_output_per_hour):
    N_TURBINES = fi.floris.farm.n_turbines
    optimized_yaw_angles = np.array([yaw_angles] * len(wind_directions)).reshape((len(wind_directions), 1, N_TURBINES))
    #optimized_yaw_angles = np.array([yaw_angles] * len(wind_directions)).reshape((len(wind_directions), 1, N_TURBINES))
    fi.calculate_wake(yaw_angles=optimized_yaw_angles)
    optimized_turbine_powers = fi.get_turbine_powers() / 1e6  # Convert from W to MW
    optimized_average_power_output_per_hour = np.mean(optimized_turbine_powers, axis=(0, 2))
    optimized_total_AEP = np.sum(optimized_average_power_output_per_hour) * total_hours_in_year / len(wind_directions)
    power_difference_percentage = 100 * (optimized_average_power_output_per_hour - baseline_power_output_per_hour) / baseline_power_output_per_hour
    return optimized_total_AEP * 0.1, np.mean(power_difference_percentage)  # Scale down AEP for simplification

# Calculate AEP for each optimization method and store results in a DataFrame
# 결과를 저장할 리스트 초기화
results_list = []

# 모든 최적화 방법에 대한 AEP 계산 및 결과 저장
for method, yaw_angles in {
    "Geometric Optimization": yaw_angles_geo, # Bayesian 추가
    "SciPy Optimization": yaw_angles_scipy, 
    "Serial Refine Optimization": yaw_angles_sr, 
    "Bayesian Optimization": yaw_angles_bayesian, 
}.items():
    optimized_total_AEP, power_difference_percentage = calculate_aep_and_power_difference(
        fi, wind_directions, yaw_angles, total_hours_in_year, average_power_output_per_hour_baseline)
    results_list.append({
        "Optimization Method": method,
        "Total AEP (GWh)": optimized_total_AEP,
        "Power Difference (%)": power_difference_percentage
    }) 
# 기준선 AEP
fi.calculate_wake()
baseline_turbine_powers = fi.get_turbine_powers() / 1e6  # Convert from W to MW
baseline_AEP = np.mean(baseline_turbine_powers) * 8760 / len(np.arange(150.0, 180.0, 5.0)) * 0.1 #  #*0.1 # /    

results_list = []

# 기준선 AEP 
results_list.append({
    "Optimization Method": "Baseline",
    "Total AEP (GWh)": baseline_AEP,
    "Power Difference (%)": 0.000001  # 기준선이므로 증감율은 0%
})

# Geometric Optimization
geometric_AEP = calculate_aep_and_power_difference(
    fi, wind_directions, yaw_angles_geo, 8760, baseline_turbine_powers
)[0]
results_list.append({
    "Optimization Method": "Geometric Optimization",
    "Total AEP (GWh)": geometric_AEP,
    "Power Difference (%)": (geometric_AEP - baseline_AEP) / baseline_AEP * 100
})

# SciPy Optimization
scipy_AEP = calculate_aep_and_power_difference(
    fi, wind_directions, yaw_angles_scipy, 8760, baseline_turbine_powers
)[0]
results_list.append({
    "Optimization Method": "SciPy Optimization",
    "Total AEP (GWh)": scipy_AEP,
    "Power Difference (%)": (scipy_AEP - baseline_AEP) / baseline_AEP * 100
})

# Serial Refine Optimization
serial_refine_AEP = calculate_aep_and_power_difference(
    fi, wind_directions, yaw_angles_sr, 8760, baseline_turbine_powers
)[0]
results_list.append({
    "Optimization Method": "Serial Refine Optimization",
    "Total AEP (GWh)": serial_refine_AEP,
    "Power Difference (%)": (serial_refine_AEP - baseline_AEP) / baseline_AEP * 100
})

# Bayesian Optimization
bayesian_AEP, _ = calculate_aep_and_power_difference(
    fi, wind_directions, yaw_angles_bayesian,  8760, baseline_turbine_powers #total_hours_in_year,
)

# Add the results to the DataFrame.
results_list.append({
    "Optimization Method": "Bayesian Optimization",
    "Total AEP (GWh)": bayesian_AEP,
    "Power Difference (%)": (bayesian_AEP - baseline_AEP) / baseline_AEP * 100
})


# 결과
results_df = pd.DataFrame(results_list)

#  재배치
desired_order = ["Baseline", "Geometric Optimization", "SciPy Optimization", "Serial Refine Optimization", "Bayesian Optimization"]
results_df['Optimization Method'] = pd.Categorical(
    results_df['Optimization Method'], 
    categories=desired_order, 
    ordered=True
)
results_df.sort_values('Optimization Method', inplace=True)


print("\n")
print("======================================================================================================================")
print(results_df)
print("\n")
print("======================================================================================================================")

fi.calculate_wake()
baseline_turbine_powers = fi.get_turbine_powers() / 1e6  # Convert from W to MW
baseline_average_power_output_per_hour = np.mean(baseline_turbine_powers, axis=(0, 2))
total_hours_in_year = 8760  # Total hours in a year

# Assuming wind directions are evenly distributed over the year
scaled_baseline_total_AEP = np.sum(baseline_average_power_output_per_hour) * total_hours_in_year / len(wind_directions) * 0.1  #0.1 #  #Scale down for simplification#######################

# Now 'scaled_baseline_total_AEP' is defined and can be used in the DataFrame and plotting


# Create a DataFrame to store the AEP results in the requested order
aep_results = {
    "Method": ["Baseline", "Geometric Optimization", "SciPy Optimization", "Serial Refine Optimization", "Bayesian Optimization"],
    "AEP (GWh)": [baseline_AEP, geometric_AEP, scipy_AEP, serial_refine_AEP, bayesian_AEP]
}

df_aep_results = pd.DataFrame(aep_results)

# Save to CSV file
csv_file_path = "AEP_Results_Optimizations_V3.csv"  # Adjust path as needed
df_aep_results.to_csv(csv_file_path, index=False)

print(f"AEP results saved to {csv_file_path}")

# Now for plotting the results, we include the Bayesian Optimization as well
aep_values = [
    baseline_AEP, geometric_AEP, scipy_AEP, serial_refine_AEP, bayesian_AEP
]
methods_with_baseline = ['Baseline', 'Geometric', 'SciPy', 'Serial Refine', 'Bayesian']

# Create figure and bar plot for AEP values
plt.figure(figsize=(12, 6))
bars = plt.bar(methods_with_baseline, aep_values, color=['cyan', 'blue', 'green', 'red', 'orange'])
plt.ylabel('AEP (GWh)')
plt.title('AEP Values by Optimization Method')

# Set y-axis limits to zoom in on the data
min_aep = min(aep_values) * 0.9  # 10% less than the smallest value for lower limit
max_aep = max(aep_values) * 1.1  # 10% more than the largest value for upper limit
plt.ylim(min_aep, max_aep)

# Set y-axis ticks for more detail
ytick_values = np.linspace(min_aep, max_aep, num=10)  # Creates 10 evenly spaced ticks
plt.yticks(ytick_values)

# Annotate bars with AEP values
for bar, value in zip(bars, aep_values):
    plt.text(bar.get_x() + bar.get_width() / 2, bar.get_height(), f'{value:.2f}', ha='center', va='bottom')

# Here, it is assumed that the bayesian_AEP variable contains the AEP result of Bayesian Optimization.

# Create a DataFrame including AEP results.
df_aep_results = pd.DataFrame({
    'Method': ['Baseline', 'Geometric Optimization', 'Serial Refine Optimization', 'SciPy Optimization', 'Bayesian Optimization'],
    'AEP (GWh)': [baseline_AEP, geometric_AEP, scipy_AEP, serial_refine_AEP, bayesian_AEP]
})

# Check the results.
print(df_aep_results)

# Create a bar plot using this DataFrame.
plt.figure(figsize=(12, 6))
bars = plt.bar(df_aep_results['Method'], df_aep_results['AEP (GWh)'], color=['cyan', 'blue', 'green', 'red', 'orange'])
plt.ylabel('AEP (GWh)')
plt.title('AEP Values by Optimization Method')

# Set y-axis limits.
min_aep = df_aep_results['AEP (GWh)'].min() * 0.9
max_aep = df_aep_results['AEP (GWh)'].max() * 1.1
plt.ylim(min_aep, max_aep)

# Set y-axis ticks.
ytick_values = np.linspace(min_aep, max_aep, num=10)
plt.yticks(ytick_values)

# Display AEP values on the bars.
for bar, value in zip(bars, df_aep_results['AEP (GWh)']):
    plt.text(bar.get_x() + bar.get_width() / 2, bar.get_height(), f'{value:.2f}', ha='center', va='bottom')

# Get the baseline AEP value.
baseline_AEP = df_aep_results[df_aep_results['Method'] == 'Baseline']['AEP (GWh)'].values[0]

# Calculate the increase/decrease rate compared to other optimization methods.
df_aep_results['AEP Increase (%)'] = (df_aep_results['AEP (GWh)'] - baseline_AEP) / baseline_AEP * 100

# Set the baseline (Baseline) to 0% increase/decrease.
df_aep_results.loc[df_aep_results['Method'] == 'Baseline', 'AEP Increase (%)'] = 0.0001

# Display the increase/decrease rate.
plt.figure(figsize=(12, 6))
bars = plt.bar(df_aep_results['Method'], df_aep_results['AEP Increase (%)'], color=['grey', 'blue', 'green', 'red', 'orange'])
plt.ylabel('AEP Increase (%)')
plt.title('AEP Increase Rate by Optimization Method Compared to Baseline')

# Display the increase/decrease rate on the bars.
for bar, value in zip(bars, df_aep_results['AEP Increase (%)']):
    yval = bar.get_height()
    plt.text(bar.get_x() + bar.get_width() / 2, yval, f'{value:.2f}%', ha='center', va='bottom' if value < 0 else 'top')

plt.tight_layout()
plt.show()

fi = FlorisInterface("examples\inputs\gch.yaml")

#######################################################################################################
D = 40.0
F = D * 3   # Factor based on turbine diameter
layout_x = [882.51 * F, 883.71 * F, 884.82 * F, 886.34 * F, 888.48 * F, 888.07 * F, 886.76 * F, 885.32 * F, 887.92 * F, 889.53 * F, 890.96 * F, 892.56 * F, 890.73 * F]
layout_y = [1758.3 * F, 1755.87 * F, 1753.26 * F, 1751.23 * F, 1752.14 * F, 1755.57 * F, 1758.35 * F, 1760.55 * F, 1761.95 * F, 1759.4 * F, 1757.2 * F, 1754.7 * F, 1753.35 * F]

wind_directions =  np.arange(150.0, 180.0, 5.0)  #[270]# [280.0]#     #####################################################################
wind_speeds = [8.0]   #np.arange(7.0, 8.0, 0.1) #

# Pass the new data to FlorisInterface
fi.reinitialize( 
    layout_x=layout_x,
    layout_y=layout_y,
    wind_directions=wind_directions,
    wind_speeds=wind_speeds
)

Please refer to the code above for the layout. I used the method of exporting the optimize values of FLORIS and other optimization algorithms to a CSV file, and then importing them to recalculate the AEP.

[donghyeoncheon]

Lastly, it is displayed in a prompt window and plotted for comparison of AEP production and the rate of increase or decrease. and In the layout, D and F are defined for adjusting the Diameter because the coordinate values do not differ significantly during the scaling process of the actual power plant layout.