Merge pull request #5 from byuflowlab/master

iea cs 3 and 4 analyses are running and appear correct, but cs4 shoul…

test/example_opt_3_ieacs3.jl

@@ -100,7 +100,7 @@ function wind_farm_opt(x)
 end
 
 # import model set with wind farm and related details
-include("./model_sets/model_set_7.jl")
+include("./model_sets/model_set_7_ieacs3.jl")
 
 # scale objective to be between 0 and 1
 obj_scale = 1E-11

test/example_opt_4_ieacs4.jl

@@ -0,0 +1,212 @@
+using Snopt
+using DelimitedFiles 
+using PyPlot
+import ForwardDiff
+
+# set up boundary constraint wrapper function
+function boundary_wrapper(x, params)
+    # include relevant globals
+    params.boundary_center
+    params.boundary_radius
+
+    # get number of turbines
+    nturbines = Int(length(x)/2)
+
+    # extract x and y locations of turbines from design variables vector
+    turbine_x = x[1:nturbines]
+    turbine_y = x[nturbines+1:end]
+
+    # get and return boundary distances
+    return ff.circle_boundary(boundary_center, boundary_radius, turbine_x, turbine_y)
+end
+
+# set up spacing constraint wrapper function
+function spacing_wrapper(x, params)
+    # include relevant globals
+    params.rotor_diameter
+
+    # get number of turbines
+    nturbines = Int(length(x)/2)
+
+    # extract x and y locations of turbines from design variables vector
+    turbine_x = x[1:nturbines]
+    turbine_y = x[nturbines+1:end]
+
+    # get and return spacing distances
+    return 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x,turbine_y)
+end
+
+# set up objective wrapper function
+function aep_wrapper(x, params)
+    # include relevant globals
+    params.turbine_z
+    params.rotor_diameter
+    params.hub_height
+    params.turbine_yaw
+    params.ct_model
+    params.generator_efficiency
+    params.cut_in_speed
+    params.cut_out_speed
+    params.rated_speed
+    params.rated_power
+    params.windresource
+    params.power_model
+    params.model_set
+    params.rotor_points_y
+    params.rotor_points_z
+    params.obj_scale
+
+    # get number of turbines
+    nturbines = Int(length(x)/2)
+
+    # extract x and y locations of turbines from design variables vector
+    turbine_x = x[1:nturbines] 
+    turbine_y = x[nturbines+1:end]
+
+    # calculate AEP
+    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
+                hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
+                cut_out_speed, rated_speed, rated_power, windresource, power_model, model_set,
+                rotor_sample_points_y=rotor_points_y,rotor_sample_points_z=rotor_points_z, hours_per_year=365.0*24.0)
+
+    # return the objective as an array
+    return [AEP]
+end
+
+# set up optimization problem wrapper function
+function wind_farm_opt(x)
+
+    # calculate spacing constraint value and jacobian
+    spacing_con = spacing_wrapper(x)
+    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)
+
+    # calculate boundary constraint and jacobian
+    boundary_con = boundary_wrapper(x)
+    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)
+
+    # combine constaint values and jacobians into overall constaint value and jacobian arrays
+    c = [spacing_con; boundary_con]
+    dcdx = [ds_dx; db_dx]
+
+    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
+    AEP = -aep_wrapper(x)[1]
+    dAEP_dx = -ForwardDiff.jacobian(aep_wrapper,x)
+
+    # set fail flag to false
+    fail = false
+
+    # return objective, constraint, and jacobian values
+    return AEP, c, dAEP_dx, dcdx, fail
+end
+
+# import model set with wind farm and related details
+include("./model_sets/model_set_7_ieacs4.jl")
+
+# scale objective to be between 0 and 1
+obj_scale = 1E-11
+
+# set wind farm boundary parameters
+boundary_center = [0.0,0.0]
+boundary_radius = 3000.0
+
+# set globals for use in wrapper functions
+struct params_struct{}
+    model_set
+    rotor_points_y
+    rotor_points_z
+    turbine_z
+    ambient_ti
+    rotor_diameter
+    boundary_center
+    boundary_radius
+    obj_scale
+    hub_height
+    turbine_yaw
+    ct_model
+    generator_efficiency
+    cut_in_speed
+    cut_out_speed
+    rated_speed
+    rated_power
+    windresource
+    power_model
+end
+
+params = params_struct(model_set, rotor_points_y, rotor_points_z, turbine_z, ambient_ti, 
+    rotor_diameter, boundary_center, boundary_radius, obj_scale, hub_height, turbine_yaw, 
+    ct_model, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed, rated_power, 
+    windresource, power_model)
+
+# initialize design variable array
+x = [copy(turbine_x);copy(turbine_y)]
+
+state_aeps = ff.calculate_state_aeps(turbine_x, turbine_y, turbine_z, rotor_diameter,
+                hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
+                cut_out_speed, rated_speed, rated_power, windresource, power_model, model_set;
+                rotor_sample_points_y=[0.0], rotor_sample_points_z=[0.0], hours_per_year=365.0*24.0)
+
+dir_aep = zeros(20)
+for i in 1:20
+    for j in 1:20
+        dir_aep[i] += state_aeps[(i-1)*20 + j]
+    end
+end
+
+# report initial objective value
+println("nturbines: ", nturbines)
+println("rotor diameter: ", rotor_diameter[1])
+println("starting AEP value (MWh): ", aep_wrapper(x, params)[1]*1E5)
+# println("frequencies ", windprobabilities)
+println("Directional AEP at start: ", dir_aep.*1E-6)
+continue
+# add initial turbine location to plot
+for i = 1:length(turbine_x)
+    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
+end
+
+# set general lower and upper bounds for design variables
+lb = zeros(length(x)) .- boundary_radius
+ub = zeros(length(x)) .+ boundary_radius
+
+# set up options for SNOPT
+options = Dict{String, Any}()
+options["Derivative option"] = 1
+options["Verify level"] = 3
+options["Major optimality tolerance"] = 1e-5
+options["Major iteration limit"] = 1e6
+options["Summary file"] = "summary.out"
+options["Print file"] = "print.out"
+
+# generate wrapper function surrogates
+spacing_wrapper(x) = spacing_wrapper(x, params)
+aep_wrapper(x) = aep_wrapper(x, params)
+boundary_wrapper(x) = boundary_wrapper(x, params)
+
+# run and time optimization
+t1 = time()
+xopt, fopt, info = snopt(obj_func, x, lb, ub, options)
+t2 = time()
+clkt = t2-t2
+
+# print optimization results
+println("Finished in : ", clkt, " (s)")
+println("info: ", info)
+println("end objective value: ", aep_wrapper(xopt))
+
+# extract final turbine locations
+turbine_x = copy(xopt[1:nturbines])
+turbine_y = copy(xopt[nturbines+1:end])
+
+# add final turbine locations to plot
+for i = 1:length(turbine_x)
+    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="--")) 
+end
+
+# add wind farm boundary to plot
+plt.gcf().gca().add_artist(plt.Circle((boundary_center[1],boundary_center[2]), boundary_radius, fill=false,color="C2"))
+
+# set up and show plot
+axis("square")
+xlim(-boundary_radius-200,boundary_radius+200)
+ylim(-boundary_radius-200,boundary_radius+200)
+plt.show()

