Repository: NREL/floris

Issue:
Bug report?: optimize yaw with wind direction uncertainty does not work in with ParFlorisModel
<!--
IMPORTANT NOTES

Thank you for taking the time to report a bug. If you aren't certain whether an issue
is a bug, please first open a Discussion. Before submitting, please reread your
description to ensure that other readers can reasonably understand the issue
you're facing and the impact on your workflow or results.

This form is written in GitHub's Markdown format. For a reference on this type
of syntax, see GitHub's documentation:
https://docs.github.com/en/get-started/writing-on-github/getting-started-with-writing-and-formatting-on-github/basic-writing-and-formatting-syntax

This template contains guidance for your submission within the < ! - -, - - > blocks.
These are comments in HTML syntax and will not appear in the submission.
Be sure to use the "Preview" feature on GitHub to ensure your submission is formatted as intended.

When including code snippets, please paste the text itself and wrap the code block with
ticks (see the other character on the tilde ~ key in a US keyboard) to format it as code.
For example, Python code should be wrapped in ticks like this:
```python
def a_func():
    return 1

a = 1
b = a_func()
print(a + b)
```
This is preferred over screen shots since it is searchable and others can copy/paste
the text to run it.
-->

# optimize yaw with wind direction uncertainty does not work in with ParFlorisModel
<!--
A clear and concise description of the bug including what happened
and what you expected to happen.
-->

It seems that it is not possible to run a parallel yaw optimization with wind direction uncertainty.

This is unexpected to me, as the UncertainFlorisModel allows for parallelization (`floris/examples/examples_uncertain/003_uncertain_model_with_parallelization.py`). 

I put a question mark into the title, as I'm not sure it is indeed a bug or just a missing feature (in that case I think a warning should be displayed for the given example below).

## How to reproduce
<!--
Describe how another person with no context can recreate this issue.
It is typically very helpful to reduce the problem as much as possible
and share a minimum working example. See the note above on including
code as text rather than screenshots.
-->


```
from time import perf_counter as timerpc

import numpy as np

from floris import (
    ParFlorisModel,
    TimeSeries,
    WindRose,
)
from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR


# When using parallel optimization it is important the "root" script include this
# if __name__ == "__main__": block to avoid problems
if __name__ == "__main__":

    # Load the wind rose from csv
    wind_rose = WindRose.read_csv_long(
        "../inputs/wind_rose.csv", wd_col="wd", ws_col="ws", freq_col="freq_val",
        ti_col_or_value=0.06
    )

    # Load FLORIS as a parallel model
    max_workers = 16
    pfmodel = ParFlorisModel(
        "../inputs/gch.yaml",
        max_workers=max_workers,
        n_wind_condition_splits=max_workers,
        interface="pathos",
        print_timings=True,
    )

    # transform the model into the uncertain model
    from floris import UncertainFlorisModel
    pfmodel = UncertainFlorisModel(pfmodel)


    # Specify wind farm layout and update in the floris object
    N = 2  # number of turbines per row and per column
    X, Y = np.meshgrid(
        5.0 * pfmodel.core.farm.rotor_diameters_sorted[0][0] * np.arange(0, N, 1),
        5.0 * pfmodel.core.farm.rotor_diameters_sorted[0][0] * np.arange(0, N, 1),
    )
    pfmodel.set(layout_x=X.flatten(), layout_y=Y.flatten())

    # Get the number of turbines
    n_turbines = len(pfmodel.layout_x)

    # Optimize the yaw angles.  This could be done for every wind direction and wind speed
    # but in practice it is much faster to optimize only for one speed and infer the rest
    # using a rule of thumb
    time_series = TimeSeries(
        wind_directions=wind_rose.wind_directions, wind_speeds=8.0, turbulence_intensities=0.06
    )
    pfmodel.set(wind_data=time_series)

    start_time = timerpc()
    yaw_opt = YawOptimizationSR(
        fmodel=pfmodel,
        minimum_yaw_angle=0.0,  # Allowable yaw angles lower bound
        maximum_yaw_angle=20.0,  # Allowable yaw angles upper bound
        Ny_passes=[5, 4],
        exclude_downstream_turbines=True,
    )
    df_opt = yaw_opt.optimize()
    end_time = timerpc()
    t_tot = end_time - start_time
    print("Optimization finished in {:.2f} seconds.".format(t_tot))
```


