NREL / floris
Showing 54 of 98 files from the diff.
Other files ignored by Codecov
docs/index.md has changed.
README.md has changed.
tests/conftest.py has changed.

@@ -14,33 +14,20 @@
Loading
14 14
15 15
from __future__ import annotations
16 16
17 -
import copy
18 -
from typing import Any, Tuple
19 17
from pathlib import Path
20 -
from itertools import repeat, product
21 -
from multiprocessing import cpu_count
22 -
from multiprocessing.pool import Pool
23 18
24 19
import numpy as np
25 20
import pandas as pd
26 -
import numpy.typing as npt
27 -
import matplotlib.pyplot as plt
28 -
from scipy.stats import norm
29 21
from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator
30 -
from numpy.lib.arraysetops import unique
31 22
32 -
from floris.utilities import Vec3
33 23
from floris.type_dec import NDArrayFloat
34 -
from floris.simulation import Farm, Floris, FlowField, WakeModelManager, farm, floris, flow_field
24 +
from floris.simulation import Floris
35 25
from floris.logging_manager import LoggerBase
36 -
from floris.tools.cut_plane import get_plane_from_flow_data
37 -
# from floris.tools.flow_data import FlowData
38 -
from floris.simulation.turbine import Ct, power, axial_induction, average_velocity
39 -
from floris.tools.interface_utilities import get_params, set_params, show_params
40 -
from floris.tools.cut_plane import CutPlane, change_resolution, get_plane_from_flow_data
41 26
42 -
# from .visualization import visualize_cut_plane
43 -
# from .layout_functions import visualize_layout, build_turbine_loc
27 +
from floris.simulation import State
28 +
29 +
from floris.tools.cut_plane import CutPlane
30 +
from floris.simulation.turbine import Ct, power, axial_induction, average_velocity
44 31
45 32
46 33
class FlorisInterface(LoggerBase):
@@ -70,7 +57,7 @@
Loading
70 57
            self.floris = Floris.from_dict(self.configuration)
71 58
72 59
        else:
73 -
            raise TypeError("The Floris `configuration` must of type 'dict', 'str', or 'Path'.")
60 +
            raise TypeError("The Floris `configuration` must be of type 'dict', 'str', or 'Path'.")
74 61
75 62
        # Store the heterogeneous map for use after reinitailization
76 63
        self.het_map = het_map
@@ -84,22 +71,32 @@
Loading
84 71
85 72
        # Make a check on reference height and provide a helpful warning
86 73
        unique_heights = np.unique(self.floris.farm.hub_heights)
87 -
        if ((len(unique_heights) == 1) and (self.floris.flow_field.reference_wind_height!=unique_heights[0])):
88 -
            err_msg = 'The only unique hub-height is not the equal to the specified reference wind height.  If this was unintended use -1 as the reference hub height to indicate use of hub-height as reference wind height.'
74 +
        if (len(unique_heights) == 1) and (self.floris.flow_field.reference_wind_height != unique_heights[0]):
75 +
            err_msg = "The only unique hub-height is not the equal to the specified reference wind height.  If this was unintended use -1 as the reference hub height to indicate use of hub-height as reference wind height."
89 76
            self.logger.warning(err_msg, stack_info=True)
90 77
78 +
        # Check the turbine_grid_points is reasonable
79 +
        if self.floris.solver["type"] == "turbine_grid":
80 +
            if self.floris.solver["turbine_grid_points"] > 3:
81 +
                self.logger.error(f"turbine_grid_points value is {self.floris.solver['turbine_grid_points']} which is larger than the recommended value of less than or equal to 3. High amounts of turbine grid points reduce the computational performance but have a small change on accuracy.")
82 +
                raise ValueError("turbine_grid_points must be less than or equal to 3.")
83 +
91 84
    def assign_hub_height_to_ref_height(self):
92 85
93 86
        # Confirm can do this operation
94 87
        unique_heights = np.unique(self.floris.farm.hub_heights)
95 -
        if (len(unique_heights) > 1):
96 -
            raise ValueError("To assign hub heights to reference height, can not have more than one specified height. Current length is {}.".format(len(unique_heights)))
88 +
        if len(unique_heights) > 1:
89 +
            raise ValueError(
90 +
                "To assign hub heights to reference height, can not have more than one specified height. Current length is {}.".format(
91 +
                    len(unique_heights)
92 +
                )
93 +
            )
97 94
98 95
        self.floris.flow_field.reference_wind_height = unique_heights[0]
99 96
100 97
    def copy(self):
101 98
        """Create an independent copy of the current FlorisInterface object"""
102 -
        return FlorisInterface(self.floris.as_dict())
99 +
        return FlorisInterface(self.floris.as_dict(), het_map=self.het_map)
103 100
104 101
    def calculate_wake(
105 102
        self,
@@ -128,9 +125,6 @@
Loading
128 125
129 126
        # TODO decide where to handle this sign issue
130 127
        if (yaw_angles is not None) and not (np.all(yaw_angles==0.)):
131 -
            if self.floris.wake.model_strings["velocity_model"] == "turbopark":
132 -
                # TODO: Implement wake steering for the TurbOPark model
133 -
                raise ValueError("Non-zero yaw angles given and for TurbOPark model; wake steering with this model is not yet implemented.")
134 128
            self.floris.farm.yaw_angles = yaw_angles
135 129
136 130
        # Initialize solution space
@@ -157,9 +151,6 @@
Loading
157 151
158 152
        # TODO decide where to handle this sign issue
159 153
        if (yaw_angles is not None) and not (np.all(yaw_angles==0.)):
160 -
            if self.floris.wake.model_strings["velocity_model"] == "turbopark":
161 -
                # TODO: Implement wake steering for the TurbOPark model
162 -
                raise ValueError("Non-zero yaw angles given and for TurbOPark model; wake steering with this model is not yet implemented.")
163 154
            self.floris.farm.yaw_angles = yaw_angles
164 155
165 156
        # Initialize solution space
@@ -180,12 +171,15 @@
Loading
180 171
        # turbulence_kinetic_energy=None,
181 172
        air_density: float | None = None,
182 173
        # wake: WakeModelManager = None,
183 -
        layout: Tuple[list[float], list[float]] | Tuple[NDArrayFloat, NDArrayFloat] | None = None,
174 +
        layout_x: list[float] | NDArrayFloat | None = None,
175 +
        layout_y: list[float] | NDArrayFloat | None = None,
184 176
        turbine_type: list | None = None,
185 177
        # turbine_id: list[str] | None = None,
186 178
        # wtg_id: list[str] | None = None,
187 179
        # with_resolution: float | None = None,
188 -
        solver_settings: dict | None = None
180 +
        solver_settings: dict | None = None,
181 +
        time_series: bool | None = False,
182 +
        layout: tuple[list[float], list[float]] | tuple[NDArrayFloat, NDArrayFloat] | None = None,
189 183
    ):
190 184
        # Export the floris object recursively as a dictionary
191 185
        floris_dict = self.floris.as_dict()
@@ -212,11 +206,22 @@
Loading
212 206
213 207
        ## Farm
214 208
        if layout is not None:
215 -
            farm_dict["layout_x"] = layout[0]
216 -
            farm_dict["layout_y"] = layout[1]
209 +
            msg = "Use the `layout_x` and `layout_y` parameters in place of `layout` because the `layout` parameter will be deprecated in 3.3."
210 +
            self.logger.warning(msg)
211 +
            layout_x = layout[0]
212 +
            layout_y = layout[1]
213 +
        if layout_x is not None:
214 +
            farm_dict["layout_x"] = layout_x
215 +
        if layout_y is not None:
216 +
            farm_dict["layout_y"] = layout_y
217 217
        if turbine_type is not None:
218 218
            farm_dict["turbine_type"] = turbine_type
219 219
220 +
        if time_series:
221 +
            flow_field_dict["time_series"] = True
222 +
        else:
223 +
            flow_field_dict["time_series"] = False
224 +
220 225
        ## Wake
221 226
        # if wake is not None:
222 227
        #     self.floris.wake = wake
@@ -300,7 +305,7 @@
Loading
300 305
        # Subset to plane
301 306
        # TODO: Seems sloppy as need more than one plane in the z-direction for GCH
302 307
        if planar_coordinate is not None:
303 -
            df = df[np.isclose(df.x3, planar_coordinate)] # , atol=0.1, rtol=0.0)]
308 +
            df = df[np.isclose(df.x3, planar_coordinate)]  # , atol=0.1, rtol=0.0)]
304 309
305 310
        # Drop duplicates
306 311
        # TODO is this still needed now that we setup a grid for just this plane?
@@ -342,7 +347,7 @@
Loading
342 347
            :py:class:`~.tools.cut_plane.CutPlane`: containing values
343 348
            of x, y, u, v, w
344 349
        """
345 -
        #TODO update docstring
350 +
        # TODO update docstring
346 351
        if wd is None:
347 352
            wd = self.floris.flow_field.wind_directions
348 353
        if ws is None:
@@ -361,9 +366,7 @@
Loading
361 366
            "flow_field_grid_points": [x_resolution, y_resolution],
362 367
            "flow_field_bounds": [x_bounds, y_bounds],
363 368
        }
364 -
        self.reinitialize(
365 -
            wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings
366 -
        )
369 +
        self.reinitialize(wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings)
367 370
368 371
        # TODO this has to be done here as it seems to be lost with reinitialize
369 372
        if yaw_angles is not None:
@@ -441,9 +444,7 @@
Loading
441 444
            "flow_field_grid_points": [y_resolution, z_resolution],
442 445
            "flow_field_bounds": [y_bounds, z_bounds],
443 446
        }
444 -
        self.reinitialize(
445 -
            wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings
446 -
        )
447 +
        self.reinitialize(wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings)
447 448
448 449
        # TODO this has to be done here as it seems to be lost with reinitialize
449 450
        if yaw_angles is not None:
@@ -502,7 +503,7 @@
Loading
502 503
            :py:class:`~.tools.cut_plane.CutPlane`: containing values
503 504
            of x, y, u, v, w
504 505
        """
505 -
        #TODO update docstring
506 +
        # TODO update docstring
506 507
        if wd is None:
507 508
            wd = self.floris.flow_field.wind_directions
508 509
        if ws is None:
@@ -521,9 +522,7 @@
Loading
521 522
            "flow_field_grid_points": [x_resolution, z_resolution],
522 523
            "flow_field_bounds": [x_bounds, z_bounds],
523 524
        }
524 -
        self.reinitialize(
525 -
            wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings
526 -
        )
525 +
        self.reinitialize(wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings)
527 526
528 527
        # TODO this has to be done here as it seems to be lost with reinitialize
529 528
        if yaw_angles is not None:
@@ -553,10 +552,14 @@
Loading
553 552
554 553
    def check_wind_condition_for_viz(self, wd=None, ws=None):
555 554
        if len(wd) > 1 or len(wd) < 1:
556 -
            raise ValueError("Wind direction input must be of length 1 for visualization. Current length is {}.".format(len(wd)))
555 +
            raise ValueError(
556 +
                "Wind direction input must be of length 1 for visualization. Current length is {}.".format(len(wd))
557 +
            )
557 558
558 559
        if len(ws) > 1 or len(ws) < 1:
559 -
            raise ValueError("Wind speed input must be of length 1 for visualization. Current length is {}.".format(len(ws)))
560 +
            raise ValueError(
561 +
                "Wind speed input must be of length 1 for visualization. Current length is {}.".format(len(ws))
562 +
            )
560 563
561 564
    def get_turbine_powers(self) -> NDArrayFloat:
562 565
        """Calculates the power at each turbine in the windfarm.
@@ -564,8 +567,14 @@
Loading
564 567
        Returns:
565 568
            NDArrayFloat: [description]
566 569
        """
570 +
571 +
        # Confirm calculate wake has been run
572 +
        if self.floris.state is not State.USED:
573 +
            raise RuntimeError(f"Can't run function `FlorisInterface.get_turbine_powers` without first running `FlorisInterface.calculate_wake`.")
574 +
567 575
        turbine_powers = power(
568 576
            air_density=self.floris.flow_field.air_density,
577 +
            ref_density_cp_ct=self.floris.farm.ref_density_cp_cts,
569 578
            velocities=self.floris.flow_field.u,
570 579
            yaw_angle=self.floris.farm.yaw_angles,
571 580
            pP=self.floris.farm.pPs,
@@ -598,8 +607,12 @@
Loading
598 607
        )
599 608
        return turbine_avg_vels
600 609
610 +
    def get_turbine_TIs(self) -> NDArrayFloat:
611 +
        return self.floris.flow_field.turbulence_intensity_field
612 +
601 613
    def get_farm_power(
602 614
        self,
615 +
        turbine_weights=None,
603 616
        use_turbulence_correction=False,
604 617
    ):
605 618
        """
@@ -610,6 +623,19 @@
Loading
610 623
        original wind direction and yaw angles.
611 624
612 625
        Args:
626 +
            turbine_weights (NDArrayFloat | list[float] | None, optional):
627 +
                weighing terms that allow the user to emphasize power at
628 +
                particular turbines and/or completely ignore the power 
629 +
                from other turbines. This is useful when, for example, you are
630 +
                modeling multiple wind farms in a single floris object. If you
631 +
                only want to calculate the power production for one of those
632 +
                farms and include the wake effects of the neighboring farms,
633 +
                you can set the turbine_weights for the neighboring farms'
634 +
                turbines to 0.0. The array of turbine powers from floris
635 +
                is multiplied with this array in the calculation of the
636 +
                objective function. If None, this  is an array with all values
637 +
                1.0 and with shape equal to (n_wind_directions, n_wind_speeds,
638 +
                n_turbines). Defaults to None.
613 639
            use_turbulence_correction: (bool, optional): When *True* uses a
614 640
                turbulence parameter to adjust power output calculations.
615 641
                Defaults to *False*.
@@ -624,7 +650,34 @@
Loading
624 650
        # for turbine in self.floris.farm.turbines:
625 651
        #     turbine.use_turbulence_correction = use_turbulence_correction
626 652
653 +
        # Confirm calculate wake has been run
654 +
        if self.floris.state is not State.USED:
655 +
            raise RuntimeError(f"Can't run function `FlorisInterface.get_turbine_powers` without running `FlorisInterface.calculate_wake`.")
656 +
657 +
        if turbine_weights is None:
658 +
            # Default to equal weighing of all turbines when turbine_weights is None
659 +
            turbine_weights = np.ones(
660 +
                (
661 +
                    self.floris.flow_field.n_wind_directions,
662 +
                    self.floris.flow_field.n_wind_speeds,
663 +
                    self.floris.farm.n_turbines
664 +
                )
665 +
            )
666 +
        elif len(np.shape(turbine_weights)) == 1:
667 +
            # Deal with situation when 1D array is provided
668 +
            turbine_weights = np.tile(
669 +
                turbine_weights,
670 +
                (
671 +
                    self.floris.flow_field.n_wind_directions,
672 +
                    self.floris.flow_field.n_wind_speeds,
673 +
                    1
674 +
                )
675 +
            )
676 +
677 +
        # Calculate all turbine powers and apply weights
627 678
        turbine_powers = self.get_turbine_powers()
679 +
        turbine_powers = np.multiply(turbine_weights, turbine_powers)
680 +
628 681
        return np.sum(turbine_powers, axis=2)
629 682
630 683
    def get_farm_AEP(
@@ -633,6 +686,7 @@
Loading
633 686
        cut_in_wind_speed=0.001,
634 687
        cut_out_wind_speed=None,
635 688
        yaw_angles=None,
689 +
        turbine_weights=None,
636 690
        no_wake=False,
637 691
    ) -> float:
638 692
        """
@@ -646,7 +700,7 @@
Loading
646 700
                up to 1.0 and are used to weigh the wind farm power for every
647 701
                condition in calculating the wind farm's AEP.
648 702
            cut_in_wind_speed (float, optional): Wind speed in m/s below which
649 -
                any calculations are ignored and the wind farm is known to 
703 +
                any calculations are ignored and the wind farm is known to
650 704
                produce 0.0 W of power. Note that to prevent problems with the
651 705
                wake models at negative / zero wind speeds, this variable must
652 706
                always have a positive value. Defaults to 0.001 [m/s].
@@ -658,13 +712,26 @@
Loading
658 712
                The relative turbine yaw angles in degrees. If None is
659 713
                specified, will assume that the turbine yaw angles are all
660 714
                zero degrees for all conditions. Defaults to None.
715 +
            turbine_weights (NDArrayFloat | list[float] | None, optional):
716 +
                weighing terms that allow the user to emphasize power at
717 +
                particular turbines and/or completely ignore the power 
718 +
                from other turbines. This is useful when, for example, you are
719 +
                modeling multiple wind farms in a single floris object. If you
720 +
                only want to calculate the power production for one of those
721 +
                farms and include the wake effects of the neighboring farms,
722 +
                you can set the turbine_weights for the neighboring farms'
723 +
                turbines to 0.0. The array of turbine powers from floris
724 +
                is multiplied with this array in the calculation of the
725 +
                objective function. If None, this  is an array with all values
726 +
                1.0 and with shape equal to (n_wind_directions, n_wind_speeds,
727 +
                n_turbines). Defaults to None.
661 728
            no_wake: (bool, optional): When *True* updates the turbine
662 729
                quantities without calculating the wake or adding the wake to
663 730
                the flow field. This can be useful when quantifying the loss
664 731
                in AEP due to wakes. Defaults to *False*.
665 732
666 733
        Returns:
667 -
            float: 
734 +
            float:
668 735
                The Annual Energy Production (AEP) for the wind farm in
669 736
                watt-hours.
670 737
        """
@@ -676,30 +743,22 @@
Loading
676 743
            & (len(np.shape(freq)) == 2)
677 744
        ):
678 745
            raise UserWarning(
679 -
                "'freq' should be a two-dimensional array with dimensions"
680 -
                + " (n_wind_directions, n_wind_speeds)."
746 +
                "'freq' should be a two-dimensional array with dimensions" + " (n_wind_directions, n_wind_speeds)."
681 747
            )
682 748
683 749
        # Check if frequency vector sums to 1.0. If not, raise a warning
684 750
        if np.abs(np.sum(freq) - 1.0) > 0.001:
685 -
            self.logger.warning(
686 -
                "WARNING: The frequency array provided to get_farm_AEP() "
687 -
                + "does not sum to 1.0. "
688 -
            )
751 +
            self.logger.warning("WARNING: The frequency array provided to get_farm_AEP() " + "does not sum to 1.0. ")
689 752
690 753
        # Copy the full wind speed array from the floris object and initialize
691 754
        # the the farm_power variable as an empty array.
692 755
        wind_speeds = np.array(self.floris.flow_field.wind_speeds, copy=True)
693 -
        farm_power = np.zeros(
694 -
            (self.floris.flow_field.n_wind_directions, len(wind_speeds))
695 -
        )
756 +
        farm_power = np.zeros((self.floris.flow_field.n_wind_directions, len(wind_speeds)))
696 757
697 758
        # Determine which wind speeds we must evaluate in floris
698 -
        conditions_to_evaluate = (wind_speeds >= cut_in_wind_speed)
759 +
        conditions_to_evaluate = wind_speeds >= cut_in_wind_speed
699 760
        if cut_out_wind_speed is not None:
700 -
            conditions_to_evaluate = conditions_to_evaluate & (
701 -
                wind_speeds < cut_out_wind_speed
702 -
            )
761 +
            conditions_to_evaluate = conditions_to_evaluate & (wind_speeds < cut_out_wind_speed)
703 762
704 763
        # Evaluate the conditions in floris
705 764
        if np.any(conditions_to_evaluate):
@@ -712,7 +771,9 @@
Loading
712 771
                self.calculate_no_wake(yaw_angles=yaw_angles_subset)
713 772
            else:
714 773
                self.calculate_wake(yaw_angles=yaw_angles_subset)
715 -
            farm_power[:, conditions_to_evaluate] = self.get_farm_power()
774 +
            farm_power[:, conditions_to_evaluate] = (
775 +
                self.get_farm_power(turbine_weights=turbine_weights)
776 +
            )
716 777
717 778
        # Finally, calculate AEP in GWh
718 779
        aep = np.sum(np.multiply(freq, farm_power) * 365 * 24)
@@ -722,6 +783,74 @@
Loading
722 783
723 784
        return aep
724 785
786 +
    def get_farm_AEP_wind_rose_class(
787 +
        self,
788 +
        wind_rose,
789 +
        cut_in_wind_speed=0.001,
790 +
        cut_out_wind_speed=None,
791 +
        yaw_angles=None,
792 +
        no_wake=False,
793 +
    ) -> float:
794 +
        """
