add optional ivcurve calculation to singlediode (#223)

* add optional ivcurve calculation to singlediode

* only include i v keys if needed

* update test for singlediode p_mp fix

* add a moderate point test sd series

Author: wholmgren

docs/sphinx/source/whatsnew/v0.4.0.txt

@@ -51,6 +51,8 @@ Enhancements
 * Add solarposition.nrel_earthsun_distance function and option to
   calculate extraterrestrial radiation using the NREL solar position
   algorithm. (:issue:`211`, :issue:`215`)
+* pvsystem.singlediode can now calculate IV curves if a user supplies
+  an ivcurve_pnts keyword argument. (:issue:`83`)
 * Includes SAM data files in the distribution. (:issue:`52`)
 
 

pvlib/pvsystem.py

@@ -413,7 +413,8 @@ def sapm_effective_irradiance(self, poa_direct, poa_diffuse,
             self.module_parameters, reference_irradiance=reference_irradiance)
 
     def singlediode(self, photocurrent, saturation_current,
-                    resistance_series, resistance_shunt, nNsVth):
+                    resistance_series, resistance_shunt, nNsVth,
+                    ivcurve_pnts=None):
         """Wrapper around the :py:func:`singlediode` function.
 
         Parameters
@@ -425,7 +426,8 @@ def singlediode(self, photocurrent, saturation_current,
         See pvsystem.singlediode for details
         """
         return singlediode(photocurrent, saturation_current,
-                           resistance_series, resistance_shunt, nNsVth)
+                           resistance_series, resistance_shunt, nNsVth,
+                           ivcurve_pnts=ivcurve_pnts)
 
     def i_from_v(self, resistance_shunt, resistance_series, nNsVth, voltage,
                  saturation_current, photocurrent):
