add optional ivcurve calculation to singlediode

wholmgren · wholmgren · commit c69164a2a80d · 2016-07-22T08:45:15.000-07:00
diff --git a/docs/sphinx/source/whatsnew/v0.4.0.txt b/docs/sphinx/source/whatsnew/v0.4.0.txt
@@ -46,6 +46,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`)
 
 
 Bug fixes
diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
@@ -412,7 +412,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
@@ -424,7 +425,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):
@@ -1526,7 +1528,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.
 
@@ -1573,13 +1575,18 @@ 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.
+    If photocurrent is a Series and ivcurve_pnts is None, a DataFrame
+    with the columns described below. All columns have the same number
+    of rows as the largest input DataFrame.
 
-    If ``photocurrent`` is a scalar, a dict with the following keys.
+    If photocurrent is a scalar or ivcurve_pnts is not None, an
+    OrderedDict with the following keys.
 
     * i_sc -  short circuit current in amperes.
     * v_oc -  open circuit voltage in volts.
@@ -1588,6 +1595,8 @@ def singlediode(photocurrent, saturation_current, resistance_series,
     * 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)``.
+    * i - None or iv curve current.
+    * v - None or iv curve voltage.
 
     Notes
     -----
@@ -1641,37 +1650,32 @@ 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
+    # 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
+    else:
+        ivcurve_v = None
+        ivcurve_i = None
 
-    try:
-        dfout = pd.DataFrame(dfout, index=photocurrent.index)
-    except AttributeError:
-        pass
+    out = OrderedDict()
+    out['i_sc'] = i_sc
+    out['i_mp'] = i_mp
+    out['v_oc'] = v_oc
+    out['v_mp'] = v_mp
+    out['p_mp'] = p_mp
+    out['i_x'] = i_x
+    out['i_xx'] = i_xx
+    out['i'] = ivcurve_i
+    out['v'] = ivcurve_v
+
+    if isinstance(photocurrent, pd.Series) and not ivcurve_pnts:
+        out = pd.DataFrame(out, index=photocurrent.index)
 
-    return dfout
+    return out
 
 
 # Created April,2014
@@ -1872,11 +1876,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))) /
diff --git a/pvlib/test/test_pvsystem.py b/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,66 @@ 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-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=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,  6.00675648])),
+             ('i_mp', array([       nan,  5.6129056])),
+             ('v_oc', array([         nan,  10.29530483])),
+             ('v_mp', array([        nan,  7.25364707])),
+             ('p_mp', array([         nan,  40.71403625])),
+             ('i_x', array([        nan,  5.74622046])),
+             ('i_xx', array([        nan,  4.97138154])),
+             ('i',
+              array([[        nan,         nan,         nan],
+       [ 6.00726296,  5.74622046,  0.        ]])),
+             ('v',
+              array([[         nan,          nan,          nan],
+       [  0.        ,   5.14765242,  10.29530483]]))])
+
+    for k, v in out.items():
+        assert_allclose(expected[k], v, atol=3)
+
+
 def test_scale_voltage_current_power(sam_data):
     data = pd.DataFrame(
         np.array([[2, 1.5, 10, 8, 12, 0.5, 1.5]]),
