From 12ba8ee00dbe1cae3d4b4b4d8c0098f9144156e3 Mon Sep 17 00:00:00 2001
From: ajonesr <105310033+ajonesr@users.noreply.github.com>
Date: Mon, 18 Dec 2023 15:19:06 -0600
Subject: [PATCH] Add functions to fit and convert IAM models (#1827)

* Adding functions, docstrings draft

* Updating imports, adding docstrings

* Updated rst file

* Make linter edits

* Docstring experiment

* Added tests

* More linter edits

* Docstrings and linter edits

* Docstrings and linter

* LINTER

* Docstrings edit

* Added more tests

* Annihilate spaces

* Spacing

* Changed default weight function

* Silence numpy warning

* Updating tests to work with new default

* Forgot a comment

* Return dict contains scalars now, instead of arrays

* Adding option to not fix n

* Adding straggler tests

* Removing examples specific to old default weight function

* Linter nitpicks

* Update docstrings

* Experimenting with example

* Adjusting figure size

* Edit gallery example

* Fixing bounds

* Linter

* Example experimentation

* exact ashrae intercept

* editing docstrings mostly

* whatsnew

* fix errors

* remove test for weight function size

* editing

* simplify weight function

* improve martin_ruiz to physical, generalize tests

* fix examples, split convert and fit examples

* linter, improve coverage

* spacing

* fix reverse order test

* improve examples

* print parameters

* whatsnew

* remove v0.10.2 whatsnew

* Revert "remove v0.10.2 whatsnew"

This reverts commit ed357311b584392714dd8d86b36e64d9efb46349.

* put v0.10.2.rst right again

* require scipy>=1.5.0

* linter

* linter

* suggestions from review

* add reference

* edits to examples

* add note to convert

* edit note on convert

* edit both notes

* polish the notes

* sum not Sum

* edits

* remove test for scipy

* edits from review

* its not it's

* change internal linspace to one degree intervals

* use linspace(0, 90, 91)

---------

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
---
 .../reflections/plot_convert_iam_models.py    | 162 ++++++++
 .../reflections/plot_fit_iam_models.py        |  99 +++++
 .../source/reference/pv_modeling/iam.rst      |   2 +
 docs/sphinx/source/whatsnew/v0.10.3.rst       |   4 +
 pvlib/iam.py                                  | 364 +++++++++++++++++-
 pvlib/tests/conftest.py                       |   1 +
 pvlib/tests/test_iam.py                       | 164 +++++++-
 7 files changed, 791 insertions(+), 5 deletions(-)
 create mode 100644 docs/examples/reflections/plot_convert_iam_models.py
 create mode 100644 docs/examples/reflections/plot_fit_iam_models.py

diff --git a/docs/examples/reflections/plot_convert_iam_models.py b/docs/examples/reflections/plot_convert_iam_models.py
new file mode 100644
index 0000000000..6b0ec78ab3
--- /dev/null
+++ b/docs/examples/reflections/plot_convert_iam_models.py
@@ -0,0 +1,162 @@
+
+"""
+IAM Model Conversion
+====================
+
+Illustrates how to convert from one IAM model to a different model using
+:py:func:`~pvlib.iam.convert`.
+
+"""
+
+# %%
+# An incidence angle modifier (IAM) model quantifies the fraction of direct
+# irradiance that is reflected away from a module's surface. Three popular
+# IAM models are Martin-Ruiz :py:func:`~pvlib.iam.martin_ruiz`, physical
+# :py:func:`~pvlib.iam.physical`, and ASHRAE :py:func:`~pvlib.iam.ashrae`.
+# Each model requires one or more parameters.
+#
+# Here, we show how to use
+# :py:func:`~pvlib.iam.convert` to estimate parameters for a desired target
+# IAM model from a source IAM model. Model conversion uses a weight
+# function that can assign more influence to some AOI values than others.
+# We illustrate how to provide a custom weight function to
+# :py:func:`~pvlib.iam.convert`.
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pvlib.tools import cosd
+from pvlib.iam import (ashrae, martin_ruiz, physical, convert)
+
+# %%
+# Converting from one IAM model to another model
+# ----------------------------------------------
+#
+# Here we'll show how to convert from the Martin-Ruiz model to the
+# physical and the ASHRAE models.
+
+# Compute IAM values using the martin_ruiz model.
+aoi = np.linspace(0, 90, 100)
+martin_ruiz_params = {'a_r': 0.16}
+martin_ruiz_iam = martin_ruiz(aoi, **martin_ruiz_params)
+
+# Get parameters for the physical model and compute IAM using these parameters.
+physical_params = convert('martin_ruiz', martin_ruiz_params, 'physical')
+physical_iam = physical(aoi, **physical_params)
+
+# Get parameters for the ASHRAE model and compute IAM using these parameters.
+ashrae_params = convert('martin_ruiz', martin_ruiz_params, 'ashrae')
+ashrae_iam = ashrae(aoi, **ashrae_params)
+
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 5), sharey=True)
+
+# Plot each model's IAM vs. angle-of-incidence (AOI).
+ax1.plot(aoi, martin_ruiz_iam, label='Martin-Ruiz')
+ax1.plot(aoi, physical_iam, label='physical')
+ax1.set_xlabel('AOI (degrees)')
+ax1.set_title('Convert from Martin-Ruiz to physical')
+ax1.legend()
+
+ax2.plot(aoi, martin_ruiz_iam, label='Martin-Ruiz')
+ax2.plot(aoi, ashrae_iam, label='ASHRAE')
+ax2.set_xlabel('AOI (degrees)')
+ax2.set_title('Convert from Martin-Ruiz to ASHRAE')
+ax2.legend()
+
+ax1.set_ylabel('IAM')
+plt.show()
+
+
+# %%
+# The weight function
+# -------------------
+# :py:func:`pvlib.iam.convert` uses a weight function when computing residuals
+# between the two models. The default weight
+# function is :math:`1 - \sin(aoi)`. We can instead pass a custom weight
+# function to :py:func:`pvlib.iam.convert`.
+#
+# In some cases, the choice of weight function has a minimal effect on the
+# returned model parameters. This is especially true when converting between
+# the Martin-Ruiz and physical models, because the curves described by these
+# models can match quite closely. However, when conversion involves the ASHRAE
+# model, the choice of weight function can have a meaningful effect on the
+# returned parameters for the target model.
+#
+# Here we'll show examples of both of these cases, starting with an example
+# where the choice of weight function does not have much impact. In doing
+# so, we'll show how to pass in a custom weight function of our choice.
+
+# Compute IAM using the Martin-Ruiz model.
+aoi = np.linspace(0, 90, 100)
+martin_ruiz_params = {'a_r': 0.16}
+martin_ruiz_iam = martin_ruiz(aoi, **martin_ruiz_params)
+
+# Get parameters for the physical model ...
+
+# ... using the default weight function.
+physical_params_default = convert('martin_ruiz', martin_ruiz_params,
+                                  'physical')
+physical_iam_default = physical(aoi, **physical_params_default)
+
+
+# ... using a custom weight function. The weight function must take ``aoi``
+# as its argument and return a vector of the same length as ``aoi``.
+def weight_function(aoi):
+    return cosd(aoi)
+
+
+physical_params_custom = convert('martin_ruiz', martin_ruiz_params, 'physical',
+                                 weight=weight_function)
+physical_iam_custom = physical(aoi, **physical_params_custom)
+
+# Plot IAM vs AOI.
+plt.plot(aoi, martin_ruiz_iam, label='Martin-Ruiz')
+plt.plot(aoi, physical_iam_default, label='Default weight function')
+plt.plot(aoi, physical_iam_custom, label='Custom weight function')
+plt.xlabel('AOI (degrees)')
+plt.ylabel('IAM')
+plt.title('Martin-Ruiz to physical')
+plt.legend()
+plt.show()
+
+# %%
+# For this choice of source and target models, the weight function has little
+# effect on the target model's parameters.
+#
+# Now we'll look at an example where the weight function does affect the
+# output.
+
+# Get parameters for the ASHRAE model ...
+
+# ... using the default weight function.
+ashrae_params_default = convert('martin_ruiz', martin_ruiz_params, 'ashrae')
+ashrae_iam_default = ashrae(aoi, **ashrae_params_default)
+
+# ... using the custom weight function
+ashrae_params_custom = convert('martin_ruiz', martin_ruiz_params, 'ashrae',
+                               weight=weight_function)
+ashrae_iam_custom = ashrae(aoi, **ashrae_params_custom)
+
+# Plot IAM vs AOI.
+plt.plot(aoi, martin_ruiz_iam, label='Martin-Ruiz')
+plt.plot(aoi, ashrae_iam_default, label='Default weight function')
+plt.plot(aoi, ashrae_iam_custom, label='Custom weight function')
+plt.xlabel('AOI (degrees)')
+plt.ylabel('IAM')
+plt.title('Martin-Ruiz to ASHRAE')
+plt.legend()
+plt.show()
+
+# %%
+# In this case, each of the two ASHRAE looks quite different.
+# Finding the right weight function and parameters in such cases will require
+# knowing where you want the target model to be more accurate. The default
+# weight function was chosen because it yielded IAM models that produce
+# similar annual insolation for a simulated PV system.
+
+# %%
+# Reference
+# ---------
+# .. [1] Jones, A. R., Hansen, C. W., Anderson, K. S. Parameter estimation
+#    for incidence angle modifier models for photovoltaic modules. Sandia
+#    report SAND2023-13944 (2023).
diff --git a/docs/examples/reflections/plot_fit_iam_models.py b/docs/examples/reflections/plot_fit_iam_models.py
new file mode 100644
index 0000000000..6feb84423a
--- /dev/null
+++ b/docs/examples/reflections/plot_fit_iam_models.py
@@ -0,0 +1,99 @@
+
+"""
+IAM Model Fitting
+================================
+
+Illustrates how to fit an IAM model to data using :py:func:`~pvlib.iam.fit`.
+
+"""
+
+# %%
+# An incidence angle modifier (IAM) model quantifies the fraction of direct
+# irradiance is that is reflected away from a module's surface. Three popular
+# IAM models are Martin-Ruiz :py:func:`~pvlib.iam.martin_ruiz`, physical
+# :py:func:`~pvlib.iam.physical`, and ASHRAE :py:func:`~pvlib.iam.ashrae`.
+# Each model requires one or more parameters.
+#
+# Here, we show how to use
+# :py:func:`~pvlib.iam.fit` to estimate a model's parameters from data.
+#
+# Model fitting require a weight function that can assign
+# more influence to some AOI values than others. We illustrate how to provide
+# a custom weight function to :py:func:`~pvlib.iam.fit`.
+
+import numpy as np
+from random import uniform
+import matplotlib.pyplot as plt
+
+from pvlib.tools import cosd
+from pvlib.iam import (martin_ruiz, physical, fit)
+
+
+# %%
+# Fitting an IAM model to data
+# ----------------------------
+#
+# Here, we'll show how to fit an IAM model to data.
+# We'll generate some data by perturbing output from the Martin-Ruiz model to
+# mimic measured data and then we'll fit the physical model to the perturbed
+# data.
+
+# Create some IAM data.
+aoi = np.linspace(0, 85, 10)
+params = {'a_r': 0.16}
+iam = martin_ruiz(aoi, **params)
+data = iam * np.array([uniform(0.98, 1.02) for _ in range(len(iam))])
+
+# Get parameters for the physical model by fitting to the perturbed data.
+physical_params = fit(aoi, data, 'physical')
+
+# Compute IAM with the fitted physical model parameters.
+physical_iam = physical(aoi, **physical_params)
+
+# Plot IAM vs. AOI
+plt.scatter(aoi, data, c='darkorange', label='Data')
+plt.plot(aoi, physical_iam, label='physical')
+plt.xlabel('AOI (degrees)')
+plt.ylabel('IAM')
+plt.title('Fitting the physical model to data')
+plt.legend()
+plt.show()
+
+
+# %%
+# The weight function
+# -------------------
+# :py:func:`pvlib.iam.fit` uses a weight function when computing residuals
+# between the model and data. The default weight
+# function is :math:`1 - \sin(aoi)`. We can instead pass a custom weight
+# function to :py:func:`pvlib.iam.fit`.
+#
+
+# Define a custom weight function.  The weight function must take ``aoi``
+# as its argument and return a vector of the same length as ``aoi``.
+def weight_function(aoi):
+    return cosd(aoi)
+
+
+physical_params_custom = fit(aoi, data, 'physical', weight=weight_function)
+
+physical_iam_custom = physical(aoi, **physical_params_custom)
+
+# Plot IAM vs AOI.
+fig, ax = plt.subplots(2, 1, figsize=(5, 8))
+ax[0].plot(aoi, data, '.', label='Data (from Martin-Ruiz model)')
+ax[0].plot(aoi, physical_iam, label='With default weight function')
+ax[0].plot(aoi, physical_iam_custom, label='With custom weight function')
+ax[0].set_xlabel('AOI (degrees)')
+ax[0].set_ylabel('IAM')
+ax[0].legend()
+
+ax[1].plot(aoi, physical_iam_custom - physical_iam, label='Custom - default')
+ax[1].set_xlabel('AOI (degrees)')
+ax[1].set_ylabel('Diff. in IAM')
+ax[1].legend()
+plt.tight_layout()
+plt.show()
+
+print("Parameters with default weights: " + str(physical_params))
+print("Parameters with custom weights: " + str(physical_params_custom))
diff --git a/docs/sphinx/source/reference/pv_modeling/iam.rst b/docs/sphinx/source/reference/pv_modeling/iam.rst
index 1871f9b4a2..e12a0c519e 100644
--- a/docs/sphinx/source/reference/pv_modeling/iam.rst
+++ b/docs/sphinx/source/reference/pv_modeling/iam.rst
@@ -17,3 +17,5 @@ Incident angle modifiers
    iam.marion_integrate
    iam.schlick
    iam.schlick_diffuse