795 +
        Estimate annual energy production (AEP) for distributions of wind speed, wind
796 +
        direction, frequency of occurrence, and yaw offset.
797 +
798 +
        Args:
799 +
            wind_rose (wind_rose): An object of the wind rose class
800 +
            cut_in_wind_speed (float, optional): Wind speed in m/s below which
801 +
                any calculations are ignored and the wind farm is known to 
802 +
                produce 0.0 W of power. Note that to prevent problems with the
803 +
                wake models at negative / zero wind speeds, this variable must
804 +
                always have a positive value. Defaults to 0.001 [m/s].
805 +
            cut_out_wind_speed (float, optional): Wind speed above which the
806 +
                wind farm is known to produce 0.0 W of power. If None is
807 +
                specified, will assume that the wind farm does not cut out
808 +
                at high wind speeds. Defaults to None.
809 +
            yaw_angles (NDArrayFloat | list[float] | None, optional):
810 +
                The relative turbine yaw angles in degrees. If None is
811 +
                specified, will assume that the turbine yaw angles are all
812 +
                zero degrees for all conditions. Defaults to None.
813 +
            no_wake: (bool, optional): When *True* updates the turbine
814 +
                quantities without calculating the wake or adding the wake to
815 +
                the flow field. This can be useful when quantifying the loss
816 +
                in AEP due to wakes. Defaults to *False*.
817 +
818 +
        Returns:
819 +
            float: 
820 +
                The Annual Energy Production (AEP) for the wind farm in
821 +
                watt-hours.
822 +
        """
823 +
824 +
        # Hold the starting values of wind speed and direction
825 +
        wind_speeds = np.array(self.floris.flow_field.wind_speeds, copy=True)
826 +
        wind_directions = np.array(self.floris.flow_field.wind_directions, copy=True)
827 +
828 +
        # Now set FLORIS wind speed and wind direction
829 +
        # over to those values in the wind rose class
830 +
        wind_speeds_wind_rose = wind_rose.df.ws.unique()
831 +
        wind_directions_wind_rose = wind_rose.df.wd.unique()
832 +
        self.reinitialize(wind_speeds=wind_speeds_wind_rose, wind_directions=wind_directions_wind_rose)
833 +
834 +
        # Build the frequency matrix from wind rose
835 +
        freq = wind_rose.df.set_index(['wd','ws']).unstack().values
836 +
837 +
        # Now compute aep
838 +
        aep = self.get_farm_AEP(
839 +
            freq,
840 +
            cut_in_wind_speed=cut_in_wind_speed,
841 +
            cut_out_wind_speed=cut_out_wind_speed,
842 +
            yaw_angles=yaw_angles,
843 +
            no_wake=no_wake)
844 +
845 +
846 +
        # Reset the FLORIS object to the original wind speed and directions
847 +
        self.reinitialize(wind_speeds=wind_speeds, wind_directions=wind_directions)
848 +
        
849 +
850 +
        return aep
851 +
852 +
853 +
725 854
    @property
726 855
    def layout_x(self):
727 856
        """
@@ -742,7 +871,6 @@
Loading
742 871
        """
743 872
        return self.floris.farm.layout_y
744 873
745 -
746 874
    def get_turbine_layout(self, z=False):
747 875
        """
748 876
        Get turbine layout
@@ -778,759 +906,6 @@
Loading
778 906
779 907
    return [in_region, out_region]
780 908
781 -
# def global_calc_one_AEP_case(FlorisInterface, wd, ws, freq, yaw=None):
782 -
#     return FlorisInterface._calc_one_AEP_case(wd, ws, freq, yaw)
783 -
784 -
DEFAULT_UNCERTAINTY = {"std_wd": 4.95, "std_yaw": 1.75, "pmf_res": 1.0, "pdf_cutoff": 0.995}
785 -
786 -
787 -
def _generate_uncertainty_parameters(unc_options: dict, unc_pmfs: dict) -> dict:
788 -
    """Generates the uncertainty parameters for `FlorisInterface.get_farm_power` and
789 -
    `FlorisInterface.get_turbine_power` for more details.
790 -
791 -
    Args:
792 -
        unc_options (dict): See `FlorisInterface.get_farm_power` or `FlorisInterface.get_turbine_power`.
793 -
        unc_pmfs (dict): See `FlorisInterface.get_farm_power` or `FlorisInterface.get_turbine_power`.
794 -
795 -
    Returns:
796 -
        dict: [description]
797 -
    """
798 -
    if (unc_options is None) & (unc_pmfs is None):
799 -
        unc_options = DEFAULT_UNCERTAINTY
800 -
801 -
    if unc_pmfs is not None:
802 -
        return unc_pmfs
803 -
804 -
    wd_unc = np.zeros(1)
805 -
    wd_unc_pmf = np.ones(1)
806 -
    yaw_unc = np.zeros(1)
807 -
    yaw_unc_pmf = np.ones(1)
808 -
809 -
    # create normally distributed wd and yaw uncertaitny pmfs if appropriate
810 -
    if unc_options["std_wd"] > 0:
811 -
        wd_bnd = int(np.ceil(norm.ppf(unc_options["pdf_cutoff"], scale=unc_options["std_wd"]) / unc_options["pmf_res"]))
812 -
        bound = wd_bnd * unc_options["pmf_res"]
813 -
        wd_unc = np.linspace(-1 * bound, bound, 2 * wd_bnd + 1)
814 -
        wd_unc_pmf = norm.pdf(wd_unc, scale=unc_options["std_wd"])
815 -
        wd_unc_pmf /= np.sum(wd_unc_pmf)  # normalize so sum = 1.0
816 -
817 -
    if unc_options["std_yaw"] > 0:
818 -
        yaw_bnd = int(
819 -
            np.ceil(norm.ppf(unc_options["pdf_cutoff"], scale=unc_options["std_yaw"]) / unc_options["pmf_res"])
820 -
        )
821 -
        bound = yaw_bnd * unc_options["pmf_res"]
822 -
        yaw_unc = np.linspace(-1 * bound, bound, 2 * yaw_bnd + 1)
823 -
        yaw_unc_pmf = norm.pdf(yaw_unc, scale=unc_options["std_yaw"])
824 -
        yaw_unc_pmf /= np.sum(yaw_unc_pmf)  # normalize so sum = 1.0
825 -
826 -
    unc_pmfs = {
827 -
        "wd_unc": wd_unc,
828 -
        "wd_unc_pmf": wd_unc_pmf,
829 -
        "yaw_unc": yaw_unc,
830 -
        "yaw_unc_pmf": yaw_unc_pmf,
831 -
    }
832 -
    return unc_pmfs
833 -
834 -
835 -
# def correct_for_all_combinations(
836 -
#     wd: NDArrayFloat,
837 -
#     ws: NDArrayFloat,
838 -
#     freq: NDArrayFloat,
839 -
#     yaw: NDArrayFloat | None = None,
840 -
# ) -> tuple[NDArrayFloat]:
841 -
#     """Computes the probabilities for the complete windrose from the desired wind
842 -
#     direction and wind speed combinations and their associated probabilities so that
843 -
#     any undesired combinations are filled with a 0.0 probability.
844 -
845 -
#     Args:
846 -
#         wd (NDArrayFloat): List or array of wind direction values.
847 -
#         ws (NDArrayFloat): List or array of wind speed values.
848 -
#         freq (NDArrayFloat): Frequencies corresponding to wind
849 -
#             speeds and directions in wind rose with dimensions
850 -
#             (N wind directions x N wind speeds).
851 -
#         yaw (NDArrayFloat | None): The corresponding yaw angles for each of the wind
852 -
#             direction and wind speed combinations, or None. Defaults to None.
853 -
854 -
#     Returns:
855 -
#         NDArrayFloat, NDArrayFloat, NDArrayFloat: The unique wind directions, wind
856 -
#             speeds, and the associated probability of their combination combinations in
857 -
#             an array of shape (N wind directions x N wind speeds).
858 -
#     """
859 -
860 -
#     combos_to_compute = np.array(list(zip(wd, ws, freq)))
861 -
862 -
#     unique_wd = wd.unique()
863 -
#     unique_ws = ws.unique()
864 -
#     all_combos = np.array(list(product(unique_wd, unique_ws)), dtype=float)
865 -
#     all_combos = np.hstack((all_combos, np.zeros((all_combos.shape[0], 1), dtype=float)))
866 -
#     expanded_yaw = np.array([None] * all_combos.shape[0]).reshape(unique_wd.size, unique_ws.size)
867 -
868 -
#     ix_match = [np.where((all_combos[:, :2] == combo[:2]).all(1))[0][0] for combo in combos_to_compute]
869 -
#     all_combos[ix_match, 2] = combos_to_compute[:, 2]
870 -
#     if yaw is not None:
871 -
#         expanded_yaw[ix_match] = yaw
872 -
#     freq = all_combos.T[2].reshape((unique_wd.size, unique_ws.size))
873 -
#     return unique_wd, unique_ws, freq
874 -
875 -
876 -
    # def get_set_of_points(self, x_points, y_points, z_points):
877 -
    #     """
878 -
    #     Calculates velocity values through the
879 -
    #     :py:meth:`~.FlowField.calculate_wake` method at points specified by
880 -
    #     inputs.
881 -
882 -
    #     Args:
883 -
    #         x_points (float): X-locations to get velocity values at.
884 -
    #         y_points (float): Y-locations to get velocity values at.
885 -
    #         z_points (float): Z-locations to get velocity values at.
886 -
887 -
    #     Returns:
888 -
    #         :py:class:`pandas.DataFrame`: containing values of x, y, z, u, v, w
889 -
    #     """
890 -
    #     # Get a copy for the flow field so don't change underlying grid points
891 -
    #     flow_field = copy.deepcopy(self.floris.flow_field)
892 -
893 -
    #     if hasattr(self.floris.wake.velocity_model, "requires_resolution"):
894 -
    #         if self.floris.velocity_model.requires_resolution:
895 -
896 -
    #             # If this is a gridded model, must extract from full flow field
897 -
    #             self.logger.info(
898 -
    #                 "Model identified as %s requires use of underlying grid print"
899 -
    #                 % self.floris.wake.velocity_model.model_string
900 -
    #             )
901 -
    #             self.logger.warning("FUNCTION NOT AVAILABLE CURRENTLY")
902 -
903 -
    #     # Set up points matrix
904 -
    #     points = np.row_stack((x_points, y_points, z_points))
905 -
906 -
    #     # TODO: Calculate wake inputs need to be mapped
907 -
    #     raise_error = True
908 -
    #     if raise_error:
909 -
    #         raise NotImplementedError("Additional point calculation is not yet supported!")
910 -
    #     # Recalculate wake with these points
911 -
    #     flow_field.calculate_wake(points=points)
912 -
913 -
    #     # Get results vectors
914 -
    #     x_flat = flow_field.x.flatten()
915 -
    #     y_flat = flow_field.y.flatten()
916 -
    #     z_flat = flow_field.z.flatten()
917 -
    #     u_flat = flow_field.u.flatten()
918 -
    #     v_flat = flow_field.v.flatten()
919 -
    #     w_flat = flow_field.w.flatten()
920 -
921 -
    #     df = pd.DataFrame(
922 -
    #         {
923 -
    #             "x": x_flat,
924 -
    #             "y": y_flat,
925 -
    #             "z": z_flat,
926 -
    #             "u": u_flat,
927 -
    #             "v": v_flat,
928 -
    #             "w": w_flat,
929 -
    #         }
930 -
    #     )
931 -
932 -
    #     # Subset to points requests
933 -
    #     df = df[df.x.isin(x_points)]
934 -
    #     df = df[df.y.isin(y_points)]
935 -
    #     df = df[df.z.isin(z_points)]
936 -
937 -
    #     # Drop duplicates
938 -
    #     df = df.drop_duplicates()
939 -
940 -
    #     # Return the dataframe
941 -
    #     return df
942 -
    
943 -
    # def get_flow_data(self, resolution=None, grid_spacing=10, velocity_deficit=False):
944 -
    #     """
945 -
    #     Generate :py:class:`~.tools.flow_data.FlowData` object corresponding to
946 -
    #     active FLORIS instance.
947 -
948 -
    #     Velocity and wake models requiring calculation on a grid implement a
949 -
    #     discretized domain at resolution **grid_spacing**. This is distinct
950 -
    #     from the resolution of the returned flow field domain.
951 -
952 -
    #     Args:
953 -
    #         resolution (float, optional): Resolution of output data.
954 -
    #             Only used for wake models that require spatial
955 -
    #             resolution (e.g. curl). Defaults to None.
956 -
    #         grid_spacing (int, optional): Resolution of grid used for
957 -
    #             simulation. Model results may be sensitive to resolution.
958 -
    #             Defaults to 10.
959 -
    #         velocity_deficit (bool, optional): When *True*, normalizes velocity
960 -
    #             with respect to initial flow field velocity to show relative
961 -
    #             velocity deficit (%). Defaults to *False*.
962 -
963 -
    #     Returns:
964 -
    #         :py:class:`~.tools.flow_data.FlowData`: FlowData object
965 -
    #     """
966 -
967 -
    #     if resolution is None:
968 -
    #         if not self.floris.wake.velocity_model.requires_resolution:
969 -
    #             self.logger.info("Assuming grid with spacing %d" % grid_spacing)
970 -
    #             (
971 -
    #                 xmin,
972 -
    #                 xmax,
973 -
    #                 ymin,
974 -
    #                 ymax,
975 -
    #                 zmin,
976 -
    #                 zmax,
977 -
    #             ) = self.floris.flow_field.domain_bounds  # TODO: No grid attribute within FlowField
978 -
    #             resolution = Vec3(
979 -
    #                 1 + (xmax - xmin) / grid_spacing,
980 -
    #                 1 + (ymax - ymin) / grid_spacing,
981 -
    #                 1 + (zmax - zmin) / grid_spacing,
982 -
    #             )
983 -
    #         else:
984 -
    #             self.logger.info("Assuming model resolution")
985 -
    #             resolution = self.floris.wake.velocity_model.model_grid_resolution
986 -
987 -
    #     # Get a copy for the flow field so don't change underlying grid points
988 -
    #     flow_field = copy.deepcopy(self.floris.flow_field)
989 -
990 -
    #     if (
991 -
    #         flow_field.wake.velocity_model.requires_resolution
992 -
    #         and flow_field.wake.velocity_model.model_grid_resolution != resolution
993 -
    #     ):
994 -
    #         self.logger.warning(
995 -
    #             "WARNING: The current wake velocity model contains a "
996 -
    #             + "required grid resolution; the Resolution given to "
997 -
    #             + "FlorisInterface.get_flow_field is ignored."
998 -
    #         )
999 -
    #         resolution = flow_field.wake.velocity_model.model_grid_resolution
1000 -
    #     flow_field.reinitialize(with_resolution=resolution)  # TODO: Not implemented
1001 -
    #     self.logger.info(resolution)
1002 -
    #     # print(resolution)
1003 -
    #     flow_field.steady_state_atmospheric_condition()
1004 -
1005 -
    #     order = "f"
1006 -
    #     x = flow_field.x.flatten(order=order)
1007 -
    #     y = flow_field.y.flatten(order=order)
1008 -
    #     z = flow_field.z.flatten(order=order)
1009 -
1010 -
    #     u = flow_field.u.flatten(order=order)
1011 -
    #     v = flow_field.v.flatten(order=order)
1012 -
    #     w = flow_field.w.flatten(order=order)
1013 -
1014 -
    #     # find percent velocity deficit
1015 -
    #     if velocity_deficit:
1016 -
    #         u = abs(u - flow_field.u_initial.flatten(order=order)) / flow_field.u_initial.flatten(order=order) * 100
1017 -
    #         v = abs(v - flow_field.v_initial.flatten(order=order)) / flow_field.v_initial.flatten(order=order) * 100
1018 -
    #         w = abs(w - flow_field.w_initial.flatten(order=order)) / flow_field.w_initial.flatten(order=order) * 100
1019 -
1020 -
    #     # Determine spacing, dimensions and origin
1021 -
    #     unique_x = np.sort(np.unique(x))
1022 -
    #     unique_y = np.sort(np.unique(y))
1023 -
    #     unique_z = np.sort(np.unique(z))
1024 -
    #     spacing = Vec3(
1025 -
    #         unique_x[1] - unique_x[0],
1026 -
    #         unique_y[1] - unique_y[0],
1027 -
    #         unique_z[1] - unique_z[0],
1028 -
    #     )
1029 -
    #     dimensions = Vec3(len(unique_x), len(unique_y), len(unique_z))
1030 -
    #     origin = Vec3(0.0, 0.0, 0.0)
1031 -
    #     return FlowData(x, y, z, u, v, w, spacing=spacing, dimensions=dimensions, origin=origin)
1032 -
1033 -
1034 -
    # def get_turbine_power(
1035 -
    #     self,
1036 -
    #     include_unc=False,
1037 -
    #     unc_pmfs=None,
1038 -
    #     unc_options=None,
1039 -
    #     no_wake=False,
1040 -
    #     use_turbulence_correction=False,
1041 -
    # ):
1042 -
    #     """
1043 -
    #     Report power from each wind turbine.
1044 -
1045 -
    #     Args:
1046 -
    #         include_unc (bool): If *True*, uncertainty in wind direction
1047 -
    #             and/or yaw position is included when determining turbine
1048 -
    #             powers. Defaults to *False*.
1049 -
    #         unc_pmfs (dictionary, optional): A dictionary containing optional
1050 -
    #             probability mass functions describing the distribution of wind
1051 -
    #             direction and yaw position deviations when wind direction and/or
1052 -
    #             yaw position uncertainty is included in the power calculations.
1053 -
    #             Contains the following key-value pairs:
1054 -
1055 -
    #             -   **wd_unc** (*np.array*): Wind direction deviations from the
1056 -
    #                 original wind direction.
1057 -
    #             -   **wd_unc_pmf** (*np.array*): Probability of each wind
1058 -
    #                 direction deviation in **wd_unc** occuring.
1059 -
    #             -   **yaw_unc** (*np.array*): Yaw angle deviations from the
1060 -
    #                 original yaw angles.
1061 -
    #             -   **yaw_unc_pmf** (*np.array*): Probability of each yaw angle
1062 -
    #                 deviation in **yaw_unc** occuring.
