diff --git a/docs/api.rst b/docs/api.rst index abdf0377..3d0ed60c 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -241,6 +241,7 @@ power or plane of array irradiance measurements. :toctree: generated/ system.infer_orientation_daily_peak + system.infer_orientation_fit_pvwatts Metrics ======= diff --git a/pvanalytics/system.py b/pvanalytics/system.py index b3e2ccf3..94b83a6d 100644 --- a/pvanalytics/system.py +++ b/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,195 @@ 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, + power_ac, + solar_zenith, solar_azimuth, + temperature, + wind_speed, + temperature_coefficient, + temperature_model_parameters): + """Return the residuals between a system with parameters given in + `system_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 + power_ac : Series + Measured AC power under clear sky conditions. + 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` + 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_coefficient : float + Temperature coefficient of DC power. [1/C] + temperature_model_parameters : dict + Parameters for the cell temperature model. + + Returns + ------- + Series + Difference between `power_ac` and the PVWatts output with the + given parameters. + + Notes + ------ + Uses the defaults in :py:func:`pvlib.irradiance.get_total_irradiance` to + calculated plane-of-array irradiance, i.e., the isotropic model for sky + diffuse irradiance, and the Perez irradiance transposition model. + """ + tilt = system_params[0] + azimuth = system_params[1] + dc_capacity = system_params[2] + dc_inverter_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_inverter_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 infer_orientation_fit_pvwatts(power_ac, ghi, dhi, dni, + solar_zenith, solar_azimuth, + temperature=25, wind_speed=0, + temperature_coefficient=-0.004, + temperature_model_parameters=None): + """Get the tilt and azimuth that give PVWatts output that most closely + fits the data in `power_ac`. + + Input data `power_ac`, `ghi`, `dhi`, `dni` should reflect clear-sky + conditions. + + 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. + + Of these four parameters, only tilt and azimuth are returned. While, DC + capacity and the DC input limit are calculated, their values may not be + accurate. While its value is not returned, because the DC input limit is + used as a free variable for the optimization process, this function + can operate on `power_ac` data that includes inverter clipping. + + All parameters passed as a Series must have the same index and must not + contain any undefined values (i.e. NaNs). If any input contains NaNs a + ValueError is raised. + + Parameters + ---------- + power_ac : Series + AC power from the system in clear sky conditions. + ghi : Series + Clear sky GHI with same index as `power_ac`. [W/m^2] + dhi : Series + Clear sky DHI with same index as `power_ac`. [W/m^2] + dni : Series + Clear sky DNI with same index as `power_ac`. [W/m^2] + 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. + + Raises + ------ + ValueError + If any input passed as a Series contains undefined values (i.e. NaNs). + """ + if power_ac.hasnans: + raise ValueError("power_ac must not contain undefined values") + if ghi.hasnans or dhi.hasnans or dni.hasnans: + raise ValueError("ghi, dhi, and dni must not contain undefined values") + if isinstance(temperature, pd.Series) and temperature.hasnans: + raise ValueError("temperature must not contain undefined values") + if isinstance(wind_speed, pd.Series) and wind_speed.hasnans: + raise ValueError("wind_speed must not contain undefined values") + 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 diff --git a/pvanalytics/tests/conftest.py b/pvanalytics/tests/conftest.py index c13c5ba0..fe2841e9 100644 --- a/pvanalytics/tests/conftest.py +++ b/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 + } diff --git a/pvanalytics/tests/features/test_orientation.py b/pvanalytics/tests/features/test_orientation.py index b85d19a3..dd80564d 100644 --- a/pvanalytics/tests/features/test_orientation.py +++ b/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( diff --git a/pvanalytics/tests/test_system.py b/pvanalytics/tests/test_system.py index 33d07364..2a7cad2e 100644 --- a/pvanalytics/tests/test_system.py +++ b/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,225 @@ 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.infer_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) + + +def test_orientation_fit_pvwatts_missing_data(naive_times): + tilt = 30 + azimuth = 100 + system_location = location.Location(35, -106) + local_time = naive_times.tz_localize('MST') + clearsky = system_location.get_clearsky( + local_time, model='simplified_solis' + ) + clearsky.loc['3/1/2020':'3/15/2020'] = np.nan + solar_position = system_location.get_solarposition(clearsky.index) + solar_position.loc['3/1/2020':'3/15/2020'] = np.nan + 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) + solar_position.dropna(inplace=True) + with pytest.raises(ValueError, + match=".* must not contain undefined values"): + system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + **clearsky + ) + pac.dropna(inplace=True) + with pytest.raises(ValueError, + match=".* must not contain undefined values"): + system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + **clearsky + ) + clearsky.dropna(inplace=True) + tilt_out, azimuth_out, rsquared = system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + **clearsky + ) + assert rsquared > 0.9 + assert tilt_out == pytest.approx(tilt, abs=10) + assert azimuth_out == pytest.approx(azimuth, abs=10) + + +def test_orientation_fit_pvwatts_temp_wind_as_series(naive_times): + tilt = 30 + azimuth = 100 + system_location = location.Location(35, -106) + local_time = naive_times.tz_localize('MST') + clearsky = system_location.get_clearsky( + local_time, model='simplified_solis' + ) + solar_position = system_location.get_solarposition(clearsky.index) + poa = irradiance.get_total_irradiance( + tilt, azimuth, + solar_position['zenith'], + solar_position['azimuth'], + **clearsky + ) + temp_cell = pvlib.temperature.sapm_cell( + poa['poa_global'], + 25, 1, + **pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS[ + 'sapm' + ][ + 'open_rack_glass_glass' + ] + ) + temperature = pd.Series(25, index=clearsky.index) + wind_speed = pd.Series(1, index=clearsky.index) + temperature_missing = temperature.copy() + temperature_missing.loc['4/5/2020':'4/10/2020'] = np.nan + wind_speed_missing = wind_speed.copy() + wind_speed_missing.loc['5/5/2020':'5/15/2020'] = np.nan + pdc = pvsystem.pvwatts_dc(poa['poa_global'], temp_cell, 100, -0.002) + pac = pvsystem.inverter.pvwatts(pdc, 120) + with pytest.raises(ValueError, + match=".* must not contain undefined values"): + system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + temperature=temperature_missing, + wind_speed=wind_speed_missing, + **clearsky + ) + with pytest.raises(ValueError, + match="temperature must not contain undefined values"): + system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + temperature=temperature_missing, + wind_speed=wind_speed, + **clearsky + ) + with pytest.raises(ValueError, + match="wind_speed must not contain undefined values"): + system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + temperature=temperature, + wind_speed=wind_speed_missing, + **clearsky + ) + # ValueError if indices don't match + with pytest.raises(ValueError): + system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + temperature=temperature_missing.dropna(), + wind_speed=wind_speed_missing.dropna(), + **clearsky + ) + tilt_out, azimuth_out, rsquared = system.infer_orientation_fit_pvwatts( + pac, + solar_zenith=solar_position['zenith'], + solar_azimuth=solar_position['azimuth'], + **clearsky + ) + assert rsquared > 0.9 + assert tilt_out == pytest.approx(tilt, abs=10) + assert azimuth_out == pytest.approx(azimuth, abs=10) diff --git a/requirements.txt b/requirements.txt index 8b8e54a0..6981aca0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ numpy>=1.16.0 pandas>=0.25.0,<1.1.0 -pvlib>=0.7.0 +pvlib>=0.8.0 scipy>=1.3.0 diff --git a/setup.py b/setup.py index 7bf1af55..36417b84 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,7 @@ INSTALL_REQUIRES = [ 'numpy >= 1.16.0', 'pandas >= 0.25.1', - 'pvlib >= 0.7.0', + 'pvlib >= 0.8.0', 'scipy >= 1.3.0' ]