Implement initial storage support

docs/sphinx/source/user_guide/index.rst

@@ -15,6 +15,7 @@ User Guide
    clearsky
    bifacial
    forecasts
+   storage
    comparison_pvlib_matlab
    variables_style_rules
    singlediode

pvlib/battery.py

@@ -0,0 +1,265 @@
+"""
+This module contains functions for modeling batteries.
+"""
+from pandas import DataFrame
+
+try:
+    from PySAM.BatteryStateful import default as sam_default
+    from PySAM.BatteryStateful import new as sam_new
+    from PySAM.BatteryTools import battery_model_sizing as sam_sizing
+except ImportError:  # pragma: no cover
+
+    def missing_nrel_pysam(*args, **kwargs):
+        raise ImportError(
+            "NREL's PySAM package required! (`pip install nrel-pysam`)"
+        )
+
+    sam_default = missing_nrel_pysam
+    sam_new = missing_nrel_pysam
+    sam_sizing = missing_nrel_pysam
+
+
+def offset_to_hours(offset):
+    """
+    Convert a Pandas offset into hours.
+
+    Parameters
+    ----------
+    offset : pd.tseries.offsets.BaseOffset
+        The input offset to convert.
+
+    Returns
+    -------
+    numeric
+        The resulting period, in hours.
+    """
+    if offset.name == "H":
+        return offset.n
+    if offset.name == "T":
+        return offset.n / 60
+    raise ValueError("Unsupported offset {}".format(offset))
+
+
+def power_to_energy(power):
+    """
+    Converts a power series to an energy series.
+
+    Assuming Watts as the input power unit, the output energy unit will be
+    Watt Hours.
+
+    Parameters
+    ----------
+    power : Series
+        The input power series. [W]
+
+    Returns
+    -------
+    The converted energy Series. [Wh]
+    """
+    return power * offset_to_hours(power.index.freq)
+
+
+def fit_boc(model):
+    """
+    Determine the BOC model matching the given characteristics.
+
+    Parameters
+    ----------
+    datasheet : dict
+        The datasheet parameters of the battery.
+
+    Returns
+    -------
+    dict
+        The BOC parameters.
+
+    Notes
+    -----
+    This function does not really perform a fitting procedure. Instead, it just
+    calculates the model parameters that match the provided information from a
+    datasheet.
+    """
+    params = {
+        "soc_percent": 50,
+    }
+    model.update(params)
+    return model
+
+
+def boc(model, dispatch):
+    """
+    Run a battery simulation with a provided dispatch series. Positive power
+    represents the power provided by the battery (i.e.: discharging) while
+    negative power represents power provided to the battery (i.e.: charging).
+
+    The provided dispatch series is the goal/target power, but the battery may
+    not be able to provide or store that much energy given its characteristics.
+    This function will calculate how much power will actually flow from/into
+    the battery.
+
+    Uses a simple "bag of Coulombs" model.
+
+    Parameters
+    ----------
+    model : dict
+        The initial BOC parameters.
+    dispatch : Series
+        The target power series. [W]
+
+    Returns
+    -------
+    export : dict
+        The final BOC parameters.
+    results : DataFrame
+        The resulting:
+
+        - Power flow. [W]
+        - SOC. [%]
+    """
+    min_soc = model.get("min_soc_percent", 10)
+    max_soc = model.get("max_soc_percent", 90)
+    factor = offset_to_hours(dispatch.index.freq)
+
+    states = []
+    current_energy = model["dc_energy_wh"] * model["soc_percent"] / 100
+    max_energy = model["dc_energy_wh"] * max_soc / 100
+    min_energy = model["dc_energy_wh"] * min_soc / 100
+
+    dispatch = dispatch.copy()
+    discharge_efficiency = model.get("discharge_efficiency", 1.0)
+    charge_efficiency = model.get("charge_efficiency", 1.0)
+    dispatch.loc[dispatch < 0] *= charge_efficiency
+    dispatch.loc[dispatch > 0] /= discharge_efficiency
+
+    for power in dispatch:
+        if power > 0:
+            power = min(power, model["dc_max_power_w"])
+            energy = power * factor
+            available = current_energy - min_energy
+            energy = min(energy, available)
+            power = energy / factor * discharge_efficiency
+        else:
+            power = max(power, -model["dc_max_power_w"])
+            energy = power * factor
+            available = current_energy - max_energy
+            energy = max(energy, available)
+            power = energy / factor / charge_efficiency
+        current_energy -= energy
+        soc = current_energy / model["dc_energy_wh"] * 100
+        states.append((power, soc))
+
+    results = DataFrame(states, index=dispatch.index, columns=["Power", "SOC"])
+
+    final_state = model.copy()
+    final_state["soc_percent"] = results.iloc[-1]["SOC"]
+
+    return (final_state, results)
+
+
+def fit_sam(datasheet):
+    """
+    Determine the SAM BatteryStateful model matching the given characteristics.
+
+    Parameters
+    ----------
+    datasheet : dict
+        The datasheet parameters of the battery.
+
+    Returns
+    -------
+    dict
+        The SAM BatteryStateful parameters.
+
+    Notes
+    -----
+    This function does not really perform a fitting procedure. Instead, it just
+    calculates the model parameters that match the provided information from a
+    datasheet.
+    """
+    chemistry = {
+        "LFP": "LFPGraphite",
+    }
+    model = sam_default(chemistry[datasheet["chemistry"]])
+    sam_sizing(
+        model=model,
+        desired_power=datasheet["dc_max_power_w"] / 1000,
+        desired_capacity=datasheet["dc_energy_wh"] / 1000,
+        desired_voltage=datasheet["dc_nominal_voltage"],
+    )
+    model.ParamsCell.initial_SOC = 50
+    model.ParamsCell.minimum_SOC = datasheet.get("min_soc_percent", 10)
+    model.ParamsCell.maximum_SOC = datasheet.get("max_soc_percent", 90)
+    export = model.export()
+    del export["Controls"]
+    result = {}
+    result["sam"] = export
+    result["charge_efficiency"] = datasheet.get("charge_efficiency", 1.0)
+    result["discharge_efficiency"] = datasheet.get("discharge_efficiency", 1.0)
+    return result
+
+
+def sam(model, power):
+    """
+    Run a battery simulation with a provided dispatch series. Positive power
+    represents the power provided by the battery (i.e.: discharging) while
+    negative power represents power provided to the battery (i.e.: charging).
+
+    The provided dispatch series is the goal/target power, but the battery may
+    not be able to provide or store that much energy given its characteristics.
+    This function will calculate how much power will actually flow from/into
+    the battery.
+
+    Uses SAM's BatteryStateful model.
+
+    Parameters
+    ----------
+    model : dict
+        The initial SAM BatteryStateful parameters.
+    dispatch : Series
+        The target dispatch power series. [W]
+
+    Returns
+    -------
+    export : dict
+        The final SAM BatteryStateful parameters.
+    results : DataFrame
+        The resulting:
+
+        - Power flow. [W]
+        - SOC. [%]
+    """
+    battery = sam_new()
+    battery.ParamsCell.assign(model["sam"].get("ParamsCell", {}))
+    battery.ParamsPack.assign(model["sam"].get("ParamsPack", {}))
+    battery.Controls.replace(
+        {
+            "control_mode": 1,
+            "dt_hr": offset_to_hours(power.index.freq),
+            "input_power": 0,
+        }
+    )
+    battery.setup()
+    battery.StateCell.assign(model["sam"].get("StateCell", {}))
+    battery.StatePack.assign(model["sam"].get("StatePack", {}))
+
+    battery_dispatch = power.copy()
+    discharge_efficiency = model.get("discharge_efficiency", 1.0)
+    charge_efficiency = model.get("charge_efficiency", 1.0)
+    battery_dispatch.loc[power < 0] *= charge_efficiency
+    battery_dispatch.loc[power > 0] /= discharge_efficiency
+
+    states = []
+    for p in battery_dispatch:
+        battery.Controls.input_power = p / 1000
+        battery.execute(0)
+        states.append((battery.StatePack.P * 1000, battery.StatePack.SOC))
+
+    results = DataFrame(states, index=power.index, columns=["Power", "SOC"])
+    results.loc[results["Power"] < 0, "Power"] /= charge_efficiency
+    results.loc[results["Power"] > 0, "Power"] *= discharge_efficiency
+    export = battery.export()
+    del export["Controls"]
+
+    state = model.copy()
+    state["sam"] = export
+    return (state, results)