1063 -
1064 -
    #             Defaults to None, in which case default PMFs are calculated
1065 -
    #             using values provided in **unc_options**.
1066 -
    #         unc_options (dictionary, optional): A dictionary containing values
1067 -
    #             used to create normally-distributed, zero-mean probability mass
1068 -
    #             functions describing the distribution of wind direction and yaw
1069 -
    #             position deviations when wind direction and/or yaw position
1070 -
    #             uncertainty is included. This argument is only used when
1071 -
    #             **unc_pmfs** is None and contains the following key-value pairs:
1072 -
1073 -
    #             -   **std_wd** (*float*): A float containing the standard
1074 -
    #                 deviation of the wind direction deviations from the
1075 -
    #                 original wind direction.
1076 -
    #             -   **std_yaw** (*float*): A float containing the standard
1077 -
    #                 deviation of the yaw angle deviations from the original yaw
1078 -
    #                 angles.
1079 -
    #             -   **pmf_res** (*float*): A float containing the resolution in
1080 -
    #                 degrees of the wind direction and yaw angle PMFs.
1081 -
    #             -   **pdf_cutoff** (*float*): A float containing the cumulative
1082 -
    #                 distribution function value at which the tails of the
1083 -
    #                 PMFs are truncated.
1084 -
1085 -
    #             Defaults to None. Initializes to {'std_wd': 4.95, 'std_yaw': 1.
1086 -
    #             75, 'pmf_res': 1.0, 'pdf_cutoff': 0.995}.
1087 -
    #         no_wake: (bool, optional): When *True* updates the turbine
1088 -
    #             quantities without calculating the wake or adding the
1089 -
    #             wake to the flow field. Defaults to *False*.
1090 -
    #         use_turbulence_correction: (bool, optional): When *True* uses a
1091 -
    #             turbulence parameter to adjust power output calculations.
1092 -
    #             Defaults to *False*.
1093 -
1094 -
    #     Returns:
1095 -
    #         np.array: Power produced by each wind turbine.
1096 -
    #     """
1097 -
    #     # TODO: Turbulence correction used in the power calculation, but may not be in
1098 -
    #     # the model yet
1099 -
    #     # TODO: Turbines need a switch for using turbulence correction
1100 -
    #     # TODO: Uncomment out the following two lines once the above are resolved
1101 -
    #     # for turbine in self.floris.farm.turbines:
1102 -
    #     #     turbine.use_turbulence_correction = use_turbulence_correction
1103 -
1104 -
    #     if include_unc:
1105 -
    #         unc_pmfs = _generate_uncertainty_parameters(unc_options, unc_pmfs)
1106 -
1107 -
    #         mean_farm_power = np.zeros(self.floris.farm.n_turbines)
1108 -
    #         wd_orig = self.floris.flow_field.wind_directions  # TODO: same comment as in get_farm_power
1109 -
1110 -
    #         yaw_angles = self.get_yaw_angles()
1111 -
    #         self.reinitialize(wind_direction=wd_orig[0] + unc_pmfs["wd_unc"])
1112 -
    #         for i, delta_yaw in enumerate(unc_pmfs["yaw_unc"]):
1113 -
    #             self.calculate_wake(
1114 -
    #                 yaw_angles=list(np.array(yaw_angles) + delta_yaw),
1115 -
    #                 no_wake=no_wake,
1116 -
    #             )
1117 -
    #             mean_farm_power += unc_pmfs["wd_unc_pmf"] * unc_pmfs["yaw_unc_pmf"][i] * self._get_turbine_powers()
1118 -
1119 -
    #         # reinitialize with original values
1120 -
    #         self.reinitialize(wind_direction=wd_orig)
1121 -
    #         self.calculate_wake(yaw_angles=yaw_angles, no_wake=no_wake)
1122 -
    #         return mean_farm_power
1123 -
1124 -
    #     return self._get_turbine_powers()
1125 -
1126 -
    # def get_power_curve(self, wind_speeds):
1127 -
    #     """
1128 -
    #     Return the power curve given a set of wind speeds
1129 -
1130 -
    #     Args:
1131 -
    #         wind_speeds (np.array): array of wind speeds to get power curve
1132 -
    #     """
1133 -
1134 -
    #     # TODO: Why is this done? Should we expand for evenutal multiple turbines types
1135 -
    #     # or just allow a filter on the turbine index?
1136 -
    #     # Temporarily set the farm to a single turbine
1137 -
    #     saved_layout_x = self.layout_x
1138 -
    #     saved_layout_y = self.layout_y
1139 -
1140 -
    #     self.reinitialize(wind_speed=wind_speeds, layout_array=([0], [0]))
1141 -
    #     self.calculate_wake()
1142 -
    #     turbine_power = self._get_turbine_powers()
1143 -
1144 -
    #     # Set it back
1145 -
    #     self.reinitialize(layout_array=(saved_layout_x, saved_layout_y))
1146 -
1147 -
    #     return turbine_power
1148 -
1149 -
    # def get_farm_power_for_yaw_angle(
1150 -
    #     self,
1151 -
    #     yaw_angles,
1152 -
    #     include_unc=False,
1153 -
    #     unc_pmfs=None,
1154 -
    #     unc_options=None,
1155 -
    #     no_wake=False,
1156 -
    # ):
1157 -
    #     """
1158 -
    #     Assign yaw angles to turbines, calculate wake, and report farm power.
1159 -
1160 -
    #     Args:
1161 -
    #         yaw_angles (np.array): Yaw to apply to each turbine.
1162 -
    #         include_unc (bool, optional): When *True*, includes wind direction
1163 -
    #             uncertainty in estimate of wind farm power. Defaults to *False*.
1164 -
    #         unc_pmfs (dictionary, optional): A dictionary containing optional
1165 -
    #             probability mass functions describing the distribution of wind
1166 -
    #             direction and yaw position deviations when wind direction and/or
1167 -
    #             yaw position uncertainty is included in the power calculations.
1168 -
    #             Contains the following key-value pairs:
1169 -
1170 -
    #             -   **wd_unc** (*np.array*): Wind direction deviations from the
1171 -
    #                 original wind direction.
1172 -
    #             -   **wd_unc_pmf** (*np.array*): Probability of each wind
1173 -
    #                 direction deviation in **wd_unc** occuring.
1174 -
    #             -   **yaw_unc** (*np.array*): Yaw angle deviations from the
1175 -
    #                 original yaw angles.
1176 -
    #             -   **yaw_unc_pmf** (*np.array*): Probability of each yaw angle
1177 -
    #                 deviation in **yaw_unc** occuring.
1178 -
1179 -
    #             Defaults to None, in which case default PMFs are calculated
1180 -
    #             using values provided in **unc_options**.
1181 -
    #         unc_options (dictionary, optional): A dictionary containing values
1182 -
    #             used to create normally-distributed, zero-mean probability mass
1183 -
    #             functions describing the distribution of wind direction and yaw
1184 -
    #             position deviations when wind direction and/or yaw position
1185 -
    #             uncertainty is included. This argument is only used when
1186 -
    #             **unc_pmfs** is None and contains the following key-value pairs:
1187 -
1188 -
    #             -   **std_wd** (*float*): A float containing the standard
1189 -
    #                 deviation of the wind direction deviations from the
1190 -
    #                 original wind direction.
1191 -
    #             -   **std_yaw** (*float*): A float containing the standard
1192 -
    #                 deviation of the yaw angle deviations from the original yaw
1193 -
    #                 angles.
1194 -
    #             -   **pmf_res** (*float*): A float containing the resolution in
1195 -
    #                 degrees of the wind direction and yaw angle PMFs.
1196 -
    #             -   **pdf_cutoff** (*float*): A float containing the cumulative
1197 -
    #                 distribution function value at which the tails of the
1198 -
    #                 PMFs are truncated.
1199 -
1200 -
    #             Defaults to None. Initializes to {'std_wd': 4.95, 'std_yaw': 1.
1201 -
    #             75, 'pmf_res': 1.0, 'pdf_cutoff': 0.995}.
1202 -
    #         no_wake: (bool, optional): When *True* updates the turbine
1203 -
    #             quantities without calculating the wake or adding the
1204 -
    #             wake to the flow field. Defaults to *False*.
1205 -
1206 -
    #     Returns:
1207 -
    #         float: Wind plant power. #TODO negative? in kW?
1208 -
    #     """
1209 -
1210 -
    #     self.calculate_wake(yaw_angles=yaw_angles, no_wake=no_wake)
1211 -
1212 -
    #     return self.get_farm_power(include_unc=include_unc, unc_pmfs=unc_pmfs, unc_options=unc_options)
1213 -
1214 -
    # def copy_and_update_turbine_map(
1215 -
    #     self, base_turbine_id: str, update_parameters: dict, new_id: str | None = None
1216 -
    # ) -> dict:
1217 -
    #     """Creates a new copy of an existing turbine and updates the parameters based on
1218 -
    #     user input. This function is a helper to make the v2 -> v3 transition easier.
1219 -
1220 -
    #     Args:
1221 -
    #         base_turbine_id (str): The base turbine's ID in `floris.farm.turbine_id`.
1222 -
    #         update_parameters (dict): A dictionary of the turbine parameters to update
1223 -
    #             and their new valies.
1224 -
    #         new_id (str, optional): The new `turbine_id`, if `None` a unique
1225 -
    #             identifier will be appended to the end. Defaults to None.
1226 -
1227 -
    #     Returns:
1228 -
    #         dict: A turbine mapping that can be passed directly to `change_turbine`.
1229 -
    #     """
1230 -
    #     if new_id is None:
1231 -
    #         new_id = f"{base_turbine_id}_copy{self.unique_copy_id}"
1232 -
    #         self.unique_copy_id += 1
1233 -
1234 -
    #     turbine = {new_id: self.floris.turbine[base_turbine_id]._asdict()}
1235 -
    #     turbine[new_id].update(update_parameters)
1236 -
    #     return turbine
1237 -
1238 -
    # def change_turbine(
1239 -
    #     self,
1240 -
    #     turbine_indices: list[int],
1241 -
    #     new_turbine_map: dict[str, dict[str, Any]],
1242 -
    #     update_specified_wind_height: bool = False,
1243 -
    # ):
1244 -
    #     """
1245 -
    #     Change turbine properties for specified turbines.
1246 -
1247 -
    #     Args:
1248 -
    #         turbine_indices (list[int]): List of turbine indices to change.
1249 -
    #         new_turbine_map (dict[str, dict[str, Any]]): New dictionary of turbine
1250 -
    #             parameters to create the new turbines for each of `turbine_indices`.
1251 -
    #         update_specified_wind_height (bool, optional): When *True*, update specified
1252 -
    #             wind height to match new hub_height. Defaults to *False*.
1253 -
    #     """
1254 -
    #     new_turbine = True
1255 -
    #     new_turbine_id = [*new_turbine_map][0]
1256 -
    #     if new_turbine_id in self.floris.farm.turbine_map:
1257 -
    #         new_turbine = False
1258 -
    #         self.logger.info(f"Turbines {turbine_indices} will be re-mapped to the definition for: {new_turbine_id}")
1259 -
1260 -
    #     self.floris.farm.turbine_id = [
1261 -
    #         new_turbine_id if i in turbine_indices else t_id for i, t_id in enumerate(self.floris.farm.turbine_id)
1262 -
    #     ]
1263 -
    #     if new_turbine:
1264 -
    #         self.logger.info(f"Turbines {turbine_indices} have been mapped to the new definition for: {new_turbine_id}")
1265 -
1266 -
    #     # Update the turbine mapping if a new turbine was provided, then regenerate the
1267 -
    #     # farm arrays for the turbine farm
1268 -
    #     if new_turbine:
1269 -
    #         turbine_map = self.floris.farm._asdict()["turbine_map"]
1270 -
    #         turbine_map.update(new_turbine_map)
1271 -
    #         self.floris.farm.turbine_map = turbine_map
1272 -
    #     self.floris.farm.generate_farm_points()
1273 -
1274 -
    #     new_hub_height = new_turbine_map[new_turbine_id]["hub_height"]
1275 -
    #     changed_hub_height = new_hub_height != self.floris.flow_field.reference_wind_height
1276 -
1277 -
    #     # Alert user if changing hub-height and not specified wind height
1278 -
    #     if changed_hub_height and not update_specified_wind_height:
1279 -
    #         self.logger.info("Note, updating hub height but not updating " + "the specfied_wind_height")
1280 -
1281 -
    #     if changed_hub_height and update_specified_wind_height:
1282 -
    #         self.logger.info(f"Note, specfied_wind_height changed to hub-height: {new_hub_height}")
1283 -
    #         self.reinitialize(specified_wind_height=new_hub_height)
1284 -
1285 -
    #     # Finish by re-initalizing the flow field
1286 -
    #     self.reinitialize()
1287 -
1288 -
    # def set_use_points_on_perimeter(self, use_points_on_perimeter=False):
1289 -
    #     """
1290 -
    #     Set whether to use the points on the rotor diameter (perimeter) when
1291 -
    #     calculating flow field and wake.
1292 -
1293 -
    #     Args:
1294 -
    #         use_points_on_perimeter (bool): When *True*, use points at rotor
1295 -
    #             perimeter in wake and flow calculations. Defaults to *False*.
1296 -
    #     """
1297 -
    #     for turbine in self.floris.farm.turbines:
1298 -
    #         turbine.use_points_on_perimeter = use_points_on_perimeter
1299 -
    #         turbine.initialize_turbine()
1300 -
1301 -
    # def set_gch(self, enable=True):
1302 -
    #     """
1303 -
    #     Enable or disable Gauss-Curl Hybrid (GCH) functions
1304 -
    #     :py:meth:`~.GaussianModel.calculate_VW`,
1305 -
    #     :py:meth:`~.GaussianModel.yaw_added_recovery_correction`, and
1306 -
    #     :py:attr:`~.VelocityDeflection.use_secondary_steering`.
1307 -
1308 -
    #     Args:
1309 -
    #         enable (bool, optional): Flag whether or not to implement flow
1310 -
    #             corrections from GCH model. Defaults to *True*.
1311 -
    #     """
1312 -
    #     self.set_gch_yaw_added_recovery(enable)
1313 -
    #     self.set_gch_secondary_steering(enable)
1314 -
1315 -
    # def set_gch_yaw_added_recovery(self, enable=True):
1316 -
    #     """
1317 -
    #     Enable or Disable yaw-added recovery (YAR) from the Gauss-Curl Hybrid
1318 -
    #     (GCH) model and the control state of
1319 -
    #     :py:meth:`~.GaussianModel.calculate_VW_velocities` and
1320 -
    #     :py:meth:`~.GaussianModel.yaw_added_recovery_correction`.
1321 -
1322 -
    #     Args:
1323 -
    #         enable (bool, optional): Flag whether or not to implement yaw-added
1324 -
    #             recovery from GCH model. Defaults to *True*.
1325 -
    #     """
1326 -
    #     model_params = self.get_model_parameters()
1327 -
    #     use_secondary_steering = model_params["Wake Deflection Parameters"]["use_secondary_steering"]
1328 -
1329 -
    #     if enable:
1330 -
    #         model_params["Wake Velocity Parameters"]["use_yaw_added_recovery"] = True
1331 -
1332 -
    #         # If enabling be sure calc vw is on
1333 -
    #         model_params["Wake Velocity Parameters"]["calculate_VW_velocities"] = True
1334 -
1335 -
    #     if not enable:
1336 -
    #         model_params["Wake Velocity Parameters"]["use_yaw_added_recovery"] = False
1337 -
1338 -
    #         # If secondary steering is also off, disable calculate_VW_velocities
1339 -
    #         if not use_secondary_steering:
1340 -
    #             model_params["Wake Velocity Parameters"]["calculate_VW_velocities"] = False
1341 -
1342 -
    #     self.set_model_parameters(model_params)
1343 -
    #     self.reinitialize()
1344 -
1345 -
    # def set_gch_secondary_steering(self, enable=True):
1346 -
    #     """
1347 -
    #     Enable or Disable secondary steering (SS) from the Gauss-Curl Hybrid
1348 -
    #     (GCH) model and the control state of
1349 -
    #     :py:meth:`~.GaussianModel.calculate_VW_velocities` and
1350 -
    #     :py:attr:`~.VelocityDeflection.use_secondary_steering`.
1351 -
1352 -
    #     Args:
1353 -
    #         enable (bool, optional): Flag whether or not to implement secondary
1354 -
    #         steering from GCH model. Defaults to *True*.
1355 -
    #     """
1356 -
    #     model_params = self.get_model_parameters()
1357 -
    #     use_yaw_added_recovery = model_params["Wake Velocity Parameters"]["use_yaw_added_recovery"]
1358 -
1359 -
    #     if enable:
1360 -
    #         model_params["Wake Deflection Parameters"]["use_secondary_steering"] = True
1361 -
1362 -
    #         # If enabling be sure calc vw is on
1363 -
    #         model_params["Wake Velocity Parameters"]["calculate_VW_velocities"] = True
1364 -
1365 -
    #     if not enable:
1366 -
    #         model_params["Wake Deflection Parameters"]["use_secondary_steering"] = False
1367 -
1368 -
    #         # If yar is also off, disable calculate_VW_velocities
1369 -
    #         if not use_yaw_added_recovery:
1370 -
    #             model_params["Wake Velocity Parameters"]["calculate_VW_velocities"] = False
1371 -
1372 -
    #     self.set_model_parameters(model_params)
1373 -
    #     self.reinitialize()
1374 -
1375 -
    # def show_model_parameters(
1376 -
    #     self,
1377 -
    #     params=None,
1378 -
    #     verbose=False,
1379 -
    #     wake_velocity_model=True,
1380 -
    #     wake_deflection_model=True,
1381 -
    #     turbulence_model=False,
1382 -
    # ):
1383 -
    #     """
1384 -
    #     Helper function to print the current wake model parameters and values.
1385 -
    #     Shortcut to :py:meth:`~.tools.interface_utilities.show_params`.
1386 -
1387 -
    #     Args:
1388 -
    #         params (list, optional): Specific model parameters to be returned,
1389 -
    #             supplied as a list of strings. If None, then returns all
1390 -
    #             parameters. Defaults to None.
1391 -
    #         verbose (bool, optional): If set to *True*, will return the
1392 -
    #             docstrings for each parameter. Defaults to *False*.
1393 -
    #         wake_velocity_model (bool, optional): If set to *True*, will return
1394 -
    #             parameters from the wake_velocity model. If set to *False*, will
1395 -
    #             exclude parameters from the wake velocity model. Defaults to
1396 -
    #             *True*.
1397 -
    #         wake_deflection_model (bool, optional): If set to *True*, will
1398 -
    #             return parameters from the wake deflection model. If set to
1399 -
    #             *False*, will exclude parameters from the wake deflection
1400 -
    #             model. Defaults to *True*.
1401 -
    #         turbulence_model (bool, optional): If set to *True*, will return
1402 -
    #             parameters from the wake turbulence model. If set to *False*,
1403 -
    #             will exclude parameters from the wake turbulence model.
1404 -
    #             Defaults to *True*.
1405 -
    #     """
1406 -
    #     show_params(
1407 -
    #         self.floris.wake,
1408 -
    #         params,
1409 -
    #         verbose,
1410 -
    #         wake_velocity_model,
1411 -
    #         wake_deflection_model,
1412 -
    #         turbulence_model,
1413 -
    #     )
1414 -
1415 -
    # def get_model_parameters(
1416 -
    #     self,
1417 -
    #     params=None,
