Fix numerical overflow in v_from_i for some conditions (#226)

* implement log lambertw calculation in v_from_i

* add test

* more robust

* add note to whatsnew

* add issue note to test

* move the .real

* remove scipy version in appveyor config

Author: wholmgren

appveyor.yml

@@ -27,7 +27,7 @@ install:
     - cmd: conda info -a
 
     # install depenencies
-    - cmd: conda create -n test_env --yes --quiet python=%PYTHON_VERSION% pip numpy scipy=0.16.0 pandas nose pytest pytz ephem numba
+    - cmd: conda create -n test_env --yes --quiet python=%PYTHON_VERSION% pip numpy scipy pandas nose pytest pytz ephem numba
     - cmd: activate test_env
     - cmd: python --version
     - cmd: conda list

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

@@ -65,6 +65,8 @@ Bug fixes
   rather than requiring that the solution for all points be within the
   error tolerance. Results in test scenarios changed by 1-10%.
   (:issue:`221`)
+* Fixed a numerical overflow error in pvsystem.singlediode's v_oc
+  determination for some combinations of parameters. (:issue:`225`)
 * dirint function yielded the wrong results for non-sea-level pressures.
   Fixed. (:issue:`212`)
 * Fixed a bug in the day angle calculation used by the 'spencer' and 'asce'

pvlib/pvsystem.py

@@ -1827,14 +1827,29 @@ def v_from_i(resistance_shunt, resistance_series, nNsVth, current,
     IL = photocurrent
     I = current
 
-    argW = I0 * Rsh / nNsVth * np.exp(Rsh *(-I + IL + I0) / nNsVth)
-    lambertwterm = lambertw(argW)
+    argW = I0 * Rsh / nNsVth * np.exp(Rsh * (-I + IL + I0) / nNsVth)
+    lambertwterm = lambertw(argW).real
 
-    # Eqn. 3 in Jain and Kapoor, 2004
+    # Calculate using log(argW) in case argW is really big
+    logargW = (np.log(I0) + np.log(Rsh) - np.log(nNsVth) +
+               Rsh * (-I + IL + I0) / nNsVth)
+
+    # Three iterations of Newton-Raphson method to solve
+    # w+log(w)=logargW. The initial guess is w=logargW. Where direct
+    # evaluation (above) results in NaN from overflow, 3 iterations
+    # of Newton's method gives approximately 8 digits of precision.
+    w = logargW
+    for i in range(0, 3):
+        w = w * (1 - np.log(w) + logargW) / (1 + w)
+    lambertwterm_log = w
 
+    lambertwterm = np.where(np.isfinite(lambertwterm), lambertwterm,
+                            lambertwterm_log)
+
+    # Eqn. 3 in Jain and Kapoor, 2004
     V = -I*(Rs + Rsh) + IL*Rsh - nNsVth*lambertwterm + I0*Rsh
 
-    return V.real
+    return V
 
 
 def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
@@ -1898,12 +1913,12 @@ def i_from_v(resistance_shunt, resistance_series, nNsVth, voltage,
     argW = (Rs*I0*Rsh *
             np.exp(Rsh*(Rs*(IL+I0)+V) / (nNsVth*(Rs+Rsh))) /
             (nNsVth*(Rs + Rsh)))
-    lambertwterm = lambertw(argW)
+    lambertwterm = lambertw(argW).real
 
     # Eqn. 4 in Jain and Kapoor, 2004
     I = -V/(Rs + Rsh) - (nNsVth/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh)
 
-    return I.real
+    return I
 
 
 def snlinverter(v_dc, p_dc, inverter):

pvlib/test/test_pvsystem.py

@@ -355,26 +355,34 @@ def test_PVSystem_calcparams_desoto(cec_module_params):
 @requires_scipy
 def test_v_from_i():
     output = pvsystem.v_from_i(20, .1, .5, 3, 6e-7, 7)
-    assert_allclose(7.5049875193450521, output, 5)
+    assert_allclose(7.5049875193450521, output, atol=1e-5)
 
 
 @requires_scipy
 def test_v_from_i_big():
     output = pvsystem.v_from_i(500, 10, 4.06, 0, 6e-10, 1.2)
-    assert_allclose(86.320000493521079, output, 5)
+    assert_allclose(86.320000493521079, output, atol=1e-5)
+
+
+@requires_scipy
+def test_v_from_i_bigger():
+    # 1000 W/m^2 on a Canadian Solar 220M with 20 C ambient temp
+    # github issue 225
+    output = pvsystem.v_from_i(190, 1.065, 2.89, 0, 7.05196029e-08, 10.491262)
+    assert_allclose(54.303958833791455, output, atol=1e-5)
 
 
 @requires_scipy
 def test_i_from_v():
     output = pvsystem.i_from_v(20, .1, .5, 40, 6e-7, 7)
-    assert_allclose(-299.746389916, output, 5)
+    assert_allclose(-299.746389916, output, atol=1e-5)
 
 
 @requires_scipy
 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)
+    assert_allclose(-299.746389916, output, atol=1e-5)
 
 
 @requires_scipy