I’m currently trying to develop a code that optimizes some parameters within a boundary to maximize another parameter. I need to insert 58 points inside a given boundary arranged based on 5 parameters: dx (distance between points in a row), dy (distance between rows), theta (angle of the rows), b (offset), and s (distance between turbines on the boundary).
To achieve this, I’ve developed an algorithm that works well for square boundaries. However, it seems to encounter issues with different geometries. After some arrangements are made, I encounter the error: NotImplementedError: Sub-geometries may have coordinate sequences, but multi-part geometries do not.
Is it possible that the problem arises from the intersection of the line in my code with the boundary occurring more than twice? Here attached is some output generated before the simulation gives the error.
What I did is generate the objective function, then calculate the function for each layout generated (generated with the bounds parameters)
def obj_fun(ga_instance, solution, solution_idx):
dx,dy,theta,b,s = solution
if np.abs(b) > dx: # penalty function to enforce the constraint on the value of b (logically must be smaller than dx)
print('z')
return 1e-10
# Extract layout coordinates from the parameters
x_coords, y_coords = fun.five_params2coords([dx,dy,theta,b,s], n_turbs, rect_boundaries)
if len(x_coords) != n_turbs: # penalty function to enforce the constraint on the number of turbines --> not necessary but makes the optimization faster. It can however become an issue when too few good solutions are found.
return 1e-10
# Reinitialize the FLORIS interface with the new layout
fi.reinitialize(
layout_x = x_coords,
layout_y = y_coords,
)
# Calculate AEP
aep = fi.get_farm_AEP( # the following parameters are defined outside the function
freq = freq,
cut_in_wind_speed = cut_in_ws,
cut_out_wind_speed = cut_out_ws,
)
plt.figure(figsize=(6, 4))
plt.scatter(x_coords, y_coords, color='blue')
plt.title('Optimized Turbine Layout')
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.grid(True)
plt.show()
return aep
bounds = [range(int(3.5*D),20*D+10,10), range(int(7.5*D),15*D+10,10), range(-90,91,1), range(-4*D,4*D+10,10), range(int(8*D),20*D+10,10)]
n_var = 5 #number of variables
n_gen = 50
n_parents = 2 # usually between 2 and 5
sol_per_pop = 50 # usually between 10 and 100
n_proc = 8 # number of processors
ga_instance = pygad.GA(
fitness_func = obj_fun, # objective function
num_generations = n_gen, # iterations
num_parents_mating = n_parents,
sol_per_pop = sol_per_pop, # number of solutions per population
num_genes = n_var, # number of variables to optimize
gene_space = bounds, # bounds for the parameters --> the notation with range is necessary --> if you use a list, only the two extreme values are considered
# notice that the step is 10 --> this discretizes the variables --> might be useful for the 64 dimensions problem
parent_selection_type = "sss",
keep_elitism = 1, # means keep only the best solution of the previous generation
crossover_type = "uniform", # alternative to 'single_point' --> seems more logical
mutation_type = "adaptive", # mutation adapts to the fitness of the solution --> if good solution, low probability of mutation.
mutation_percent_genes = [60,20], # the first is the probability of mutation of the bad solutions, the second of the good ones
#mutation_type = "random", # Alternative to adaptive. This however causes good solution to get worse
#mutation_num_genes = 1, # genes = variables --> after mating, 1 gene is changed (used only with random mutation type)
#mutation_by_replacement = True, # to avoid that variables go out of bounds after mutation
suppress_warnings=True,
parallel_processing=['process', n_proc], # to parallelize the code
)
if __name__ == "__main__": # fundamental for parallelization
start = time.time()
ga_instance.run()
solution, solution_fitness, solution_idx = ga_instance.best_solution()
print(f"nParameters of the best solution : {solution}")
print(f'AEP = {solution_fitness/1e9:.2f} GWh')
end = time.time()
elapsed_time_seconds = end - start
elapsed_minutes = int(elapsed_time_seconds // 60)
elapsed_seconds = int(elapsed_time_seconds % 60)
print(f'Time elapsed for parallel optimization with {n_proc} cores : {elapsed_minutes} min {elapsed_seconds} secn')
dx_opt, dy_opt, theta_opt, b_opt, s_opt = solution
# Save the results in a txt file
data = {
'dx opt [m]': dx_opt,
'dy opt [m]': dy_opt,
'theta opt [°]': theta_opt,
'b opt [m]': b_opt,
's opt [m]': s_opt,
'opt AEP [GWh]': solution_fitness/1e9,
'time [sec]': elapsed_time_seconds,
'algorithm': 'Genetic algorithm with pyGAD',
'settings': {
'num_generations': n_gen,
'num_parents_mating': n_parents,
'sol_per_pop': sol_per_pop,
'bounds': bounds,
'num CPU': n_proc,
}
}
# Store the optimized parameters in a txt file
file_path = r"C:FLORISENI_siteWork_folder[Thesis][PowerFloatingWindFarm]backup_codicippt_e_immaginioptimization_resultsprova_bottom_fixed_farm_5param_opt_layout.txt.txt"
with open(file_path, 'w') as file:
json.dump(data, file, indent=2)
ga_instance.plot_fitness()
# Obtain the optimized turbine coordinates
x_coords_opt, y_coords_opt = fun.five_params2coords([dx_opt, dy_opt, theta_opt, b_opt, s_opt], n_turbs, rect_boundaries)
# Plot the optimized turbine layout
plt.figure(figsize=(8, 6))
plt.scatter(x_coords_opt, y_coords_opt, color='blue')
plt.title('Optimized Turbine Layout')
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.grid(True)
plt.show()
And here below the function for realizing the grid:
# Function that finds the intersection between a line and a polygon --> used below
def line_polygon_intersection(angle, point, polygon, D=240):
# Angle in degrees, point as a tuple (x, y), polygon as a shapely Polygon object
# Output is a list of tuples (x, y) containing the intersection points or None if there is no intersection
# Calculate the direction of the line based on the angle
segment_length_x = 10*(5*D) # To be tuned depending on the size of the polygon (line in Python means a finite segment) so that the segments intesects the boundary
segment_length_y = segment_length_x * np.tan(np.radians(angle))
# Place the segment so that the midpoint is at the given point
x_mid, y_mid = point
x_start = x_mid - 0.5 * segment_length_x
y_start = y_mid - 0.5 * segment_length_y
x_end = x_mid + 0.5 * segment_length_x
y_end = y_mid + 0.5 * segment_length_y
# Define the line using the calculated starting and ending points
line = LineString([(x_start, y_start), (x_end, y_end)])
# Plot the polygon and the line
x_bound, y_bound = polygon.exterior.xy
plt.plot(x_bound, y_bound, color='b')
plt.plot(*line.xy, color='r')
plt.gca().set_aspect('equal', adjustable='box') # Equal aspect ratio
plt.show() # Display the plot
# Check if the line intersects the polygon
if line.intersects(polygon):
intersection = line.intersection(polygon)
if hasattr(intersection, 'coords'):
return [tuple(coord) for coord in intersection.coords]
elif hasattr(intersection, 'geoms'): # Check if it's a multipart geometry
coords = []
for geom in intersection.geoms:
if hasattr(geom, 'coords'):
coords.extend([tuple(coord) for coord in geom.coords])
else:
print('Unexpected intersection geometry:', geom)
return coords
else:
print('Unexpected intersection geometry:', intersection)
return None
else:
return None # No intersection
# Function that extracts the coordinates of the turbines from the 4 parameters --> used below
def params2coords(params, n_turbs_max, rect_boundaries):
# Function that takes as input the 4 parameters that define the regular layout and returns
# the corresponding coordinates of the turbines
# ATTENTION: works only for a rectangular boundary (it shouldn't be hard to generalize it though)
# STEPS OF THE INITIALIZATION
# 1) extract boundary from file and define the polygon
# 2) define a rotating frame of reference in the center of the polygon
# 3) check if the straight line from that point and with given slope, intersects the polygon, if not, move to the next iteration remebering this fact
# 4) if there is intersection, start from the 'start_point' and intialize turbines in that row with logic 0 1 -1 2 -2 3 -3...
# 5) Once the first row is initialized, the second one is initialized with the same logic but after moving by dy (rotated) towards up or down (logic 0 1 -1 2 -2 3 -3...)
# 6) This goes on until the number of turbines is reached or the rows don't intersect the polygon anymore (both from bottom and top)
# Parameters extraction
dx,dy,theta,b = params
theta_rad = np.deg2rad(theta)
# Extract rectangular boundary coordinates from the json file
x_rect, y_rect = zip(*rect_boundaries)
# Define values that will matter when initializing the single turbines
x_rect_max = max(x_rect)
x_rect_min = min(x_rect)
y_rect_max = max(y_rect)
y_rect_min = min(y_rect)
# Transform the boundaries in a polygon so that geometrical operations can be performed easily
polygon = Polygon(rect_boundaries)
# Find the center of the polygon
x0, y0 = polygon.centroid.xy[0][0], polygon.centroid.xy[1][0]
# Start_point is the point from which the initialization starts AT EACH NEW ROW
start_x = x0
start_y = y0
# Initialize the lists of coordinates
x_coords =[]
y_coords =[]
# I define 4 boolean variable that will be used to stop the initialization
# 'upper_row_no_inters' is turned True when the upper straight line does not intersect the polygon anymore
# 'bottom_row_no_inters' is turned True when the lower straight line does not intersect the polygon anymore
# 'left_point_out' is turned True when the points on a row exceed the left boundary of the polygon
# 'right_point_out' is turned True when the points on a row exceed the right boundary of the polygon
upper_row_no_inters = False
bottom_row_no_inters = False
n_turbs = 0
is_n_max_reached = False
m = 1 #Starting from 1 instead than zero allows to generate the desired sequence: 0, 1, -1, 2, -2, 3, -3.....
tol=1e-8
while upper_row_no_inters*bottom_row_no_inters==False and is_n_max_reached==False:
left_point_out = False
right_point_out = False
n = 1 #Starting from 1 instead than zero allows to generate the desired sequence: 0, 1, -1, 2, -2, 3, -3.....
# Compute the starting point for the next row and add the offset b # ATTENZIONE CHE PER IL PRIMO CASO TUTTO CIò CHE AGGIUNGO A X0 DEVE RISULTARE NULLO
start_x = x0 + ((-1)**m)*(m//2) * dy * np.cos(theta_rad+np.pi/2) + ((-1)**m)*(m//2) * b * np.cos(theta_rad)
start_y = y0 + ((-1)**m)*(m//2) * dy * np.sin(theta_rad+np.pi/2) + ((-1)**m)*(m//2) * b * np.sin(theta_rad)
start_point = (start_x,start_y)
# Check that the straight line with slope theta and passing through start_point is intersecting the polygon
# if not, don't even enter the while loop
intersection_points = line_polygon_intersection(theta, start_point, polygon) # function created by me
if intersection_points == None:
if m%2!=0:
upper_row_no_inters = True
else:
bottom_row_no_inters = True
"""
if theta >= 0:
if start_x < x_rect_min or start_y > y_rect_max:
upper_row_no_inters = True
elif start_x > x_rect_max or start_y < y_rect_min:
bottom_row_no_inters = True
else:
print('A') # This appears only if some cases were not considered
elif theta < 0:
if start_x > x_rect_max or start_y > y_rect_max:
upper_row_no_inters = True
elif start_x < x_rect_min or start_y < y_rect_min:
bottom_row_no_inters = True
else:
print('B')
else:
print('C')
"""
else:
# If the starting point is outside of the polygon, one direction of initializing is unuseful for sure.
max_x_inters_point = max(x[0] for x in intersection_points)
min_x_inters_point = min(x[0] for x in intersection_points)
while right_point_out*left_point_out == False and is_n_max_reached==False:
# Compute the new point
new_x = start_x + ((-1) ** n) * (n // 2) * dx * np.cos(theta_rad)
new_y = start_y + ((-1) ** n) * (n // 2) * dx * np.sin(theta_rad)
new_point = Point(new_x,new_y)
# Check if the new point is inside the polygon
#if polygon.contains(new_point): # is True also if the point is on the border but it has issues with floating points
if polygon.distance(new_point) < tol:
x_coords = x_coords + [new_x]
y_coords = y_coords + [new_y]
n_turbs = n_turbs + 1
# Check if the max number of turbines is reached
if n_turbs == n_turbs_max:
is_n_max_reached = True
# if the new point is outside the polygon: put either right_point_out or left_point_out to True
else:
if new_x > max_x_inters_point:
right_point_out = True
elif new_x < min_x_inters_point:
left_point_out = True
"""
if theta >= 0:
if new_x > x_rect_max or new_y > y_rect_max:
right_point_out = True
elif new_x < x_rect_min or new_y < y_rect_min:
left_point_out = True
elif new_x == x_rect_max or new_y == y_rect_max or new_x == x_rect_min or new_y == y_rect_min:
print('coordinata sul bordo')
else :
print('D') # in questo caso non capisco cosa stia succedendo
elif theta < 0:
if new_x > x_rect_max or new_y < y_rect_min:
right_point_out = True
elif new_x < x_rect_min or new_y > y_rect_max:
left_point_out = True
elif new_x == x_rect_max or new_y == y_rect_max or new_x == x_rect_min or new_y == y_rect_min:
print('coordinata sul bordo')
else:
print('E')
else:
print('F')
"""
n = n+1
m = m+1
#if n_turbs < n_turbs_max:
# print(f'Only {n_turbs} turbines fit in the area with these 4 parameters: [dx,dy,theta,b]={dx,dy,theta,b}' )
#elif n_turbs == n_turbs_max:
# print(f'Valid guess: 32 turbs fit with these 4 param: [dx,dy,theta,b]={dx,dy,theta,b}')
#else:
# print(f'Too many turbines ({n_turbs}) fit in the area with these 4 parameters: [dx,dy,theta,b]={dx,dy,theta,b}' )
return x_coords,y_coords,
# Function that extracts the Boundary Grid (regular layout function of 5 parameters) --> uses the function above
def five_params2coords(five_params, n_turbs, rect_boundaries):
# function that receives as input the 5 parameters, the number of turbines and the boundaries of the area
# and returns the coordinates of the turbines in a regular layout function of those 5 params
# It uses the function params2coords to built the pattern on a shrunk version of the provided area
# --> the distance between the outer and inner area is defined by the parameter 'min_dist_between_turb'
D = 240 # diameter of the IEA 15MW
min_dist_between_turb = 3.5*D
[dx,dy,theta,b,s] = five_params
# work in shapely since it is easier with geometries
polygon = Polygon(rect_boundaries)
Perimeter = LineString(rect_boundaries)
# compute the perimeter of the area
perimeter_length = polygon.length
# compute the number of turbines on the border
n_turbs_border = math.floor(perimeter_length/s)
# Make sure that turbines on the border are lower than the total number of turbines
if n_turbs_border >= n_turbs:
print(f'Too many turbines ({n_turbs_border}) on the border. Increase the value of s ({s}).')
x_coords = [] # decide to leave the array empty so that the AEP should be zero (?)
y_coords = []
return x_coords, y_coords
# Otherwise initialize the turbines on the border
coords_turbs_border = []
for i in range(n_turbs_border):
point = Perimeter.interpolate(i * s)
point_tuple = (point.x, point.y)
coords_turbs_border.append(point_tuple)
# compute the number of turbines inside the area
n_turbs_inside = n_turbs - n_turbs_border
# Create a buffer around the polygon of a fixed values from each border of the polygon --> chosen value = s
distance = -min_dist_between_turb
inner_polygon = polygon.buffer(distance, join_style = 'mitre')
# extract the coordinates of the inner area
inner_rect_boundaries = list(inner_polygon.exterior.coords)
# initialize the turbines inside the area with the params2coords function
x_coords_inside,y_coords_inside = params2coords([dx,dy,theta,b], n_turbs_inside, inner_rect_boundaries)
coords_turbs_inside = list(zip(x_coords_inside, y_coords_inside))
# merge the two lists of coordinates
coords_turbs = coords_turbs_border + coords_turbs_inside
x_coords, y_coords = zip(*coords_turbs)
if len(x_coords) != n_turbs:
print(f'Too low number of turbines ({len(x_coords)}) fit in the area with these 5 parameters: [dx,dy,theta,b,s]={dx,dy,theta,b,s}')
else:
print(f'Valid guess: {len(x_coords)} fit in the area with these 5 parameters: [dx,dy,theta,b,s]={dx,dy,theta,b,s}')
return x_coords, y_coords