Add function to fit Sandia inverter model (#1011)

* WIP: fitting inverter model

* correct calculated data, improve test

* move code to pvlib modules, rename data file

* fix test, stickler

* sort out np.polyfit vs. np.polynomial.polynomial.polyfit

* correct file name

* api, whatsnew

* docstring improvements

* revisions from review

* correct use of np.roots

* use pytest.approx instead of converting dicts to np.array

* correct kwarg for pytest.approx

* back to solve_quad instead of np.roots, lint

* change to Series inputs

* correct test, Series to array_like

Author: cwhanse

docs/sphinx/source/api.rst

@@ -265,6 +265,15 @@ Low-level functions for solving the single diode equation.
    singlediode.bishop88_v_from_i
    singlediode.bishop88_mpp
 
+Functions for fitting diode models
+
+.. autosummary::
+   :toctree: generated/
+
+    ivtools.fit_sde_sandia
+    ivtools.fit_sdm_cec_sam
+    ivtools.fit_sdm_desoto
+
 Inverter models (DC to AC conversion)
 -------------------------------------
 
@@ -275,6 +284,14 @@ Inverter models (DC to AC conversion)
    inverter.adr
    inverter.pvwatts
 
+Functions for fitting inverter models
+
+.. autosummary::
+   :toctree: generated/
+
+   inverter.fit_sandia
+
+
 PV System Models
 ----------------
 
@@ -311,16 +328,6 @@ PVWatts model
    inverter.pvwatts
    pvsystem.pvwatts_losses
 
-Functions for fitting diode models
-----------------------------------
-
-.. autosummary::
-   :toctree: generated/
-
-    ivtools.fit_sde_sandia
-    ivtools.fit_sdm_cec_sam
-    ivtools.fit_sdm_desoto
-
 Other
 -----
 

docs/sphinx/source/whatsnew/v0.8.0.rst

@@ -38,6 +38,8 @@ Enhancements
 * Add :py:func:`pvlib.iam.marion_diffuse` and
   :py:func:`pvlib.iam.marion_integrate` to calculate IAM values for
   diffuse irradiance. (:pull:`984`)
+* Add :py:func:`pvlib.inverter.fit_sandia` that fits the Sandia inverter model
+  to a set of inverter efficiency curves. (:pull:`1011`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/data/inverter_fit_snl_meas.csv

@@ -0,0 +1,127 @@
+fraction_of_rated_power,dc_voltage_level,ac_power,dc_voltage,efficiency
+0.1,Vmin,32800,660.5,0.95814
+0.2,Vmin,73000,660.9,0.9755
+0.3,Vmin,107500,660.73,0.97787
+0.5,Vmin,168100,660.1,0.97998
+0.75,Vmin,235467,660.27,0.97785
+1,Vmin,318067,660.03,0.97258
+0.1,Vnom,32800,740.1,0.95441
+0.2,Vnom,72900,740.2,0.96985
+0.3,Vnom,107600,740.13,0.97611
+0.5,Vnom,167500,740.57,0.97554
+0.75,Vnom,234967,741.87,0.97429
+1,Vnom,317267,737.7,0.97261
+0.1,Vmax,32800,959.07,0.94165
+0.2,Vmax,71600,959.43,0.95979
+0.3,Vmax,107300,959.1,0.96551
+0.5,Vmax,166700,959.5,0.96787
+0.75,Vmax,234767,958.8,0.96612
+1,Vmax,317467,957,0.96358
+0.1,Vmin,32800,660.77,0.95721
+0.2,Vmin,73000,660.77,0.97247
+0.3,Vmin,107500,660.47,0.97668
+0.5,Vmin,168100,660.23,0.98018
+0.75,Vmin,235333.3333,660.3,0.97716
+1,Vmin,317466.6667,659.8,0.97184
+0.1,Vnom,32800,740.27,0.95534
+0.2,Vnom,72900,740.27,0.97071
+0.3,Vnom,107600,740.2,0.97523
+0.5,Vnom,167500,740.8,0.97592
+0.75,Vnom,234966.6667,741.67,0.97429
+1,Vnom,317300,737.97,0.97252
+0.1,Vmax,32800,959.23,0.93718
+0.2,Vmax,71600,959.4,0.96107
+0.3,Vmax,107300,959.27,0.96638
+0.5,Vmax,166700,959.57,0.96825
+0.75,Vmax,234733.3333,959.17,0.96731
+1,Vmax,317466.6667,957.07,0.96241
+0.1,Vmin,32800,660.57,0.95814
+0.2,Vmin,73000,660.67,0.97333
+0.3,Vmin,107500,660.5,0.97609
+0.5,Vmin,168100,660.1,0.97884
+0.75,Vmin,235066.6667,660.3,0.97781
+1,Vmin,316900,659.27,0.97209
+0.1,Vnom,32800,740.17,0.95441
+0.2,Vnom,72900,740.27,0.97028
+0.3,Vnom,107600,740.23,0.97464
+0.5,Vnom,167500,740.3,0.97573
+0.75,Vnom,235133.3333,742.13,0.97417
+1,Vnom,317300,737.9,0.97252
+0.1,Vmax,32800,959.2,0.93626
+0.2,Vmax,71600,959.43,0.95979
+0.3,Vmax,107300,959.2,0.96493
+0.5,Vmax,166700,959.5,0.96806
+0.75,Vmax,234833.3333,958.97,0.96573
+1,Vmax,317400,956.87,0.96279
+0.1,Vmin,32800,660.63,0.95627
+0.2,Vmin,73000,660.9,0.97377
+0.3,Vmin,107500,661.07,0.97846
+0.5,Vmin,168100,660.13,0.97827
+0.75,Vmin,235200,660.43,0.97701
+1,Vmin,316933.3333,660.07,0.97308
+0.1,Vnom,32800,740.27,0.95441
+0.2,Vnom,72900,740.37,0.96985
+0.3,Vnom,107600,740.27,0.97464
+0.5,Vnom,167500,740.53,0.97592
+0.75,Vnom,234800,742.13,0.97374
+1,Vnom,317300,737.73,0.97202
+0.1,Vmax,32800,959.2,0.93271
+0.2,Vmax,71600,959.27,0.95594
+0.3,Vmax,107300,959.2,0.96783
+0.5,Vmax,166700,959.47,0.96806
+0.75,Vmax,234700,958.67,0.96505
+1,Vmax,317433.3333,956.8,0.96299
+0.1,Vmin,32800,660.67,0.95534
+0.2,Vmin,73000,660.8,0.9755
+0.3,Vmin,107500,661.23,0.97905
+0.5,Vmin,168100,660.33,0.97941
+0.75,Vmin,236566.6667,660.43,0.97741
+1,Vmin,317866.6667,659.53,0.97366
+0.1,Vnom,32800,740.13,0.95627
+0.2,Vnom,72900,740.37,0.97071
+0.3,Vnom,107600,740.4,0.97523
+0.5,Vnom,167500,740.57,0.97649
+0.75,Vnom,234733.3333,741.83,0.97413
+1,Vnom,317333.3333,737.77,0.97222
+0.1,Vmax,32800,959.03,0.9336
+0.2,Vmax,71600,959.33,0.96108
+0.3,Vmax,107300,959.2,0.96464
+0.5,Vmax,166700,959.57,0.96975
+0.75,Vmax,234700,958.83,0.96584
+1,Vmax,317400,956.7,0.96338
+0.1,Vmin,32800,660.43,0.95349
+0.2,Vmin,73000,660.83,0.97247
+0.3,Vmin,107500,660.47,0.97668
+0.5,Vmin,168100,660.27,0.97941
+0.75,Vmin,236167,660.57,0.97657
+1,Vmin,317833,660.47,0.97177
+0.1,Vnom,32800,740.2,0.95534
+0.2,Vnom,72900,740.3,0.96985
+0.3,Vnom,107600,740.33,0.97434
+0.5,Vnom,167500,740.53,0.9763
+0.75,Vnom,234833,741.93,0.97468
+1,Vnom,317333,737.73,0.97242
+0.1,Vmax,32800,959.03,0.93626
+0.2,Vmax,71600,959.37,0.95936
+0.3,Vmax,107300,959.23,0.96464
+0.5,Vmax,166700,959.5,0.96731
+0.75,Vmax,235267,958.67,0.96592
+1,Vmax,317400,957.07,0.96269
+0.1,Vmin,32800,660.73,0.95627
+0.2,Vmin,73000,660.57,0.97204
+0.3,Vmin,107500,660.97,0.97787
+0.5,Vmin,168100,660,0.97865
+0.75,Vmin,236200,659.77,0.97778
+1,Vmin,317200,659.2,0.97221
+0.1,Vnom,32800,740.2,0.95257
+0.2,Vnom,72900,740.33,0.97071
+0.3,Vnom,107600,740.43,0.97464
+0.5,Vnom,167500,740.63,0.97592
+0.75,Vnom,235100,741.87,0.97458
+1,Vnom,317333.3333,737.83,0.97232
+0.1,Vmax,32800,959,0.93182
+0.2,Vmax,71600,959.4,0.9585
+0.3,Vmax,107300,959.07,0.96783
+0.5,Vmax,166700,959.7,0.96806
+0.75,Vmax,235466.6667,958.77,0.96569
+1,Vmax,317400,956.6,0.96308

pvlib/data/inverter_fit_snl_sim.csv

@@ -0,0 +1,19 @@
+fraction_of_rated_power,efficiency,dc_voltage_level,dc_voltage,dc_power,ac_power,efficiency
+0.1,0.892146066,Vmin,220,112.0892685,100,0.892146067
+0.1,0.876414009,Vnom,240,114.1013254,100,0.876414009
+0.1,0.861227164,Vmax,260,116.1133835,100,0.861227164
+0.2,0.925255801,Vmin,220,216.15644,200,0.925255801
+0.2,0.916673906,Vnom,240,218.1800951,200,0.916673906
+0.2,0.908249736,Vmax,260,220.2037524,200,0.908249736
+0.3,0.936909957,Vmin,220,320.2015283,300,0.936909957
+0.3,0.93099374,Vnom,240,322.2363236,300,0.93099374
+0.3,0.925151763,Vmax,260,324.2711217,300,0.925151763
+0.5,0.946565413,Vmin,220,528.2255121,500,0.946565413
+0.5,0.94289593,Vnom,240,530.2812159,500,0.94289593
+0.5,0.939254781,Vmax,260,532.336923,500,0.939254781
+0.75,0.951617818,Vmin,220,788.1315225,750,0.951617818
+0.75,0.949113828,Vnom,240,790.2108027,750,0.949113828
+0.75,0.946622992,Vmax,260,792.2900736,750,0.946622992
+1,0.954289529,Vmin,220,1047.900002,1000,0.95428953
+1,0.952380952,Vnom,240,1050,1000,0.952380952
+1,0.950479992,Vmax,260,1052.1,1000,0.950479992

pvlib/inverter.py

@@ -1,12 +1,21 @@
 # -*- coding: utf-8 -*-
 """
-This module contains functions for inverter modeling, primarily conversion of
-DC to AC power.
+This module contains functions for inverter modeling and for fitting inverter
+models to data.
+
+Inverter models calculate AC power output from DC input. Model parameters
+should be passed as a single dict.
+
+Functions for estimating parameters for inverter models should follow the
+naming pattern 'fit_<model name>', e.g., fit_sandia.
+
 """
 
 import numpy as np
 import pandas as pd
 
+from numpy.polynomial.polynomial import polyfit  # different than np.polyfit
+
 
 def sandia(v_dc, p_dc, inverter):
     r'''
@@ -176,7 +185,7 @@ def adr(v_dc, p_dc, inverter, vtol=0.10):
 
     References
     ----------
-    .. [1] Driesse, A. "Beyond the Curves: Modeling the Electrical Efficiency
+    .. [1] A. Driesse, "Beyond the Curves: Modeling the Electrical Efficiency
        of Photovoltaic Inverters", 33rd IEEE Photovoltaic Specialist
        Conference (PVSC), June 2008
 
@@ -285,8 +294,7 @@ def pvwatts(pdc, pdc0, eta_inv_nom=0.96, eta_inv_ref=0.9637):
     References
     ----------
     .. [1] A. P. Dobos, "PVWatts Version 5 Manual,"
-           http://pvwatts.nrel.gov/downloads/pvwattsv5.pdf
-           (2014).
+       http://pvwatts.nrel.gov/downloads/pvwattsv5.pdf (2014).
     """
 
     pac0 = eta_inv_nom * pdc0
@@ -306,3 +314,110 @@ def pvwatts(pdc, pdc0, eta_inv_nom=0.96, eta_inv_ref=0.9637):
     power_ac = np.maximum(0, power_ac)     # GH 541
 
     return power_ac
+
+
+def fit_sandia(ac_power, dc_power, dc_voltage, dc_voltage_level, p_ac_0, p_nt):
+    r'''
+    Determine parameters for the Sandia inverter model.
+
+    Parameters
+    ----------
+    ac_power : array_like
+        AC power output at each data point [W].
+    dc_power : array_like
+        DC power input at each data point [W].
+    dc_voltage : array_like
+        DC input voltage at each data point [V].
+    dc_voltage_level : array_like
+        DC input voltage level at each data point. Values must be 'Vmin',
+        'Vnom' or 'Vmax'.
+    p_ac_0 : float
+        Rated AC power of the inverter [W].
+    p_nt : float
+        Night tare, i.e., power consumed while inverter is not delivering
+        AC power. [W]
+
+    Returns
+    -------
+    dict
+        A set of parameters for the Sandia inverter model [1]_. See
+        :py:func:`pvlib.inverter.sandia` for a description of keys and values.
+
+    See Also
+    --------
+    pvlib.inverter.sandia
+
+    Notes
+    -----
+    The fitting procedure to estimate parameters is described at [2]_.
+    A data point is a pair of values (dc_power, ac_power). Typically, inverter
+    performance is measured or described at three DC input voltage levels,
+    denoted 'Vmin', 'Vnom' and 'Vmax' and at each level, inverter efficiency
+    is determined at various output power levels. For example,
+    the CEC inverter test protocol [3]_ specifies measurement of input DC
+    power that delivers AC output power of 0.1, 0.2, 0.3, 0.5, 0.75 and 1.0 of
+    the inverter's AC power rating.
+
+    References
+    ----------
+    .. [1] D. King, S. Gonzalez, G. Galbraith, W. Boyson, "Performance Model
+       for Grid-Connected Photovoltaic Inverters", SAND2007-5036, Sandia
+       National Laboratories.
+    .. [2] Sandia Inverter Model page, PV Performance Modeling Collaborative
+       https://pvpmc.sandia.gov/modeling-steps/dc-to-ac-conversion/sandia-inverter-model/
+    .. [3] W. Bower, et al., "Performance Test Protocol for Evaluating
+       Inverters Used in Grid-Connected Photovoltaic Systems", available at
+       https://www.energy.ca.gov/sites/default/files/2020-06/2004-11-22_Sandia_Test_Protocol_ada.pdf
+    '''  # noqa: E501
+
+    voltage_levels = ['Vmin', 'Vnom', 'Vmax']
+
+    # average dc input voltage at each voltage level
+    v_d = np.array(
+        [dc_voltage[dc_voltage_level == 'Vmin'].mean(),
+         dc_voltage[dc_voltage_level == 'Vnom'].mean(),
+         dc_voltage[dc_voltage_level == 'Vmax'].mean()])
+    v_nom = v_d[1]  # model parameter
+    # independent variable for regressions, x_d
+    x_d = v_d - v_nom
+
+    # empty dataframe to contain intermediate variables
+    coeffs = pd.DataFrame(index=voltage_levels,
+                          columns=['a', 'b', 'c', 'p_dc', 'p_s0'], data=np.nan)
+
+    def solve_quad(a, b, c):
+        return (-b + (b**2 - 4 * a * c)**.5) / (2 * a)
+
+    # [2] STEP 3E, fit a line to (DC voltage, model_coefficient)
+    def extract_c(x_d, add):
+        beta0, beta1 = polyfit(x_d, add, 1)
+        c = beta1 / beta0
+        return beta0, beta1, c
+
+    for d in voltage_levels:
+        x = dc_power[dc_voltage_level == d]
+        y = ac_power[dc_voltage_level == d]
+        # [2] STEP 3B
+        # fit a quadratic to (DC power, AC power)
+        c, b, a = polyfit(x, y, 2)
+
+        # [2] STEP 3D, solve for p_dc and p_s0
+        p_dc = solve_quad(a, b, (c - p_ac_0))
+        p_s0 = solve_quad(a, b, c)
+
+        # Add values to dataframe at index d
+        coeffs['a'][d] = a
+        coeffs['p_dc'][d] = p_dc
+        coeffs['p_s0'][d] = p_s0
+
+    b_dc0, b_dc1, c1 = extract_c(x_d, coeffs['p_dc'])
+    b_s0, b_s1, c2 = extract_c(x_d, coeffs['p_s0'])
+    b_c0, b_c1, c3 = extract_c(x_d, coeffs['a'])
+
+    p_dc0 = b_dc0
+    p_s0 = b_s0
+    c0 = b_c0
+
+    # prepare dict and return
+    return {'Paco': p_ac_0, 'Pdco': p_dc0, 'Vdco': v_nom, 'Pso': p_s0,
+            'C0': c0, 'C1': c1, 'C2': c2, 'C3': c3, 'Pnt': p_nt}

pvlib/tests/test_inverter.py

@@ -6,8 +6,10 @@
 from conftest import assert_series_equal
 from numpy.testing import assert_allclose
 
+from conftest import needs_numpy_1_10, DATA_DIR
+import pytest
+
 from pvlib import inverter
-from conftest import needs_numpy_1_10
 
 
 def test_adr(adr_inverter_parameters):
@@ -131,3 +133,26 @@ def test_pvwatts_series():
     expected = pd.Series(np.array([np.nan, 0., 47.608436, 95.]))
     out = inverter.pvwatts(pdc, pdc0, 0.95)
     assert_series_equal(expected, out)
+
+
+INVERTER_TEST_MEAS = DATA_DIR / 'inverter_fit_snl_meas.csv'
+INVERTER_TEST_SIM = DATA_DIR / 'inverter_fit_snl_sim.csv'
+
+
+@pytest.mark.parametrize('infilen, expected', [
+    (INVERTER_TEST_MEAS, {'Paco': 333000., 'Pdco': 343251., 'Vdco': 740.,
+                          'Pso': 1427.746, 'C0': -5.768e-08, 'C1': 3.596e-05,
+                          'C2': 1.038e-03, 'C3': 2.978e-05, 'Pnt': 1.}),
+    (INVERTER_TEST_SIM,  {'Paco': 1000., 'Pdco': 1050., 'Vdco': 240.,
+                          'Pso': 10., 'C0': 1e-6, 'C1': 1e-4, 'C2': 1e-2,
+                          'C3': 1e-3, 'Pnt': 1.}),
+])
+def test_fit_sandia(infilen, expected):
+    curves = pd.read_csv(infilen)
+    dc_power = curves['ac_power'] / curves['efficiency']
+    result = inverter.fit_sandia(ac_power=curves['ac_power'],
+                                 dc_power=dc_power,
+                                 dc_voltage=curves['dc_voltage'],
+                                 dc_voltage_level=curves['dc_voltage_level'],
+                                 p_ac_0=expected['Paco'], p_nt=expected['Pnt'])
+    assert expected == pytest.approx(result, rel=1e-3)