8. Continuous farm optimization¶
8.1. Introduction¶
This example runs a tidal farm optimization based on a continuous turbine farm representation. As outputs, one obtains an optimal turbine density function, from which one can derive the optimal number of turbines to be deployed, and their approximate layout within the farm.
It shows how to
set up a time-dependent shallow water model with a continuous farm representation;
define a callback function that is called after every optimisation iteration;
run the optimization
extract the optimal turbine density function and compute the optimal number of turbines.
The equations to be solved are the shallow water equations
where
\(u\) is the velocity,
\(\eta\) is the free-sufface displacement,
\(H=\eta + h\) is the total water depth where \(h\) is the water depth at rest,
\(c_b\) is the (quadratic) natural bottom friction coefficient,
\(\nu\) is the viscosity coefficient,
\(g\) is the gravitational constant.
8.2. Implementation¶
We begin with importing the OpenTidalFarm module.
import argparse
from opentidalfarm import *
from model_turbine import ModelTurbine
from vorticity_solver import VorticitySolver
import time
model_turbine = ModelTurbine()
print model_turbine
# Read the command line arguments
parser = argparse.ArgumentParser()
parser.add_argument('--turbines', required=True, type=int, help='number of turbines')
parser.add_argument('--optimize', action='store_true', help='Optimise instead of just simulate')
parser.add_argument('--withcuts', action='store_true', help='with cut in/out speeds')
parser.add_argument('--cost', type=float, default=0., help='the cost coefficient')
args = parser.parse_args()
Next we get the default parameters of a shallow water problem and configure it to our needs.
prob_params = SWProblem.default_parameters()
First we define the computational domain. We load a previously generated mesh
(using `Gmsh`
and converted with `dolfin-convert`
):
domain = FileDomain("mesh/headland.xml")
Once the domain is created we attach it to the problem parameters:
prob_params.domain = domain
turbine = SmearedTurbine()
V = FunctionSpace(domain.mesh, "DG", 0)
farm = Farm(domain, turbine, function_space=V)
# Sub domain for inflow (right)
class FarmDomain(SubDomain):
def inside(self, x, on_boundary):
return (9200 <= x[0] <= 10800 and
2500 <= x[1] <= 3500)
farm_domain = FarmDomain()
domains = MeshFunction("size_t", domain.mesh, domain.mesh.topology().dim())
domains.set_all(0)
farm_domain.mark(domains, 1)
site_dx = Measure("dx")(subdomain_data=domains)
farm.site_dx = site_dx(1)
#plot(domains, interactive=True)
prob_params.tidal_farm = farm
Next we specify the boundary conditions. The parameter t in the
dolfin.Expression
is special in OpenTidalFarm as it will be
automatically updated to the current timelevel during the solve.
tidal_amplitude = 5.
tidal_period = 12.42*60*60
H = 40
bcs = BoundaryConditionSet()
eta_channel = "amp*sin(omega*t + omega/pow(g*H, 0.5)*x[0])"
eta_expr = Expression(eta_channel, t=Constant(0), amp=tidal_amplitude,
omega=2*pi/tidal_period, g=9.81, H=H, degree=3)
bcs.add_bc("eta", eta_expr, facet_id=1, bctype="strong_dirichlet")
bcs.add_bc("eta", eta_expr, facet_id=2, bctype="strong_dirichlet")
# Apply a strong no-slip boundary condition. This can be changed to
# free slip (weakly enforced), by leaving out the Constant((0, 0))
# argument and changing bctype to "free_slip"
bcs.add_bc("u", Constant((0, 0)), facet_id=3, bctype="strong_dirichlet")
Again we attach boundary conditions to the problem parameters:
prob_params.bcs = bcs
The other parameters are straight forward:
# Equation settings
nu = Constant(60)
prob_params.viscosity = nu
prob_params.depth = Constant(H)
prob_params.friction = Constant(0.0025)
# Temporal settings
prob_params.theta = Constant(0.6)
prob_params.start_time = Constant(0)
prob_params.finish_time = Constant(2*tidal_period)
prob_params.dt = Constant(tidal_period/100)
prob_params.functional_final_time_only = False
# The initial condition consists of three components: u_x, u_y and eta
# Note that we do not set all components to zero, as some components of the
# Jacobian of the quadratic friction term is non-differentiable.
prob_params.initial_condition = Expression(("1e-7", "0", eta_channel), t=Constant(0),
amp=tidal_amplitude, omega=2*pi/tidal_period, g=9.81, H=H, degree=3)
Here we only set the necessary options. A full option list with its current values can be viewed with:
print prob_params
Once the parameter have been set, we create the shallow water problem:
problem = SWProblem(prob_params)
Next we create a shallow water solver. Here we choose to solve the shallow water equations in its fully coupled form. Again, we first ask for the default parameters, adjust them to our needs and then create the solver object.
sol_params = CoupledSWSolver.default_parameters()
sol_params.dump_period = 1
sol_params.output_dir = "output_{}_turbines_optimize_{}_cutinout_{}_cost_{}".format(args.turbines,
args.optimize, args.withcuts, args.cost)
sol_params.cache_forward_state = False
solver = CoupledSWSolver(problem, sol_params)
V = solver.function_space.extract_sub_space([0]).collapse()
Q = solver.function_space.extract_sub_space([1]).collapse()
base_u = Function(V, name="base_u")
base_u_tmp = Function(V, name="base_u_tmp")
# Define the functional
if args.withcuts:
power_functional = PowerFunctional(problem, cut_in_speed=1.0, cut_out_speed=3.)
else:
power_functional = PowerFunctional(problem)
cost_functional = args.cost * CostFunctional(problem)
functional = power_functional - cost_functional
# Define the control
control = TurbineFarmControl(farm)
# Set up the reduced functional
rf_params = ReducedFunctional.default_parameters()
rf_params.automatic_scaling = None
if args.optimize:
rf_params.save_checkpoints = True
rf_params.load_checkpoints = True
rf = ReducedFunctional(functional, control, solver, rf_params)
As always, we can print all options of the ReducedFunctional
with:
print rf_params
Now we can define the constraints for the controls and start the optimisation.
init_tf = model_turbine.maximum_smeared_friction/1000*args.turbines
farm.friction_function.assign(Constant(init_tf))
# Comment this for only forward modelling
if args.optimize:
maximize(rf, bounds=[0, model_turbine.maximum_smeared_friction],
method="L-BFGS-B", options={'maxiter': 15})
# Recompute the energy for the optimal farm array and store the results
sol_h5 = HDF5File(mpi_comm_world(), "{}/solution.h5".format(sol_params.output_dir), "w")
vort_solver = VorticitySolver(V)
def callback(s):
print "*** Storing timestep to solution.h5 ***"
u = project(s["u"], V)
eta = project(s["eta"], Q)
vort = vort_solver.solve(u)
sol_h5.write(u, "u_{}".format(float(s["time"])))
sol_h5.write(eta, "eta_{}".format(float(s["time"])))
sol_h5.write(vort, "vorticity_{}".format(float(s["time"])))
sol_h5.write(farm.friction_function, "turbine_friction_{}".format(float(s["time"])))
sol_h5.flush()
total_friction = assemble(farm.friction_function*farm.site_dx(1))
num_turbines = total_friction/model_turbine.friction
print "Estimated number of turbines: ", float(num_turbines)
# Recompute the cost, but this time with the power functional only
j = rf(farm.control_array)
# Power functional only
solver.parameters.callback = callback
rf_params.save_checkpoints = False
rf_params.load_checkpoints = False
rf = ReducedFunctional(power_functional, control, solver, rf_params)
energy = rf(farm.control_array)
# We are done with the solution.h5 file. Close it.
sol_h5.close()
# Save optimal friction as xdmf
optimal_turbine_friction_file = XDMFFile(mpi_comm_world(),
sol_params.output_dir+"/optimal_turbine_friction.xdmf")
optimal_turbine_friction_file.write(farm.friction_function)
# Compute the total turbine friction
total_friction = assemble(farm.friction_function*farm.site_dx(1))
# Compute the total cost
cost = float((prob_params.finish_time-prob_params.start_time) * args.cost * total_friction)
# Compute the site area
site_area = assemble(Constant(1)*farm.site_dx(1, domain=domain.mesh))
avg_power = energy/1e6/float(prob_params.finish_time-prob_params.start_time)
num_turbines = total_friction/model_turbine.friction
print "="*40
print "Site area (m^2): ", site_area
print "Cost coefficient: {}".format(args.cost)
print "Total energy (MWh): %e." % (energy/1e6/60/60)
print "Average power (MW): %e." % avg_power
print "Total cost: %e." % cost
print "Maximum smeared turbine friction: %e." % model_turbine.maximum_smeared_friction
print "Total turbine friction: %e." % total_friction
print "Average smeared turbine friction: %e." % (total_friction / site_area)
print "Average power / total friction: %e." % (avg_power / total_friction)
print "Friction per discrete turbine: {}".format(model_turbine.friction)
print "Estimated number of discrete turbines: {}".format(num_turbines)
8.3. How to run the example¶
The example code can be found in examples/headland-optimization/
in the
OpenTidalFarm
source tree, and executed as follows:
$ python headland-optimization.py --turbines N [--optimize]
where N is the number of turbines to deploy and –optimize is an optional flag to switch from simulation (default) to optimization.
The results can be viewed with [Paraview](http://www.paraview.org/).