Sg/tilted bem dev

examples/MITWindfarm_quickstart.ipynb

@@ -20,7 +20,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 25,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -38,7 +38,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 26,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -71,7 +71,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 27,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -108,7 +108,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -129,7 +129,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 29,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -144,14 +144,14 @@
     "\n",
     "Analogously to how we are currently building a `WindFarm` object so that we can solve for a `WindFarmSolution`, a `Rotor` has a `RotorSolution`. So, if you are only interested in a single rotor solution, you can now solve for that rotor's solution, rather than building an entire wind farm, which will solve for all of the rotor's solutions together in tandem.\n",
     "\n",
-    "For the `AD` and `UniformAD` models, you need to provide a `x`, `y`, and `z` arguments that specify the location you want the solution for, measured in distance from the centroid of the rotor. You also need to provide `windfield` model, as well as `Ctprime` and `yaw` values. We will later refer to `Ctprime` and `yaw` as \"setpoints\".\n",
+    "For the `AD` and `UniformAD` models, you need to provide a `x`, `y`, and `z` arguments that specify the location you want the solution for, measured in distance from the centroid of the rotor. You also need to provide `windfield` model, as well as `Ctprime`, `yaw`, and `tilt` values. We will later refer to `Ctprime`, `yaw`, and `tilt` as \"setpoints\".\n",
     "\n",