1418 -
    #     wake_velocity_model=True,
1419 -
    #     wake_deflection_model=True,
1420 -
    #     turbulence_model=True,
1421 -
    # ):
1422 -
    #     """
1423 -
    #     Helper function to return the current wake model parameters and values.
1424 -
    #     Shortcut to :py:meth:`~.tools.interface_utilities.get_params`.
1425 -
1426 -
    #     Args:
1427 -
    #         params (list, optional): Specific model parameters to be returned,
1428 -
    #             supplied as a list of strings. If None, then returns all
1429 -
    #             parameters. Defaults to None.
1430 -
    #         wake_velocity_model (bool, optional): If set to *True*, will return
1431 -
    #             parameters from the wake_velocity model. If set to *False*, will
1432 -
    #             exclude parameters from the wake velocity model. Defaults to
1433 -
    #             *True*.
1434 -
    #         wake_deflection_model (bool, optional): If set to *True*, will
1435 -
    #             return parameters from the wake deflection model. If set to
1436 -
    #             *False*, will exclude parameters from the wake deflection
1437 -
    #             model. Defaults to *True*.
1438 -
    #         turbulence_model ([type], optional): If set to *True*, will return
1439 -
    #             parameters from the wake turbulence model. If set to *False*,
1440 -
    #             will exclude parameters from the wake turbulence model.
1441 -
    #             Defaults to *True*.
1442 -
1443 -
    #     Returns:
1444 -
    #         dict: Dictionary containing model parameters and their values.
1445 -
    #     """
1446 -
    #     model_params = get_params(
1447 -
    #         self.floris.wake, params, wake_velocity_model, wake_deflection_model, turbulence_model
1448 -
    #     )
1449 -
1450 -
    #     return model_params
1451 -
1452 -
    # def set_model_parameters(self, params, verbose=True):
1453 -
    #     """
1454 -
    #     Helper function to set current wake model parameters.
1455 -
    #     Shortcut to :py:meth:`~.tools.interface_utilities.set_params`.
1456 -
1457 -
    #     Args:
1458 -
    #         params (dict): Specific model parameters to be set, supplied as a
1459 -
    #             dictionary of key:value pairs.
1460 -
    #         verbose (bool, optional): If set to *True*, will print information
1461 -
    #             about each model parameter that is changed. Defaults to *True*.
1462 -
    #     """
1463 -
    #     self.floris.wake = set_params(self.floris.wake, params, verbose)
1464 -
1465 -
1466 -
1467 -
1468 -
1469 -
1470 -
    # def vis_layout(
1471 -
    #     self,
1472 -
    #     ax=None,
1473 -
    #     show_wake_lines=False,
1474 -
    #     limit_dist=None,
1475 -
    #     turbine_face_north=False,
1476 -
    #     one_index_turbine=False,
1477 -
    #     black_and_white=False,
1478 -
    # ):
1479 -
    #     """
1480 -
    #     Visualize the layout of the wind farm in the floris instance.
1481 -
    #     Shortcut to :py:meth:`~.tools.layout_functions.visualize_layout`.
1482 -
1483 -
    #     Args:
1484 -
    #         ax (:py:class:`matplotlib.pyplot.axes`, optional):
1485 -
    #             Figure axes. Defaults to None.
1486 -
    #         show_wake_lines (bool, optional): Flag to control plotting of
1487 -
    #             wake boundaries. Defaults to False.
1488 -
    #         limit_dist (float, optional): Downstream limit to plot wakes.
1489 -
    #             Defaults to None.
1490 -
    #         turbine_face_north (bool, optional): Force orientation of wind
1491 -
    #             turbines. Defaults to False.
1492 -
    #         one_index_turbine (bool, optional): If *True*, 1st turbine is
1493 -
    #             turbine 1.
1494 -
    #     """
1495 -
    #     for i, turbine in enumerate(self.floris.farm.turbines):
1496 -
    #         D = turbine.rotor_diameter
1497 -
    #         break
1498 -
    #     layout_x, layout_y = self.get_turbine_layout()
1499 -
1500 -
    #     turbineLoc = build_turbine_loc(layout_x, layout_y)
1501 -
1502 -
    #     # Show visualize the turbine layout
1503 -
    #     visualize_layout(
1504 -
    #         turbineLoc,
1505 -
    #         D,
1506 -
    #         ax=ax,
1507 -
    #         show_wake_lines=show_wake_lines,
1508 -
    #         limit_dist=limit_dist,
1509 -
    #         turbine_face_north=turbine_face_north,
1510 -
    #         one_index_turbine=one_index_turbine,
1511 -
    #         black_and_white=black_and_white,
1512 -
    #     )
1513 -
1514 -
    # def show_flow_field(self, ax=None):
1515 -
    #     """
1516 -
    #     Shortcut method to
1517 -
    #     :py:meth:`~.tools.visualization.visualize_cut_plane`.
1518 -
1519 -
    #     Args:
1520 -
    #         ax (:py:class:`matplotlib.pyplot.axes` optional):
1521 -
    #             Figure axes. Defaults to None.
1522 -
    #     """
1523 -
    #     # Get horizontal plane at default height (hub-height)
1524 -
    #     hor_plane = self.get_hor_plane()
1525 -
1526 -
    #     # Plot and show
1527 -
    #     if ax is None:
1528 -
    #         fig, ax = plt.subplots()
1529 -
    #     visualize_cut_plane(hor_plane, ax=ax)
1530 -
    #     plt.show()
1531 -
1532 -
1533 -
1534 909
    ## Functionality removed in v3
1535 910
1536 911
    def set_rotor_diameter(self, rotor_diameter):

@@ -79,6 +79,7 @@
Loading
79 79
80 80
def power(
81 81
    air_density: float,
82 +
    ref_density_cp_ct: float,
82 83
    velocities: NDArrayFloat,
83 84
    yaw_angle: NDArrayFloat,
84 85
    pP: float,
@@ -91,6 +92,7 @@
Loading
91 92
92 93
    Args:
93 94
        air_density (NDArrayFloat[wd, ws, turbines]): The air density value(s) at each turbine.
95 +
        ref_density_cp_cts (NDArrayFloat[wd, ws, turbines]): The reference density for each turbine
94 96
        velocities (NDArrayFloat[wd, ws, turbines, grid1, grid2]): The velocity field at a turbine.
95 97
        pP (NDArrayFloat[wd, ws, turbines]): The pP value(s) of the cosine exponent relating
96 98
            the yaw misalignment angle to power for each turbine.
@@ -134,7 +136,7 @@
Loading
134 136
135 137
    # Compute the yaw effective velocity
136 138
    pW = pP / 3.0  # Convert from pP to w
137 -
    yaw_effective_velocity = ((air_density/1.225)**(1/3)) * average_velocity(velocities) * cosd(yaw_angle) ** pW
139 +
    yaw_effective_velocity = ((air_density/ref_density_cp_ct)**(1/3)) * average_velocity(velocities) * cosd(yaw_angle) ** pW
138 140
139 141
    # Loop over each turbine type given to get thrust coefficient for all turbines
140 142
    p = np.zeros(np.shape(yaw_effective_velocity))
@@ -145,7 +147,7 @@
Loading
145 147
        # type to the main thrust coefficient array
146 148
        p += power_interp[turb_type](yaw_effective_velocity) * np.array(turbine_type_map == turb_type)
147 149
148 -
    return p * 1.225
150 +
    return p * ref_density_cp_ct
149 151
150 152
151 153
def Ct(
@@ -317,6 +319,8 @@
Loading
317 319
            tilt angle to power.
318 320
        generator_efficiency (:py:obj: float): The generator
319 321
            efficiency factor used to scale the power production.
322 +
        ref_density_cp_ct (:py:obj: float): The density at which the provided
323 +
            cp and ct is defined
320 324
        power_thrust_table (PowerThrustTable): A dictionary containing the
321 325
            following key-value pairs:
322 326
@@ -343,8 +347,11 @@
Loading
343 347
    pT: float
344 348
    TSR: float
345 349
    generator_efficiency: float
350 +
    ref_density_cp_ct: float
346 351
    power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict)
347 352
353 +
354 +
348 355
    # rloc: float = float_attrib()  # TODO: goes here or on the Grid?
349 356
    # use_points_on_perimeter: bool = bool_attrib()
350 357
@@ -355,6 +362,7 @@
Loading
355 362
    fCt_interp: interp1d = field(init=False)
356 363
    power_interp: interp1d = field(init=False)
357 364
365 +
358 366
    # For the following parameters, use default values if not user-specified
359 367
    # self.rloc = float(input_dictionary["rloc"]) if "rloc" in input_dictionary else 0.5
360 368
    # if "use_points_on_perimeter" in input_dictionary:

@@ -18,20 +18,29 @@
Loading
18 18
"""
19 19
20 20
from abc import ABC, abstractmethod
21 +
from enum import Enum
21 22
from typing import Any, Dict, Final
22 23
23 24
import attrs
24 -
from attrs import define
25 25
26 26
from floris.type_dec import FromDictMixin
27 27
from floris.logging_manager import LoggerBase
28 28
29 29
30 +
class State(Enum):
31 +
    UNINITIALIZED = 0
32 +
    INITIALIZED = 1
33 +
    USED = 2
34 +
35 +
30 36
class BaseClass(LoggerBase, FromDictMixin):
31 37
    """
32 38
    BaseClass object class. This class does the logging and MixIn class inheritance.
33 39
    """
34 40
41 +
    state = State.UNINITIALIZED
42 +
43 +
35 44
    @classmethod
36 45
    def get_model_defaults(cls) -> Dict[str, Any]:
37 46
        """Produces a dictionary of the keyword arguments and their defaults.
@@ -63,15 +72,6 @@
Loading
63 72
64 73
    NUM_EPS: Final[float] = 0.001  # This is a numerical epsilon to prevent divide by zeros
65 74
66 -
    @property
67 -
    def model_string(self):
68 -
        return self.model_string
69 -
70 -
    @model_string.setter
71 -
    @abstractmethod
72 -
    def model_string(self, string):
73 -
        raise NotImplementedError("BaseModel.model_string")
74 -
75 75
    @abstractmethod
76 76
    def prepare_function() -> dict:
77 77
        raise NotImplementedError("BaseModel.prepare_function")

@@ -1, +1, @@
Loading
1 -
from . import other, scipy, pyoptsparse, yaw_optimization
1 +
from . import other, legacy, yaw_optimization, layout_optimization

@@ -0,0 +1,183 @@
Loading
1 +
# Copyright 2022 NREL
2 +
3 +
# Licensed under the Apache License, Version 2.0 (the "License"); you may not
4 +
# use this file except in compliance with the License. You may obtain a copy of
5 +
# the License at http://www.apache.org/licenses/LICENSE-2.0
6 +
7 +
# Unless required by applicable law or agreed to in writing, software
8 +
# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
9 +
# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
10 +
# License for the specific language governing permissions and limitations under
11 +
# the License.
12 +
13 +
# See https://floris.readthedocs.io for documentation
14 +
15 +
16 +
import numpy as np
17 +
import matplotlib.pyplot as plt
18 +
from shapely.geometry import Point
19 +
from scipy.spatial.distance import cdist
20 +
21 +
from .layout_optimization_base import LayoutOptimization
22 +
23 +
class LayoutOptimizationPyOptSparse(LayoutOptimization):
24 +
    def __init__(
25 +
        self,
26 +
        fi,
27 +
        boundaries,
28 +
        min_dist=None,
29 +
        freq=None,
30 +
        solver=None,
31 +
        optOptions=None,
32 +
        timeLimit=None,
33 +
        storeHistory='hist.hist',
34 +
        hotStart=None
35 +
    ):
36 +
        super().__init__(fi, boundaries, min_dist=min_dist, freq=freq)
37 +
38 +
        self.x0 = self._norm(self.fi.layout_x, self.xmin, self.xmax)
39 +
        self.y0 = self._norm(self.fi.layout_y, self.ymin, self.ymax)
40 +
41 +
        self.storeHistory = storeHistory
42 +
        self.timeLimit = timeLimit
43 +
        self.hotStart = hotStart
44 +
45 +
        try:
46 +
            import pyoptsparse
47 +
        except ImportError:
48 +
            err_msg = (
49 +
                "It appears you do not have pyOptSparse installed. "
50 +
                + "Please refer to https://pyoptsparse.readthedocs.io/ for "
51 +
                + "guidance on how to properly install the module."
52 +
            )
53 +
            self.logger.error(err_msg, stack_info=True)
54 +
            raise ImportError(err_msg)
55 +
56 +
        # Insantiate ptOptSparse optimization object with name and objective function
57 +
        self.optProb = pyoptsparse.Optimization('layout', self._obj_func)
58 +
59 +
        self.optProb = self.add_var_group(self.optProb)
60 +
        self.optProb = self.add_con_group(self.optProb)
61 +
        self.optProb.addObj("obj")
62 +
63 +
        if solver is not None:
64 +
            self.solver = solver
65 +
            print("Setting up optimization with user's choice of solver: ", self.solver)
66 +
        else:
67 +
            self.solver = "SLSQP"
68 +
            print("Setting up optimization with default solver: SLSQP.")
69 +
        if optOptions is not None:
70 +
            self.optOptions = optOptions
71 +
        else:
72 +
            if self.solver == "SNOPT":
73 +
                self.optOptions = {"Major optimality tolerance": 1e-7}
74 +
            else:
75 +
                self.optOptions = {}
76 +
77 +
        exec("self.opt = pyoptsparse." + self.solver + "(options=self.optOptions)")
78 +
79 +
    def _optimize(self):
80 +
        if hasattr(self, "_sens"):
81 +
            self.sol = self.opt(self.optProb, sens=self._sens)
82 +
        else:
83 +
            if self.timeLimit is not None:
84 +
                self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, timeLimit=self.timeLimit, hotStart=self.hotStart)
85 +
            else:
86 +
                self.sol = self.opt(self.optProb, sens="CDR", storeHistory=self.storeHistory, hotStart=self.hotStart)
87 +
        return self.sol
88 +
89 +
    def _obj_func(self, varDict):
90 +
        # Parse the variable dictionary
91 +
        self.parse_opt_vars(varDict)
92 +
93 +
        # Update turbine map with turbince locations
94 +
        self.fi.reinitialize(layout_x = self.x, layout_y = self.y)
95 +
96 +
        # Compute the objective function
97 +
        funcs = {}
98 +
        funcs["obj"] = (
99 +
            -1 * self.fi.get_farm_AEP(self.freq) / self.initial_AEP
100 +
        )
101 +
102 +
        # Compute constraints, if any are defined for the optimization
103 +
        funcs = self.compute_cons(funcs, self.x, self.y)
104 +
105 +
        fail = False
106 +
        return funcs, fail
107 +
108 +
    # Optionally, the user can supply the optimization with gradients
109 +
    # def _sens(self, varDict, funcs):
110 +
    #     funcsSens = {}
111 +
    #     fail = False
112 +
    #     return funcsSens, fail
113 +
114 +
    def parse_opt_vars(self, varDict):
115 +
        self.x = self._unnorm(varDict["x"], self.xmin, self.xmax)
116 +
        self.y = self._unnorm(varDict["y"], self.ymin, self.ymax)
117 +
118 +
    def parse_sol_vars(self, sol):
119 +
        self.x = list(self._unnorm(sol.getDVs()["x"], self.xmin, self.xmax))[0]
120 +
        self.y = list(self._unnorm(sol.getDVs()["y"], self.ymin, self.ymax))[1]
121 +
122 +
    def add_var_group(self, optProb):
123 +
        optProb.addVarGroup(
124 +
            "x", self.nturbs, varType="c", lower=0.0, upper=1.0, value=self.x0
125 +
        )
126 +
        optProb.addVarGroup(
127 +
            "y", self.nturbs, varType="c", lower=0.0, upper=1.0, value=self.y0
128 +
        )
129 +
130 +
        return optProb
131 +
132 +
    def add_con_group(self, optProb):
133 +
        optProb.addConGroup("boundary_con", self.nturbs, upper=0.0)
134 +
        optProb.addConGroup("spacing_con", 1, upper=0.0)
135 +
136 +
        return optProb
137 +
138 +
    def compute_cons(self, funcs, x, y):
139 +
        funcs["boundary_con"] = self.distance_from_boundaries(x, y)
140 +
        funcs["spacing_con"] = self.space_constraint(x, y)
141 +
142 +
        return funcs
143 +
144 +
    def space_constraint(self, x, y, rho=500):
145 +
        # Calculate distances between turbines
146 +
        locs = np.vstack((x, y)).T
147 +
        distances = cdist(locs, locs)
148 +
        arange = np.arange(distances.shape[0])
149 +
        distances[arange, arange] = 1e10
150 +
        dist = np.min(distances, axis=0)
151 +
152 +
        g = 1 - np.array(dist) / self.min_dist
153 +
154 +
        # Following code copied from OpenMDAO KSComp().
155 +
        # Constraint is satisfied when KS_constraint <= 0
156 +
        g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis]
157 +
        g_diff = g - g_max
158 +
        exponents = np.exp(rho * g_diff)
159 +
        summation = np.sum(exponents, axis=-1)[:, np.newaxis]
160 +
        KS_constraint = g_max + 1.0 / rho * np.log(summation)
161 +
162 +
        return KS_constraint[0][0]
163 +
164 +
    def distance_from_boundaries(self, x, y):
165 +
        boundary_con = np.zeros(self.nturbs)
166 +
        for i in range(self.nturbs):
167 +
            loc = Point(x[i], y[i])
168 +
            boundary_con[i] = loc.distance(self._boundary_line)
169 +
            if self._boundary_polygon.contains(loc)==True:
170 +
                boundary_con[i] *= -1.0
171 +
172 +
        return boundary_con
173 +
174 +
    def _get_initial_and_final_locs(self):
175 +
        x_initial = self._unnorm(self.x0, self.xmin, self.xmax)
176 +
        y_initial = self._unnorm(self.y0, self.ymin, self.ymax)
177 +
        x_opt, y_opt = self.get_optimized_locs()
178 +
        return x_initial, y_initial, x_opt, y_opt
179 +
180 +
    def get_optimized_locs(self):
181 +
        x_opt = self._unnorm(self.sol.getDVs()["x"], self.xmin, self.xmax)
182 +
        y_opt = self._unnorm(self.sol.getDVs()["y"], self.ymin, self.ymax)
183 +
        return x_opt, y_opt

@@ -1, +1, @@
Loading
1 -
3.1.1
1 +
3.2

@@ -13,16 +13,16 @@
Loading
13 13
14 14
15 15
import copy
16 +
16 17
import numpy as np
17 18
from scipy.stats import norm
18 19
19 20
from floris.tools import FlorisInterface
20 -
from floris.logging_manager import LoggerBase
21 21
from floris.utilities import wrap_360
22 +
from floris.logging_manager import LoggerBase
22 23
23 24
24 25
class UncertaintyInterface(LoggerBase):
25 -
26 26
    def __init__(
27 27
        self,
28 28
        configuration,
@@ -77,7 +77,7 @@
Loading
77 77
            will essentially come down to a Gaussian smoothing of FLORIS
78 78
            solutions over the wind directions. This calculation can therefore
79 79
            be really fast, since it does not require additional calculations
80 -
            compared to a non-uncertainty FLORIS evaluation. 
80 +
            compared to a non-uncertainty FLORIS evaluation.
81 81
            When fix_yaw_in_relative_frame=False, the yaw angles are fixed in
82 82
            the absolute (compass) reference frame, meaning that for each
83 83
            probablistic wind direction evaluation, our probablistic (relative)
@@ -118,6 +118,9 @@
Loading
118 118
            fix_yaw_in_relative_frame=fix_yaw_in_relative_frame,
119 119
        )