pvlib/powerflow.py

@@ -0,0 +1,207 @@
+"""
+This module contains functions for simulating power flow.
+"""
+import numpy as np
+from pandas import DataFrame
+
+from pvlib.inverter import _sandia_eff
+
+
+def self_consumption(generation, load):
+    """
+    Calculate the power flow for a self-consumption use case. It assumes the
+    system is connected to the grid.
+
+    Parameters
+    ----------
+    generation : Series
+        The AC generation profile. [W]
+    load : Series
+        The load profile. [W]
+
+    Returns
+    -------
+    DataFrame
+        The resulting power flow provided by the system and the grid into the
+        system, grid and load. [W]
+    """
+    df = DataFrame(index=generation.index)
+    df["Grid to system"] = -generation.loc[generation < 0]
+    df["Grid to system"] = df["Grid to system"].fillna(0.0)
+    df["Generation"] = generation.loc[generation > 0]
+    df["Generation"] = df["Generation"].fillna(0.0)
+    df["Load"] = load
+    df["System to load"] = df[["Generation", "Load"]].min(axis=1, skipna=False)
+    df.loc[df["System to load"] < 0, "System to load"] = 0.0
+    df["System to grid"] = df["Generation"] - df["System to load"]
+    df["Grid to load"] = df["Load"] - df["System to load"]
+    df["Grid"] = df[["Grid to system", "Grid to load"]].sum(
+        axis=1, skipna=False
+    )
+    return df
+
+
+def self_consumption_ac_battery(df, dispatch, battery, model):
+    """
+    Calculate the power flow for a self-consumption use case with an
+    AC-connected battery and a custom dispatch series. It assumes the system is
+    connected to the grid.
+
+    Parameters
+    ----------
+    df : DataFrame
+        The self-consumption power flow solution. [W]
+    dispatch : Series
+        The dispatch series to use.
+    battery : dict
+        The battery parameters.
+    model : str
+        The battery model to use.
+
+    Returns
+    -------
+    DataFrame
+        The resulting power flow provided by the system, the grid and the
+        battery into the system, grid, battery and load. [W]
+    """
+    final_state, results = model(battery, dispatch)
+    df = df.copy()
+    df["System to battery"] = -results["Power"]
+    df.loc[df["System to battery"] < 0, "System to battery"] = 0.0
+    df["System to battery"] = df[["System to battery", "System to grid"]].min(
+        axis=1
+    )
+    df["System to grid"] -= df["System to battery"]
+    df["Battery to load"] = results["Power"]
+    df.loc[df["Battery to load"] < 0, "Battery to load"] = 0.0
+    df["Battery to load"] = df[["Battery to load", "Grid to load"]].min(axis=1)
+    df["Grid to load"] -= df["Battery to load"]
+    df["Grid"] = df[["Grid to system", "Grid to load"]].sum(
+        axis=1, skipna=False
+    )
+    return final_state, df
+
+
+def self_consumption_dc_battery(dc_solution, load):
+    """
+    Calculate the power flow for a self-consumption use case with a
+    DC-connected battery. It assumes the system is connected to the grid.
+
+    Parameters
+    ----------
+    dc_solution : DataFrame
+        The DC-connected inverter power flow solution. [W]
+    load : Series
+        The load profile. [W]
+
+    Returns
+    -------
+    DataFrame
+        The resulting power flow provided by the system, the grid and the
+        battery into the system, grid, battery and load. [W]
+    """
+    df = self_consumption(dc_solution["AC power"], load)
+    df["Battery"] = df["Generation"] * dc_solution["Battery factor"]
+    df["Battery to load"] = df[["Battery", "System to load"]].min(axis=1)
+    df["Battery to grid"] = df["Battery"] - df["Battery to load"]
+    df["PV to battery"] = -dc_solution["Battery power flow"]
+    df.loc[df["PV to battery"] < 0, "PV to battery"] = 0.0
+    df["PV to load"] = df["System to load"] - df["Battery to load"]
+    return df
+
+
+def multi_dc_battery(
+    v_dc, p_dc, inverter, battery_dispatch, battery_parameters, battery_model
+):
+    """
+    Calculate the power flow for a self-consumption use case with a
+    DC-connected battery. It assumes the system is connected to the grid.
+
+    Parameters
+    ----------
+    v_dc : numeric
+        DC voltage input to the inverter. [V]
+    p_dc : numeric
+        DC power input to the inverter. [W]
+    inverter : dict
+        Inverter parameters.
+    battery_dispatch : Series
+        Battery power dispatch series. [W]
+    battery_parameters : dict
+        Battery parameters.
+    battery_model : str
+        Battery model.
+
+    Returns
+    -------
+    DataFrame
+        The resulting inverter power flow.
+    """
+    dispatch = battery_dispatch.copy()
+
+    # Limit charging to the available DC power
+    power_dc = sum(p_dc)
+    max_charging = -power_dc
+    charging_mask = dispatch < 0
+    dispatch[charging_mask] = np.max([dispatch, max_charging], axis=0)[
+        charging_mask
+    ]
+
+    # Limit discharging to the inverter's maximum output power (approximately)
+    # Note this can revert the dispatch and charge when there is too much DC
+    # power (prevents clipping)
+    max_discharging = inverter['Paco'] - power_dc
+    discharging_mask = dispatch > 0
+    dispatch[discharging_mask] = np.min([dispatch, max_discharging], axis=0)[
+        discharging_mask
+    ]
+
+    # Calculate the actual battery power flow
+    final_state, battery_flow = battery_model(battery_parameters, dispatch)
+    charge = -battery_flow['Power'].copy()
+    charge.loc[charge < 0] = 0
+    discharge = battery_flow['Power'].copy()
+    discharge.loc[discharge < 0] = 0
+
+    # Adjust the DC power
+    ratios = [sum(power) / sum(power_dc) for power in p_dc]
+    adjusted_p_dc = [
+        power - ratio * charge for (power, ratio) in zip(p_dc, ratios)
+    ]
+    final_dc_power = sum(adjusted_p_dc) + discharge
+
+    # PV-contributed AC power
+    pv_ac_power = 0.0 * final_dc_power
+    for vdc, pdc in zip(v_dc, adjusted_p_dc):
+        array_contribution = (
+            pdc / final_dc_power * _sandia_eff(vdc, final_dc_power, inverter)
+        )
+        array_contribution[np.isnan(array_contribution)] = 0.0
+        pv_ac_power += array_contribution
+
+    # Battery-contributed AC power
+    vdc = inverter["Vdcmax"] / 2
+    pdc = discharge
+    battery_ac_power = (
+        pdc / final_dc_power * _sandia_eff(vdc, final_dc_power, inverter)
+    )
+    battery_ac_power[np.isnan(battery_ac_power)] = 0.0
+
+    # Total AC power
+    total_ac_power = pv_ac_power + battery_ac_power
+
+    # Limit output power (Sandia limits)
+    clipping = total_ac_power - inverter["Paco"]
+    clipping[clipping < 0] = 0
+    limited_ac_power = total_ac_power - clipping
+    battery_factor = battery_ac_power / limited_ac_power
+    min_ac_power = -1.0 * abs(inverter["Pnt"])
+    below_limit = final_dc_power < inverter["Pso"]
+    limited_ac_power[below_limit] = min_ac_power
+
+    result = DataFrame(index=dispatch.index)
+    result["Battery power flow"] = battery_flow["Power"]
+    result["AC power"] = limited_ac_power
+    result["Clipping"] = clipping
+    result["Battery factor"] = battery_factor
+    return result

