Merge ead8058 into aa04b16

docs/api.rst

@@ -241,6 +241,7 @@ power or plane of array irradiance measurements.
    :toctree: generated/
 
    system.infer_orientation_daily_peak
+   system.orientation_fit_pvwatts
 
 Metrics
 =======

pvanalytics/system.py

@@ -2,8 +2,10 @@
 import enum
 import warnings
 import numpy as np
+import scipy
 import pandas as pd
 import pvlib
+from pvlib.temperature import TEMPERATURE_MODEL_PARAMETERS
 from pvanalytics.util import _fit, _group
 
 
@@ -401,3 +403,160 @@ def infer_orientation_daily_peak(power_or_poa, sunny, tilts,
                 best_azimuth = azimuth
                 best_tilt = tilt
     return best_azimuth, best_tilt
+
+
+def _power_residuals_from_clearsky(system_params,
+                                   ghi, dhi, dni,
+                                   solar_zenith, solar_azimuth,
+                                   power_ac, temperature,
+                                   wind_speed,
+                                   temperature_coefficient,
+                                   temperature_model_parameters):
+    # Return the residuals between a system with parameters given in `params`
+    # and the data in `power_ac`.
+    #
+    # Parameters
+    # ----------
+    # system_params : array-like
+    #     array of four floats: tilt, azimuth, DC capacity, and inverter
+    #     DC input limit.
+    # ghi : Series
+    #     Clear sky GHI
+    # dhi : Series
+    #     Clear sky DHI
+    # dni : Series
+    #     Clear sky DNI
+    # solar_zenith : Series
+    #     Solar zenith at the same times as data in `power_ac`
+    # solar_azimuth : Series
+    #     Solar azimuth at the same times as data in `power_ac`
+    # power_ac : Series
+    #     Measured AC power under clear sky conditions.
+    # temperature : float or Series
+    #     Air temperature at which to model the hypothetical system. If a
+    #     float then a constant temperature is used. If a Series, must have
+    #     the same index as `power_ac`. [C]
+    # wind_speed : float or Series
+    #     Wind speed. If a float then a constant wind speed is used. If a
+    #     Series, must have the same index as `power_ac`. [m/s]
+    # temperature_model_parameters : dict
+    #     Parameters fot the cell temperature model.
+    #
+    # Returns
+    # -------
+    # Series
+    #     Difference between `power_ac` and the PVWatts output with the
+    #     given parameters.
+    tilt = system_params[0]
+    azimuth = system_params[1]
+    dc_capacity = system_params[2]
+    dc_limit = system_params[3]
+    poa = pvlib.irradiance.get_total_irradiance(
+        tilt, azimuth,
+        solar_zenith,
+        solar_azimuth,
+        dni, ghi, dhi
+    )
+    temp_cell = pvlib.temperature.sapm_cell(
+        poa['poa_global'],
+        temperature,
+        wind_speed,
+        **temperature_model_parameters
+    )
+    pdc = pvlib.pvsystem.pvwatts_dc(
+        poa['poa_global'],
+        temp_cell,
+        dc_capacity,
+        temperature_coefficient
+    )
+    return power_ac - pvlib.inverter.pvwatts(pdc, dc_limit)
+
+
+def _rsquared(data, residuals):
+    # Calculate the coefficient of determination from `residuals`
+    model = data + residuals
+    _, _, r, _, _ = scipy.stats.linregress(model, data)
+    return r*r
+
+
+def orientation_fit_pvwatts(power_ac, ghi, dhi, dni,
+                            solar_zenith, solar_azimuth,
+                            temperature=25, wind_speed=0,
+                            temperature_coefficient=-0.002,
+                            temperature_model_parameters=None):
+    """Get the tilt and azimuth that give pvwatts output that most closely
+    fits the data in `power_ac`.
+
+    Uses non-linear least squares to optimize over four free variables
+    to find the values that result in the best fit between power modeled
+    using PVWatts and `power_ac`. The four free variables are
+
+    - surface tilt
+    - surface azimuth
+    - the DC capacity of the system
+    - the DC input limit of the inverter.
+
+    Parameters
+    ----------
+    power_ac : Series
+        AC power from the system in clear sky conditions.
+    ghi : Series
+        Clear sky GHI with same index as `power_ac`.
+    dhi : Series
+        Clear sky DHI with same index as `power_ac`.
+    dni : Series
+        Clear sky DNI with same index as `power_ac`.
+    solar_zenith : Series
+        Solar zenith. [degrees]
+    solar_azimuth : Series
+        Solar azimuth. [degrees]
+    temperature : float or Series, default 25
+        Air temperature at which to model the hypothetical system. If a
+        float then a constant temperature is used. If a Series, must have
+        the same index as `power_ac`. [C]
+    wind_speed : float or Series, default 0
+        Wind speed. If a float then a constant wind speed is used. If a
+        Series, must have the same index as `power_ac`. [m/s]
+    temperature_model_parameters : dict, optional
+        Parameters fot the cell temperature model. If not specified
+        ``pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS['sapm'][
+        'open_rack_glass_glass']`` is used. See
+        :py:func:`pvlib.temperature.sapm_cell` for more information.
+
+    Returns
+    -------
+    surface_tilt : float
+        Tilt of the array. [degrees]
+    surface_azimuth : float
+        Azimuth of the array. [degrees]
+    r_squared : float
+        :math:`r^2` value for the fit at the returned orientation.
+
+    """
+    if temperature_model_parameters is None:
+        temperature_model_parameters = \
+            TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
+    initial_tilt = 45
+    initial_azimuth = 180
+    initial_dc_capacity = power_ac.max()
+    initial_dc_limit = power_ac.max() * 1.5
+    fit_result = scipy.optimize.least_squares(
+        _power_residuals_from_clearsky,
+        [initial_tilt, initial_azimuth, initial_dc_capacity, initial_dc_limit],
+        bounds=([0, 0, power_ac.max()*0.5, power_ac.max()*0.5],
+                [90, 360, power_ac.max()*2, power_ac.max()*3]),
+        kwargs={
+            'ghi': ghi,
+            'dhi': dhi,
+            'dni': dni,
+            'solar_zenith': solar_zenith,
+            'solar_azimuth': solar_azimuth,
+            'power_ac': power_ac,
+            'temperature': temperature,
+            'temperature_coefficient': temperature_coefficient,
+            'wind_speed': wind_speed,
+            'temperature_model_parameters': temperature_model_parameters
+        }
+    )
+    r_squared = _rsquared(power_ac, fit_result.fun)
+    return fit_result.x[0], fit_result.x[1], r_squared

pvanalytics/tests/conftest.py

@@ -2,7 +2,8 @@
 import pytest
 import numpy as np
 import pandas as pd
-from pvlib import location
+from pvlib import location, pvsystem
+from pvlib.temperature import TEMPERATURE_MODEL_PARAMETERS
 
 
 def pytest_addoption(parser):
@@ -99,3 +100,20 @@ def clearsky_year(one_year_hourly, albuquerque):
 def solarposition_year(one_year_hourly, albuquerque):
     """One year of solar position data in albuquerque"""
     return albuquerque.get_solarposition(one_year_hourly)
+
+
+@pytest.fixture(scope='module')
+def system_parameters():
+    """System parameters for generating simulated power data."""
+    sandia_modules = pvsystem.retrieve_sam('SandiaMod')
+    sapm_inverters = pvsystem.retrieve_sam('cecinverter')
+    module = sandia_modules['Canadian_Solar_CS5P_220M___2009_']
+    inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208__208V_']
+    temperature_model_parameters = (
+        TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
+    )
+    return {
+        'module_parameters': module,
+        'inverter_parameters': inverter,
+        'temperature_model_parameters': temperature_model_parameters
+    }

pvanalytics/tests/features/test_orientation.py

@@ -1,28 +1,10 @@
 import pytest
 from pandas.util.testing import assert_series_equal
 import pandas as pd
-from pvlib import pvsystem, tracking, modelchain, irradiance
-from pvlib.temperature import TEMPERATURE_MODEL_PARAMETERS
+from pvlib import tracking, modelchain, irradiance
 from pvanalytics.features import orientation
 
 
-@pytest.fixture(scope='module')
-def system_parameters():
-    """System parameters for generating simulated power data."""
-    sandia_modules = pvsystem.retrieve_sam('SandiaMod')
-    sapm_inverters = pvsystem.retrieve_sam('cecinverter')
-    module = sandia_modules['Canadian_Solar_CS5P_220M___2009_']
-    inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208__208V_']
-    temperature_model_parameters = (
-        TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
-    )
-    return {
-        'module_parameters': module,
-        'inverter_parameters': inverter,
-        'temperature_model_parameters': temperature_model_parameters
-    }
-
-
 def test_clearsky_ghi_fixed(clearsky, solarposition):
     """Identify every day as fixed, since clearsky GHI is sunny."""
     assert orientation.fixed_nrel(

pvanalytics/tests/test_system.py

@@ -1,8 +1,9 @@
-"""Tests for funcitons that identify system characteristics."""
+"""Tests for system parameter identification functions."""
 import pytest
-import numpy as np
 import pandas as pd
-from pvlib import pvsystem, tracking, modelchain, irradiance
+import numpy as np
+import pvlib
+from pvlib import location, pvsystem, tracking, modelchain, irradiance
 from pvlib.temperature import TEMPERATURE_MODEL_PARAMETERS
 from pvanalytics import system
 
@@ -51,9 +52,9 @@ def summer_ghi(summer_clearsky):
 @pytest.fixture
 def summer_power_fixed(summer_clearsky, albuquerque, system_parameters):
     """Simulated power from a FIXED PVSystem in Albuquerque, NM."""
-    system = pvsystem.PVSystem(**system_parameters)
+    pv_system = pvsystem.PVSystem(**system_parameters)
     mc = modelchain.ModelChain(
-        system,
+        pv_system,
         albuquerque,
         orientation_strategy='south_at_latitude_tilt'
     )
@@ -64,9 +65,9 @@ def summer_power_fixed(summer_clearsky, albuquerque, system_parameters):
 @pytest.fixture
 def summer_power_tracking(summer_clearsky, albuquerque, system_parameters):
     """Simulated power for a pvlib SingleAxisTracker PVSystem in Albuquerque"""
-    system = tracking.SingleAxisTracker(**system_parameters)
+    pv_system = tracking.SingleAxisTracker(**system_parameters)
     mc = modelchain.ModelChain(
-        system,
+        pv_system,
         albuquerque,
         orientation_strategy='south_at_latitude_tilt'
     )
@@ -406,3 +407,85 @@ def test_orientation_with_gaps(clearsky_year, solarposition_year):
     )
     assert azimuth == 180
     assert tilt == 15
+
+
+@pytest.fixture(scope='module')
+def naive_times():
+    """One year at 1-hour intervals"""
+    return pd.date_range(
+        start='2020',
+        end='2021',
+        freq='H'
+    )
+
+
+@pytest.fixture(scope='module',
+                params=[(35, -106, 'Etc/GMT+7'),
+                        (50, 10, 'Etc/GMT-1'),
+                        (-37, 144, 'Etc/GMT-10')],
+                ids=['Albuquerque', 'Berlin', 'Melbourne'])
+def system_location(request):
+    """Location of the system."""
+    return location.Location(
+        request.param[0], request.param[1], tz=request.param[2]
+    )
+
+
+@pytest.fixture(scope='module',
+                params=[(0, 180), (30, 180), (30, 90), (30, 270), (30, 0)],
+                ids=['South-0', 'South-30', 'East-30', 'West-30', 'North-30'])
+def system_power(request, system_location, naive_times):
+    tilt = request.param[0]
+    azimuth = request.param[1]
+    local_time = naive_times.tz_localize(system_location.tz)
+    clearsky = system_location.get_clearsky(
+        local_time, model='simplified_solis'
+    )
+    solar_position = system_location.get_solarposition(local_time)
+    poa = irradiance.get_total_irradiance(
+        tilt, azimuth,
+        solar_position['zenith'],
+        solar_position['azimuth'],
+        **clearsky
+    )
+    temp_cell = pvlib.temperature.sapm_cell(
+        poa['poa_global'],
+        25, 0,
+        **pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS[
+            'sapm'
+        ][
+            'open_rack_glass_glass'
+        ]
+    )
+    pdc = pvsystem.pvwatts_dc(poa['poa_global'], temp_cell, 100, -0.002)
+    pac = pvsystem.inverter.pvwatts(pdc, 120)
+    return {
+        'location': system_location,
+        'tilt': tilt,
+        'azimuth': azimuth,
+        'clearsky': clearsky,
+        'solar_position': solar_position,
+        'ac': pac
+    }
+
+
+@pytest.mark.slow
+def test_orientation_fit_pvwatts(system_power):
+    day_mask = system_power['ac'] > 0
+    tilt, azimuth, rsquared = system.orientation_fit_pvwatts(
+        system_power['ac'][day_mask],
+        solar_zenith=system_power['solar_position']['zenith'][day_mask],
+        solar_azimuth=system_power['solar_position']['azimuth'][day_mask],
+        **system_power['clearsky'][day_mask]
+    )
+    assert rsquared > 0.9
+    assert tilt == pytest.approx(system_power['tilt'], abs=10)
+    if system_power['tilt'] == 0:
+        # Any azimuth will give the same results at tilt 0.
+        return
+    if system_power['azimuth'] == 0:
+        # 0 degrees equals 360 degrees.
+        assert (azimuth == pytest.approx(0, abs=10)
+                or azimuth == pytest.approx(360, abs=10))
+    else:
+        assert azimuth == pytest.approx(system_power['azimuth'], abs=10)