Source code for opentidalfarm.optimisation_helpers

import os.path
import numpy
import dolfin
from dolfin import Constant, log, INFO
from .helpers import function_eval
from dolfin_adjoint import InequalityConstraint, EqualityConstraint

__all__ = ["MinimumDistanceConstraints", "MinimumDistanceConstraintsLargeArrays",
    "friction_constraints", "get_domain_constraints", "position_constraints",
    "get_distance_function", "ConvexPolygonSiteConstraint", "DomainRestrictionConstraints"]

[docs]def position_constraints(config): ''' This function returns the constraints to ensure that the turbine positions remain inside the domain. ''' n = len(config.params["position"]) if n == 0: raise ValueError("You need to deploy the turbines before computing the position constraints.") lb_x = config.domain.site_x_start + config.params["turbine_x"] / 2 lb_y = config.domain.site_y_start + config.params["turbine_y"] / 2 ub_x = config.domain.site_x_end - config.params["turbine_x"] / 2 ub_y = config.domain.site_y_end - config.params["turbine_y"] / 2 if not lb_x < ub_x or not lb_y < ub_y: raise ValueError("Lower bound is larger than upper bound. Is your domain large enough?") # The control variable is ordered as [t1_x, t1_y, t2_x, t2_y, t3_x, ...] lb = n * [Constant(lb_x), Constant(lb_y)] ub = n * [Constant(ub_x), Constant(ub_y)] return lb, ub
[docs]def friction_constraints(config, lb=0.0, ub=None): ''' This function returns the constraints to ensure that the turbine friction controls remain reasonable. ''' if ub is not None and not lb < ub: raise ValueError("Lower bound is larger than upper bound") if ub is None: ub = 10 ** 12 n = len(config.params["position"]) return n * [Constant(lb)], n * [Constant(ub)]
[docs]class DomainRestrictionConstraints(InequalityConstraint): def __init__(self, config, feasible_area, attraction_center): ''' Generates the inequality constraints to enforce the turbines in the feasible area. If the turbine is outside the domain, the constraints is equal to the distance between the turbine and the attraction center. ''' self.config = config self.feasible_area = feasible_area # Compute the gradient of the feasible area fs = dolfin.FunctionSpace(feasible_area.function_space().mesh(), "DG", feasible_area.function_space().ufl_element().degree() - 1) feasible_area_grad = (dolfin.Function(fs), dolfin.Function(fs)) t = dolfin.TestFunction(fs) log(INFO, "Solving for gradient of feasible area") for i in range(2): form = dolfin.inner(feasible_area_grad[i], t) * dolfin.dx - dolfin.inner(feasible_area.dx(i), t) * dolfin.dx if dolfin.NonlinearVariationalSolver.default_parameters().has_parameter("linear_solver"): dolfin.solve(form == 0, feasible_area_grad[i], solver_parameters={"linear_solver": "cg", "preconditioner": "amg"}) else: dolfin.solve(form == 0, feasible_area_grad[i], solver_parameters={"newton_solver": {"linear_solver": "cg", "preconditioner": "amg"}}) self.feasible_area_grad = feasible_area_grad self.attraction_center = attraction_center
[docs] def length(self): m_pos = self.config.params['turbine_pos'] return len(m_pos)
[docs] def function(self, m): ieqcons = [] if len(self.config.params['controls']) == 2: # If the controls consists of the the friction and the positions, then we need to first extract the position part assert(len(m) % 3 == 0) m_pos = m[len(m) / 3:] else: m_pos = m for i in range(len(m_pos) / 2): x = m_pos[2 * i] y = m_pos[2 * i + 1] try: ieqcons.append(function_eval(self.feasible_area, (x, y))) except RuntimeError: print("Warning: a turbine is outside the domain") ieqcons.append((x - self.attraction_center[0]) ** 2 + (y - self.attraction_center[1]) ** 2) # Point is outside domain arr = -numpy.array(ieqcons) if any(arr <= 0): log(INFO, "Domain restriction inequality constraints (should be >= 0): %s" % arr) return arr
[docs] def jacobian(self, m): ieqcons = [] if len(self.config.params['controls']) == 2: # If the controls consists of the the friction and the positions, then we need to first extract the position part assert(len(m) % 3 == 0) m_pos = m[len(m) / 3:] else: m_pos = m for i in range(len(m_pos) / 2): x = m_pos[2 * i] y = m_pos[2 * i + 1] primes = numpy.zeros(len(m)) try: primes[2 * i] = function_eval(self.feasible_area_grad[0], (x, y)) primes[2 * i + 1] = function_eval(self.feasible_area_grad[1], (x, y)) except RuntimeError: primes[2 * i] = 2 * (x - self.attraction_center[0]) primes[2 * i + 1] = 2 * (y - self.attraction_center[1]) ieqcons.append(primes) return -numpy.array(ieqcons)
[docs]def get_domain_constraints(config, feasible_area, attraction_center): return DomainRestrictionConstraints(config, feasible_area, attraction_center)
[docs]def get_distance_function(config, domains): V = dolfin.FunctionSpace(config.domain.mesh, "CG", 1) v = dolfin.TestFunction(V) d = dolfin.TrialFunction(V) sol = dolfin.Function(V) s = dolfin.interpolate(Constant(1.0), V) domains_func = dolfin.Function(dolfin.FunctionSpace(config.domain.mesh, "DG", 0)) domains_func.vector().set_local(domains.array().astype(numpy.float)) def boundary(x): eps_x = config.params["turbine_x"] eps_y = config.params["turbine_y"] min_val = 1 for e_x, e_y in [(-eps_x, 0), (eps_x, 0), (0, -eps_y), (0, eps_y)]: try: min_val = min(min_val, domains_func((x[0] + e_x, x[1] + e_y))) except RuntimeError: pass return min_val == 1.0 bc = dolfin.DirichletBC(V, 0.0, boundary) # Solve the diffusion problem with a constant source term log(INFO, "Solving diffusion problem to identify feasible area ...") a = dolfin.inner(dolfin.grad(d), dolfin.grad(v)) * dolfin.dx L = dolfin.inner(s, v) * dolfin.dx dolfin.solve(a == L, sol, bc) return sol
[docs]class MinimumDistanceConstraints(InequalityConstraint): """This class implements minimum distance constraints between turbines. .. note:: This class subclasses `dolfin_adjoint.InequalityConstraint`_. The following method names must not change: * ``length(self)`` * ``function(self, m)`` * ``jacobian(self, m)`` .. _dolfin_adjoint.InequalityConstraint: http://www.dolfin-adjoint.org/en/latest/documentation/api.html#dolfin_adjoint.InequalityConstraint """ def __init__(self, turbine_positions, minimum_distance, controls): """Create MinimumDistanceConstraints :param serialized_turbines: The serialized turbine paramaterisation. :type serialized_turbines: numpy.ndarray. :param minimum_distance: The minimum distance allowed between turbines. :type minimum_distance: float. :raises: NotImplementedError """ if len(turbine_positions)==0: raise NotImplementedError("Turbines must be deployed for distance " "constraints to be used.") self._turbines = numpy.asarray(turbine_positions).flatten().tolist() self._minimum_distance = minimum_distance self._controls = controls def _sl2norm(self, x): """Calculates the squared l2norm of a vector x.""" return sum([v**2 for v in x])
[docs] def length(self): """Returns the number of constraints ``len(function(m))``.""" n_constraints = 0 for i in range(len(self._turbines)): for j in range(len(self._turbines)): if i <= j: continue n_constraints += 1 return n_constraints
[docs] def function(self, m): """Return an object which must be zero for the point to be feasible. :param m: The serialized paramaterisation of the turbines. :tpye m: numpy.ndarray. :returns: numpy.ndarray -- each entry must be zero for a poinst to be feasible. """ dolfin.log(dolfin.PROGRESS, "Calculating minimum distance constraints.") inequality_constraints = [] for i in range(len(m)/2): for j in range(len(m)/2): if i <= j: continue inequality_constraints.append(self._sl2norm([m[2*i]-m[2*j], m[2*i+1]-m[2*j+1]]) - self._minimum_distance**2) inequality_constraints = numpy.array(inequality_constraints) if any(inequality_constraints <= 0): dolfin.log(dolfin.WARNING, "Minimum distance inequality constraints (should all " "be > 0): %s" % inequality_constraints) return inequality_constraints
[docs] def jacobian(self, m): """Returns the gradient of the constraint function. Return a list of vector-like objects representing the gradient of the constraint function with respect to the parameter m. :param m: The serialized paramaterisation of the turbines. :tpye m: numpy.ndarray. :returns: numpy.ndarray -- the gradient of the constraint function with respect to each input parameter m. """ dolfin.log(dolfin.PROGRESS, "Calculating the jacobian of minimum " "distance constraints function.") inequality_constraints = [] for i in range(len(m)/2): for j in range(len(m)/2): if i <= j: continue # Need to add space for zeros for the friction if self._controls.position and self._controls.friction: prime_inequality_constraints = numpy.zeros(len(m)*3/2) friction_length = len(m)/2 else: prime_inequality_constraints = numpy.zeros(len(m)) friction_length = 0 # Provide a shorter handle p_ineq_c = prime_inequality_constraints # The control vector contains the friction coefficients first, # so we need to shift here p_ineq_c[friction_length+2*i] = 2*(m[2*i] - m[2*j]) p_ineq_c[friction_length+2*j] = -2*(m[2*i] - m[2*j]) p_ineq_c[friction_length+2*i+1] = 2*(m[2*i+1] - m[2*j+1]) p_ineq_c[friction_length+2*j+1] = -2*(m[2*i+1] - m[2*j+1]) inequality_constraints.append(p_ineq_c) return numpy.array(inequality_constraints)
[docs]class MinimumDistanceConstraintsLargeArrays(InequalityConstraint): """This class implements minimum distance constraints between turbines which is well suited for large arrays. It enforces the minimum distance with only one enquality constraint of the form: .. math:: \sum_{i,j=0, i>j}^{N-1} \min(||p_i - p_j||^2 - D^2, 0) \ge 0 where N is the number of turbines, :math:`p_i` is the position of the i'th turbine and D is the minimum distance between two turbines. .. note:: This class subclasses `dolfin_adjoint.InequalityConstraint`_. The following method names must not change: * ``length(self)`` * ``function(self, m)`` * ``jacobian(self, m)`` .. _dolfin_adjoint.InequalityConstraint: http://www.dolfin-adjoint.org/en/latest/documentation/api.html#dolfin_adjoint.InequalityConstraint """ def __init__(self, turbine_positions, minimum_distance, controls): """Create MinimumDistanceConstraints :param serialized_turbines: The serialized turbine paramaterisation. :type serialized_turbines: numpy.ndarray. :param minimum_distance: The minimum distance allowed between turbines. :type minimum_distance: float. :raises: NotImplementedError """ if len(turbine_positions)==0: raise NotImplementedError("Turbines must be deployed for distance " "constraints to be used.") self._turbines = numpy.asarray(turbine_positions).flatten().tolist() self._minimum_distance = minimum_distance self._controls = controls def _sl2norm(self, x): """Calculates the squared l2norm of a vector x.""" return sum([v**2 for v in x])
[docs] def length(self): """Returns the number of constraints ``len(function(m))``.""" return 1
def _penalty(self, x_sq): return min(0, x_sq - self._minimum_distance**2) def _dpenalty(self, x_sq): """ Returns the derivative of the penalty function """ if x_sq - self._minimum_distance**2 > 0: return 0.0 else: return 1.0
[docs] def function(self, m): """Return an object which must be >0 for the point to be feasible. :param m: The serialized paramaterisation of the turbines. :tpye m: numpy.ndarray. :returns: numpy.ndarray -- each entry must be zero for a poinst to be feasible. """ dolfin.log(dolfin.PROGRESS, "Calculating minimum distance constraints.") if self._controls.position and self._controls.friction: friction_length = len(m)/3 m = m[friction_length:] value = 0 for i in range(len(m)/2): for j in range(len(m)/2): if i <= j: continue dist_sq = self._sl2norm([m[2*i]-m[2*j], m[2*i+1]-m[2*j+1]]) value += self._penalty(dist_sq) if value <= 0: dolfin.log(dolfin.WARNING, "Minimum distance inequality constraints (should all " "be => 0): %s" % value) return numpy.array([value])
[docs] def jacobian(self, m): """Returns the gradient of the constraint function. Return a list of vector-like objects representing the gradient of the constraint function with respect to the parameter m. :param m: The serialized paramaterisation of the turbines. :tpye m: numpy.ndarray. :returns: numpy.ndarray -- the gradient of the constraint function with respect to each input parameter m. """ dolfin.log(dolfin.PROGRESS, "Calculating the jacobian of minimum " "distance constraints function.") inequality_constraints = [] if self._controls.position and self._controls.friction: friction_length = len(m)/3 m = m[friction_length:] else: friction_length = 0 p_ineq_c = numpy.zeros(friction_length + len(m)) for i in range(len(m)/2): for j in range(len(m)/2): if i <= j: continue dist_sq = self._sl2norm([m[2*i]-m[2*j], m[2*i+1]-m[2*j+1]]) dvalue = self._dpenalty(dist_sq) # The control vector contains the friction coefficients first, # so we need to shift here p_ineq_c[friction_length+2*i] += dvalue * 2*(m[2*i] - m[2*j]) p_ineq_c[friction_length+2*j] += dvalue * (-2*(m[2*i] - m[2*j])) p_ineq_c[friction_length+2*i+1] += dvalue * 2*(m[2*i+1] - m[2*j+1]) p_ineq_c[friction_length+2*j+1] += dvalue * (-2*(m[2*i+1] - m[2*j+1])) return numpy.array([p_ineq_c])
[docs]class ConvexPolygonSiteConstraint(InequalityConstraint): ''' Generates the inequality constraints for generic polygon constraints. The parameter polygon must be a list of point coordinates that describes the site edges in anti-clockwise order. ''' def __init__(self, farm, vertices): self.farm = farm V = numpy.array(vertices) assert len(V.shape) == 2 assert V.shape[1] == 2 nvertices = V.shape[0] # algorithm taken from vert2con.m c = numpy.sum(V, axis=0)/float(nvertices) V = V - numpy.tile(c, (nvertices, 1)) A = numpy.nan * numpy.zeros((nvertices, 2)) for i in range(nvertices): j = (i + 1) % nvertices F = V[[i, j], :] ones = numpy.ones(2) A[i, :] = numpy.linalg.solve(F, ones) b = numpy.ones(nvertices) + numpy.dot(A, c) self.A = A self.b = b self.nvertices = nvertices # The region inside is defined by Ax <= b. # So, to make something that should be >= 0, # our function is b - A*x.
[docs] def length(self): return len(self.config.params['turbine_pos']) * self.nvertices
[docs] def output_workspace(self): return numpy.array([0]*self.length())
[docs] def function(self, m): ieqcons = [] controlled_by = self.farm.turbine_specification.controls if (controlled_by.position and (controlled_by.friction or controlled_by.dynamic_friction)): # If the controls consists of the the friction and the positions, then we need to first extract the position part assert(len(m) % 3 == 0) m_pos = m[len(m) / 3:] else: m_pos = m for i in range(len(m_pos) / 2): pos = numpy.array([m_pos[2 * i], m_pos[2 * i + 1]]) c = self.b - numpy.dot(self.A, pos) ieqcons = ieqcons + list(c) arr = numpy.array(ieqcons) if any(arr < 0): log(INFO, "Convex site position constraints (should be >= 0): %s" % arr) return arr
[docs] def jacobian(self, m): ieqcons = [] controlled_by = self.farm.turbine_specification.controls if (controlled_by.position and (controlled_by.friction or controlled_by.dynamic_friction)): # If the controls consists of the the friction and the positions, then we need to first extract the position part assert(len(m) % 3 == 0) m_pos = m[len(m) / 3:] mf_len = len(m_pos) / 2 else: m_pos = m mf_len = 0 for i in range(len(m_pos) / 2): pos = numpy.array([m_pos[2 * i], m_pos[2 * i + 1]]) for constraint in range(self.nvertices): d = -self.A[constraint, :] prime = numpy.zeros(len(m)) prime[mf_len + 2 * i] = d[0] prime[mf_len + 2 * i + 1] = d[1] ieqcons.append(prime) arr = numpy.array(ieqcons) return arr