pvlib/tests/conftest.py

@@ -1,12 +1,16 @@
-from pathlib import Path
+import os
 import platform
 import warnings
+from functools import lru_cache
+from functools import wraps
+from importlib.resources import files
+from pathlib import Path
 
 import pandas as pd
-import os
-from pkg_resources import parse_version
 import pytest
-from functools import wraps
+from numpy import nan
+from numpy.random import uniform
+from pkg_resources import parse_version
 
 import pvlib
 from pvlib.location import Location
@@ -493,3 +497,102 @@ def sapm_module_params():
                   'IXXO': 3.18803,
                   'FD': 1}
     return parameters
+
+
+@pytest.fixture(scope="function")
+def datasheet_battery_params():
+    """
+    Define some datasheet battery parameters for testing.
+
+    The scope of the fixture is set to ``'function'`` to allow tests to modify
+    parameters if required without affecting other tests.
+    """
+    parameters = {
+        "brand": "BYD",
+        "model": "HVS 5.1",
+        "width": 0.585,
+        "height": 0.712,
+        "depth": 0.298,
+        "weight": 91,
+        "chemistry": "LFP",
+        "mode": "AC",
+        "charge_efficiency": 0.96,
+        "discharge_efficiency": 0.96,
+        "min_soc_percent": 10,
+        "max_soc_percent": 90,
+        "dc_modules": 2,
+        "dc_modules_in_series": 2,
+        "dc_energy_wh": 5120,
+        "dc_nominal_voltage": 204,
+        "dc_max_power_w": 5100,
+    }
+    return parameters
+
+
+@pytest.fixture(scope="session")
+def residential_load_profile_generator():
+    """
+    Get a sample residential hourly load profile for testing purposes.
+
+    Returns
+    -------
+    The load profile.
+    """
+    def profile_generator(index):
+        load = pd.Series(data=nan, index=index)
+        load[load.index.hour == 0] = 600
+        load[load.index.hour == 4] = 400
+        load[load.index.hour == 13] = 1100
+        load[load.index.hour == 17] = 800
+        load[load.index.hour == 21] = 1300
+        load *= uniform(low=0.6, high=1.4, size=len(load))
+        load = load.interpolate(method="spline", order=2)
+        load = load.bfill().ffill()
+        return load
+    return profile_generator
+
+
+@pytest.fixture(scope="session")
+def residential_model_chain():
+    """
+    Get a sample residential hourly generation profile for testing purposes.
+
+    Returns
+    -------
+    The generation profile.
+    """
+    name = 'Madrid'
+    latitude = 40.31672645215922
+    longitude = -3.674695061062714
+    altitude = 603
+    timezone = 'Europe/Madrid'
+    module = pvlib.pvsystem.retrieve_sam('SandiaMod')['Canadian_Solar_CS5P_220M___2009_']
+    inverter = pvlib.pvsystem.retrieve_sam('cecinverter')['Powercom__SLK_1500__240V_']
+    weather = pvlib.iotools.get_pvgis_tmy(latitude, longitude, map_variables=True)[0]
+    weather.index = pd.date_range(
+        start=weather.index[0].replace(year=2021),
+        end=weather.index[-1].replace(year=2021),
+        freq="H",
+    )
+    weather.index = weather.index.tz_convert(timezone)
+    weather.index.name = "Timestamp"
+    location = pvlib.location.Location(
+        latitude,
+        longitude,
+        name=name,
+        altitude=altitude,
+        tz=timezone,
+    )
+    mount = pvlib.pvsystem.FixedMount(surface_tilt=latitude, surface_azimuth=180)
+    temperature_model_parameters = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
+    array = pvlib.pvsystem.Array(
+        mount=mount,
+        module_parameters=module,
+        modules_per_string=16,
+        strings=1,
+        temperature_model_parameters=temperature_model_parameters,
+    )
+    system = pvlib.pvsystem.PVSystem(arrays=[array], inverter_parameters=inverter)
+    mc = pvlib.modelchain.ModelChain(system, location)
+    mc.run_model(weather)
+    return mc