-    "Similarly, for the `BEM` model, you also need to provide `x`, `y`, `z`, and `windfield` arguments. The `BEM` setpoints, however are different: `pitch`, `tsr` (tip speed ratio), and `yaw`."
+    "Similarly, for the `BEM` model, you also need to provide `x`, `y`, `z`, and `windfield` arguments. The `BEM` setpoints, however are different: `pitch`, `tsr` (tip speed ratio), `yaw`, and `tilt`. Note that within `MITRotor`, there are a few specifications on which setups support tilt. See the `MITRotor` quickstart guide for that information."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 30,
    "metadata": {},
    "outputs": [
     {
@@ -184,7 +184,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 31,
    "metadata": {},
    "outputs": [
     {
@@ -235,7 +235,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 32,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -261,7 +261,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 33,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -286,7 +286,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 34,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -320,7 +320,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 35,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -334,7 +334,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 36,
    "metadata": {},
    "outputs": [
     {
@@ -374,7 +374,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 37,
    "metadata": {},
    "outputs": [
     {
@@ -395,7 +395,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 38,
    "metadata": {},
    "outputs": [
     {
@@ -420,12 +420,12 @@
    "source": [
     "### Wind Farm `setpoints`\n",
     "\n",
-    "Setpoints are properties specific to the rotor model. For both the `AD` and `UnifiedAD` models, the `setpoints` are `Ctprime`, `yaw`, and `tilt`. For `BEM`, the `setpoints` are `pitch`, `tsr`, and `yaw` (`tilt` is not yet implemnted for BEM). A tuple of setpoints must be provided for each turbine in a wind farm, in the form of a vector of tuples."
+    "Setpoints are properties specific to the rotor model. For both the `AD` and `UnifiedAD` models, the `setpoints` are `Ctprime`, `yaw`, and `tilt`. For `BEM`, the `setpoints` are `pitch`, `tsr`, and `yaw`, and `tilt`. A tuple of setpoints must be provided for each turbine in a wind farm, in the form of a vector of tuples."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 39,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -448,7 +448,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 40,
    "metadata": {},
    "outputs": [
     {

examples/example_05_basic_windfarm.py

@@ -1,4 +1,5 @@
 from pathlib import Path
+import numpy as np
 
 import matplotlib.pyplot as plt
 from MITRotor import IEA15MW
@@ -12,10 +13,10 @@
 if __name__ == "__main__":
     windfarm = Windfarm(rotor_model=BEM(IEA15MW()))
 
-    setpoints = [
-        (0, 7, 0),
-        (0, 7, 0),
-        (0, 7, 0),
+    setpoints = [ # pitch, tsr, yaw, tilt
+        (0, 7, 0, 0),
+        (0, 7, 0, 0),
+        (0, 7, 0, 0),
     ]
 
     layout = Layout([0, 12 / 2, 24 / 2], [0, 0, 0])

examples/example_07_BEM_yaw_optimisation.py

@@ -6,43 +6,39 @@
 from mitwindfarm import Plotting
 from mitwindfarm.windfarm import Windfarm, WindfarmSolution, Layout
 from mitwindfarm.Rotor import BEM
+from MITRotor.Momentum import UnifiedMomentum
 from MITRotor.ReferenceTurbines import IEA15MW
 from rich import print
 import numpy as np
 
 FIGDIR = Path(__file__).parent.parent / "fig"
 FIGDIR.mkdir(exist_ok=True, parents=True)
 
-windfarm = Windfarm(rotor_model=BEM(IEA15MW()), TIamb=0.1)
-layout = Layout([0, 12, 24], [0.0, 0.5, 1.0])
-
+momentum_model= UnifiedMomentum(averaging = "rotor")
+windfarm = Windfarm(
+    rotor_model=BEM(IEA15MW(), momentum_model = momentum_model),
+    TIamb=0.1,
+)
+layout = Layout([0, 12, 24], [0.0, 0.5, 1.0], [0.0, 0.0, 0.0])
 
 def solve_for_setpoints(x) -> WindfarmSolution:
-    setpoints = [(x[0], x[1], x[2]), (x[3], x[4], x[5]), (x[6], x[7], x[8])]
+    setpoints = [(x[0], x[1], x[2], 0.0), (x[3], x[4], x[5], 0.0), (x[6], x[7], x[8], 0.0)]
     return windfarm(layout, setpoints)
 
-
 def objective_func(x):
     windfarm_sol = solve_for_setpoints(x)
     return -windfarm_sol.Cp
 
-
 if __name__ == "__main__":
-    x0 = [0, 7, 0, 0, 7, 0, 0, 7, 0]
+    nturbs = 3
+    angle_bounds = (-np.deg2rad(15), np.deg2rad(15))
+    tsr_bounds = (3, 10)
+    x0 = [0, 7, 0] * nturbs
+
     sol = minimize(
         objective_func,
         x0,
-        bounds=[
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (3, 10),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (3, 10),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (3, 10),
-            (-np.deg2rad(15), np.deg2rad(15)),
-        ],
+        bounds=[angle_bounds,tsr_bounds,angle_bounds] * nturbs,
     )
     print(sol)
 
@@ -53,12 +49,17 @@ def objective_func(x):
     fig, axes = plt.subplots(2, 1)
     Plotting.plot_windfarm(windfarm_sol_ref, axes[0])
     Plotting.plot_windfarm(windfarm_sol_opt, axes[1])
+    fig.suptitle("Original and Optimized Wind Farm Layout for Yaw ($z = 0$)")
+    axes[0].set_xlabel("$x/D$")
+    axes[1].set_xlabel("$x/D$")
+    axes[0].set_ylabel("$y/D$")
+    axes[1].set_ylabel("$y/D$")
 
     axes[0].set_title(f"$C_P$: {windfarm_sol_ref.Cp:2.3f}")
     axes[1].set_title(
         f"$C_P$: {windfarm_sol_opt.Cp:2.3f} ({100*(windfarm_sol_opt.Cp/windfarm_sol_ref.Cp - 1):+2.1f}%)"
     )
-
+    plt.tight_layout()
     plt.savefig(
         FIGDIR / "example_07_BEM_yaw_optimisation.png", dpi=300, bbox_inches="tight"
     )

examples/example_09_BEM_tilt_optimization.py

@@ -0,0 +1,68 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from scipy.optimize import minimize
+
+from mitwindfarm import Plotting
+from mitwindfarm.windfarm import Windfarm, WindfarmSolution, Layout
+from mitwindfarm.Rotor import BEM
+from MITRotor.ReferenceTurbines import IEA15MW
+from MITRotor.Momentum import UnifiedMomentum
+# from rich import print
+import numpy as np
+
+FIGDIR = Path(__file__).parent.parent / "fig"
+FIGDIR.mkdir(exist_ok=True, parents=True)
+
+# Note that this setup is meant to mimic the example 7 yaw optimization setup.
+# Thus we are optimizing pitch, tsr, and tilt.
+# It becomes visually apparent that optimizing for yaw vs tilt return the same thing!
+momentum_model= UnifiedMomentum(averaging = "rotor")
+windfarm = Windfarm(
+    rotor_model=BEM(IEA15MW(), momentum_model = momentum_model),
+    TIamb=0.1,
+)
+layout = Layout([0, 12, 24], [0.0, 0.0, 0.0], [2.0, 1.5, 1.0])
+
+def solve_for_setpoints(x) -> WindfarmSolution:
+    setpoints = [(x[0], x[1], 0.0, x[2]), (x[3], x[4], 0.0, x[5]), (x[6], x[7], 0.0, x[8])]
+    return windfarm(layout, setpoints)
+
+def objective_func(x):
+    windfarm_sol = solve_for_setpoints(x)
+    return -windfarm_sol.Cp
+
+if __name__ == "__main__":
+    nturbs = 3
+    # note that we limit to positive tilt so the blades don't hit to shaft
+    half_angle_bounds = (0, np.deg2rad(15))
+    tsr_bounds = (3, 10)
+    x0 = [0, 7, 0] * nturbs
+
+    sol = minimize(
+        objective_func,
+        x0,
+        bounds=[half_angle_bounds,tsr_bounds,half_angle_bounds] * nturbs,
+    )
+    print(sol)
+
+    windfarm_sol_ref = solve_for_setpoints(x0)
+    windfarm_sol_opt = solve_for_setpoints(sol.x)
+    print(windfarm_sol_opt)
+
+    fig, axes = plt.subplots(2, 1)
+    Plotting.plot_windfarm(windfarm_sol_ref, axes[0], y = 0)
+    Plotting.plot_windfarm(windfarm_sol_opt, axes[1], y = 0)
+    fig.suptitle("Original and Optimized Wind Farm Layout for Tilt ($y = 0$)")
+    axes[0].set_xlabel("$x/D$")
+    axes[1].set_xlabel("$x/D$")
+    axes[0].set_ylabel("$y/D$")
+    axes[1].set_ylabel("$y/D$")
+    axes[0].set_title(f"$C_P$: {windfarm_sol_ref.Cp:2.3f}")
+    axes[1].set_title(
+        f"$C_P$: {windfarm_sol_opt.Cp:2.3f} ({100*(windfarm_sol_opt.Cp/windfarm_sol_ref.Cp - 1):+2.1f}%)"
+    )
+
+    plt.savefig(
+        FIGDIR / "example_09_BEM_tilt_optimisation.png", dpi=300, bbox_inches="tight"
+    )

examples/example_10_BEM_yaw_tilt.py

@@ -0,0 +1,73 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from scipy.optimize import minimize
+
+from mitwindfarm import Plotting
+from mitwindfarm.windfarm import Windfarm, WindfarmSolution, Layout
+from mitwindfarm.Rotor import BEM
+from MITRotor.Momentum import UnifiedMomentum
+from MITRotor.ReferenceTurbines import IEA15MW
+from UnifiedMomentumModel.Utilities.Geometry import calc_eff_yaw
+# from rich import print
+import numpy as np
+
+FIGDIR = Path(__file__).parent.parent / "fig"
+FIGDIR.mkdir(exist_ok=True, parents=True)
+
+momentum_model= UnifiedMomentum(averaging = "rotor")
+windfarm = Windfarm(
+    rotor_model=BEM(IEA15MW(), momentum_model = momentum_model),
+    TIamb=0.1,
+)
+layout = Layout([0, 12], [0.0, 0.0], [0.0, 0.0])
+
+def solve_for_setpoints(x) -> WindfarmSolution:
+    setpoints = [(x[0], x[1], x[2], 0.0), (x[3], x[4], x[5], 0.0), (x[6], x[7], x[8], 0.0)]
+    return windfarm(layout, setpoints)
+
+def objective_func(x):
+    windfarm_sol = solve_for_setpoints(x)
+    return -windfarm_sol.Cp
+
+if __name__ == "__main__":
+    pitch_val = 0
+    tsr_val = 7
+
+    yaw_tilt_val = np.deg2rad(15)
+    eff_yaw = calc_eff_yaw(yaw_tilt_val, yaw_tilt_val)
+    yaw_val = eff_yaw
+    tilt_val = eff_yaw
+
+    aligned_setpoints = (pitch_val, tsr_val, 0, 0)
+    setpoints_yaw = [(pitch_val, tsr_val, eff_yaw, 0), aligned_setpoints]
+    setpoints_tilt = [(pitch_val, tsr_val, 0, eff_yaw), aligned_setpoints]
+    setpoints_yaw_and_tilt = [(pitch_val, tsr_val, yaw_tilt_val, yaw_tilt_val), aligned_setpoints]   
+
+    yaw_sol = windfarm(layout, setpoints_yaw)
+    tilt_sol = windfarm(layout, setpoints_tilt) 
+    yaw_and_tilt_sol = windfarm(layout, setpoints_yaw_and_tilt) 
+
+    fig, axes = plt.subplots(2, 3)
+    Plotting.plot_windfarm(yaw_sol, axes[0, 0], z = 0)
+    Plotting.plot_windfarm(tilt_sol, axes[0, 1], z = 0)
+    Plotting.plot_windfarm(yaw_and_tilt_sol, axes[0, 2], z = 0)
+    axes[0, 0].set_ylabel("$z/D$")
+
+    Plotting.plot_windfarm(yaw_sol, axes[1, 0], y = 0)
+    Plotting.plot_windfarm(tilt_sol, axes[1, 1], y = 0)
+    Plotting.plot_windfarm(yaw_and_tilt_sol, axes[1, 2], y = 0)
+    axes[1, 0].set_ylabel("$y/D$")
+    axes[1, 0].set_xlabel("$x/D$")
+    axes[1, 1].set_xlabel("$x/D$")
+    axes[1, 2].set_xlabel("$x/D$")
+
+    fig.suptitle("$C_P$ for Equivalent Turbine Setups", size = 16)
+    deg_yaw_tilt, deg_eff_yaw = np.rad2deg(yaw_tilt_val), np.rad2deg(eff_yaw)
+    axes[0, 0].set_title(f"Yaw {np.round(deg_eff_yaw, decimals=1)}$^\circ$\n$C_P$: {yaw_sol.Cp:2.3f}")
+    axes[0, 1].set_title(f"Tilt {np.round(deg_eff_yaw, decimals=1)}$^\circ$\n$C_P$: {tilt_sol.Cp:2.3f}")
+    axes[0, 2].set_title(f"Yaw {np.round(deg_yaw_tilt, decimals=1)}$^\circ$ & Tilt {np.round(deg_yaw_tilt, decimals=1)}$^\circ$\n$C_P$: {yaw_and_tilt_sol.Cp:2.3f}")
+    plt.tight_layout()
+    plt.savefig(
+        FIGDIR / "example_10_yaw_tilt_comparison.png", dpi=300, bbox_inches="tight"
+    )