120 120
121 +
        # Add a _no_wake switch to keep track of calculate_wake/calculate_no_wake
122 +
        self._no_wake = False
123 +
121 124
    # Private methods
122 125
123 126
    def _generate_pdfs_from_dict(self):
@@ -132,7 +135,9 @@
Loading
132 135
        # create normally distributed wd and yaw uncertaitny pmfs if appropriate
133 136
        unc_options = self.unc_options
134 137
        if unc_options["std_wd"] > 0:
135 -
            wd_bnd = int(np.ceil(norm.ppf(unc_options["pdf_cutoff"], scale=unc_options["std_wd"]) / unc_options["pmf_res"]))
138 +
            wd_bnd = int(
139 +
                np.ceil(norm.ppf(unc_options["pdf_cutoff"], scale=unc_options["std_wd"]) / unc_options["pmf_res"])
140 +
            )
136 141
            bound = wd_bnd * unc_options["pmf_res"]
137 142
            wd_unc = np.linspace(-1 * bound, bound, 2 * wd_bnd + 1)
138 143
            wd_unc_pmf = norm.pdf(wd_unc, scale=unc_options["std_wd"])
@@ -176,10 +181,7 @@
Loading
176 181
177 182
        # Expand wind direction and yaw angle array into the direction
178 183
        # of uncertainty over the ambient wind direction.
179 -
        wd_array_probablistic = np.vstack(
180 -
            [np.expand_dims(wd_array_nominal, axis=0) + dy
181 -
            for dy in unc_pmfs["wd_unc"]]
182 -
        )
184 +
        wd_array_probablistic = np.vstack([np.expand_dims(wd_array_nominal, axis=0) + dy for dy in unc_pmfs["wd_unc"]])
183 185
184 186
        if self.fix_yaw_in_relative_frame:
185 187
            # The relative yaw angle is fixed and always has the nominal
@@ -190,8 +192,7 @@
Loading
190 192
            # not require any additional calculations compared to the
191 193
            # non-uncertainty FLORIS evaluation.
192 194
            yaw_angles_probablistic = np.vstack(
193 -
                [np.expand_dims(yaw_angles_nominal, axis=0)
194 -
                for _ in unc_pmfs["wd_unc"]]
195 +
                [np.expand_dims(yaw_angles_nominal, axis=0) for _ in unc_pmfs["wd_unc"]]
195 196
            )
196 197
        else:
197 198
            # Fix yaw angles in the absolute (compass) reference frame,
@@ -202,8 +203,7 @@
Loading
202 203
            # it with a relative yaw angle that is 3 deg below its nominal
203 204
            # value.
204 205
            yaw_angles_probablistic = np.vstack(
205 -
                [np.expand_dims(yaw_angles_nominal, axis=0) - dy
206 -
                for dy in unc_pmfs["wd_unc"]]
206 +
                [np.expand_dims(yaw_angles_nominal, axis=0) - dy for dy in unc_pmfs["wd_unc"]]
207 207
            )
208 208
209 209
        self.wd_array_probablistic = wd_array_probablistic
@@ -223,12 +223,7 @@
Loading
223 223
        fi_unc_copy.fi = self.fi.copy()
224 224
        return fi_unc_copy
225 225
226 -
    def reinitialize_uncertainty(
227 -
        self,
228 -
        unc_options=None,
229 -
        unc_pmfs=None,
230 -
        fix_yaw_in_relative_frame=None
231 -
    ):
226 +
    def reinitialize_uncertainty(self, unc_options=None, unc_pmfs=None, fix_yaw_in_relative_frame=None):
232 227
        """Reinitialize the wind direction and yaw angle probability
233 228
        distributions used in evaluating FLORIS. Must either specify
234 229
        'unc_options', in which case distributions are calculated assuming
@@ -281,7 +276,7 @@
Loading
281 276
                will essentially come down to a Gaussian smoothing of FLORIS
282 277
                solutions over the wind directions. This calculation can therefore
283 278
                be really fast, since it does not require additional calculations
284 -
                compared to a non-uncertainty FLORIS evaluation. 
279 +
                compared to a non-uncertainty FLORIS evaluation.
285 280
                When fix_yaw_in_relative_frame=False, the yaw angles are fixed in
286 281
                the absolute (compass) reference frame, meaning that for each
287 282
                probablistic wind direction evaluation, our probablistic (relative)
@@ -300,23 +295,21 @@
Loading
300 295
                often does not perfectly know the true wind direction, and that a
301 296
                turbine often does not perfectly achieve its desired yaw angle offset.
302 297
                Defaults to fix_yaw_in_relative_frame=False.
303 -
                
298 +
304 299
        """
305 300
306 301
        # Check inputs
307 -
        if ((unc_options is not None) and (unc_pmfs is not None)):
308 -
            self.logger.error(
309 -
                "Must specify either 'unc_options' or 'unc_pmfs', not both."
310 -
            )
302 +
        if (unc_options is not None) and (unc_pmfs is not None):
303 +
            self.logger.error("Must specify either 'unc_options' or 'unc_pmfs', not both.")
311 304
312 305
        # Assign uncertainty probability distributions
313 306
        if unc_options is not None:
314 307
            self.unc_options = unc_options
315 308
            self._generate_pdfs_from_dict()
316 -
        
309 +
317 310
        if unc_pmfs is not None:
318 311
            self.unc_pmfs = unc_pmfs
319 -
        
312 +
320 313
        if fix_yaw_in_relative_frame is not None:
321 314
            self.fix_yaw_in_relative_frame = bool(fix_yaw_in_relative_frame)
322 315
@@ -330,6 +323,8 @@
Loading
330 323
        turbulence_intensity=None,
331 324
        air_density=None,
332 325
        layout=None,
326 +
        layout_x=None,
327 +
        layout_y=None,
333 328
        turbine_type=None,
334 329
        solver_settings=None,
335 330
    ):
@@ -337,6 +332,12 @@
Loading
337 332
        to directly replace a FlorisInterface object with this
338 333
        UncertaintyInterface object, this function is required."""
339 334
335 +
        if layout is not None:
336 +
            msg = "Use the `layout_x` and `layout_y` parameters in place of `layout` because the `layout` parameter will be deprecated in 3.3."
337 +
            self.logger.warning(msg)
338 +
            layout_x = layout[0]
339 +
            layout_y = layout[1]
340 +
340 341
        # Just passes arguments to the floris object
341 342
        self.fi.reinitialize(
342 343
            wind_speeds=wind_speeds,
@@ -346,7 +347,8 @@
Loading
346 347
            reference_wind_height=reference_wind_height,
347 348
            turbulence_intensity=turbulence_intensity,
348 349
            air_density=air_density,
349 -
            layout=layout,
350 +
            layout_x=layout_x,
351 +
            layout_y=layout_y,
350 352
            turbine_type=turbine_type,
351 353
            solver_settings=solver_settings,
352 354
        )
@@ -364,16 +366,26 @@
Loading
364 366
            yaw_angles: NDArrayFloat | list[float] | None = None,
365 367
        """
366 368
        self._reassign_yaw_angles(yaw_angles)
369 +
        self._no_wake = False
367 370
368 -
    def get_turbine_powers(self, no_wake=False):
369 -
        """Calculates the probability-weighted power production of each
370 -
        turbine in the wind farm.
371 +
    def calculate_no_wake(self, yaw_angles=None):
372 +
        """Replaces the 'calculate_no_wake' function in the FlorisInterface
373 +
        object. Fundamentally, this function only overwrites the nominal
374 +
        yaw angles in the FlorisInterface object. The actual wake calculations
375 +
        are performed once 'get_turbine_powers' or 'get_farm_powers' is
376 +
        called. However, to allow users to directly replace a FlorisInterface
377 +
        object with this UncertaintyInterface object, this function is
378 +
        required.
371 379
372 380
        Args:
373 -
            no_wake (bool, optional): disable the wakes in the flow model.
374 -
            This can be useful to determine the (probablistic) power
375 -
            production of the farm in the artificial scenario where there
376 -
            would never be any wake losses. Defaults to False.
381 +
            yaw_angles: NDArrayFloat | list[float] | None = None,
382 +
        """
383 +
        self._reassign_yaw_angles(yaw_angles)
384 +
        self._no_wake = True
385 +
386 +
    def get_turbine_powers(self):
387 +
        """Calculates the probability-weighted power production of each
388 +
        turbine in the wind farm.
377 389
378 390
        Returns:
379 391
            NDArrayFloat: Power production of all turbines in the wind farm.
@@ -399,9 +411,7 @@
Loading
399 411
400 412
        # Format into conventional floris format by reshaping
401 413
        wd_array_probablistic = np.reshape(self.wd_array_probablistic, -1)
402 -
        yaw_angles_probablistic = np.reshape(
403 -
            self.yaw_angles_probablistic, (-1, num_ws, num_turbines)
404 -
        )
414 +
        yaw_angles_probablistic = np.reshape(self.yaw_angles_probablistic, (-1, num_ws, num_turbines))
405 415
406 416
        # Wrap wind direction array around 360 deg
407 417
        wd_array_probablistic = wrap_360(wd_array_probablistic)
@@ -409,17 +419,14 @@
Loading
409 419
        # Find minimal set of solutions to evaluate
410 420
        wd_exp = np.tile(wd_array_probablistic, (1, num_ws, 1)).T
411 421
        _, id_unq, id_unq_rev = np.unique(
412 -
            np.append(yaw_angles_probablistic, wd_exp, axis=2),
413 -
            axis=0,
414 -
            return_index=True,
415 -
            return_inverse=True
422 +
            np.append(yaw_angles_probablistic, wd_exp, axis=2), axis=0, return_index=True, return_inverse=True
416 423
        )
417 424
        wd_array_probablistic_min = wd_array_probablistic[id_unq]
418 425
        yaw_angles_probablistic_min = yaw_angles_probablistic[id_unq, :, :]
419 426
420 427
        # Evaluate floris for minimal probablistic set
421 428
        self.fi.reinitialize(wind_directions=wd_array_probablistic_min)
422 -
        if no_wake:
429 +
        if self._no_wake:
423 430
            self.fi.calculate_no_wake(yaw_angles=yaw_angles_probablistic_min)
424 431
        else:
425 432
            self.fi.calculate_wake(yaw_angles=yaw_angles_probablistic_min)
@@ -430,37 +437,188 @@
Loading
430 437
431 438
        # Reshape solutions back to full set
432 439
        power_probablistic = turbine_powers[id_unq_rev, :]
433 -
        power_probablistic = np.reshape(
434 -
            power_probablistic, 
435 -
            (num_wd_unc, num_wd, num_ws, num_turbines)
436 -
        )
440 +
        power_probablistic = np.reshape(power_probablistic, (num_wd_unc, num_wd, num_ws, num_turbines))
437 441
438 442
        # Calculate probability weighing terms
439 443
        wd_weighing = (
440 -
            np.expand_dims(unc_pmfs["wd_unc_pmf"], axis=(1, 2, 3))
441 -
        ).repeat(num_wd, 1).repeat(num_ws, 2).repeat(num_turbines, 3)
444 +
            (np.expand_dims(unc_pmfs["wd_unc_pmf"], axis=(1, 2, 3)))
445 +
            .repeat(num_wd, 1)
446 +
            .repeat(num_ws, 2)
447 +
            .repeat(num_turbines, 3)
448 +
        )
442 449
443 450
        # Now apply probability distribution weighing to get turbine powers
444 451
        return np.sum(wd_weighing * power_probablistic, axis=0)
445 452
446 -
    def get_farm_power(self, no_wake=False):
453 +
    def get_farm_power(self, turbine_weights=None):
447 454
        """Calculates the probability-weighted power production of the
448 455
        collective of all turbines in the farm, for each wind direction
449 456
        and wind speed specified.
450 457
451 458
        Args:
452 -
            no_wake (bool, optional): disable the wakes in the flow model.
453 -
            This can be useful to determine the (probablistic) power
454 -
            production of the farm in the artificial scenario where there
455 -
            would never be any wake losses. Defaults to False.
459 +
            turbine_weights (NDArrayFloat | list[float] | None, optional):
460 +
                weighing terms that allow the user to emphasize power at
461 +
                particular turbines and/or completely ignore the power 
462 +
                from other turbines. This is useful when, for example, you are
463 +
                modeling multiple wind farms in a single floris object. If you
464 +
                only want to calculate the power production for one of those
465 +
                farms and include the wake effects of the neighboring farms,
466 +
                you can set the turbine_weights for the neighboring farms'
467 +
                turbines to 0.0. The array of turbine powers from floris
468 +
                is multiplied with this array in the calculation of the
469 +
                objective function. If None, this  is an array with all values
470 +
                1.0 and with shape equal to (n_wind_directions, n_wind_speeds,
471 +
                n_turbines). Defaults to None.
456 472
457 473
        Returns:
458 474
            NDArrayFloat: Expectation of power production of the wind farm.
459 475
            This array has the shape (num_wind_directions, num_wind_speeds).
460 476
        """
461 -
        turbine_powers = self.get_turbine_powers(no_wake=no_wake)
477 +
478 +
        if turbine_weights is None:
479 +
            # Default to equal weighing of all turbines when turbine_weights is None
480 +
            turbine_weights = np.ones(
481 +
                (
482 +
                    self.floris.flow_field.n_wind_directions,
483 +
                    self.floris.flow_field.n_wind_speeds,
484 +
                    self.floris.farm.n_turbines
485 +
                )
486 +
            )
487 +
        elif len(np.shape(turbine_weights)) == 1:
488 +
            # Deal with situation when 1D array is provided
489 +
            turbine_weights = np.tile(
490 +
                turbine_weights,
491 +
                (
492 +
                    self.floris.flow_field.n_wind_directions,
493 +
                    self.floris.flow_field.n_wind_speeds,
494 +
                    1
495 +
                )
496 +
            )
497 +
498 +
        # Calculate all turbine powers and apply weights
499 +
        turbine_powers = self.get_turbine_powers()
500 +
        turbine_powers = np.multiply(turbine_weights, turbine_powers)
501 +
462 502
        return np.sum(turbine_powers, axis=2)
463 503
504 +
    def get_farm_AEP(
505 +
        self,
506 +
        freq,
507 +
        cut_in_wind_speed=0.001,
508 +
        cut_out_wind_speed=None,
509 +
        yaw_angles=None,
510 +
        turbine_weights=None,
511 +
        no_wake=False,
512 +
    ) -> float:
513 +
        """
514 +
        Estimate annual energy production (AEP) for distributions of wind speed, wind
515 +
        direction, frequency of occurrence, and yaw offset.
516 +
517 +
        Args:
518 +
            freq (NDArrayFloat): NumPy array with shape (n_wind_directions,
519 +
                n_wind_speeds) with the frequencies of each wind direction and
520 +
                wind speed combination. These frequencies should typically sum
521 +
                up to 1.0 and are used to weigh the wind farm power for every
522 +
                condition in calculating the wind farm's AEP.
523 +
            cut_in_wind_speed (float, optional): Wind speed in m/s below which
524 +
                any calculations are ignored and the wind farm is known to
525 +
                produce 0.0 W of power. Note that to prevent problems with the
526 +
                wake models at negative / zero wind speeds, this variable must
527 +
                always have a positive value. Defaults to 0.001 [m/s].
528 +
            cut_out_wind_speed (float, optional): Wind speed above which the
529 +
                wind farm is known to produce 0.0 W of power. If None is
530 +
                specified, will assume that the wind farm does not cut out
531 +
                at high wind speeds. Defaults to None.
532 +
            yaw_angles (NDArrayFloat | list[float] | None, optional):
533 +
                The relative turbine yaw angles in degrees. If None is
534 +
                specified, will assume that the turbine yaw angles are all
535 +
                zero degrees for all conditions. Defaults to None.
536 +
            turbine_weights (NDArrayFloat | list[float] | None, optional):
537 +
                weighing terms that allow the user to emphasize power at
538 +
                particular turbines and/or completely ignore the power 
539 +
                from other turbines. This is useful when, for example, you are
540 +
                modeling multiple wind farms in a single floris object. If you
541 +
                only want to calculate the power production for one of those
542 +
                farms and include the wake effects of the neighboring farms,
543 +
                you can set the turbine_weights for the neighboring farms'
544 +
                turbines to 0.0. The array of turbine powers from floris
545 +
                is multiplied with this array in the calculation of the
546 +
                objective function. If None, this  is an array with all values
547 +
                1.0 and with shape equal to (n_wind_directions, n_wind_speeds,
548 +
                n_turbines). Defaults to None.
549 +
            no_wake: (bool, optional): When *True* updates the turbine
550 +
                quantities without calculating the wake or adding the wake to
551 +
                the flow field. This can be useful when quantifying the loss
552 +
                in AEP due to wakes. Defaults to *False*.
553 +
554 +
        Returns:
555 +
            float:
556 +
                The Annual Energy Production (AEP) for the wind farm in
557 +
                watt-hours.
558 +
        """
559 +
560 +
        # Verify dimensions of the variable "freq"
561 +
        if not (
562 +
            (np.shape(freq)[0] == self.floris.flow_field.n_wind_directions)
563 +
            & (np.shape(freq)[1] == self.floris.flow_field.n_wind_speeds)
564 +
            & (len(np.shape(freq)) == 2)
565 +
        ):
566 +
            raise UserWarning(
567 +
                "'freq' should be a two-dimensional array with dimensions (n_wind_directions, n_wind_speeds)."
568 +
            )
569 +
570 +
        # Check if frequency vector sums to 1.0. If not, raise a warning
571 +
        if np.abs(np.sum(freq) - 1.0) > 0.001:
572 +
            self.logger.warning("WARNING: The frequency array provided to get_farm_AEP() does not sum to 1.0. ")
573 +
574 +
        # Copy the full wind speed array from the floris object and initialize
575 +
        # the the farm_power variable as an empty array.
576 +
        wind_speeds = np.array(self.fi.floris.flow_field.wind_speeds, copy=True)
577 +
        farm_power = np.zeros((self.fi.floris.flow_field.n_wind_directions, len(wind_speeds)))
578 +
579 +
        # Determine which wind speeds we must evaluate in floris
580 +
        conditions_to_evaluate = wind_speeds >= cut_in_wind_speed
581 +
        if cut_out_wind_speed is not None:
582 +
            conditions_to_evaluate = conditions_to_evaluate & (wind_speeds < cut_out_wind_speed)
583 +
584 +
        # Evaluate the conditions in floris
585 +
        if np.any(conditions_to_evaluate):
586 +
            wind_speeds_subset = wind_speeds[conditions_to_evaluate]
587 +
            yaw_angles_subset = None
588 +
            if yaw_angles is not None:
589 +
                yaw_angles_subset = yaw_angles[:, conditions_to_evaluate]
590 +
            self.reinitialize(wind_speeds=wind_speeds_subset)
591 +
            if no_wake:
592 +
                self.calculate_no_wake(yaw_angles=yaw_angles_subset)
593 +
            else:
594 +
                self.calculate_wake(yaw_angles=yaw_angles_subset)
595 +
            farm_power[:, conditions_to_evaluate] = (
596 +
                self.get_farm_power(turbine_weights=turbine_weights)
597 +
            )
598 +
599 +
        # Finally, calculate AEP in GWh