The previous snipped is identical to `floris/examples/examples_control_optimization/005_optimize_yaw_aep_parallel.py` with exception of three added lines:
```
# transform the model into the uncertain model
from floris import UncertainFlorisModel
pfmodel = UncertainFlorisModel(pfmodel)
```

## Relevant output
<!--
Include any output, plots, or other means of communication here to add context to the problem.
 -->

Without the added lines, the computation is parallelized as indicated by the output line:
`` 
"Total time spent for parallel calculation (16 workers): ..."
`` 

With the added lines the computations are not parallelized anymore, and the output line doesn't show up. 


## Floris version
<!--
Share your floris version and how you installed it. You can print the floris version from
a Python REPL or script with these commands:
```python
import floris
print(floris.__version__)
```
 --> 4.4.1

## System Information

<!-- Add your information here. -->
 - OS: Windows 11 Pro
 - Python version: Python 3.12.7
 - Library versions
```
attrs==25.3.0
coloredlogs==15.0.1
contourpy==1.3.2
cycler==0.12.1
dill==0.4.0
-e git+https://github.com/NREL/floris.git@0a250744b2bbf5d3163949c193f53efdedde1085#egg=floris
fonttools==4.57.0
humanfriendly==10.0
kiwisolver==1.4.8
matplotlib==3.10.1
multiprocess==0.70.18
numexpr==2.10.2
numpy==2.2.5
packaging==25.0
pandas==2.2.3
pathos==0.3.4
pillow==11.2.1
pox==0.3.6
ppft==1.7.7
pyparsing==3.2.3
pyreadline3==3.5.4
python-dateutil==2.9.0.post0
pytz==2025.2
PyYAML==6.0.2
scipy==1.15.2
shapely==2.1.0
six==1.17.0
tzdata==2025.2
```

Discussion:
@misi9170: Hi @gogannes ,

I believe this is indeed a bug. I'm going to dig into this some more and hope to get back to you soon.

Misha@misi9170: Hi @gogannes , thank you for finding and reporting this bug. It looks like this is in fact a problem not just with yaw optimization, but with using the `UncertainFlorisModel` with an underlying `ParFlorisModel` generally. I have now created a bug fix over on #1107 if you'd like to check it out, and I hope to get that merged into develop and then released onto main in the coming days. Again, thank you for bringing this to our attention!

As you suggest, it could be considered a "missing feature", in that the results are not affected. However, our intention was for uncertain + parallel to work, and I do consider it a bug that that is not enabled.

Misha@gogannes: Dear Misha, 

thank you so much for looking into this immediatley and even fixing the issue!
I tested the bug fix and it seemst work for me. 

While testing I found out that there's another likely related issue:
Using the `set_param()` method for an `UncertainFlorisModel` object, seems to break the parallel computing (in a sense that it is not parallelized anymore).

I think you should be able to reproduce this in the original example above by replacing the three lines with 
```
# transform the model into the uncertain model
from floris import UncertainFlorisModel
pfmodel = UncertainFlorisModel(pfmodel)
pfmodel.set_param(["solver", "turbine_grid_points"], 3)
```

Note that it doesn't matter whether the `set_param()` method is calld on the  `UncertainFlorisModel` or `ParFlorisModel`, in both cases it breaks the parallel computing.