pvlib/tests/test_battery.py

@@ -0,0 +1,347 @@
+from itertools import product
+
+import pytest
+from pandas import Series
+from pandas import date_range
+from pvlib.tests.conftest import requires_pysam
+from pytest import approx
+from pytest import mark
+from pytest import raises
+
+from pvlib.battery import boc
+from pvlib.battery import fit_boc
+from pvlib.battery import fit_sam
+from pvlib.battery import power_to_energy
+from pvlib.battery import sam
+
+all_models = mark.parametrize(
+    "fit,run",
+    [
+        (fit_boc, boc),
+        pytest.param(fit_sam, sam, marks=requires_pysam),
+    ],
+)
+
+all_efficiencies = mark.parametrize(
+    "charge_efficiency,discharge_efficiency",
+    list(product([1.0, 0.98, 0.95], repeat=2)),
+)
+
+
+@mark.parametrize(
+    "power_value,frequency,energy_value",
+    [
+        (1000, "H", 1000),
+        (1000, "2H", 2000),
+        (1000, "15T", 250),
+        (1000, "60T", 1000),
+    ],
+)
+def test_power_to_energy(power_value, frequency, energy_value):
+    """
+    The function should be able to convert power to energy for different power
+    series' frequencies.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=10,
+        freq=frequency,
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(power_value, index=index)
+    energy = power_to_energy(power)
+    assert approx(energy) == energy_value
+
+
+def test_power_to_energy_unsupported_frequency():
+    """
+    When the power series' frequency is unsupported, the function raises an
+    exception.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=10,
+        freq="1M",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(1000, index=index)
+    with raises(ValueError, match=r"Unsupported offset"):
+        power_to_energy(power)
+
+
+def test_fit_boc(datasheet_battery_params):
+    """
+    The function returns a dictionary with the BOC model parameters.
+    """
+    model = fit_boc(datasheet_battery_params)
+    assert model["soc_percent"] == 50.0
+
+
+@requires_pysam
+def test_fit_sam(datasheet_battery_params):
+    """
+    The function returns a dictionary with a `"sam"` key that must be
+    assignable to a SAM BatteryStateful model. Parameters like the nominal
+    voltage and the battery energy must also be properly inherited.
+    """
+    from PySAM.BatteryStateful import new
+
+    model = fit_sam(datasheet_battery_params)
+    battery = new()
+    battery.assign(model["sam"])
+    assert approx(battery.value("nominal_voltage")) == 204
+    assert approx(battery.value("nominal_energy")) == 5.12
+
+
+@requires_pysam
+def test_fit_sam_controls(datasheet_battery_params):
+    """
+    The controls group should not be exported as part of the battery state when
+    creating a new SAM battery.
+    """
+    model = fit_sam(datasheet_battery_params)
+    assert "Controls" not in set(model.keys())
+
+
+@requires_pysam
+def test_sam_controls(datasheet_battery_params):
+    """
+    The controls group should not be exported as part of the battery state when
+    running a simulation with the SAM model.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=100,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(1000.0, index=index)
+    state = fit_sam(datasheet_battery_params)
+    state, _ = sam(state, power)
+    assert "Controls" not in set(state.keys())
+
+
+@all_models
+def test_model_return_index(datasheet_battery_params, fit, run):
+    """
+    The returned index must match the index of the given input power series.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=100,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(1000.0, index=index)
+    state = fit(datasheet_battery_params)
+    _, result = run(state, power)
+    assert all(result.index == power.index)
+
+
+@all_models
+def test_model_offset_valueerror(datasheet_battery_params, fit, run):
+    """
+    When the power series' frequency is unsupported, the models must raise an
+    exception.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=10,
+        freq="1M",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(1000, index=index)
+    state = fit(datasheet_battery_params)
+    with raises(ValueError, match=r"Unsupported offset"):
+        run(state, power)
+
+
+@all_models
+@all_efficiencies
+def test_model_dispatch_power(
+    datasheet_battery_params,
+    fit,
+    run,
+    charge_efficiency,
+    discharge_efficiency,
+):
+    """
+    The dispatch power series represents the power flow as seen from the
+    outside. As long as the battery is capable to provide sufficient power and
+    sufficient energy, the resulting power flow should match the provided
+    dispatch series independently of the charging/discharging efficiencies.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=10,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    dispatch = Series(100.0, index=index)
+    state = fit(datasheet_battery_params)
+    _, result = run(state, dispatch)
+    assert approx(result["Power"], rel=0.005) == dispatch
+
+
+@all_models
+@all_efficiencies
+def test_model_soc_value(
+    datasheet_battery_params,
+    fit,
+    run,
+    charge_efficiency,
+    discharge_efficiency,
+):
+    """
+    The SOC should be updated according to the power flow and battery capacity.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=20,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    step_percent = 1.5
+    step_power = step_percent / 100 * datasheet_battery_params["dc_energy_wh"]
+    power = Series(step_power, index=index)
+
+    state = fit(datasheet_battery_params)
+    _, result = run(state, power)
+    assert (
+        approx(result["SOC"].diff().iloc[-10:].mean(), rel=0.05)
+        == -step_percent / datasheet_battery_params["discharge_efficiency"]
+    )
+
+
+@all_models
+def test_model_charge_convergence(datasheet_battery_params, fit, run):
+    """
+    Charging only should converge into almost no power flow and maximum SOC.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=100,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(-1000, index=index)
+    state = fit(datasheet_battery_params)
+    _, result = run(state, power)
+    assert result["Power"].iloc[-1] == approx(0, abs=0.01)
+    assert result["SOC"].iloc[-1] == approx(90, rel=0.01)
+
+
+@all_models
+def test_model_discharge_convergence(datasheet_battery_params, fit, run):
+    """
+    Discharging only should converge into almost no power flow and minimum SOC.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=100,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(1000, index=index)
+    state = fit(datasheet_battery_params)
+    _, result = run(state, power)
+    assert result["Power"].iloc[-1] == approx(0, abs=0.01)
+    assert result["SOC"].iloc[-1] == approx(10, rel=0.01)
+
+
+@all_models
+def test_model_chain(datasheet_battery_params, fit, run):
+    """
+    The returning state must be reusable. Simulating continuously for ``2n``
+    steps should be the same as splitting the simulation in 2 for ``n`` steps.
+    """
+    index = date_range(
+        start="2022-01-01",
+        periods=100,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    power = Series(2000.0, index=index)
+    power.iloc[::2] = -2000.0
+    half_length = int(len(power) / 2)
+
+    continuous_state = fit(datasheet_battery_params)
+    continuous_state, continuous_power = run(continuous_state, power)
+
+    split_state = fit(datasheet_battery_params)
+    split_state, split_power_0 = run(split_state, power[:half_length])
+    split_state, split_power_1 = run(split_state, power[half_length:])
+    split_power = split_power_0.append(split_power_1)
+
+    assert split_state == continuous_state
+    assert approx(split_power) == continuous_power
+
+
+@all_models
+def test_model_equivalent_periods(datasheet_battery_params, fit, run):
+    """
+    The results of a simulation with a 1-hour period should match those of a
+    simulation with a 60-minutes period.
+    """
+    battery = fit(datasheet_battery_params)
+    hourly_index = date_range(
+        start="2022-01-01",
+        periods=50,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    minutely_index = date_range(
+        start="2022-01-01",
+        periods=50,
+        freq="60T",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+
+    _, hourly = run(battery, Series(20.0, index=hourly_index))
+    _, minutely = run(battery, Series(20.0, index=minutely_index))
+
+    assert approx(hourly) == minutely
+
+
+@all_models
+def test_model_equivalent_power_timespan(datasheet_battery_params, fit, run):
+    """
+    Simulating with the same constant input power over the same time span but
+    with different frequency should yield similar results.
+    """
+    battery = fit(datasheet_battery_params)
+    half_index = date_range(
+        start="2022-01-01 00:00",
+        end="2022-01-02 00:00",
+        freq="30T",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+    double_index = date_range(
+        start="2022-01-01 00:00",
+        end="2022-01-02 00:00",
+        freq="60T",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+
+    _, half = run(battery, Series(20.0, index=half_index))
+    _, double = run(battery, Series(20.0, index=double_index))
+
+    assert approx(half.iloc[-1]["SOC"], rel=0.001) == double.iloc[-1]["SOC"]
+    assert (
+        approx(power_to_energy(half["Power"]).sum(), rel=0.001)
+        == power_to_energy(double["Power"]).sum()
+    )

pvlib/tests/test_powerflow.py

@@ -0,0 +1,363 @@
+from numpy import array
+from numpy import nan
+from numpy.random import uniform
+from pandas import Series
+from pandas import Timestamp
+from pandas import date_range
+from pytest import approx
+from pytest import mark
+
+from pvlib.battery import boc
+from pvlib.battery import fit_boc
+from pvlib.powerflow import multi_dc_battery
+from pvlib.powerflow import self_consumption
+from pvlib.powerflow import self_consumption_ac_battery
+
+
+def gen_hourly_index(periods):
+    return date_range(
+        "2022-01-01",
+        periods=periods,
+        freq="1H",
+        tz="Europe/Madrid",
+        closed="left",
+    )
+
+
+@mark.parametrize(
+    "generation,load,flow",
+    [
+        (
+            42,
+            20,
+            {
+                "Generation": 42,
+                "Load": 20,
+                "System to load": 20,
+                "System to grid": 22,
+                "Grid to load": 0,
+                "Grid to system": 0,
+                "Grid": 0,
+            },
+        ),
+        (
+            42,
+            42,
+            {
+                "Generation": 42,
+                "Load": 42,
+                "System to load": 42,
+                "System to grid": 0,
+                "Grid to load": 0,
+                "Grid to system": 0,
+                "Grid": 0,
+            },
+        ),
+        (
+            42,
+            50,
+            {
+                "Generation": 42,
+                "Load": 50,
+                "System to load": 42,
+                "System to grid": 0,
+                "Grid to load": 8,
+                "Grid to system": 0,
+                "Grid": 8,
+            },
+        ),
+        (
+            -3,
+            0,
+            {
+                "Generation": 0,
+                "Load": 0,
+                "System to load": 0,
+                "System to grid": 0,
+                "Grid to load": 0,
+                "Grid to system": 3,
+                "Grid": 3,
+            },
+        ),
+        (
+            -3,
+            42,
+            {
+                "Generation": 0,
+                "Load": 42,
+                "System to load": 0,
+                "System to grid": 0,
+                "Grid to load": 42,
+                "Grid to system": 3,
+                "Grid": 45,
+            },
+        ),
+    ],
+    ids=[
+        "Positive generation with lower load",
+        "Positive generation with same load",
+        "Positive generation with higher load",
+        "Negative generation with zero load",
+        "Negative generation with positive load",
+    ],
+)
+def test_self_consumption(generation, load, flow):
+    """
+    Check multiple conditions with well-known cases:
+
+    - Excess generation must flow into grid
+    - Load must be fed with the system when possible, otherwise from grid
+    - Grid must provide energy for the system when required (i.e.: night hours)
+    - Negative values from the input generation are removed from the output
+      generation and added to the grid-to-system flow
+    """
+    result = (
+        self_consumption(
+            generation=Series([generation]),
+            load=Series([load]),
+        )
+        .iloc[0]
+        .to_dict()
+    )
+    assert approx(result) == flow
+
+
+def test_self_consumption_sum():
+    """
+    The sum of the flows with respect to the system, load and grid must be
+    balanced.
+    """
+    flow = self_consumption(
+        generation=Series(uniform(0, 1, 1000)),
+        load=Series(uniform(0, 1, 1000)),
+    )
+    assert (
+        approx(flow["Generation"])
+        == flow["System to load"] + flow["System to grid"]
+    )
+    assert (
+        approx(flow["Grid"]) == flow["Grid to load"] + flow["Grid to system"]
+    )
+    assert (
+        approx(flow["Load"]) == flow["System to load"] + flow["Grid to load"]
+    )
+    assert (
+        approx(flow["Load"] + flow["Grid to system"])
+        == flow["System to load"] + flow["Grid"]
+    )
+
+
+def test_self_consumption_ac_battery_sum(datasheet_battery_params):
+    """
+    The sum of the flows with respect to the system, load, grid and battery
+    must be balanced.
+    """
+    self_consumption_flow = self_consumption(
+        generation=Series(uniform(0, 1, 1000), index=gen_hourly_index(1000)),
+        load=Series(uniform(0, 1, 1000), index=gen_hourly_index(1000)),
+    )
+    dispatch = (
+        self_consumption_flow["Grid to load"]
+        - self_consumption_flow["System to grid"]
+    )
+    _, flow = self_consumption_ac_battery(
+        self_consumption_flow,
+        dispatch=dispatch,
+        battery=fit_boc(datasheet_battery_params),
+        model=boc,
+    )
+    assert (
+        approx(flow["Generation"])
+        == flow["System to load"]
+        + flow["System to battery"]
+        + flow["System to grid"]
+    )
+    assert (
+        approx(flow["Grid"]) == flow["Grid to load"] + flow["Grid to system"]
+    )
+    assert (
+        approx(flow["Load"])
+        == flow["System to load"]
+        + flow["Battery to load"]
+        + flow["Grid to load"]
+    )
+
+
+@mark.parametrize(
+    "charge_efficiency,discharge_efficiency,efficiency",
+    [
+        (1.0, 1.0, 1.0),
+        (0.97, 1.0, 0.97),
+        (1.0, 0.95, 0.95),
+        (0.97, 0.95, 0.97 * 0.95),
+    ],
+)
+def test_self_consumption_ac_battery_losses(
+    datasheet_battery_params,
+    residential_model_chain,
+    residential_load_profile_generator,
+    charge_efficiency,
+    discharge_efficiency,
+    efficiency,
+):
+    """
+    AC-DC conversion losses must be taken into account.
+
+    With the BOC model these losses are easy to track if we setup a simulation
+    in which we make sure to begin and end with an "empty" battery.
+    """
+    datasheet_battery_params["charge_efficiency"] = charge_efficiency
+    datasheet_battery_params["discharge_efficiency"] = discharge_efficiency
+    battery = fit_boc(datasheet_battery_params)
+    generation = residential_model_chain.results.ac.copy()
+    generation.iloc[:1000] = 0.0
+    generation.iloc[-1000:] = 0.0
+    load = residential_load_profile_generator(generation.index)
+    self_consumption_flow = self_consumption(
+        generation=generation,
+        load=load,
+    )
+    dispatch = (
+        self_consumption_flow["Grid to load"]
+        - self_consumption_flow["System to grid"]
+    )
+    _, lossy = self_consumption_ac_battery(
+        self_consumption_flow,
+        dispatch=dispatch,
+        battery=battery,
+        model=boc,
+    )
+    lossy = lossy.iloc[1000:]
+    assert lossy["Battery to load"].sum() / lossy[
+        "System to battery"
+    ].sum() == approx(efficiency)
+
+
+def test_self_consumption_nan_load():
+    """
+    When the load is unknown (NaN), the calculated flow to load should also be
+    unknown.
+    """
+    flow = self_consumption(
+        generation=Series([1, -2, 3, -4]),
+        load=Series([nan, nan, nan, nan]),
+    )
+    assert flow["System to load"].isna().all()
+    assert flow["Grid to load"].isna().all()
+
+
+@mark.parametrize(
+    "inputs,outputs",
+    [
+        (
+            {"pv_power": 800, "dispatch": -400},
+            {
+                "Battery power flow": -400,
+                "AC power": 400,
+                "Clipping": 0,
+                "Battery factor": 0,
+            },
+        ),
+        (
+            {"pv_power": 200, "dispatch": -600},
+            {
+                "Battery power flow": -200,
+                "AC power": 0,
+                "Clipping": 0,
+                "Battery factor": nan,
+            },
+        ),
+        (
+            {"pv_power": 1200, "dispatch": 400},
+            {
+                "Battery power flow": -200,
+                "AC power": 1000,
+                "Clipping": 0,
+                "Battery factor": 0,
+            },
+        ),
+        (
+            {"pv_power": 2000, "dispatch": 400},
+            {
+                "Battery power flow": -850,
+                "AC power": 1000,
+                "Clipping": 150,
+                "Battery factor": 0,
+            },
+        ),
+        (
+            {"pv_power": 100, "dispatch": 400},
+            {
+                "Battery power flow": 400,
+                "AC power": 500,
+                "Clipping": 0,
+                "Battery factor": 0.8,
+            },
+        ),
+        (
+            {"pv_power": 400, "dispatch": 1000},
+            {
+                "Battery power flow": 600,
+                "AC power": 1000,
+                "Clipping": 0,
+                "Battery factor": 0.6,
+            },
+        ),
+    ],
+    ids=[
+        "Charging is prioritized over AC conversion while enough PV power is available",
+        "Charging is limited by the available PV power",
+        "Clipping forces battery to charge, even when dispatch is set to discharge",
+        "Clipping cannot be avoided if the battery is unable to handle too much input power",
+        "Battery discharge can be combined with PV power to provide higher AC output power",
+        "Battery discharge is limited to the inverter's nominal power",
+    ],
+)
+def test_multi_dc_battery(inputs, outputs, datasheet_battery_params):
+    """
+    Test well-known cases of a multi-input (PV) and DC-connected battery
+    inverter.
+
+    - Assume an ideal inverter with 100 % DC-AC conversion efficiency
+    - Assume an ideal battery with "infinite" capacity and 100 % efficiency
+    - The inverter must try to follow the custom dispatch series as close as
+      possible
+    - Battery can only charge from PV
+    - The inverter is smart enough to charge the battery in order to avoid
+      clipping losses while still maintaining the MPP tracking
+    """
+    datasheet_battery_params.update(
+        {
+            "dc_energy_wh": 100000,
+            "dc_max_power_w": 850,
+            "charge_efficiency": 1.0,
+            "discharge_efficiency": 1.0,
+        }
+    )
+    inverter = {
+        'Vac': '240',
+        'Pso': 0.0,
+        'Paco': 1000.0,
+        'Pdco': 1000.0,
+        'Vdco': 325.0,
+        'C0': 0.0,
+        'C1': 0.0,
+        'C2': 0.0,
+        'C3': 0.0,
+        'Pnt': 0.5,
+        'Vdcmax': 600.0,
+        'Idcmax': 12.0,
+        'Mppt_low': 100.0,
+        'Mppt_high': 600.0,
+    }
+    dispatch = array([inputs["dispatch"]])
+    dispatch = Series(data=dispatch, index=gen_hourly_index(len(dispatch)))
+    result = multi_dc_battery(
+        v_dc=[array([400] * len(dispatch))],
+        p_dc=[array([inputs["pv_power"]])],
+        inverter=inverter,
+        battery_dispatch=dispatch,
+        battery_parameters=fit_boc(datasheet_battery_params),
+        battery_model=boc,
+    )
+    assert approx(result.iloc[0].to_dict(), nan_ok=True) == outputs

setup.py

@@ -47,13 +47,13 @@
                  'requests-mock', 'pytest-timeout', 'pytest-rerunfailures',
                  'pytest-remotedata']
 EXTRAS_REQUIRE = {
-    'optional': ['cython', 'ephem', 'netcdf4', 'nrel-pysam', 'numba',
+    'optional': ['cython', 'ephem', 'netcdf4', 'nrel-pysam >= 3.0.1', 'numba',
                  'pvfactors', 'siphon', 'statsmodels',
                  'cftime >= 1.1.1'],
     'doc': ['ipython', 'matplotlib', 'sphinx == 4.5.0',
             'pydata-sphinx-theme == 0.8.1', 'sphinx-gallery',
             'docutils == 0.15.2', 'pillow', 'netcdf4', 'siphon',
-            'sphinx-toggleprompt >= 0.0.5', 'pvfactors'],
+            'sphinx-toggleprompt >= 0.0.5', 'pvfactors', 'nrel-pysam >= 3.0.1'],
     'test': TESTS_REQUIRE
 }
 EXTRAS_REQUIRE['all'] = sorted(set(sum(EXTRAS_REQUIRE.values(), [])))