600 +
        aep = np.sum(np.multiply(freq, farm_power) * 365 * 24)
601 +
602 +
        # Reset the FLORIS object to the full wind speed array
603 +
        self.reinitialize(wind_speeds=wind_speeds)
604 +
605 +
        return aep
606 +
607 +
    def assign_hub_height_to_ref_height(self):
608 +
        return self.fi.assign_hub_height_to_ref_height()
609 +
610 +
    def get_turbine_layout(self, z=False):
611 +
        return self.fi.get_turbine_layout(z=z)
612 +
613 +
    def get_turbine_Cts(self):
614 +
        return self.fi.get_turbine_Cts()
615 +
616 +
    def get_turbine_ais(self):
617 +
        return self.fi.get_turbine_ais()
618 +
619 +
    def get_turbine_average_velocities(self):
620 +
        return self.fi.get_turbine_average_velocities()
621 +
464 622
    # Define getter functions that just pass information from FlorisInterface
465 623
    @property
466 624
    def floris(self):

@@ -188,6 +188,7 @@
Loading
188 188
        "rotor_diameter": tp["rotor_diameter"],
189 189
        "TSR": tp["TSR"],
190 190
        "power_thrust_table": tp["power_thrust_table"],
191 +
        "ref_density_cp_ct": 1.225 # This was implicit in the former input file
191 192
    }
192 193
193 194
    return dict_floris, dict_turbine

@@ -61,7 +61,7 @@
Loading
61 61
        self.parse_opt_vars(varDict)
62 62
63 63
        # Update turbine map with turbince locations
64 -
        self.fi.reinitialize(layout=[self.x, self.y])
64 +
        self.fi.reinitialize(layout_x=self.x, layout_y=self.y)
65 65
        self.fi.calculate_wake()
66 66
67 67
        # Compute the objective function
68 68
imilarity index 100%
69 69
ename from floris/tools/optimization/pyoptsparse/optimization.py
70 70
ename to floris/tools/optimization/legacy/pyoptsparse/optimization.py
71 71
imilarity index 100%
72 72
ename from floris/tools/optimization/pyoptsparse/power_density.py
73 73
ename to floris/tools/optimization/legacy/pyoptsparse/power_density.py
74 74
imilarity index 100%
75 75
ename from floris/tools/optimization/pyoptsparse/yaw.py
76 76
ename to floris/tools/optimization/legacy/pyoptsparse/yaw.py
77 77
imilarity index 100%
78 78
ename from floris/tools/optimization/scipy/__init__.py
79 79
ename to floris/tools/optimization/legacy/scipy/__init__.py
80 80
imilarity index 100%
81 81
ename from floris/tools/optimization/scipy/base_COE.py
82 82
ename to floris/tools/optimization/legacy/scipy/base_COE.py
83 83
imilarity index 100%
84 84
ename from floris/tools/optimization/scipy/cluster_turbines.py
85 85
ename to floris/tools/optimization/legacy/scipy/cluster_turbines.py
86 86
imilarity index 100%
87 87
ename from floris/tools/optimization/scipy/derive_downstream_turbines.py
88 88
ename to floris/tools/optimization/legacy/scipy/derive_downstream_turbines.py
89 89
imilarity index 100%
90 90
ename from floris/tools/optimization/scipy/layout.py
91 91
ename to floris/tools/optimization/legacy/scipy/layout.py
92 92
imilarity index 100%
93 93
ename from floris/tools/optimization/scipy/layout_height.py
94 94
ename to floris/tools/optimization/legacy/scipy/layout_height.py
95 95
imilarity index 100%
96 96
ename from floris/tools/optimization/scipy/optimization.py
97 97
ename to floris/tools/optimization/legacy/scipy/optimization.py
98 98
imilarity index 100%
99 99
ename from floris/tools/optimization/scipy/power_density.py
100 100
ename to floris/tools/optimization/legacy/scipy/power_density.py
101 101
imilarity index 100%
102 102
ename from floris/tools/optimization/scipy/power_density_1D.py
103 103
ename to floris/tools/optimization/legacy/scipy/power_density_1D.py
104 104
imilarity index 100%
105 105
ename from floris/tools/optimization/scipy/yaw.py
106 106
ename to floris/tools/optimization/legacy/scipy/yaw.py
107 107
imilarity index 100%
108 108
ename from floris/tools/optimization/scipy/yaw_clustered.py
109 109
ename to floris/tools/optimization/legacy/scipy/yaw_clustered.py
110 110
imilarity index 100%
111 111
ename from floris/tools/optimization/scipy/yaw_wind_rose.py
112 112
ename to floris/tools/optimization/legacy/scipy/yaw_wind_rose.py
113 113
imilarity index 100%
114 114
ename from floris/tools/optimization/scipy/yaw_wind_rose_clustered.py
115 115
ename to floris/tools/optimization/legacy/scipy/yaw_wind_rose_clustered.py
116 116
imilarity index 100%
117 117
ename from floris/tools/optimization/scipy/yaw_wind_rose_parallel.py
118 118
ename to floris/tools/optimization/legacy/scipy/yaw_wind_rose_parallel.py
119 119
imilarity index 100%
120 120
ename from floris/tools/optimization/scipy/yaw_wind_rose_parallel_clustered.py
121 121
ename to floris/tools/optimization/legacy/scipy/yaw_wind_rose_parallel_clustered.py

@@ -18,11 +18,13 @@
Loading
18 18
from scipy import integrate
19 19
from scipy.interpolate import RegularGridInterpolator
20 20
import scipy.io
21 -
import os
22 21
23 22
from floris.simulation import BaseModel
23 +
from floris.simulation import Farm
24 24
from floris.simulation import FlowField
25 25
from floris.simulation import Grid
26 +
from floris.simulation import Turbine
27 +
from floris.utilities import cosd, sind, tand
26 28
27 29
28 30
@define
@@ -36,7 +38,6 @@
Loading
36 38
    A: float = field(default=0.04)
37 39
    sigma_max_rel: float = field(default=4.0)
38 40
    overlap_gauss_interp: RegularGridInterpolator = field(init=False)
39 -
    model_string = "turbopark"
40 41
41 42
    def __attrs_post_init__(self) -> None:
42 43
        lookup_table_matlab_file = Path(__file__).parent / "turbopark_lookup_table.mat"
@@ -71,6 +72,7 @@
Loading
71 72
        rotor_diameter_i: np.ndarray,
72 73
        rotor_diameters: np.ndarray,
73 74
        i: int,
75 +
        deflection_field: np.ndarray,
74 76
        # enforces the use of the below as keyword arguments and adherence to the
75 77
        # unpacking of the results from prepare_function()
76 78
        *,
@@ -88,8 +90,8 @@
Loading
88 90
        x_dist = (x_i - x) * downstream_mask / rotor_diameters
89 91
90 92
        # Radial distance between turbine i and the centerlines of wakes from all real/image turbines
91 -
        r_dist = np.sqrt((y_i - y) ** 2 + (z_i - z) ** 2)
92 -
        r_dist_image = np.sqrt((y_i - y) ** 2 + (z_i - (-z)) ** 2)
93 +
        r_dist = np.sqrt((y_i - (y + deflection_field)) ** 2 + (z_i - z) ** 2)
94 +
        r_dist_image = np.sqrt((y_i - (y + deflection_field)) ** 2 + (z_i - (-z)) ** 2)
93 95
94 96
        Cts[:,:,i:,:,:] = 0.00001
95 97

@@ -13,9 +13,7 @@
Loading
13 13
from typing import Any, Dict
14 14
15 15
from attrs import define, field
16 -
import numexpr as ne
17 16
import numpy as np
18 -
from numpy import newaxis as na
19 17
from scipy.special import gamma
20 18
21 19
from floris.simulation import BaseModel
@@ -23,7 +21,7 @@
Loading
23 21
from floris.simulation import FlowField
24 22
from floris.simulation import Grid
25 23
from floris.simulation import Turbine
26 -
from floris.utilities import cosd, sind, tand, pshape
24 +
from floris.utilities import cosd, sind, tand
27 25
28 26
29 27
@define
@@ -38,8 +36,6 @@
Loading
38 36
    c_f: float = field(default=2.41)
39 37
    alpha_mod: float = field(default=1.0)
40 38
41 -
    model_string = "cumulative_gauss_curl"