I guess this is not intended and maybe it can be fix together with the original bug (unfortunatley I don't have the time right now to look into this deeper myself, therefore I can just report).

Best and thanks again,
Johannes@misi9170: Thanks Johannes! I'll look into it.

Update: I see what's happening here. When `set_param()` is called, the `print_timings` flag is not copied over (as well as some other important `ParFlorisModel` instantiation arguments, such as the `interface` and `max_workers`). So, the underlying model after `set_param()` does remain a `ParFlorisModel` and will run in parallel, but won't retain the other parallel settings (and so didn't produce the expected printouts). I'll fix this over on #1107.@gogannes: Hey Misha,

I tested the bugfix (https://github.com/NREL/floris/pull/1107) again and now it works as expected :-)

Thank you so much!

Best,
Johannes@misi9170: Hi Johannes, glad to hear it! I'll try to get #1107 merged in soon and make a patch release for FLORIS including the bug fix.

Solution patch:
diff --git a/examples/examples_uncertain/003_uncertain_model_with_parallelization.py b/examples/examples_uncertain/003_uncertain_model_with_parallelization.py
index 37f785c77..64f48a92a 100644
--- a/examples/examples_uncertain/003_uncertain_model_with_parallelization.py
+++ b/examples/examples_uncertain/003_uncertain_model_with_parallelization.py
@@ -14,72 +14,73 @@
 from floris.par_floris_model import ParFlorisModel
 
 
-# Following the refactoring of ParFlorisModel, the UncertainFlorisModel can be
-# parallelized by passing the ParFlorisModel as the model to be run.  This example
-# demonstrates the usage and shows that the result obtained from the UncertainFlorisModel
-# with and without parallelization is the same.  The results are compared to the nominal
-# results.
-
-# Instantiate a FlorisModel and ParallelFlorisModel using the GCH model
-fmodel = FlorisModel("../inputs/gch.yaml")
-pfmodel = ParFlorisModel("../inputs/gch.yaml")
-
-# Use the above model to declare a serial and parallel UncertainFlorisModel
-ufmodel = UncertainFlorisModel(fmodel)
-pufmodel = UncertainFlorisModel(pfmodel)
-
-
-# Define an inflow where wind direction is swept while
-# wind speed and turbulence intensity are held constant
-wind_directions = np.arange(240.0, 300.0, 1.0)
-time_series = TimeSeries(
-    wind_directions=wind_directions,
-    wind_speeds=8.0,
-    turbulence_intensities=0.06,
-)
-
-# Define a two turbine farm and apply the inflow
-D = 126.0
-layout_x = np.array([0, D * 6])
-layout_y = [0, 0]
-
-# Apply to fmodel, ufmodel, and pufmodel
-fmodel.set(
-    layout_x=layout_x,
-    layout_y=layout_y,
-    wind_data=time_series,
-)
-
-ufmodel.set(
-    layout_x=layout_x,
-    layout_y=layout_y,
-    wind_data=time_series,
-)
-
-pufmodel.set(
-    layout_x=layout_x,
-    layout_y=layout_y,
-    wind_data=time_series,
-)
-
-# Run the models
-fmodel.run()
-ufmodel.run()
-pufmodel.run()
-
-# Collect the farm power results from each model
-farm_powers_nom = fmodel.get_farm_power() / 1e3
-farm_powers_unc = ufmodel.get_farm_power() / 1e3
-farm_powers_punc = pufmodel.get_farm_power() / 1e3
-
-# Compare the results
-fig, ax = plt.subplots()
-ax.plot(wind_directions, farm_powers_nom.flatten(), 'k-', label="Nominal power")
-ax.plot(wind_directions, farm_powers_unc.flatten(), 'bs-', label="Uncertain power")
-ax.plot(wind_directions, farm_powers_punc.flatten(), 'r.--', label="Parallel uncertain power")
-ax.grid(True)
-ax.legend()
-ax.set_xlabel("Wind Direction (deg)")
-ax.set_ylabel("Power (kW)")
-
-plt.show()
+if __name__ == "__main__":
+    # Following the refactoring of ParFlorisModel, the UncertainFlorisModel can be
+    # parallelized by passing the ParFlorisModel as the model to be run.  This example
+    # demonstrates the usage and shows that the result obtained from the UncertainFlorisModel
+    # with and without parallelization is the same.  The results are compared to the nominal
+    # results.
+
+    # Instantiate a FlorisModel and ParallelFlorisModel using the GCH model
+    fmodel = FlorisModel("../inputs/gch.yaml")
+    pfmodel = ParFlorisModel("../inputs/gch.yaml")
+
+    # Use the above model to declare a serial and parallel UncertainFlorisModel
+    ufmodel = UncertainFlorisModel(fmodel)
+    pufmodel = UncertainFlorisModel(pfmodel)
+
+
+    # Define an inflow where wind direction is swept while
+    # wind speed and turbulence intensity are held constant
+    wind_directions = np.arange(240.0, 300.0, 1.0)
+    time_series = TimeSeries(
+        wind_directions=wind_directions,
+        wind_speeds=8.0,
+        turbulence_intensities=0.06,
+    )
+
+    # Define a two turbine farm and apply the inflow
+    D = 126.0
+    layout_x = np.array([0, D * 6])
+    layout_y = [0, 0]
+
+    # Apply to fmodel, ufmodel, and pufmodel
+    fmodel.set(
+        layout_x=layout_x,
+        layout_y=layout_y,
+        wind_data=time_series,
+    )
+
+    ufmodel.set(
+        layout_x=layout_x,
+        layout_y=layout_y,
+        wind_data=time_series,
+    )
+
+    pufmodel.set(
+        layout_x=layout_x,
+        layout_y=layout_y,
+        wind_data=time_series,
+    )
+
+    # Run the models
+    fmodel.run()
+    ufmodel.run()
+    pufmodel.run()
+
+    # Collect the farm power results from each model
+    farm_powers_nom = fmodel.get_farm_power() / 1e3
+    farm_powers_unc = ufmodel.get_farm_power() / 1e3
+    farm_powers_punc = pufmodel.get_farm_power() / 1e3
+
+    # Compare the results
+    fig, ax = plt.subplots()
+    ax.plot(wind_directions, farm_powers_nom.flatten(), 'k-', label="Nominal power")
+    ax.plot(wind_directions, farm_powers_unc.flatten(), 'bs-', label="Uncertain power")
+    ax.plot(wind_directions, farm_powers_punc.flatten(), 'r.--', label="Parallel uncertain power")
+    ax.grid(True)
+    ax.legend()
+    ax.set_xlabel("Wind Direction (deg)")
+    ax.set_ylabel("Power (kW)")
+
+    plt.show()
diff --git a/floris/floris_model.py b/floris/floris_model.py
index 4bb926688..3547afe0f 100644
--- a/floris/floris_model.py
+++ b/floris/floris_model.py
@@ -1576,8 +1576,13 @@ def set_operation_model(self, operation_model: str | List[str]):
         )
 
     def copy(self):
-        """Create an independent copy of the current FlorisModel object"""
-        return FlorisModel(self.core.as_dict())
+        """Create an independent copy of the current FlorisModel object
+
+        When creating the copy, this method uses self.__class__(), rather than FlorisModel()
+        directly, so that subclasses of FlorisModel can inherit this method and return
+        instantiations of their own class, rather than the FlorisModel class.
+        """
+        return self.__class__(self.core.as_dict(), **self.secondary_init_kwargs)
 
     def get_param(
         self,
@@ -1616,7 +1621,7 @@ def set_param(
         """
         fm_dict_mod = self.core.as_dict()
         nested_set(fm_dict_mod, param, value, param_idx)
-        self.__init__(fm_dict_mod)
+        self.__init__(fm_dict_mod, **self.secondary_init_kwargs)
 
     def get_turbine_layout(self, z=False):
         """
@@ -1662,6 +1667,14 @@ def print_dict(self) -> None:
 
     ### Properties
 
+    @property
+    def secondary_init_kwargs(self):
+        """
+        FlorisModel takes only the configuration argument. This method is a placeholder for
+        subclasses of FlorisModel.
+        """
+        return {}
+
     @property
     def layout_x(self):
         """
diff --git a/floris/optimization/yaw_optimization/yaw_optimization_base.py b/floris/optimization/yaw_optimization/yaw_optimization_base.py
index 5dab7fe30..4b5c6a927 100644
--- a/floris/optimization/yaw_optimization/yaw_optimization_base.py
+++ b/floris/optimization/yaw_optimization/yaw_optimization_base.py
@@ -231,7 +231,10 @@ def _reduce_control_problem(self):
 
         # Initialize subset variables as full set
         self.fmodel_subset = copy.deepcopy(self.fmodel)
-        self.fmodel_subset._wind_data = None # Accessing private attribute!
+        if hasattr(self.fmodel_subset, "_wind_data"): # Normal FlorisModel, ParFlorisModel
+            self.fmodel_subset._wind_data = None # Accessing private attribute!
+        elif hasattr(self.fmodel_subset, "fmodel_expanded"): # UncertainFlorisModel
+            self.fmodel_subset.fmodel_unexpanded._wind_data = None # Accessing private attribute!
         n_findex_subset = copy.deepcopy(self.fmodel.core.flow_field.n_findex)
         minimum_yaw_angle_subset = copy.deepcopy(self.minimum_yaw_angle)
         maximum_yaw_angle_subset = copy.deepcopy(self.maximum_yaw_angle)
diff --git a/floris/par_floris_model.py b/floris/par_floris_model.py
index 8c3e5dce8..1efa4263f 100644
--- a/floris/par_floris_model.py
+++ b/floris/par_floris_model.py
@@ -350,6 +350,19 @@ def _get_turbine_powers(self):
         else:
             return super()._get_turbine_powers()
 
+    @property
+    def secondary_init_kwargs(self):
+        """
+        ParFlorisModel secondary keyword arguments (after configuration).
+        """
+        return {
+            "interface": self.interface,
+            "max_workers": self.max_workers,
+            "n_wind_condition_splits": self.n_wind_condition_splits,
+            "return_turbine_powers_only": self.return_turbine_powers_only,
+            "print_timings": self.print_timings
+        }
+
     @property
     def fmodel(self):
         """
diff --git a/floris/uncertain_floris_model.py b/floris/uncertain_floris_model.py
index ab5b404ad..5fcc03f0b 100644
--- a/floris/uncertain_floris_model.py
+++ b/floris/uncertain_floris_model.py
@@ -109,6 +109,15 @@ def __init__(
         # Instantiate the un-expanded FlorisModel
         if isinstance(configuration, (FlorisModel, ParFlorisModel)):
             self.fmodel_unexpanded = configuration.copy()
+            # Copy over any control setpoints, wind data, if not already done.
+            self.fmodel_unexpanded.set(
+                yaw_angles=configuration.core.farm.yaw_angles,
+                power_setpoints=configuration.core.farm.power_setpoints,
+                awc_modes=configuration.core.farm.awc_modes,
+                awc_amplitudes=configuration.core.farm.awc_amplitudes,
+                awc_frequencies=configuration.core.farm.awc_frequencies,
+                wind_data=configuration.wind_data,
+            )
         elif isinstance(configuration, (dict, str, Path)):
             self.fmodel_unexpanded = FlorisModel(configuration)
         else:
@@ -985,20 +994,14 @@ def set_operation_model(self, operation_model: str | List[str]):
         )
 
     def copy(self):
-        """Create an independent copy of the current UncertainFlorisModel object"""
-        return UncertainFlorisModel(
-            self.fmodel_unexpanded.core.as_dict(),
-            wd_resolution=self.wd_resolution,
-            ws_resolution=self.ws_resolution,
-            ti_resolution=self.ti_resolution,
-            yaw_resolution=self.yaw_resolution,
-            power_setpoint_resolution=self.power_setpoint_resolution,
-            awc_amplitude_resolution=self.awc_amplitude_resolution,
-            wd_std=self.wd_std,
-            wd_sample_points=self.wd_sample_points,
-            fix_yaw_to_nominal_direction=self.fix_yaw_to_nominal_direction,
-            verbose=self.verbose,
-        )
+        """Create an independent copy of the current UncertainFlorisModel object
+
+        When creating the copy, this method uses self.__class__(), rather than
+        UncertainFlorisModel() directly, so that subclasses of UncertainFlorisModel can inherit this
+        method and return instantiations of their own class, rather than the UncertainFlorisModel
+        class.
+        """
+        return self.__class__(self.fmodel_unexpanded.copy(), **self.secondary_init_kwargs)
 
     def get_param(self, param: List[str], param_idx: Optional[int] = None) -> Any:
         """Get a parameter from a FlorisModel object.
@@ -1028,9 +1031,27 @@ def set_param(self, param: List[str], value: Any, param_idx: Optional[int] = Non
         """
         fm_dict_mod = self.fmodel_unexpanded.core.as_dict()
         nested_set(fm_dict_mod, param, value, param_idx)
-        self.fmodel_unexpanded.__init__(fm_dict_mod)
+        self.fmodel_unexpanded.__init__(fm_dict_mod, **self.fmodel_unexpanded.secondary_init_kwargs)
         self.set()
 
+    @property
+    def secondary_init_kwargs(self):
+        """
+        UncertainFlorisModel secondary keyword arguments (after configuration).
+        """
+        return {
+            "wd_resolution": self.wd_resolution,
+            "ws_resolution": self.ws_resolution,
+            "ti_resolution": self.ti_resolution,
+            "yaw_resolution": self.yaw_resolution,
+            "power_setpoint_resolution": self.power_setpoint_resolution,
+            "awc_amplitude_resolution": self.awc_amplitude_resolution,
+            "wd_std": self.wd_std,
+            "wd_sample_points": self.wd_sample_points,
+            "fix_yaw_to_nominal_direction": self.fix_yaw_to_nominal_direction,
+            "verbose": self.verbose,
+        }
+
     @property
     def layout_x(self):
         """
@@ -1212,3 +1233,18 @@ def __init__(
         self.yaw_resolution = yaw_resolution
         self.power_setpoint_resolution = power_setpoint_resolution
         self.awc_amplitude_resolution = awc_amplitude_resolution
+
+    @property
+    def secondary_init_kwargs(self):
+        """
+        ApproxFlorisModel secondary keyword arguments (after configuration).
+        """
+        return {
+            "wd_resolution": self.wd_resolution,
+            "ws_resolution": self.ws_resolution,
+            "ti_resolution": self.ti_resolution,
+            "yaw_resolution": self.yaw_resolution,
+            "power_setpoint_resolution": self.power_setpoint_resolution,
+            "awc_amplitude_resolution": self.awc_amplitude_resolution,
+            "verbose": self.verbose,
+        }


Test patch:
diff --git a/tests/par_floris_model_unit_test.py b/tests/par_floris_model_unit_test.py
index 9f5669a3c..0c4d80b07 100644
--- a/tests/par_floris_model_unit_test.py
+++ b/tests/par_floris_model_unit_test.py
@@ -312,3 +312,16 @@ def test_sample_flow_at_points(sample_inputs_fixture):
         pfmodel = ParFlorisModel(fmodel, max_workers=2, interface=interface)
         ws_test = pfmodel.sample_flow_at_points(x_test, y_test, z_test)
         assert np.allclose(ws_base, ws_test)
+
+def test_copy(sample_inputs_fixture):
+    """
+    Check that the ParFlorisModel copies correctly as a ParFlorisModel.
+    """
+
+    sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL
+    sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL
+
+    pfmodel = ParFlorisModel(sample_inputs_fixture.core, max_workers=2)
+    pfmodel_copy = pfmodel.copy()
+    assert isinstance(pfmodel_copy, ParFlorisModel)
+    assert pfmodel_copy.max_workers == 2
diff --git a/tests/uncertain_floris_model_integration_test.py b/tests/uncertain_floris_model_integration_test.py
index 02859e611..c471758fb 100644
--- a/tests/uncertain_floris_model_integration_test.py
+++ b/tests/uncertain_floris_model_integration_test.py
@@ -271,6 +271,32 @@ def test_get_powers_with_wind_data():
 
     assert np.allclose(farm_power_weighted, ufmodel.get_turbine_powers()[:, :, :-1].sum(axis=2))
 
+def test_AEP_with_wind_data():
+    wind_speeds = np.array([8.0, 10.0])
+    wind_directions = np.array([270.0, 280.0])
+    frequencies = np.array([[0.25, 0.25], [0.1, 0.4]])
+    wind_rose = WindRose(
+        wind_directions=np.unique(wind_directions),
+        wind_speeds=np.unique(wind_speeds),
+        freq_table=frequencies,
+        ti_table=0.06,
+    )
+
+    # Set wind_data on UncertainFlorisModel directly
+    ufmodel = UncertainFlorisModel(configuration=YAML_INPUT)
+    ufmodel.set(wind_data=wind_rose)
+    ufmodel.run()
+    aep_1 = ufmodel.get_farm_AEP()
+
+    # Set wind_data on FlorisModel and then set on UncertainFlorisModel
+    fmodel = FlorisModel(configuration=YAML_INPUT)
+    fmodel.set(wind_data=wind_rose)
+    ufmodel = UncertainFlorisModel(fmodel)
+    ufmodel.run()
+    aep_2 = ufmodel.get_farm_AEP()
+
+    # Check AEPs match
+    assert np.allclose(aep_1, aep_2)
 
 def test_approx_floris_model():
     afmodel = ApproxFlorisModel(configuration=YAML_INPUT, wd_resolution=1.0)
@@ -305,6 +331,9 @@ def test_approx_floris_model():
     np.testing.assert_almost_equal(power[0], power[1])
     assert not np.allclose(power[2], power[3])
 
+    # Test copy method
+    assert isinstance(afmodel.copy(), ApproxFlorisModel)
+
 
 def test_expected_farm_power_regression():
     ufmodel = UncertainFlorisModel(
@@ -466,6 +495,8 @@ def test_parallel_uncertain_model():
 
     ufmodel = UncertainFlorisModel(FlorisModel(configuration=YAML_INPUT))
     pufmodel = UncertainFlorisModel(ParFlorisModel(configuration=YAML_INPUT))
+    assert isinstance(ufmodel.fmodel_expanded, FlorisModel)
+    assert isinstance(pufmodel.fmodel_expanded, ParFlorisModel)
 
     # Run the models and compare outputs
     ufmodel.run()
@@ -474,3 +505,19 @@ def test_parallel_uncertain_model():
     powers_punc = pufmodel.get_turbine_powers()
 
     assert np.allclose(powers_unc, powers_punc)
+
+def test_copy():
+    """
+    Check that the UncertainFlorisModel copy method works as expected for
+    both FlorisModel and ParFlorisModel.
+    """
+    ufmodel = UncertainFlorisModel(FlorisModel(configuration=YAML_INPUT))
+
+    ufmodel_copy = ufmodel.copy()
+    assert isinstance(ufmodel_copy, UncertainFlorisModel)
+    assert isinstance(ufmodel_copy.fmodel_expanded, FlorisModel)
+
+    pufmodel = UncertainFlorisModel(ParFlorisModel(configuration=YAML_INPUT))
+    pufmodel_copy = pufmodel.copy()
+    assert isinstance(pufmodel_copy, UncertainFlorisModel)
+    assert isinstance(pufmodel_copy.fmodel_expanded, ParFlorisModel)