test/model_sets/model_set_7_ieacs4.jl

@@ -0,0 +1,76 @@
+import FlowFarm; const ff = FlowFarm
+
+# based on IEA case study 3
+
+# set initial turbine x and y locations
+layout_file_name = "./inputfiles/iea37-ex-opt4.yaml"
+turbine_x, turbine_y, fname_turb, fname_wr = ff.get_turb_loc_YAML(layout_file_name)
+
+# calculate the number of turbines
+nturbines = length(turbine_x)
+
+# set turbine base heights
+turbine_z = zeros(nturbines)
+
+# set turbine yaw values
+turbine_yaw = zeros(nturbines)
+
+# set turbine design parameters
+turbine_file_name = string("./inputfiles/",fname_turb)
+turb_ci, turb_co, rated_ws, rated_pwr, turb_diam, turb_hub_height = ff.get_turb_atrbt_YAML(turbine_file_name)
+
+rotor_diameter = zeros(nturbines) .+ turb_diam # m
+hub_height = zeros(nturbines) .+ turb_hub_height   # m
+cut_in_speed = zeros(nturbines) .+ turb_ci  # m/s
+cut_out_speed = zeros(nturbines) .+ turb_co  # m/s
+rated_speed = zeros(nturbines) .+ rated_ws # m/s
+rated_power = zeros(nturbines) .+ rated_pwr # W
+generator_efficiency = zeros(nturbines) .+ 1.0
+
+# rotor swept area sample points (normalized by rotor radius)
+rotor_points_y = [0.0]
+rotor_points_z = [0.0]
+
+# set flow parameters
+windrose_file_name = string("./inputfiles/",fname_wr)
+winddirections, windspeeds, windprobabilities, ambient_ti = ff.get_wind_rose_YAML(windrose_file_name)
+nstates = length(winddirections)
+winddirections *= pi/180.0
+
+air_density = 1.1716  # kg/m^3
+shearexponent = 0.15
+ambient_tis = zeros(nstates) .+ ambient_ti
+measurementheight = zeros(nstates) .+ turb_hub_height
+
+# initialize power model
+power_model = ff.PowerModelPowerCurveCubic()
+
+# load thrust curve
+ct = 4.0*(1.0/3.0)*(1.0 - 1.0/3.0)
+
+# initialize thurst model
+ct_model1 = ff.ThrustModelConstantCt(ct)
+ct_model = Vector{typeof(ct_model1)}(undef, nturbines)
+for i = 1:nturbines
+    ct_model[i] = ct_model1
+end
+
+# initialize wind shear model
+wind_shear_model = ff.PowerLawWindShear(shearexponent)
+
+# get sorted indecies 
+sorted_turbine_index = sortperm(turbine_x)
+
+# initialize the wind resource definition
+windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheight, air_density, ambient_tis, wind_shear_model)
+
+# set up wake and related models
+k = 0.0324555
+wakedeficitmodel = ff.GaussSimple(k)
+
+wakedeflectionmodel = ff.JiminezYawDeflection()
+wakecombinationmodel = ff.SumOfSquaresFreestreamSuperposition()
+localtimodel = ff.LocalTIModelNoLocalTI()
+
+# initialize model set
+model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)