@@ -1533,7 +1535,7 @@ def sapm_effective_irradiance(poa_direct, poa_diffuse, airmass_absolute, aoi,
 
 
 def singlediode(photocurrent, saturation_current, resistance_series,
-                resistance_shunt, nNsVth):
+                resistance_shunt, nNsVth, ivcurve_pnts=None):
     r'''
     Solve the single-diode model to obtain a photovoltaic IV curve.
 
@@ -1555,23 +1557,23 @@ def singlediode(photocurrent, saturation_current, resistance_series,
 
     Parameters
     ----------
-    photocurrent : float or Series
+    photocurrent : numeric
         Light-generated current (photocurrent) in amperes under desired
         IV curve conditions. Often abbreviated ``I_L``.
 
-    saturation_current : float or Series
+    saturation_current : numeric
         Diode saturation current in amperes under desired IV curve
         conditions. Often abbreviated ``I_0``.
 
-    resistance_series : float or Series
+    resistance_series : numeric
         Series resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rs``.
 
-    resistance_shunt : float or Series
+    resistance_shunt : numeric
         Shunt resistance in ohms under desired IV curve conditions.
         Often abbreviated ``Rsh``.
 
-    nNsVth : float or Series
+    nNsVth : numeric
         The product of three components. 1) The usual diode ideal factor
         (n), 2) the number of cells in series (Ns), and 3) the cell
         thermal voltage under the desired IV curve conditions (Vth). The
@@ -1580,21 +1582,35 @@ def singlediode(photocurrent, saturation_current, resistance_series,
         temp_cell is the temperature of the p-n junction in Kelvin, and
         q is the charge of an electron (coulombs).
 
+    ivcurve_pnts : None or int
+        Number of points in the desired IV curve. If None or 0, no
+        IV curves will be produced.
+
     Returns
     -------
-    If ``photocurrent`` is a Series, a DataFrame with the following
-    columns. All columns have the same number of rows as the largest
-    input DataFrame.
+    OrderedDict or DataFrame
+
+    The returned dict-like object always contains the keys/columns:
+
+        * i_sc - short circuit current in amperes.
+        * v_oc - open circuit voltage in volts.
+        * i_mp - current at maximum power point in amperes.
+        * v_mp - voltage at maximum power point in volts.
+        * p_mp - power at maximum power point in watts.
+        * i_x - current, in amperes, at ``v = 0.5*v_oc``.
+        * i_xx - current, in amperes, at ``V = 0.5*(v_oc+v_mp)``.
+
+    If ivcurve_pnts is greater than 0, the output dictionary will also
+    include the keys:
 
-    If ``photocurrent`` is a scalar, a dict with the following keys.
+        * i - IV curve current in amperes.
+        * v - IV curve voltage in volts.
 
-    * i_sc -  short circuit current in amperes.
-    * v_oc -  open circuit voltage in volts.
-    * i_mp -  current at maximum power point in amperes.
-    * v_mp -  voltage at maximum power point in volts.
-    * p_mp -  power at maximum power point in watts.
-    * i_x -  current, in amperes, at ``v = 0.5*v_oc``.
-    * i_xx -  current, in amperes, at ``V = 0.5*(v_oc+v_mp)``.
+    The output will be an OrderedDict if photocurrent is a scalar,
+    array, or ivcurve_pnts is not None.
+
+    The output will be a DataFrame if photocurrent is a Series and
+    ivcurve_pnts is None.
 
     Notes
     -----
@@ -1648,37 +1664,29 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     i_xx = i_from_v(resistance_shunt, resistance_series, nNsVth,
                     0.5*(v_oc+v_mp), saturation_current, photocurrent)
 
-    # @wholmgren: need to move this stuff to a different function
-#     If the user says they want a curve of with number of points equal to
-#     NumPoints (must be >=2), then create a voltage array where voltage is
-#     zero in the first column, and Voc in the last column. Number of columns
-#     must equal NumPoints. Each row represents the voltage for one IV curve.
-#     Then create a current array where current is Isc in the first column, and
-#     zero in the last column, and each row represents the current in one IV
-#     curve. Thus the nth (V,I) point of curve m would be found as follows:
-#     (Result.V(m,n),Result.I(m,n)).
-#     if NumPoints >= 2
-#        s = ones(1,NumPoints);  # shaping DataFrame to shape the column
-#                                # DataFrame parameters into 2-D matrices
-#        Result.V = (Voc)*(0:1/(NumPoints-1):1);
-#        Result.I = I_from_V(Rsh*s, Rs*s, nNsVth*s, Result.V, I0*s, IL*s);
-#     end
-
-    dfout = {}
-    dfout['i_sc'] = i_sc
-    dfout['i_mp'] = i_mp
-    dfout['v_oc'] = v_oc
-    dfout['v_mp'] = v_mp
-    dfout['p_mp'] = p_mp
-    dfout['i_x'] = i_x
-    dfout['i_xx'] = i_xx
-
-    try:
-        dfout = pd.DataFrame(dfout, index=photocurrent.index)
-    except AttributeError:
-        pass
+    out = OrderedDict()
+    out['i_sc'] = i_sc
+    out['v_oc'] = v_oc
+    out['i_mp'] = i_mp
+    out['v_mp'] = v_mp
+    out['p_mp'] = p_mp
+    out['i_x'] = i_x
+    out['i_xx'] = i_xx
+
+    # create ivcurve
+    if ivcurve_pnts:
+        ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *
+                     np.linspace(0, 1, ivcurve_pnts))
+        ivcurve_i = i_from_v(
+            resistance_shunt, resistance_series, nNsVth, ivcurve_v.T,
+            saturation_current, photocurrent).T
+        out['v'] = ivcurve_v
+        out['i'] = ivcurve_i
+
+    if isinstance(photocurrent, pd.Series) and not ivcurve_pnts:
+        out = pd.DataFrame(out, index=photocurrent.index)
 
-    return dfout
+    return out
 
 
 # Created April,2014
@@ -1879,11 +1887,13 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     except ImportError:
         raise ImportError('This function requires scipy')
 
-    Rsh = resistance_shunt
-    Rs = resistance_series
-    I0 = saturation_current
-    IL = photocurrent
-    V = voltage
+    # asarray turns Series into arrays so that we don't have to worry
+    # about multidimensional broadcasting failing
+    Rsh = np.asarray(resistance_shunt)
+    Rs = np.asarray(resistance_series)
+    I0 = np.asarray(saturation_current)
+    IL = np.asarray(photocurrent)
+    V = np.asarray(voltage)
 
     argW = (Rs*I0*Rsh *
             np.exp(Rsh*(Rs*(IL+I0)+V) / (nNsVth*(Rs+Rsh))) /

pvlib/test/test_pvsystem.py

@@ -4,7 +4,7 @@
 from collections import OrderedDict
 
 import numpy as np
-from numpy import nan
+from numpy import nan, array
 import pandas as pd
 
 import pytest
@@ -150,6 +150,14 @@ def sapm_module_params(sam_data):
     return module_parameters
 
 
+@pytest.fixture(scope="session")
+def cec_module_params(sam_data):
+    modules = sam_data['cecmod']
+    module = 'Example_Module'
+    module_parameters = modules[module]
+    return module_parameters
+
+
 def test_sapm(sapm_module_params):
 
     times = pd.DatetimeIndex(start='2015-01-01', periods=5, freq='12H')
@@ -307,17 +315,15 @@ def test_PVSystem_sapm_effective_irradiance(sapm_module_params):
         aoi, reference_irradiance=reference_irradiance)
 
 
-def test_calcparams_desoto(sam_data):
-    module = 'Example_Module'
-    module_parameters = sam_data['cecmod'][module]
+def test_calcparams_desoto(cec_module_params):
     times = pd.DatetimeIndex(start='2015-01-01', periods=2, freq='12H')
     poa_data = pd.Series([0, 800], index=times)
 
     IL, I0, Rs, Rsh, nNsVth = pvsystem.calcparams_desoto(
                                   poa_data,
                                   temp_cell=25,
-                                  alpha_isc=module_parameters['alpha_sc'],
-                                  module_parameters=module_parameters,
+                                  alpha_isc=cec_module_params['alpha_sc'],
+                                  module_parameters=cec_module_params,
                                   EgRef=1.121,
                                   dEgdT=-0.0002677)
 
@@ -328,13 +334,11 @@ def test_calcparams_desoto(sam_data):
     assert_allclose(nNsVth, 0.473)
 
 
-def test_PVSystem_calcparams_desoto(sam_data):
-    module = 'Example_Module'
-    module_parameters = sam_data['cecmod'][module].copy()
+def test_PVSystem_calcparams_desoto(cec_module_params):
+    module_parameters = cec_module_params.copy()
     module_parameters['EgRef'] = 1.121
     module_parameters['dEgdT'] = -0.0002677
-    system = pvsystem.PVSystem(module=module,
-                               module_parameters=module_parameters)
+    system = pvsystem.PVSystem(module_parameters=module_parameters)
     times = pd.DatetimeIndex(start='2015-01-01', periods=2, freq='12H')
     poa_data = pd.Series([0, 800], index=times)
     temp_cell = 25
@@ -367,26 +371,21 @@ def test_i_from_v():
 
 
 @requires_scipy
-def test_PVSystem_i_from_v(sam_data):
-    module = 'Example_Module'
-    module_parameters = sam_data['cecmod'][module]
-    system = pvsystem.PVSystem(module=module,
-                               module_parameters=module_parameters)
+def test_PVSystem_i_from_v():
+    system = pvsystem.PVSystem()
     output = system.i_from_v(20, .1, .5, 40, 6e-7, 7)
     assert_allclose(-299.746389916, output, 5)
 
 
 @requires_scipy
-def test_singlediode_series(sam_data):
-    module = 'Example_Module'
-    module_parameters = sam_data['cecmod'][module]
+def test_singlediode_series(cec_module_params):
     times = pd.DatetimeIndex(start='2015-01-01', periods=2, freq='12H')
     poa_data = pd.Series([0, 800], index=times)
     IL, I0, Rs, Rsh, nNsVth = pvsystem.calcparams_desoto(
                                          poa_data,
                                          temp_cell=25,
-                                         alpha_isc=module_parameters['alpha_sc'],
-                                         module_parameters=module_parameters,
+                                         alpha_isc=cec_module_params['alpha_sc'],
+                                         module_parameters=cec_module_params,
                                          EgRef=1.121,
                                          dEgdT=-0.0002677)
     out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth)
@@ -424,31 +423,68 @@ def test_singlediode_floats(sam_data):
                 'p_mp': 38.194165464983037,
                 'i_x': 6.7556075876880621,
                 'i_sc': 6.9646747613963198,
-                'v_mp': 6.221535886625464}
+                'v_mp': 6.221535886625464,
+                'i': None,
+                'v': None}
     assert isinstance(out, dict)
     for k, v in out.items():
-        assert_allclose(expected[k], v, atol=3)
+        if k in ['i', 'v']:
+            assert v is None
+        else:
+            assert_allclose(expected[k], v, atol=3)
 
 
 @requires_scipy
-def test_PVSystem_singlediode_floats(sam_data):
-    module = 'Example_Module'
-    module_parameters = sam_data['cecmod'][module]
-    system = pvsystem.PVSystem(module=module,
-                               module_parameters=module_parameters)
-    out = system.singlediode(7, 6e-7, .1, 20, .5)
+def test_singlediode_floats_ivcurve():
+    out = pvsystem.singlediode(7, 6e-7, .1, 20, .5, ivcurve_pnts=3)
     expected = {'i_xx': 4.2685798754011426,
                 'i_mp': 6.1390251797935704,
                 'v_oc': 8.1063001465863085,
                 'p_mp': 38.194165464983037,
                 'i_x': 6.7556075876880621,
                 'i_sc': 6.9646747613963198,
-                'v_mp': 6.221535886625464}
+                'v_mp': 6.221535886625464,
+                'i': np.array([6.965172e+00,   6.755882e+00,   2.575717e-14]),
+                'v': np.array([0.     ,  4.05315,  8.1063])}
     assert isinstance(out, dict)
     for k, v in out.items():
         assert_allclose(expected[k], v, atol=3)
 
 
+@requires_scipy
+def test_singlediode_series_ivcurve(cec_module_params):
+    times = pd.DatetimeIndex(start='2015-06-01', periods=3, freq='6H')
+    poa_data = pd.Series([0, 400, 800], index=times)
+    IL, I0, Rs, Rsh, nNsVth = pvsystem.calcparams_desoto(
+                                         poa_data,
+                                         temp_cell=25,
+                                         alpha_isc=cec_module_params['alpha_sc'],
+                                         module_parameters=cec_module_params,
+                                         EgRef=1.121,
+                                         dEgdT=-0.0002677)
+
+    out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth, ivcurve_pnts=3)
+
+    expected = OrderedDict([('i_sc', array([        nan,  3.01054475,  6.00675648])),
+             ('v_oc', array([         nan,   9.96886962,  10.29530483])),
+             ('i_mp', array([        nan,  2.65191983,  5.28594672])),
+             ('v_mp', array([        nan,  8.33392491,  8.4159707 ])),
+             ('p_mp', array([         nan,  22.10090078,  44.48637274])),
+             ('i_x', array([        nan,  2.88414114,  5.74622046])),
+             ('i_xx', array([        nan,  2.04340914,  3.90007956])),
+             ('v',
+              array([[         nan,          nan,          nan],
+       [  0.        ,   4.98443481,   9.96886962],
+       [  0.        ,   5.14765242,  10.29530483]])),
+             ('i',
+              array([[             nan,              nan,              nan],
+       [  3.01079860e+00,   2.88414114e+00,   3.10862447e-14],
+       [  6.00726296e+00,   5.74622046e+00,   0.00000000e+00]]))])
+
+    for k, v in out.items():
+        assert_allclose(expected[k], v, atol=1e-2)
+
+
 def test_scale_voltage_current_power(sam_data):
     data = pd.DataFrame(
         np.array([[2, 1.5, 10, 8, 12, 0.5, 1.5]]),