+   iam.convert
+   iam.fit
diff --git a/docs/sphinx/source/whatsnew/v0.10.3.rst b/docs/sphinx/source/whatsnew/v0.10.3.rst
index 4d0cc774c8..30fd46648f 100644
--- a/docs/sphinx/source/whatsnew/v0.10.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.10.3.rst
@@ -15,6 +15,8 @@ Enhancements
 * :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa` now include
   shaded fraction in returned variables. (:pull:`1871`)
+* Added :py:func:`~pvlib.iam.convert` and :py:func:`~pvlib.iam.fit` that
+  convert between IAM models, and that fit an IAM model to data. (:issue:`1824`, :pull:`1827`)
 
 Bug fixes
 ~~~~~~~~~
@@ -50,6 +52,8 @@ Contributors
 * Miguel Sánchez de León Peque (:ghuser:`Peque`)
 * Will Hobbs (:ghuser:`williamhobbs`)
 * Anton Driesse (:ghuser:`adriesse`)
+* Abigail Jones (:ghuser:`ajonesr`)
+* Cliff Hansen (:ghuser:`cwhanse`)
 * Gilles Fischer (:ghuser:`GillesFischerV`)
 * Adam R. Jensen (:ghusuer:`AdamRJensen`)
 * :ghuser:`matsuobasho`
diff --git a/pvlib/iam.py b/pvlib/iam.py
index 4f32d352ee..83b8955e2f 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -11,7 +11,8 @@
 import numpy as np
 import pandas as pd
 import functools
-from pvlib.tools import cosd, sind
+from scipy.optimize import minimize
+from pvlib.tools import cosd, sind, acosd
 
 # a dict of required parameter names for each IAM model
 # keys are the function names for the IAM models
@@ -856,7 +857,7 @@ def schlick_diffuse(surface_tilt):
     Unlike the Fresnel reflection factor itself, Schlick's approximation can
     be integrated analytically to derive a closed-form equation for diffuse
     IAM factors for the portions of the sky and ground visible
-    from a tilted surface if isotropic distributions are assumed.  
+    from a tilted surface if isotropic distributions are assumed.
     This function implements the integration of the
     Schlick approximation provided by Xie et al. [2]_.
 
@@ -940,3 +941,362 @@ def schlick_diffuse(surface_tilt):
         cug = pd.Series(cug, surface_tilt.index)
 
     return cuk, cug
+
+
+def _get_model(model_name):
+    # check that model is implemented
+    model_dict = {'ashrae': ashrae, 'martin_ruiz': martin_ruiz,
+                  'physical': physical}
+    try:
+        model = model_dict[model_name]
+    except KeyError:
+        raise NotImplementedError(f"The {model_name} model has not been "
+                                  "implemented")
+
+    return model
+
+
+def _check_params(model_name, params):
+    # check that the parameters passed in with the model
+    # belong to the model
+    exp_params = _IAM_MODEL_PARAMS[model_name]
+    if set(params.keys()) != exp_params:
+        raise ValueError(f"The {model_name} model was expecting to be passed "
+                         "{', '.join(list(exp_params))}, but "
+                         "was handed {', '.join(list(params.keys()))}")
+
+
+def _sin_weight(aoi):
+    return 1 - sind(aoi)
+
+
+def _residual(aoi, source_iam, target, target_params,
+              weight=_sin_weight):
+    # computes a sum of weighted differences between the source model
+    # and target model, using the provided weight function
+
+    weights = weight(aoi)
+
+    # if aoi contains values outside of interval (0, 90), annihilate
+    # the associated weights (we don't want IAM values from AOI outside
+    # of (0, 90) to affect the fit; this is a possible issue when using
+    # `iam.fit`, but not an issue when using `iam.convert`, since in
+    # that case aoi is defined internally)
+    weights = weights * np.logical_and(aoi >= 0, aoi <= 90).astype(int)
+
+    diff = np.abs(source_iam - np.nan_to_num(target(aoi, *target_params)))
+    return np.sum(diff * weights)
+
+
+def _get_ashrae_intercept(b):
+    # find x-intercept of ashrae model
+    return acosd(b / (1 + b))
+
+
+def _ashrae_to_physical(aoi, ashrae_iam, weight, fix_n, b):
+    if fix_n:
+        # the ashrae model has an x-intercept less than 90
+        # we solve for this intercept, and fix n so that the physical
+        # model will have the same x-intercept
+        intercept = _get_ashrae_intercept(b)
+        n = sind(intercept)
+
+        # with n fixed, we will optimize for L (recall that K and L always
+        # appear in the physical model as a product, so it is enough to
+        # optimize for just L, and to fix K=4)
+
+        # we will pass n to the optimizer to simplify things later on,
+        # but because we are setting (n, n) as the bounds, the optimizer
+        # will leave n fixed
+        bounds = [(1e-6, 0.08), (n, n)]
+        guess = [0.002, n]
+
+    else:
+        # we don't fix n, so physical won't have same x-intercept as ashrae
+        # the fit will be worse, but the parameters returned for the physical
+        # model will be more realistic
+        bounds = [(1e-6, 0.08), (0.8, 2)]  # L, n
+        guess = [0.002, 1.0]
+
+    def residual_function(target_params):
+        L, n = target_params
+        return _residual(aoi, ashrae_iam, physical, [n, 4, L], weight)
+
+    return residual_function, guess, bounds
+
+
+def _martin_ruiz_to_physical(aoi, martin_ruiz_iam, weight, a_r):
+    # we will optimize for both n and L (recall that K and L always
+    # appear in the physical model as a product, so it is enough to
+    # optimize for just L, and to fix K=4)
+    # set lower bound for n at 1.0 so that x-intercept will be at 90
+    # order for Powell's method depends on a_r value
+    bounds = [(1e-6, 0.08), (1.05, 2)]  # L, n
+    guess = [0.002, 1.1]  # L, n
+    # get better results if we reverse order to n, L at high a_r
+    if a_r > 0.22:
+        bounds.reverse()
+        guess.reverse()
+
+    # the product of K and L is more important in determining an initial
+    # guess for the location of the minimum, so we pass L in first
+    def residual_function(target_params):
+        # unpack target_params for either search order
+        if target_params[0] < target_params[1]:
+            # L will always be less than n
+            L, n = target_params
+        else:
+            n, L = target_params
+        return _residual(aoi, martin_ruiz_iam, physical, [n, 4, L], weight)
+
+    return residual_function, guess, bounds
+
+
+def _minimize(residual_function, guess, bounds, xtol):
+    if xtol is not None:
+        options = {'xtol': xtol}
+    else:
+        options = None
+    with np.errstate(invalid='ignore'):
+        optimize_result = minimize(residual_function, guess, method="powell",
+                                   bounds=bounds, options=options)
+
+    if not optimize_result.success:
+        try:
+            message = "Optimizer exited unsuccessfully:" \
+                      + optimize_result.message
+        except AttributeError:
+            message = "Optimizer exited unsuccessfully: \
+                       No message explaining the failure was returned. \
+                       If you would like to see this message, please \
+                       update your scipy version (try version 1.8.0 \
+                       or beyond)."
+        raise RuntimeError(message)
+
+    return optimize_result
+
+
+def _process_return(target_name, optimize_result):
+    if target_name == "ashrae":
+        target_params = {'b': optimize_result.x.item()}
+
+    elif target_name == "martin_ruiz":
+        target_params = {'a_r': optimize_result.x.item()}
+
+    elif target_name == "physical":
+        L, n = optimize_result.x
+        # have to unpack order because search order may be different
+        if L > n:
+            L, n = n, L
+        target_params = {'n': n, 'K': 4, 'L': L}
+
+    return target_params
+
+
+def convert(source_name, source_params, target_name, weight=_sin_weight,
+            fix_n=True, xtol=None):
+    """
+    Convert a source IAM model to a target IAM model.
+
+    Parameters
+    ----------
+    source_name : str
+        Name of the source model. Must be ``'ashrae'``, ``'martin_ruiz'``, or
+        ``'physical'``.
+
+    source_params : dict
+        A dictionary of parameters for the source model.
+
+            If source model is ``'ashrae'``, the dictionary must contain
+            the key ``'b'``.
+
+            If source model is ``'martin_ruiz'``, the dictionary must
+            contain the key ``'a_r'``.
+
+            If source model is ``'physical'``, the dictionary must
+            contain the keys ``'n'``, ``'K'``, and ``'L'``.
+
+    target_name : str
+        Name of the target model. Must be ``'ashrae'``, ``'martin_ruiz'``, or
+        ``'physical'``.
+
+    weight : function, optional
+        A single-argument function of AOI (degrees) that calculates weights for
+        the residuals between models. Must return a float or an array-like
+        object. The default weight function is :math:`f(aoi) = 1 - sin(aoi)`.
+
+    fix_n : bool, default True
+        A flag to determine which method is used when converting from the
+        ASHRAE model to the physical model.
+
+        When ``source_name`` is ``'ashrae'`` and ``target_name`` is
+        ``'physical'``, if `fix_n` is ``True``,
+        :py:func:`iam.convert` will fix ``n`` so that the returned physical
+        model has the same x-intercept as the inputted ASHRAE model.
+        Fixing ``n`` like this improves the fit of the conversion, but often
+        returns unrealistic values for the parameters of the physical model.
+        If more physically meaningful parameters are wanted,
+        set `fix_n` to False.
+
+    xtol : float, optional
+        Passed to scipy.optimize.minimize.
+
+    Returns
+    -------
+    dict
+        Parameters for the target model.
+
+            If target model is ``'ashrae'``, the dictionary will contain
+            the key ``'b'``.
+
+            If target model is ``'martin_ruiz'``, the dictionary will
+            contain the key ``'a_r'``.
+
+            If target model is ``'physical'``, the dictionary will
+            contain the keys ``'n'``, ``'K'``, and ``'L'``.
+
+    Note
+    ----
+    Target model parameters are determined by minimizing
+
+    .. math::
+
+        \\sum_{\\theta=0}^{90} weight \\left(\\theta \\right) \\times
+        \\| source \\left(\\theta \\right) - target \\left(\\theta \\right) \\|
+
+    The sum is over :math:`\\theta = 0, 1, 2, ..., 90`.
+
+    References
+    ----------
+    .. [1] Jones, A. R., Hansen, C. W., Anderson, K. S. Parameter estimation
+       for incidence angle modifier models for photovoltaic modules. Sandia
+       report SAND2023-13944 (2023).
+
+    See Also
+    --------
+    pvlib.iam.fit
+    pvlib.iam.ashrae
+    pvlib.iam.martin_ruiz
+    pvlib.iam.physical
+    """
+    source = _get_model(source_name)
+    target = _get_model(target_name)
+
+    aoi = np.linspace(0, 90, 91)
+    _check_params(source_name, source_params)
+    source_iam = source(aoi, **source_params)
+
+    if target_name == "physical":
+        # we can do some special set-up to improve the fit when the
+        # target model is physical
+        if source_name == "ashrae":
+            residual_function, guess, bounds = \
+                _ashrae_to_physical(aoi, source_iam, weight, fix_n,
+                                    source_params['b'])
+        elif source_name == "martin_ruiz":
+            residual_function, guess, bounds = \
+                _martin_ruiz_to_physical(aoi, source_iam, weight,
+                                         source_params['a_r'])
+
+    else:
+        # otherwise, target model is ashrae or martin_ruiz, and scipy
+        # does fine without any special set-up
+        bounds = [(1e-04, 1)]
+        guess = [1e-03]
+
+        def residual_function(target_param):
+            return _residual(aoi, source_iam, target, target_param, weight)
+
+    optimize_result = _minimize(residual_function, guess, bounds,
+                                xtol=xtol)
+
+    return _process_return(target_name, optimize_result)
+
+
+def fit(measured_aoi, measured_iam, model_name, weight=_sin_weight, xtol=None):
+    """
+    Find model parameters that best fit the data.
+
+    Parameters
+    ----------
+    measured_aoi : array-like
+        Angle of incidence values associated with the
+        measured IAM values. [degrees]
+
+    measured_iam : array-like
+        IAM values. [unitless]
+
+    model_name : str
+        Name of the model to be fit. Must be ``'ashrae'``, ``'martin_ruiz'``,
+        or ``'physical'``.
+
+    weight : function, optional
+        A single-argument function of AOI (degrees) that calculates weights for
+        the residuals between models. Must return a float or an array-like
+        object. The default weight function is :math:`f(aoi) = 1 - sin(aoi)`.
+
+    xtol : float, optional
+        Passed to scipy.optimize.minimize.
+
+    Returns
+    -------
+    dict
+        Parameters for target model.
+
+            If target model is ``'ashrae'``, the dictionary will contain
+            the key ``'b'``.
+
+            If target model is ``'martin_ruiz'``, the dictionary will
+            contain the key ``'a_r'``.
+
+            If target model is ``'physical'``, the dictionary will
+            contain the keys ``'n'``, ``'K'``, and ``'L'``.
+
+    References
+    ----------
+    .. [1] Jones, A. R., Hansen, C. W., Anderson, K. S. Parameter estimation
+       for incidence angle modifier models for photovoltaic modules. Sandia
+       report SAND2023-13944 (2023).
+
+    Note
+    ----
+    Model parameters are determined by minimizing
+
+    .. math::
+
+        \\sum_{AOI} weight \\left( AOI \\right) \\times
+        \\| IAM \\left( AOI \\right) - model \\left( AOI \\right) \\|
+
+    The sum is over ``measured_aoi`` and :math:`IAM \\left( AOI \\right)`
+    is ``measured_IAM``.
+
+    See Also
+    --------
+    pvlib.iam.convert
+    pvlib.iam.ashrae
+    pvlib.iam.martin_ruiz
+    pvlib.iam.physical
+    """
+    target = _get_model(model_name)
+
+    if model_name == "physical":
+        bounds = [(0, 0.08), (1, 2)]
+        guess = [0.002, 1+1e-08]
+
+        def residual_function(target_params):
+            L, n = target_params
+            return _residual(measured_aoi, measured_iam, target, [n, 4, L],
+                             weight)
+
+    # otherwise, target_name is martin_ruiz or ashrae
+    else:
+        bounds = [(1e-08, 1)]
+        guess = [0.05]
+
+        def residual_function(target_param):
+            return _residual(measured_aoi, measured_iam, target,
+                             target_param, weight)
+
+    optimize_result = _minimize(residual_function, guess, bounds, xtol)
+
+    return _process_return(model_name, optimize_result)
diff --git a/pvlib/tests/conftest.py b/pvlib/tests/conftest.py
index 15b0cd70e8..f579ef45f2 100644
--- a/pvlib/tests/conftest.py
+++ b/pvlib/tests/conftest.py
@@ -11,6 +11,7 @@
 import pvlib
 from pvlib.location import Location
 
+
 pvlib_base_version = Version(Version(pvlib.__version__).base_version)
 
 
diff --git a/pvlib/tests/test_iam.py b/pvlib/tests/test_iam.py
index 3d017ef7c2..f5ca231bd4 100644
--- a/pvlib/tests/test_iam.py
+++ b/pvlib/tests/test_iam.py
@@ -265,10 +265,10 @@ def test_marion_diffuse_model(mocker):
     assert physical_spy.call_count == 3
 
     for k, v in ashrae_expected.items():
-        np.testing.assert_allclose(ashrae_actual[k], v)
+        assert_allclose(ashrae_actual[k], v)
 
     for k, v in physical_expected.items():
-        np.testing.assert_allclose(physical_actual[k], v)
+        assert_allclose(physical_actual[k], v)
 
 
 def test_marion_diffuse_kwargs():
@@ -281,7 +281,7 @@ def test_marion_diffuse_kwargs():
     actual = _iam.marion_diffuse('ashrae', 20, b=0.04)
 
     for k, v in expected.items():
-        np.testing.assert_allclose(actual[k], v)
+        assert_allclose(actual[k], v)
 
 
 def test_marion_diffuse_invalid():
@@ -395,3 +395,161 @@ def test_schlick_diffuse():
     assert_series_equal(pd.Series(expected_sky, idx), actual_sky)
     assert_series_equal(pd.Series(expected_ground, idx), actual_ground,
                         rtol=1e-6)
+
+
+@pytest.mark.parametrize('source,source_params,target,expected', [
+    ('physical', {'n': 1.5, 'K': 4.5, 'L': 0.004}, 'martin_ruiz',
+     {'a_r': 0.174037}),
+    ('physical', {'n': 1.5, 'K': 4.5, 'L': 0.004}, 'ashrae',
+     {'b': 0.042896}),
+    ('ashrae', {'b': 0.15}, 'physical',
+     {'n': 0.991457, 'K': 4, 'L': 0.037813}),
+    ('ashrae', {'b': 0.15}, 'martin_ruiz', {'a_r': 0.302390}),
+    ('martin_ruiz', {'a_r': 0.15}, 'physical',
+     {'n': 1.240190, 'K': 4, 'L': 0.002791055}),
+    ('martin_ruiz', {'a_r': 0.15}, 'ashrae', {'b': 0.025458})])
+def test_convert(source, source_params, target, expected):
+    target_params = _iam.convert(source, source_params, target)
+    exp = [expected[k] for k in expected]
+    tar = [target_params[k] for k in expected]
+    assert_allclose(exp, tar, rtol=1e-05)
+
+
+@pytest.mark.parametrize('source,source_params', [
+    ('ashrae', {'b': 0.15}),
+    ('ashrae', {'b': 0.05}),
+    ('martin_ruiz', {'a_r': 0.15})])
+def test_convert_recover(source, source_params):
+    # convert isn't set up to handle both source and target = 'physical'
+    target_params = _iam.convert(source, source_params, source, xtol=1e-7)
+    exp = [source_params[k] for k in source_params]
+    tar = [target_params[k] for k in source_params]
+    assert_allclose(exp, tar, rtol=1e-05)
+
+
+def test_convert_ashrae_physical_no_fix_n():
+    # convert ashrae to physical, without fixing n
+    source_params = {'b': 0.15}
+    target_params = _iam.convert('ashrae', source_params, 'physical',
+                                 fix_n=False)
+    expected = {'n': 0.989019, 'K': 4, 'L': 0.037382}
+    exp = [expected[k] for k in expected]
+    tar = [target_params[k] for k in expected]
+    assert_allclose(exp, tar, rtol=1e-05)
+
+
+def test_convert_reverse_order_in_physical():
+    source_params = {'a_r': 0.25}
+    target_params = _iam.convert('martin_ruiz', source_params, 'physical')
+    expected = {'n': 1.691398, 'K': 4, 'L': 0.071633}
+    exp = [expected[k] for k in expected]
+    tar = [target_params[k] for k in expected]
+    assert_allclose(exp, tar, rtol=1e-05)
+
+
+def test_convert_xtol():
+    source_params = {'b': 0.15}
+    target_params = _iam.convert('ashrae', source_params, 'physical',
+                                 xtol=1e-8)
+    expected = {'n': 0.9914568914, 'K': 4, 'L': 0.0378126985}
+    exp = [expected[k] for k in expected]
+    tar = [target_params[k] for k in expected]
+    assert_allclose(exp, tar, rtol=1e-6)
+
+
+def test_convert_custom_weight_func():
+    aoi = np.linspace(0, 90, 91)
+
+    # convert physical to martin_ruiz, using custom weight function
+    source_params = {'n': 1.5, 'K': 4.5, 'L': 0.004}
+    source_iam = _iam.physical(aoi, **source_params)
+
+    # define custom weight function that takes in other arguments
+    def scaled_weight(aoi):
+        return 2. * aoi
+
+    # expected value calculated from computing residual function over
+    # a range of inputs, and taking minimum of these values
+    expected_min_res = 16.39724
+
+    actual_dict = _iam.convert('physical', source_params, 'martin_ruiz',
+                               weight=scaled_weight)
+    actual_min_res = _iam._residual(aoi, source_iam, _iam.martin_ruiz,
+                                    [actual_dict['a_r']], scaled_weight)
+
+    assert np.isclose(expected_min_res, actual_min_res, atol=1e-06)
+
+
+def test_convert_model_not_implemented():
+    with pytest.raises(NotImplementedError, match='model has not been'):
+        _iam.convert('ashrae', {'b': 0.1}, 'foo')
+
+
+def test_convert_wrong_model_parameters():
+    with pytest.raises(ValueError, match='model was expecting'):
+        _iam.convert('ashrae', {'B': 0.1}, 'physical')
+
+
+def test_convert__minimize_fails():
+    # to make scipy.optimize.minimize fail, we'll pass in a nonsense
+    # weight function that only outputs nans
+    def nan_weight(aoi):
+        return np.nan
+
+    with pytest.raises(RuntimeError, match='Optimizer exited unsuccessfully'):
+        _iam.convert('ashrae', {'b': 0.1}, 'physical', weight=nan_weight)
+
+
+def test_fit():
+    aoi = np.linspace(0, 90, 5)
+    perturb = np.array([1.2, 1.01, 0.95, 1, 0.98])
+    perturbed_iam = _iam.martin_ruiz(aoi, a_r=0.14) * perturb
+
+    expected_a_r = 0.14
+
+    actual_dict = _iam.fit(aoi, perturbed_iam, 'martin_ruiz')
+    actual_a_r = actual_dict['a_r']
+
+    assert np.isclose(expected_a_r, actual_a_r, atol=1e-04)
+
+
+def test_fit_custom_weight_func():
+    # define custom weight function that takes in other arguments
+    def scaled_weight(aoi):
+        return 2. * aoi
+
+    aoi = np.linspace(0, 90, 5)
+    perturb = np.array([1.2, 1.01, 0.95, 1, 0.98])
+    perturbed_iam = _iam.martin_ruiz(aoi, a_r=0.14) * perturb
+
+    expected_a_r = 0.14
+
+    actual_dict = _iam.fit(aoi, perturbed_iam, 'martin_ruiz',
+                           weight=scaled_weight)
+    actual_a_r = actual_dict['a_r']
+
+    assert np.isclose(expected_a_r, actual_a_r, atol=1e-04)
+
+
+def test_fit_model_not_implemented():
+    with pytest.raises(NotImplementedError, match='model has not been'):
+        _iam.fit(np.array([0, 10]), np.array([1, 0.99]), 'foo')
+
+
+def test_fit__minimize_fails():
+    # to make scipy.optimize.minimize fail, we'll pass in a nonsense
+    # weight function that only outputs nans
+    def nan_weight(aoi):
+        return np.nan
+
+    with pytest.raises(RuntimeError, match='Optimizer exited unsuccessfully'):
+        _iam.fit(np.array([0, 10]), np.array([1, 0.99]), 'physical',
+                 weight=nan_weight)
+
+
+def test__residual_zero_outside_range():
+    # check that _residual annihilates any weights that come from aoi
+    # outside of interval [0, 90] (this is important for `iam.fit`, when
+    # the `measured_aoi` contains angles outside this range
+    residual = _iam._residual(101, _iam.ashrae(101), _iam.martin_ruiz, [0.16])
+    assert residual == 0.0