42 -
43 39
    def prepare_function(
44 40
        self,
45 41
        grid: Grid,

@@ -16,7 +16,6 @@
Loading
16 16
import sys
17 17
18 18
from floris.simulation import Farm
19 -
from floris.simulation import Turbine
20 19
from floris.simulation import TurbineGrid, FlowFieldGrid
21 20
from floris.simulation import Ct, axial_induction
22 21
from floris.simulation import FlowField
@@ -134,6 +133,7 @@
Loading
134 133
            v_wake, w_wake = calculate_transverse_velocity(
135 134
                u_i,
136 135
                flow_field.u_initial_sorted,
136 +
                flow_field.dudz_initial_sorted,
137 137
                grid.x_sorted - x_i,
138 138
                grid.y_sorted - y_i,
139 139
                grid.z_sorted,
@@ -224,6 +224,7 @@
Loading
224 224
    turbine_grid_farm.construct_rotor_diameters()
225 225
    turbine_grid_farm.construct_turbine_TSRs()
226 226
    turbine_grid_farm.construc_turbine_pPs()
227 +
    turbine_grid_farm.construc_turbine_ref_density_cp_cts()
227 228
    turbine_grid_farm.construct_coordinates()
228 229
229 230
@@ -233,6 +234,7 @@
Loading
233 234
        wind_directions=turbine_grid_flow_field.wind_directions,
234 235
        wind_speeds=turbine_grid_flow_field.wind_speeds,
235 236
        grid_resolution=3,
237 +
        time_series=turbine_grid_flow_field.time_series,
236 238
    )
237 239
    turbine_grid_farm.expand_farm_properties(
238 240
        turbine_grid_flow_field.n_wind_directions, turbine_grid_flow_field.n_wind_speeds, turbine_grid.sorted_coord_indices
@@ -321,6 +323,7 @@
Loading
321 323
            v_wake, w_wake = calculate_transverse_velocity(
322 324
                u_i,
323 325
                flow_field.u_initial_sorted,
326 +
                flow_field.dudz_initial_sorted,
324 327
                flow_field_grid.x_sorted - x_i,
325 328
                flow_field_grid.y_sorted - y_i,
326 329
                flow_field_grid.z_sorted,
@@ -465,6 +468,7 @@
Loading
465 468
            v_wake, w_wake = calculate_transverse_velocity(
466 469
                u_i,
467 470
                flow_field.u_initial_sorted,
471 +
                flow_field.dudz_initial_sorted,
468 472
                grid.x_sorted - x_i,
469 473
                grid.y_sorted - y_i,
470 474
                grid.z_sorted,
@@ -551,6 +555,7 @@
Loading
551 555
    turbine_grid_farm.construct_rotor_diameters()
552 556
    turbine_grid_farm.construct_turbine_TSRs()
553 557
    turbine_grid_farm.construc_turbine_pPs()
558 +
    turbine_grid_farm.construc_turbine_ref_density_cp_cts()
554 559
    turbine_grid_farm.construct_coordinates()
555 560
556 561
    turbine_grid = TurbineGrid(
@@ -559,6 +564,7 @@
Loading
559 564
        wind_directions=turbine_grid_flow_field.wind_directions,
560 565
        wind_speeds=turbine_grid_flow_field.wind_speeds,
561 566
        grid_resolution=3,
567 +
        time_series=turbine_grid_flow_field.time_series,
562 568
    )
563 569
    turbine_grid_farm.expand_farm_properties(
564 570
        turbine_grid_flow_field.n_wind_directions, turbine_grid_flow_field.n_wind_speeds, turbine_grid.sorted_coord_indices
@@ -653,6 +659,7 @@
Loading
653 659
            v_wake, w_wake = calculate_transverse_velocity(
654 660
                u_i,
655 661
                flow_field.u_initial_sorted,
662 +
                flow_field.dudz_initial_sorted,
656 663
                flow_field_grid.x_sorted - x_i,
657 664
                flow_field_grid.y_sorted - y_i,
658 665
                flow_field_grid.z_sorted,
@@ -703,6 +710,7 @@
Loading
703 710
    w_wake = np.zeros_like(flow_field.w_initial_sorted)
704 711
    shape = (farm.n_turbines,) + np.shape(flow_field.u_initial_sorted)
705 712
    velocity_deficit = np.zeros(shape)
713 +
    deflection_field = np.zeros_like(flow_field.u_initial_sorted)
706 714
707 715
    turbine_turbulence_intensity = flow_field.turbulence_intensity * np.ones((flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1))
708 716
    ambient_turbulence_intensity = flow_field.turbulence_intensity
@@ -769,20 +777,43 @@
Loading
769 777
770 778
        # Model calculations
771 779
        # NOTE: exponential
772 -
        deflection_field = model_manager.deflection_model.function(
773 -
            x_i,
774 -
            y_i,
775 -
            effective_yaw_i,
776 -
            turbulence_intensity_i,
777 -
            ct_i,
778 -
            rotor_diameter_i,
779 -
            **deflection_model_args
780 -
        )
780 +
        if not np.all(farm.yaw_angles_sorted):
781 +
            model_manager.deflection_model.logger.warning("WARNING: Deflection with the TurbOPark model has not been fully validated. This is an initial implementation, and we advise you use at your own risk and perform a thorough examination of the results.")
782 +
            for ii in range(i):
783 +
                x_ii = np.mean(grid.x_sorted[:, :, ii:ii+1], axis=(3, 4))
784 +
                x_ii = x_ii[:, :, :, None, None]
785 +
                y_ii = np.mean(grid.y_sorted[:, :, ii:ii+1], axis=(3, 4))        
786 +
                y_ii = y_ii[:, :, :, None, None]
787 +
788 +
                yaw_ii = farm.yaw_angles_sorted[:, :, ii:ii+1, None, None]
789 +
                turbulence_intensity_ii = turbine_turbulence_intensity[:, :, ii:ii+1]
790 +
                ct_ii = Ct(
791 +
                    velocities=flow_field.u_sorted,
792 +
                    yaw_angle=farm.yaw_angles_sorted,
793 +
                    fCt=farm.turbine_fCts,
794 +
                    turbine_type_map=farm.turbine_type_map_sorted,
795 +
                    ix_filter=[ii]
796 +
                )
797 +
                ct_ii = ct_ii[:, :, 0:1, None, None]
798 +
                rotor_diameter_ii = farm.rotor_diameters_sorted[: ,:, ii:ii+1, None, None]
799 +
800 +
                deflection_field_ii = model_manager.deflection_model.function(
801 +
                    x_ii,
802 +
                    y_ii,
803 +
                    yaw_ii,
804 +
                    turbulence_intensity_ii,
805 +
                    ct_ii,
806 +
                    rotor_diameter_ii,
807 +
                    **deflection_model_args
808 +
                )
809 +
810 +
                deflection_field[:,:,ii:ii+1,:,:] = deflection_field_ii[:,:,i:i+1,:,:]
781 811
782 812
        if model_manager.enable_transverse_velocities:
783 813
            v_wake, w_wake = calculate_transverse_velocity(
784 814
                u_i,
785 815
                flow_field.u_initial_sorted,
816 +
                flow_field.dudz_initial_sorted,
786 817
                grid.x_sorted - x_i,
787 818
                grid.y_sorted - y_i,
788 819
                grid.z_sorted,
@@ -816,6 +847,7 @@
Loading
816 847
            rotor_diameter_i,
817 848
            farm.rotor_diameters_sorted[:, :, :, None, None],
818 849
            i,
850 +
            deflection_field,
819 851
            **deficit_model_args
820 852
        )
821 853

@@ -17,7 +17,6 @@
Loading
17 17
from attrs import define, field
18 18
import numpy as np
19 19
from pathlib import Path
20 -
import os
21 20
import copy
22 21
23 22
from floris.type_dec import (
@@ -26,7 +25,7 @@
Loading
26 25
    NDArrayFloat
27 26
)
28 27
from floris.utilities import Vec3, load_yaml
29 -
from floris.simulation import BaseClass
28 +
from floris.simulation import BaseClass, State
30 29
from floris.simulation import Turbine
31 30
32 31
@@ -83,10 +82,22 @@
Loading
83 82
            if type(val) is str:
84 83
                _floris_dir = Path(__file__).parent.parent
85 84
                fname = _floris_dir / "turbine_library" / f"{val}.yaml"
86 -
                if not os.path.isfile(fname):
85 +
                if not Path.is_file(fname):
87 86
                    raise ValueError("User-selected turbine definition `{}` does not exist in pre-defined turbine library.".format(val))
88 87
                self.turbine_definitions[i] = load_yaml(fname)
89 88
89 +
                # This is a temporary block of code that catches that ref_density_cp_ct is not defined
90 +
                # In the yaml file and forces it in
91 +
                # A warning is issued letting the user know in future versions defining this value explicitly
92 +
                # will be required 
93 +
                if not 'ref_density_cp_ct' in self.turbine_definitions[i]:
94 +
                    self.logger.warn("The value ref_density_cp_ct is not defined in the file: %s " % fname)
95 +
                    self.logger.warn("This value is not the simulated air density but is the density at which the cp/ct curves are defined")
96 +
                    self.logger.warn("In previous versions this was assumed to be 1.225")
97 +
                    self.logger.warn("Future versions of FLORIS will give an error if this value is not explicitly defined")
98 +
                    self.logger.warn("Currently this value is being set to the prior default value of 1.225")
99 +
                    self.turbine_definitions[i]['ref_density_cp_ct'] = 1.225
100 +
90 101
    def initialize(self, sorted_indices):
91 102
        # Sort yaw angles from most upstream to most downstream wind turbine
92 103
        self.yaw_angles_sorted = np.take_along_axis(
@@ -94,6 +105,7 @@
Loading
94 105
            sorted_indices[:, :, :, 0, 0],
95 106
            axis=2,
96 107
        )
108 +
        self.state = State.INITIALIZED
97 109
98 110
    def construct_hub_heights(self):
99 111
        self.hub_heights = np.array([turb['hub_height'] for turb in self.turbine_definitions])
@@ -107,6 +119,9 @@
Loading
107 119
    def construc_turbine_pPs(self):
108 120
        self.pPs = np.array([turb['pP'] for turb in self.turbine_definitions])
109 121
122 +
    def construc_turbine_ref_density_cp_cts(self):
123 +
        self.ref_density_cp_cts = np.array([turb['ref_density_cp_ct'] for turb in self.turbine_definitions])
124 +
110 125
    def construct_turbine_map(self):
111 126
        self.turbine_map = [Turbine.from_dict(turb) for turb in self.turbine_definitions]
112 127
@@ -149,6 +164,7 @@
Loading
149 164
        self.TSRs = np.take_along_axis(self.TSRs_sorted, unsorted_indices[:,:,:,0,0], axis=2)
150 165
        self.pPs = np.take_along_axis(self.pPs_sorted, unsorted_indices[:,:,:,0,0], axis=2)
151 166
        self.turbine_type_map = np.take_along_axis(self.turbine_type_map_sorted, unsorted_indices[:,:,:,0,0], axis=2)
167 +
        self.state.USED
152 168
153 169
    @property
154 170
    def n_turbines(self):

@@ -0,0 +1,114 @@
Loading
1 +
# Copyright 2022 NREL
2 +
3 +
# Licensed under the Apache License, Version 2.0 (the "License"); you may not
4 +
# use this file except in compliance with the License. You may obtain a copy of
5 +
# the License at http://www.apache.org/licenses/LICENSE-2.0
6 +
7 +
# Unless required by applicable law or agreed to in writing, software
8 +
# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
9 +
# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
10 +
# License for the specific language governing permissions and limitations under
11 +
# the License.
12 +
13 +
# See https://floris.readthedocs.io for documentation
14 +
15 +
import numpy as np
16 +
import matplotlib.pyplot as plt
17 +
from shapely.geometry import Polygon, LineString
18 +
19 +
from ....logging_manager import LoggerBase
20 +
21 +
class LayoutOptimization(LoggerBase):
22 +
    def __init__(self, fi, boundaries, min_dist=None, freq=None):
23 +
        self.fi = fi.copy()
24 +
        self.boundaries = boundaries
25 +
26 +
        self._boundary_polygon = Polygon(self.boundaries)
27 +
        self._boundary_line = LineString(self.boundaries)
28 +
29 +
        self.xmin = np.min([tup[0] for tup in boundaries])
30 +
        self.xmax = np.max([tup[0] for tup in boundaries])
31 +
        self.ymin = np.min([tup[1] for tup in boundaries])
32 +
        self.ymax = np.max([tup[1] for tup in boundaries])
33 +
34 +
        # If no minimum distance is provided, assume a value of 2 rotor diamters
35 +
        if min_dist is None:
36 +
            self.min_dist = 2 * self.rotor_diameter
37 +
        else:
38 +
            self.min_dist = min_dist
39 +
40 +
        # If freq is not provided, give equal weight to all wind conditions
41 +
        if freq is None:
42 +
            self.freq = np.ones((self.fi.floris.flow_field.n_wind_directions, self.fi.floris.flow_field.n_wind_speeds))
43 +
            self.freq = self.freq / self.freq.sum()
44 +
        else:
45 +
            self.freq = freq
46 +
47 +
        self.initial_AEP = fi.get_farm_AEP(self.freq)
48 +
49 +
    def __str__(self):
50 +
        return "layout"
51 +
52 +
    def _norm(self, val, x1, x2):
53 +
            return (val - x1) / (x2 - x1)
54 +
55 +
    def _unnorm(self, val, x1, x2):
56 +
        return np.array(val) * (x2 - x1) + x1
57 +
58 +
    # Public methods
59 +
60 +
    def optimize(self):
61 +
        sol = self._optimize()
62 +
        return sol
63 +
64 +
    def plot_layout_opt_results(self):
65 +
        x_initial, y_initial, x_opt, y_opt = self._get_initial_and_final_locs()
66 +
67 +
        plt.figure(figsize=(9, 6))
68 +
        fontsize = 16
69 +
        plt.plot(x_initial, y_initial, "ob")
70 +
        plt.plot(x_opt, y_opt, "or")
71 +
        # plt.title('Layout Optimization Results', fontsize=fontsize)
72 +
        plt.xlabel("x (m)", fontsize=fontsize)
73 +
        plt.ylabel("y (m)", fontsize=fontsize)
74 +
        plt.axis("equal")
75 +
        plt.grid()
76 +
        plt.tick_params(which="both", labelsize=fontsize)
77 +
        plt.legend(
78 +
            ["Old locations", "New locations"],
79 +
            loc="lower center",
80 +
            bbox_to_anchor=(0.5, 1.01),
81 +
            ncol=2,
82 +
            fontsize=fontsize,
83 +
        )
84 +
85 +
        verts = self.boundaries
86 +
        for i in range(len(verts)):
87 +
            if i == len(verts) - 1:
88 +
                plt.plot([verts[i][0], verts[0][0]], [verts[i][1], verts[0][1]], "b")
89 +
            else:
90 +
                plt.plot(
91 +
                    [verts[i][0], verts[i + 1][0]], [verts[i][1], verts[i + 1][1]], "b"
92 +
                )
93 +
        
94 +
        plt.show()
95 +
96 +
    ###########################################################################
97 +
    # Properties
98 +
    ###########################################################################
99 +
100 +
    @property
101 +
    def nturbs(self):
102 +
        """
103 +
        This property returns the number of turbines in the FLORIS
104 +
        object.
105 +
106 +
        Returns:
107 +
            nturbs (int): The number of turbines in the FLORIS object.
108 +
        """
109 +
        self._nturbs = self.fi.floris.farm.n_turbines
110 +
        return self._nturbs
111 +
112 +
    @property
113 +
    def rotor_diameter(self):
114 +
        return self.fi.floris.farm.rotor_diameters_sorted[0][0][0]

@@ -34,7 +34,8 @@
Loading
34 34
    wind_shear: float = field(converter=float)
35 35
    air_density: float = field(converter=float)
36 36
    turbulence_intensity: float = field(converter=float)
37 -
    reference_wind_height: float = field(converter=float)
37 +
    reference_wind_height: int = field(converter=int)
38 +
    time_series : bool = field(default=False)
38 39
39 40
    n_wind_speeds: int = field(init=False)
40 41
    n_wind_directions: int = field(init=False)
@@ -49,13 +50,17 @@
Loading
49 50
    v: NDArrayFloat = field(init=False, default=np.array([]))
50 51
    w: NDArrayFloat = field(init=False, default=np.array([]))
51 52
    het_map: list = field(init=False, default=None)
53 +
    dudz_initial_sorted: NDArrayFloat = field(init=False, default=np.array([]))
52 54
53 55
    turbulence_intensity_field: NDArrayFloat = field(init=False, default=np.array([]))
54 56
55 57
    @wind_speeds.validator
56 58
    def wind_speeds_validator(self, instance: attrs.Attribute, value: NDArrayFloat) -> None:
57 59
        """Using the validator method to keep the `n_wind_speeds` attribute up to date."""
58 -
        self.n_wind_speeds = value.size
60 +
        if self.time_series:
61 +
            self.n_wind_speeds = 1
62 +
        else:
63 +
            self.n_wind_speeds = value.size
59 64
60 65
    @wind_directions.validator
61 66
    def wind_directions_validator(self, instance: attrs.Attribute, value: NDArrayFloat) -> None:
@@ -74,6 +79,7 @@
Loading
74 79
        # for height, using it here to apply the shear law makes that dimension store the vertical
75 80
        # wind profile.
76 81
        wind_profile_plane = (grid.z_sorted / self.reference_wind_height) ** self.wind_shear
82 +
        dwind_profile_plane = self.wind_shear * (1 / self.reference_wind_height) ** self.wind_shear * (grid.z_sorted) ** (self.wind_shear - 1)
77 83
78 84
        # If no hetergeneous inflow defined, then set all speeds ups to 1.0
79 85
        if self.het_map is None:
@@ -93,7 +99,12 @@
Loading
93 99
        # here to do broadcasting from left to right (transposed), and then transpose back.
94 100
        # The result is an array the wind speed and wind direction dimensions on the left side
95 101
        # of the shape and the grid.template array on the right
96 -
        self.u_initial_sorted = (self.wind_speeds[None, :].T * wind_profile_plane.T).T * speed_ups
102 +
        if self.time_series:
103 +
            self.u_initial_sorted = (self.wind_speeds[:].T * wind_profile_plane.T).T * speed_ups
104 +
            self.dudz_initial_sorted = (self.wind_speeds[:].T * dwind_profile_plane.T).T * speed_ups
105 +
        else:
106 +
            self.u_initial_sorted = (self.wind_speeds[None, :].T * wind_profile_plane.T).T * speed_ups
107 +
            self.dudz_initial_sorted = (self.wind_speeds[None, :].T * dwind_profile_plane.T).T * speed_ups
97 108
        self.v_initial_sorted = np.zeros(np.shape(self.u_initial_sorted), dtype=self.u_initial_sorted.dtype)
98 109
        self.w_initial_sorted = np.zeros(np.shape(self.u_initial_sorted), dtype=self.u_initial_sorted.dtype)
99 110

@@ -0,0 +1,628 @@
Loading
1 +
# Copyright 2022 NREL
2 +
3 +
# Licensed under the Apache License, Version 2.0 (the "License"); you may not
4 +
# use this file except in compliance with the License. You may obtain a copy of
5 +
# the License at http://www.apache.org/licenses/LICENSE-2.0
6 +
7 +
# Unless required by applicable law or agreed to in writing, software
8 +
# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
9 +
# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
10 +
# License for the specific language governing permissions and limitations under
11 +
# the License.
12 +
13 +
# See https://floris.readthedocs.io for documentation
14 +
15 +
16 +
import numpy as np
17 +
import matplotlib.pyplot as plt
18 +
from shapely.geometry import Point, Polygon, LineString
19 +
from scipy.spatial.distance import cdist
20 +
21 +
from .layout_optimization_base import LayoutOptimization
22 +
23 +
class LayoutOptimizationBoundaryGrid(LayoutOptimization):
24 +
    def __init__(
25 +
        self,
26 +
        fi,
27 +
        boundaries,
28 +
        start,
29 +
        x_spacing,
30 +
        y_spacing,
31 +
        shear,
32 +
        rotation,
33 +
        center_x,
34 +
        center_y,
35 +
        boundary_setback,
36 +
        n_boundary_turbines=None,
37 +
        boundary_spacing=None,
38 +
    ):
39 +
        self.fi = fi
40 +
        
41 +
        self.boundary_x = np.array([val[0] for val in boundaries])
42 +
        self.boundary_y = np.array([val[1] for val in boundaries])
43 +
        boundary = np.zeros((len(self.boundary_x), 2))
44 +
        boundary[:, 0] = self.boundary_x[:]
45 +
        boundary[:, 1] = self.boundary_y[:]
46 +
        self._boundary_polygon = Polygon(boundary)
47 +
48 +
        self.start = start
49 +
        self.x_spacing = x_spacing
50 +
        self.y_spacing = y_spacing
51 +
        self.shear = shear
52 +
        self.rotation = rotation
53 +
        self.center_x = center_x
54 +
        self.center_y = center_y
55 +
        self.boundary_setback = boundary_setback
56 +
        self.n_boundary_turbines = n_boundary_turbines
57 +
        self.boundary_spacing = boundary_spacing
58 +
59 +
    def _discontinuous_grid(
60 +
        self,
61 +
        nrows,
62 +
        ncols,
63 +
        farm_width,
64 +
        farm_height,
65 +
        shear,
66 +
        rotation,
67 +
        center_x,
68 +
        center_y,
69 +
        shrink_boundary,
70 +
        boundary_x,
71 +
        boundary_y,
72 +
        eps=1e-3,
73 +
    ):
74 +
        """
75 +
        Map from grid design variables to turbine x and y locations. Includes integer design variables and the formulation
76 +
        results in a discontinous design space.
77 +
78 +
        TODO: shrink_boundary doesn't work well with concave boundaries, or with boundary angles less than 90 deg
79 +
80 +
        Args:
81 +
            nrows (Int): number of rows in the grid.
82 +
            ncols (Int): number of columns in the grid.
83 +
            farm_width (Float): total grid width (before shear).
84 +
            farm_height (Float): total grid height.
85 +
            shear (Float): grid shear (rad).
86 +
            rotation (Float): rotation about grid center (rad).
87 +
            center_x (Float): location of grid x center.
88 +
            center_y (Float): location of grid y center.
89 +
            shrink_boundary (Float): how much to shrink the boundary that the grid can occupy.
90 +
            boundary_x (Array(Float)): x boundary points.
91 +
            boundary_y (Array(Float)): y boundary points.
92 +
93 +
        Returns:
94 +
            grid_x (Array(Float)): turbine x locations.
95 +
            grid_y (Array(Float)): turbine y locations.
96 +
        """
97 +
        # create grid
98 +
        nrows = int(nrows)
99 +
        ncols = int(ncols)
100 +
        xlocs = np.linspace(0.0, farm_width, ncols)
101 +
        ylocs = np.linspace(0.0, farm_height, nrows)
102 +
        y_spacing = ylocs[1] - ylocs[0]
103 +
        nturbs = nrows * ncols
104 +
        grid_x = np.zeros(nturbs)
105 +
        grid_y = np.zeros(nturbs)
106 +
        turb = 0
107 +
        for i in range(nrows):
108 +
            for j in range(ncols):
109 +
                grid_x[turb] = xlocs[j] + float(i) * y_spacing * np.tan(shear)
110 +
                grid_y[turb] = ylocs[i]
111 +
                turb += 1
112 +
113 +
        # rotate
114 +
        grid_x, grid_y = (
115 +
            np.cos(rotation) * grid_x - np.sin(rotation) * grid_y,
116 +
            np.sin(rotation) * grid_x + np.cos(rotation) * grid_y,
117 +
        )
118 +
119 +
        # move center of grid
120 +
        grid_x = (grid_x - np.mean(grid_x)) + center_x
121 +
        grid_y = (grid_y - np.mean(grid_y)) + center_y
122 +
123 +
        # arrange the boundary
124 +
125 +
        # boundary = np.zeros((len(boundary_x),2))
126 +
        # boundary[:,0] = boundary_x[:]
127 +
        # boundary[:,1] = boundary_y[:]
128 +
        # poly = Polygon(boundary)
129 +
        # centroid = poly.centroid
130 +
131 +
        # boundary[:,0] = (boundary_x[:]-centroid.x)*boundary_mult + centroid.x
132 +
        # boundary[:,1] = (boundary_y[:]-centroid.y)*boundary_mult + centroid.y
133 +
        # poly = Polygon(boundary)
134 +
135 +
        boundary = np.zeros((len(boundary_x), 2))
136 +
        boundary[:, 0] = boundary_x[:]
137 +
        boundary[:, 1] = boundary_y[:]
138 +
        poly = Polygon(boundary)
139 +
140 +
        if shrink_boundary != 0.0:
141 +
            nBounds = len(boundary_x)
142 +
            for i in range(nBounds):
143 +
                point = Point(boundary_x[i] + eps, boundary_y[i])
144 +
                if poly.contains(point) is True or poly.touches(point) is True:
145 +
                    boundary[i, 0] = boundary_x[i] + shrink_boundary
146 +
                else:
147 +
                    boundary[i, 0] = boundary_x[i] - shrink_boundary
148 +
149 +
                point = Point(boundary_x[i], boundary_y[i] + eps)
150 +
                if poly.contains(point) is True or poly.touches(point) is True:
151 +
                    boundary[i, 1] = boundary_y[i] + shrink_boundary
152 +
                else:
153 +
                    boundary[i, 1] = boundary_y[i] - shrink_boundary
154 +
155 +
            poly = Polygon(boundary)
156 +
157 +
        # get rid of points outside of boundary
158 +
        index = 0
159 +
        for i in range(len(grid_x)):
160 +
            point = Point(grid_x[index], grid_y[index])
161 +
            if poly.contains(point) is False and poly.touches(point) is False:
162 +
                grid_x = np.delete(grid_x, index)
163 +
                grid_y = np.delete(grid_y, index)
164 +
            else:
165 +
                index += 1
166 +
167 +
        return grid_x, grid_y
168 +
169 +
    def _discrete_grid(
170 +
        self,
171 +
        x_spacing,
172 +
        y_spacing,
173 +
        shear,
174 +
        rotation,
175 +
        center_x,
176 +
        center_y,
177 +
        boundary_setback,
178 +
        boundary_poly
179 +
    ):
180 +
        """
181 +
        returns grid turbine layout. Assumes the turbines fill the entire plant area
182 +
183 +
        Args:
184 +
        x_spacing (Float): grid spacing in the unrotated x direction (m)
185 +
        y_spacing (Float): grid spacing in the unrotated y direction (m)
186 +
        shear (Float): grid shear (rad)
187 +
        rotation (Float): grid rotation (rad)
188 +
        center_x (Float): the x coordinate of the grid center (m)
189 +
        center_y (Float): the y coordinate of the grid center (m)
190 +
        boundary_poly (Polygon): a shapely Polygon of the wind plant boundary
191 +
192 +
        Returns
193 +
        return_x (Array(Float)): turbine x locations
194 +
        return_y (Array(Float)): turbine y locations
195 +
        """
196 +
197 +
        shrunk_poly = boundary_poly.buffer(-boundary_setback)
198 +
        if shrunk_poly.area <= 0:
199 +
            return np.array([]), np.array([])
200 +
        # create grid
201 +
        minx, miny, maxx, maxy = shrunk_poly.bounds
202 +
        width = maxx-minx
203 +
        height = maxy-miny
204 +
205 +
        center_point = Point((center_x,center_y))
206 +
        poly_to_center = center_point.distance(shrunk_poly.centroid)
207 +
208 +
        width = np.max([width,poly_to_center])
209 +
        height = np.max([height,poly_to_center])
210 +
        nrows = int(np.max([width,height])/np.min([x_spacing,y_spacing]))*2 + 1
211 +
        ncols = nrows
212 +
213 +
        xlocs = np.arange(0,ncols)*x_spacing
214 +
        ylocs = np.arange(0,nrows)*y_spacing
215 +
        row_number = np.arange(0,nrows)
216 +
217 +
        d = np.array([i for x in xlocs for i in row_number])
218 +
        layout_x = np.array([x for x in xlocs for y in ylocs]) + d*y_spacing*np.tan(shear)
219 +
        layout_y = np.array([y for x in xlocs for y in ylocs])
220 +
        
221 +
        # rotate
222 +
        rotate_x = np.cos(rotation)*layout_x - np.sin(rotation)*layout_y
223 +
        rotate_y = np.sin(rotation)*layout_x + np.cos(rotation)*layout_y
224 +
225 +
        # move center of grid
226 +
        rotate_x = (rotate_x - np.mean(rotate_x)) + center_x
227 +
        rotate_y = (rotate_y - np.mean(rotate_y)) + center_y
228 +
229 +
        # get rid of points outside of boundary polygon
230 +
        meets_constraints = np.zeros(len(rotate_x),dtype=bool)
231 +
        for i in range(len(rotate_x)):
232 +
            pt = Point(rotate_x[i],rotate_y[i])
233 +
            if shrunk_poly.contains(pt) or shrunk_poly.touches(pt):
234 +
                meets_constraints[i] = True
235 +
236 +
        # arrange final x,y points
237 +
        return_x = rotate_x[meets_constraints]
238 +
        return_y = rotate_y[meets_constraints]
239 +
    
240 +
        return return_x, return_y
241 +
242 +
    def find_lengths(self, x, y, npoints):
243 +
        length = np.zeros(len(x) - 1)
244 +
        for i in range(npoints):
245 +
            length[i] = np.sqrt((x[i + 1] - x[i]) ** 2 + (y[i + 1] - y[i]) ** 2)
246 +
        return length
247 +
248 +
    # def _place_boundary_turbines(self, n_boundary_turbs, start, boundary_x, boundary_y):
249 +
    #     """
250 +
    #     Place turbines equally spaced traversing the perimiter if the wind farm along the boundary
251 +
252 +
    #     Args:
253 +
    #     n_boundary_turbs (Int): number of turbines to be placed on the boundary
254 +
    #     start (Float): where the first turbine should be placed
255 +
    #     boundary_x (Array(Float)): x boundary points
256 +
    #     boundary_y (Array(Float)): y boundary points
257 +
258 +
    #     Returns
259 +
    #     layout_x (Array(Float)): turbine x locations
260 +
    #     layout_y (Array(Float)): turbine y locations
261 +
    #     """
262 +
263 +
    #     # check if the boundary is closed, correct if not
264 +
    #     if boundary_x[-1] != boundary_x[0] or boundary_y[-1] != boundary_y[0]:
265 +
    #         boundary_x = np.append(boundary_x, boundary_x[0])
266 +
    #         boundary_y = np.append(boundary_y, boundary_y[0])
267 +
268 +
    #     # make the boundary
269 +
    #     boundary = np.zeros((len(boundary_x), 2))
270 +
    #     boundary[:, 0] = boundary_x[:]
271 +
    #     boundary[:, 1] = boundary_y[:]
272 +
    #     poly = Polygon(boundary)
273 +
    #     perimeter = poly.length
274 +
275 +
    #     # get the flattened turbine locations
276 +
    #     spacing = perimeter / float(n_boundary_turbs)
277 +
    #     flattened_locs = np.linspace(start, perimeter + start - spacing, n_boundary_turbs)
278 +
279 +
    #     # set all of the flattened values between 0 and the perimeter
280 +
    #     for i in range(n_boundary_turbs):
281 +
    #         while flattened_locs[i] < 0.0:
282 +
    #             flattened_locs[i] += perimeter
283 +
    #         if flattened_locs[i] > perimeter:
284 +
    #             flattened_locs[i] = flattened_locs[i] % perimeter
285 +
286 +
    #     # place the turbines around the perimeter
287 +
    #     nBounds = len(boundary_x)
288 +
    #     layout_x = np.zeros(n_boundary_turbs)
289 +
    #     layout_y = np.zeros(n_boundary_turbs)
290 +
291 +
    #     lenBound = np.zeros(nBounds - 1)
292 +
    #     for i in range(nBounds - 1):
293 +
    #         lenBound[i] = Point(boundary[i]).distance(Point(boundary[i + 1]))
294 +
    #     for i in range(n_boundary_turbs):
295 +
    #         for j in range(nBounds - 1):
296 +
    #             if flattened_locs[i] < sum(lenBound[0 : j + 1]):
297 +
    #                 layout_x[i] = (
298 +
    #                     boundary_x[j]
299 +
    #                     + (boundary_x[j + 1] - boundary_x[j])
300 +
    #                     * (flattened_locs[i] - sum(lenBound[0:j]))
301 +
    #                     / lenBound[j]
302 +
    #                 )
303 +
    #                 layout_y[i] = (
304 +
    #                     boundary_y[j]
305 +
    #                     + (boundary_y[j + 1] - boundary_y[j])
306 +
    #                     * (flattened_locs[i] - sum(lenBound[0:j]))
307 +
    #                     / lenBound[j]
308 +
    #                 )
309 +
    #                 break
310 +
311 +
    #     return layout_x, layout_y
312 +
313 +
    def _place_boundary_turbines(self, start, boundary_poly, nturbs=None, spacing=None):
314 +
        xBounds, yBounds = boundary_poly.boundary.coords.xy
315 +
316 +
        if xBounds[-1] != xBounds[0]:
317 +
            xBounds = np.append(xBounds, xBounds[0])
318 +
            yBounds = np.append(yBounds, yBounds[0])
319 +
320 +
        nBounds = len(xBounds)
321 +
        lenBound = self.find_lengths(xBounds, yBounds, len(xBounds) - 1)
322 +
        circumference = sum(lenBound)
323 +
        
324 +
        if nturbs is not None and spacing is None:
325 +
            # When the number of boundary turbines is specified
326 +
            nturbs = int(nturbs)
327 +
            bound_loc = np.linspace(
328 +
                start, start + circumference - circumference / float(nturbs), nturbs
329 +
            )
330 +
        elif spacing is not None and nturbs is None:
331 +
            # When the spacing of boundary turbines is specified
332 +
            nturbs = int(np.floor(circumference / spacing))
333 +
            bound_loc = np.linspace(
334 +
                start, start + circumference - circumference / float(nturbs), nturbs
335 +
            )
336 +
        else:
337 +
            raise ValueError("Please specify either nturbs or spacing.")
338 +
339 +
        x = np.zeros(nturbs)
340 +
        y = np.zeros(nturbs)
341 +
342 +
        if spacing is None:
343 +
            # When the number of boundary turbines is specified
344 +
            for i in range(nturbs):
345 +
                if bound_loc[i] > circumference:
346 +
                    bound_loc[i] = bound_loc[i] % circumference
347 +
                while bound_loc[i] < 0.0:
348 +
                    bound_loc[i] += circumference
349 +
            for i in range(nturbs):
350 +
                done = False
351 +
                for j in range(nBounds):
352 +
                    if done == False:
353 +
                        if bound_loc[i] < sum(lenBound[0:j+1]):
354 +
                            point_x = xBounds[j] + (xBounds[j+1]-xBounds[j])*(bound_loc[i]-sum(lenBound[0:j]))/lenBound[j]
355 +
                            point_y = yBounds[j] + (yBounds[j+1]-yBounds[j])*(bound_loc[i]-sum(lenBound[0:j]))/lenBound[j]
356 +
                            done = True
357 +
                            x[i] = point_x
358 +
                            y[i] = point_y
359 +
        else:
360 +
            # When the spacing of boundary turbines is specified
361 +
            additional_space = 0.0
362 +
            end_loop = False
363 +
            for i in range(nturbs):
364 +
                done = False
365 +
                for j in range(nBounds):
366 +
                    while done == False:
367 +
                        dist = start + i*spacing + additional_space
368 +
                        if dist < sum(lenBound[0:j+1]):
369 +
                            point_x = xBounds[j] + (xBounds[j+1]-xBounds[j])*(dist -sum(lenBound[0:j]))/lenBound[j]
370 +
                            point_y = yBounds[j] + (yBounds[j+1]-yBounds[j])*(dist -sum(lenBound[0:j]))/lenBound[j]
371 +
372 +
                            # Check if turbine is too close to previous turbine
373 +
                            if i > 0:
374 +
                                # Check if turbine just placed is to close to first turbine
375 +
                                min_dist = cdist([(point_x, point_y)], [(x[0], y[0])])
376 +
                                if min_dist < spacing:
377 +
                                    # TODO: make this more robust; pass is needed if 2nd turbine is too close to the first
378 +
                                    if i == 1:
379 +
                                        pass
380 +
                                    else:
381 +
                                        end_loop = True
382 +
                                        ii = i
383 +
                                        break
384 +
385 +
                                min_dist = cdist([(point_x, point_y)], [(x[i-1], y[i-1])])
386 +
                                if min_dist < spacing:
387 +
                                    additional_space += 1.0
388 +
                                else:
389 +
                                    done = True
390 +
                                    x[i] = point_x
391 +
                                    y[i] = point_y
392 +
                            elif i == 0:
393 +
                                # If first turbine, just add initial turbine point
394 +
                                done = True
395 +
                                x[i] = point_x
396 +
                                y[i] = point_y
397 +
                            else:
398 +
                                pass
399 +
                        else:
400 +
                            break
401 +
                    if end_loop == True:
402 +
                        break
403 +
                if end_loop == True:
404 +
                    x = x[:ii]
405 +
                    y = y[:ii]
406 +
                    break
407 +
        return x, y
408 +
409 +
    def _place_boundary_turbines_with_specified_spacing(self, spacing, start, boundary_x, boundary_y):
410 +
        """
411 +
        Place turbines equally spaced traversing the perimiter if the wind farm along the boundary
412 +
413 +
        Args:
414 +
        n_boundary_turbs (Int): number of turbines to be placed on the boundary
415 +
        start (Float): where the first turbine should be placed
416 +
        boundary_x (Array(Float)): x boundary points
417 +
        boundary_y (Array(Float)): y boundary points
418 +
419 +
        Returns
420 +
        layout_x (Array(Float)): turbine x locations
421 +
        layout_y (Array(Float)): turbine y locations
422 +
        """
423 +
424 +
        # check if the boundary is closed, correct if not
425 +
        if boundary_x[-1] != boundary_x[0] or boundary_y[-1] != boundary_y[0]:
426 +
            boundary_x = np.append(boundary_x, boundary_x[0])
427 +
            boundary_y = np.append(boundary_y, boundary_y[0])
428 +
429 +
        # make the boundary
430 +
        boundary = np.zeros((len(boundary_x), 2))
431 +
        boundary[:, 0] = boundary_x[:]
432 +
        boundary[:, 1] = boundary_y[:]
433 +
        poly = Polygon(boundary)
434 +
        perimeter = poly.length
435 +
436 +
        # get the flattened turbine locations
437 +
        n_boundary_turbs = int(perimeter / float(spacing))
438 +
        flattened_locs = np.linspace(start, perimeter + start - spacing, n_boundary_turbs)
439 +
440 +
        # set all of the flattened values between 0 and the perimeter
441 +
        for i in range(n_boundary_turbs):
442 +
            while flattened_locs[i] < 0.0:
443 +
                flattened_locs[i] += perimeter
444 +
            if flattened_locs[i] > perimeter:
445 +
                flattened_locs[i] = flattened_locs[i] % perimeter
446 +
447 +
        # place the turbines around the perimeter
448 +
        nBounds = len(boundary_x)
449 +
        layout_x = np.zeros(n_boundary_turbs)
450 +
        layout_y = np.zeros(n_boundary_turbs)
451 +
452 +
        lenBound = np.zeros(nBounds - 1)
453 +
        for i in range(nBounds - 1):
454 +
            lenBound[i] = Point(boundary[i]).distance(Point(boundary[i + 1]))
455 +
        for i in range(n_boundary_turbs):
456 +
            for j in range(nBounds - 1):
457 +
                if flattened_locs[i] < sum(lenBound[0 : j + 1]):
458 +
                    layout_x[i] = (
459 +
                        boundary_x[j]
460 +
                        + (boundary_x[j + 1] - boundary_x[j])
461 +
                        * (flattened_locs[i] - sum(lenBound[0:j]))
462 +
                        / lenBound[j]
463 +
                    )
464 +
                    layout_y[i] = (
465 +
                        boundary_y[j]
466 +
                        + (boundary_y[j + 1] - boundary_y[j])
467 +
                        * (flattened_locs[i] - sum(lenBound[0:j]))
468 +
                        / lenBound[j]
469 +
                    )
470 +
                    break
471 +
472 +
        return layout_x, layout_y
473 +
474 +
    def boundary_grid(
475 +
        self,
476 +
        start,
477 +
        x_spacing,
478 +
        y_spacing,
479 +
        shear,
480 +
        rotation,
481 +
        center_x,
482 +
        center_y,
483 +
        boundary_setback,
484 +
        n_boundary_turbines=None,
485 +
        boundary_spacing=None,
486 +
    ):
487 +
        """
488 +
        Place turbines equally spaced traversing the perimiter if the wind farm along the boundary
489 +
490 +
        Args:
491 +
        n_boundary_turbs,start: boundary variables
492 +
        nrows,ncols,farm_width,farm_height,shear,rotation,center_x,center_y,shrink_boundary,eps: grid variables
493 +
        boundary_x,boundary_y: boundary points
494 +
495 +
        Returns
496 +
        layout_x (Array(Float)): turbine x locations
497 +
        layout_y (Array(Float)): turbine y locations
498 +
        """
499 +
500 +
        boundary_turbines_x, boundary_turbines_y = self._place_boundary_turbines(
501 +
            start, self._boundary_polygon, nturbs=n_boundary_turbines, spacing=boundary_spacing
502 +
        )
503 +
        # boundary_turbines_x, boundary_turbines_y = self._place_boundary_turbines_with_specified_spacing(
504 +
        #     spacing, start, boundary_x, boundary_y
505 +
        # )
506 +
507 +
        # grid_turbines_x, grid_turbines_y = self._discontinuous_grid(
508 +
        #     nrows,
509 +
        #     ncols,
510 +
        #     farm_width,
511 +
        #     farm_height,
512 +
        #     shear,
513 +
        #     rotation,
514 +
        #     center_x,
515 +
        #     center_y,
516 +
        #     shrink_boundary,
517 +
        #     boundary_x,
518 +
        #     boundary_y,
519 +
        #     eps=eps,
520 +
        # )
521 +
522 +
        grid_turbines_x, grid_turbines_y = self._discrete_grid(
523 +
            x_spacing,
524 +
            y_spacing,
525 +
            shear,
526 +
            rotation,
527 +
            center_x,
528 +
            center_y,
529 +
            boundary_setback,
530 +
            self._boundary_polygon,
531 +
        )
532 +
533 +
        layout_x = np.append(boundary_turbines_x, grid_turbines_x)
534 +
        layout_y = np.append(boundary_turbines_y, grid_turbines_y)
535 +
536 +
        return layout_x, layout_y
537 +
538 +
    def reinitialize_bg(
539 +
        self,
540 +
        n_boundary_turbines=None,
541 +
        start=None,
542 +
        x_spacing=None,
543 +
        y_spacing=None,
544 +
        shear=None,
545 +
        rotation=None,
546 +
        center_x=None,
547 +
        center_y=None,
548 +
        boundary_setback=None,
549 +
        boundary_x=None,
550 +
        boundary_y=None,
551 +
        boundary_spacing=None,
552 +
    ):
553 +
554 +
        if n_boundary_turbines is not None:
555 +
            self.n_boundary_turbines = n_boundary_turbines
556 +
        if start is not None:
557 +
            self.start = start
558 +
        if x_spacing is not None:
559 +
            self.x_spacing = x_spacing
560 +
        if y_spacing is not None:
561 +
            self.y_spacing = y_spacing
562 +
        if shear is not None:
563 +
            self.shear = shear
564 +
        if rotation is not None:
565 +
            self.rotation = rotation
566 +
        if center_x is not None:
567 +
            self.center_x = center_x
568 +
        if center_y is not None:
569 +
            self.center_y = center_y
570 +
        if boundary_setback is not None:
571 +
            self.boundary_setback = boundary_setback
572 +
        if boundary_x is not None:
573 +
            self.boundary_x = boundary_x
574 +
        if boundary_y is not None:
575 +
            self.boundary_y = boundary_y
576 +
        if boundary_spacing is not None:
577 +
            self.boundary_spacing = boundary_spacing
578 +
579 +
    def reinitialize_xy(self):
580 +
        layout_x, layout_y = self.boundary_grid(
581 +
            self.start,
582 +
            self.x_spacing,
583 +
            self.y_spacing,
584 +
            self.shear,
585 +
            self.rotation,
586 +
            self.center_x,
587 +
            self.center_y,
588 +
            self.boundary_setback,
589 +
            self.n_boundary_turbines,
590 +
            self.boundary_spacing,
591 +
        )
592 +
593 +
        self.fi.reinitialize(layout=(layout_x, layout_y))
594 +
595 +
    def plot_layout(self):
596 +
        plt.figure(figsize=(9, 6))
597 +
        fontsize = 16
598 +
599 +
        plt.plot(self.fi.layout_x, self.fi.layout_y, "ob")
600 +
        # plt.plot(locsx, locsy, "or")
601 +
602 +
        plt.xlabel("x (m)", fontsize=fontsize)
603 +
        plt.ylabel("y (m)", fontsize=fontsize)
604 +
        plt.axis("equal")
605 +
        plt.grid()
606 +
        plt.tick_params(which="both", labelsize=fontsize)
607 +
608 +
        plt.show()
609 +
610 +
    def space_constraint(self, x, y, min_dist, rho=500):
611 +
        # Calculate distances between turbines
612 +
        locs = np.vstack((x, y)).T
613 +
        distances = cdist(locs, locs)
614 +
        arange = np.arange(distances.shape[0])
615 +
        distances[arange, arange] = 1e10
616 +
        dist = np.min(distances, axis=0)
617 +
618 +
        g = 1 - np.array(dist) / min_dist
619 +
620 +
        # Following code copied from OpenMDAO KSComp().
621 +
        # Constraint is satisfied when KS_constraint <= 0
622 +
        g_max = np.max(np.atleast_2d(g), axis=-1)[:, np.newaxis]
623 +
        g_diff = g - g_max
624 +
        exponents = np.exp(rho * g_diff)
625 +
        summation = np.sum(exponents, axis=-1)[:, np.newaxis]
626 +
        KS_constraint = g_max + 1.0 / rho * np.log(summation)
627 +
628 +
        return KS_constraint[0][0], dist

@@ -0,0 +1,232 @@
Loading
1 +
# Copyright 2022 NREL
2 +
3 +
# Licensed under the Apache License, Version 2.0 (the "License"); you may not
4 +
# use this file except in compliance with the License. You may obtain a copy of
5 +
# the License at http://www.apache.org/licenses/LICENSE-2.0
6 +
7 +
# Unless required by applicable law or agreed to in writing, software
8 +
# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
9 +
# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
10 +
# License for the specific language governing permissions and limitations under
11 +
# the License.
12 +
13 +
# See https://floris.readthedocs.io for documentation
14 +
15 +
import numpy as np
16 +
import matplotlib.pyplot as plt
17 +
from scipy.optimize import minimize
18 +
from shapely.geometry import Point
19 +
from scipy.spatial.distance import cdist
20 +
21 +
from .layout_optimization_base import LayoutOptimization
22 +
23 +
class LayoutOptimizationScipy(LayoutOptimization):
24 +
    def __init__(
25 +
        self,
26 +
        fi,
27 +
        boundaries,
28 +
        freq=None,
29 +
        bnds=None,
30 +
        min_dist=None,
31 +
        solver='SLSQP',
32 +
        optOptions=None,
33 +
    ):
34 +
        """
35 +
        _summary_
36 +
37 +
        Args:
38 +
            fi (_type_): _description_
39 +
            boundaries (iterable(float, float)): Pairs of x- and y-coordinates
40 +
                that represent the boundary's vertices (m).
41 +
            freq (np.array): An array of the frequencies of occurance
42 +
                correponding to each pair of wind direction and wind speed
43 +
                values. If None, equal weight is given to each pair of wind conditions
44 +
                Defaults to None.
45 +
            bnds (iterable, optional): Bounds for the optimization
46 +
                variables (pairs of min/max values for each variable (m)). If
47 +
                none are specified, they are set to 0 and 1. Defaults to None.
48 +
            min_dist (float, optional): The minimum distance to be maintained
49 +
                between turbines during the optimization (m). If not specified,
50 +
                initializes to 2 rotor diameters. Defaults to None.
51 +
            solver (str, optional): Sets the solver used by Scipy. Defaults to 'SLSQP'.
52 +
            optOptions (dict, optional): Dicitonary for setting the
53 +
                optimization options. Defaults to None.
54 +
        """
55 +
        super().__init__(fi, boundaries, min_dist=min_dist, freq=freq)
56 +
57 +
        self.boundaries_norm = [
58 +
            [
59 +
                self._norm(val[0], self.xmin, self.xmax),
60 +
                self._norm(val[1], self.ymin, self.ymax),
61 +
            ]
62 +
            for val in self.boundaries
63 +
        ]
64 +
        self.x0 = [
65 +
            self._norm(x, self.xmin, self.xmax)
66 +
            for x in self.fi.layout_x
67 +
        ] + [
68 +
            self._norm(y, self.ymin, self.ymax)
69 +
            for y in self.fi.layout_y
70 +
        ]
71 +
        if bnds is not None:
72 +
            self.bnds = bnds
73 +
        else:
74 +
            self._set_opt_bounds()
75 +
        if solver is not None:
76 +
            self.solver = solver
77 +
        if optOptions is not None:
78 +
            self.optOptions = optOptions
79 +
        else:
80 +
            self.optOptions = {"maxiter": 100, "disp": True, "iprint": 2, "ftol": 1e-9, "eps":0.01}
81 +
82 +
        self._generate_constraints()
83 +
84 +
85 +
    # Private methods
86 +