Fix to WVM methodology to track changes in MATLAB version. (#1213)

* Add theoretical DWT test for scaling.py

* Adjust test fixtures relative to outputs from corrected MATLAB version.

* Implement fix to Python version of code. Some additional work needed to be done here to get the wavelet levels to match the MATLAB output. The only place python and MATLAB don't agree now is at the edges, where there are different assumptions about how the moving average should be computed for the beginning/end of the time series. In principle, this could actually matter for time series that aren't "much" longer than 4096 seconds.

* Update docs

* Update reference to WVM in scaling.py and formatting thereof.

* Suggested update to the documentation language.

Co-authored-by: Cliff Hansen <[email protected]>

* docstring and whatsnew adjustments

Co-authored-by: Cliff Hansen <[email protected]>

Author: jranalli

docs/sphinx/source/whatsnew/v0.9.0.rst

@@ -117,6 +117,8 @@ Bug fixes
 * Reindl model fixed to generate sky_diffuse=0 when GHI=0.
   (:issue:`1153`, :pull:`1154`)
 * Update GFS product names for GFS v16. (:issue:`1202`, :pull:`1203`)
+* Corrected methodology error in :py:func:`~pvlib.scaling.wvm`. Tracks with
+  fix in PVLib for MATLAB. (:issue:`1206`, :pull:`1213`)
 
 Testing
 ~~~~~~~
@@ -126,7 +128,6 @@ Testing
 
 Documentation
 ~~~~~~~~~~~~~
-
 * Update intro tutorial to highlight the use of historical meteorological data
   and to make the procedural and OO results match exactly. (:issue:`1116`, :pull:`1144`)
 * Add a gallery example showing how to appropriately use interval-averaged
@@ -150,3 +151,4 @@ Contributors
 * Joshua Stein (:ghuser:`jsstein`)
 * Tony Lorenzo (:ghuser:`alorenzo175`)
 * Damjan Postolovski (:ghuser:`dpostolovski`)
+* Joe Ranalli (:ghuser:`jranalli`)

pvlib/scaling.py

@@ -13,8 +13,8 @@
 def wvm(clearsky_index, positions, cloud_speed, dt=None):
     """
     Compute spatial aggregation time series smoothing on clear sky index based
-    on the Wavelet Variability model of Lave et al [1-2]. Implementation is
-    basically a port of the Matlab version of the code [3].
+    on the Wavelet Variability model of Lave et al. [1]_, [2]_. Implementation
+    is basically a port of the Matlab version of the code [3]_.
 
     Parameters
     ----------
@@ -48,16 +48,16 @@ def wvm(clearsky_index, positions, cloud_speed, dt=None):
 
     References
     ----------
-    [1] M. Lave, J. Kleissl and J.S. Stein. A Wavelet-Based Variability
-    Model (WVM) for Solar PV Power Plants. IEEE Transactions on Sustainable
-    Energy, vol. 4, no. 2, pp. 501-509, 2013.
+    .. [1] M. Lave, J. Kleissl and J.S. Stein. A Wavelet-Based Variability
+       Model (WVM) for Solar PV Power Plants. IEEE Transactions on Sustainable
+       Energy, vol. 4, no. 2, pp. 501-509, 2013.
 
-    [2] M. Lave and J. Kleissl. Cloud speed impact on solar variability
-    scaling - Application to the wavelet variability model. Solar Energy,
-    vol. 91, pp. 11-21, 2013.
+    .. [2] M. Lave and J. Kleissl. Cloud speed impact on solar variability
+       scaling - Application to the wavelet variability model. Solar Energy,
+       vol. 91, pp. 11-21, 2013.
 
-    [3] Wavelet Variability Model - Matlab Code:
-    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
+    .. [3] Wavelet Variability Model - Matlab Code:
+       https://github.com/sandialabs/wvm
     """
 
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
@@ -128,13 +128,13 @@ def latlon_to_xy(coordinates):
 
     References
     ----------
-    [1] H. Moritz. Geodetic Reference System 1980, Journal of Geodesy, vol. 74,
-    no. 1, pp 128–133, 2000.
+    .. [1] H. Moritz. Geodetic Reference System 1980, Journal of Geodesy, vol.
+       74, no. 1, pp 128–133, 2000.
 
-    [2] https://pypi.org/project/pyproj/
+    .. [2] https://pypi.org/project/pyproj/
 
-    [3] Wavelet Variability Model - Matlab Code:
-    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
+    .. [3] Wavelet Variability Model - Matlab Code:
+       https://github.com/sandialabs/wvm
     """
 
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
@@ -159,7 +159,12 @@ def latlon_to_xy(coordinates):
 
 def _compute_wavelet(clearsky_index, dt=None):
     """
-    Compute the wavelet transform on the input clear_sky time series.
+    Compute the wavelet transform on the input clear_sky time series. Uses a
+    top hat wavelet [-1,1,1,-1] shape, based on the difference of successive
+    centered moving averages. Smallest scale (filter size of 2) is a degenerate
+    case that resembles a Haar wavelet. Returns one level of approximation
+    coefficient (CAn) and n levels of detail coefficients (CD1, CD2, ...,
+    CDn-1, CDn).
 
     Parameters
     ----------
@@ -174,19 +179,20 @@ def _compute_wavelet(clearsky_index, dt=None):
     Returns
     -------
     wavelet: numeric
-        The individual wavelets for the time series
+        The individual wavelets for the time series. Format follows increasing
+        scale (decreasing frequency): [CD1, CD2, ..., CDn, CAn]
 
     tmscales: numeric
         The timescales associated with the wavelets in seconds [s]
 
     References
     ----------
-    [1] M. Lave, J. Kleissl and J.S. Stein. A Wavelet-Based Variability
-    Model (WVM) for Solar PV Power Plants. IEEE Transactions on Sustainable
-    Energy, vol. 4, no. 2, pp. 501-509, 2013.
+    .. [1] M. Lave, J. Kleissl and J.S. Stein. A Wavelet-Based Variability
+       Model (WVM) for Solar PV Power Plants. IEEE Transactions on
+       Sustainable Energy, vol. 4, no. 2, pp. 501-509, 2013.
 
-    [3] Wavelet Variability Model - Matlab Code:
-    https://pvpmc.sandia.gov/applications/wavelet-variability-model/
+    .. [2] Wavelet Variability Model - Matlab Code:
+       https://github.com/sandialabs/wvm
     """
 
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2019
@@ -209,31 +215,37 @@ def _compute_wavelet(clearsky_index, dt=None):
 
     # Compute wavelet time scales
     min_tmscale = np.ceil(np.log(dt)/np.log(2))  # Minimum wavelet timescale
-    max_tmscale = int(12 - min_tmscale)  # maximum wavelet timescale
+    max_tmscale = int(13 - min_tmscale)  # maximum wavelet timescale
 
     tmscales = np.zeros(max_tmscale)
     csi_mean = np.zeros([max_tmscale, len(cs_long)])
+    # Skip averaging for the 0th scale
+    csi_mean[0, :] = cs_long.values.flatten()
+    tmscales[0] = 1
     # Loop for all time scales we will consider
-    for i in np.arange(0, max_tmscale):
-        j = i+1
-        tmscales[i] = 2**j * dt  # Wavelet integration time scale
-        intvlen = 2**j  # Wavelet integration time series interval
+    for i in np.arange(1, max_tmscale):
+        tmscales[i] = 2**i * dt  # Wavelet integration time scale
+        intvlen = 2**i  # Wavelet integration time series interval
         # Rolling average, retains only lower frequencies than interval
+        # Produces slightly different end effects than the MATLAB version
         df = cs_long.rolling(window=intvlen, center=True, min_periods=1).mean()
         # Fill nan's in both directions
         df = df.fillna(method='bfill').fillna(method='ffill')
         # Pop values back out of the dataframe and store
         csi_mean[i, :] = df.values.flatten()
+        # Shift to account for different indexing in MATLAB moving average
+        csi_mean[i, :] = np.roll(csi_mean[i, :], -1)
+        csi_mean[i, -1] = csi_mean[i, -2]
 
-    # Calculate the wavelets by isolating the rolling mean frequency ranges
+    # Calculate detail coefficients by difference between successive averages
     wavelet_long = np.zeros(csi_mean.shape)
     for i in np.arange(0, max_tmscale-1):
         wavelet_long[i, :] = csi_mean[i, :] - csi_mean[i+1, :]
-    wavelet_long[max_tmscale-1, :] = csi_mean[max_tmscale-1, :]  # Lowest freq
+    wavelet_long[-1, :] = csi_mean[-1, :]  # Lowest freq (CAn)
 
     # Clip off the padding and just return the original time window
     wavelet = np.zeros([max_tmscale, len(vals)])
     for i in np.arange(0, max_tmscale):
-        wavelet[i, :] = wavelet_long[i, len(vals)+1: 2*len(vals)+1]
+        wavelet[i, :] = wavelet_long[i, len(vals): 2*len(vals)]
 
     return wavelet, tmscales

pvlib/tests/test_scaling.py

@@ -48,21 +48,24 @@ def positions():
 @pytest.fixture
 def expect_tmscale():
     # Expected timescales for dt = 1
-    return [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+    return [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
 
 
 @pytest.fixture
 def expect_wavelet():
     # Expected wavelet for indices 5000:5004 for clear_sky_index above (Matlab)
-    return np.array([[-0.025, 0.05, 0., -0.05, 0.025],
-                     [0.025, 0., 0., 0., -0.025],
-                     [0., 0., 0., 0., 0.]])
+    e = np.zeros([13, 5])
+    e[0, :] = np.array([0, -0.05, 0.1, -0.05, 0])
+    e[1, :] = np.array([-0.025, 0.05, 0., -0.05, 0.025])
+    e[2, :] = np.array([0.025, 0., 0., 0., -0.025])
+    e[-1, :] = np.array([1, 1, 1, 1, 1])
+    return e
 
 
 @pytest.fixture
 def expect_cs_smooth():
     # Expected smoothed clear sky index for indices 5000:5004 (Matlab)
-    return np.array([1., 1.0289, 1., 0.9711, 1.])
+    return np.array([1., 1., 1.05774, 0.94226, 1.])
 
 
 def test_latlon_to_xy_zero():
@@ -94,7 +97,7 @@ def test_compute_wavelet_series(clear_sky_index, time,
     csi_series = pd.Series(clear_sky_index, index=time)
     wavelet, tmscale = scaling._compute_wavelet(csi_series)
     assert_almost_equal(tmscale, expect_tmscale)
-    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+    assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet)
 
 
 def test_compute_wavelet_series_numindex(clear_sky_index, time,
@@ -103,21 +106,29 @@ def test_compute_wavelet_series_numindex(clear_sky_index, time,
     csi_series = pd.Series(clear_sky_index, index=dtindex)
     wavelet, tmscale = scaling._compute_wavelet(csi_series)
     assert_almost_equal(tmscale, expect_tmscale)
-    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+    assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet)
 
 
 def test_compute_wavelet_array(clear_sky_index,
                                expect_tmscale, expect_wavelet):
     wavelet, tmscale = scaling._compute_wavelet(clear_sky_index, dt)
     assert_almost_equal(tmscale, expect_tmscale)
-    assert_almost_equal(wavelet[0:3, 5000:5005], expect_wavelet)
+    assert_almost_equal(wavelet[:, 5000:5005], expect_wavelet)
 
 
 def test_compute_wavelet_array_invalid(clear_sky_index):
     with pytest.raises(ValueError):
         scaling._compute_wavelet(clear_sky_index)
 
 
+def test_compute_wavelet_dwttheory(clear_sky_index, time,
+                                   expect_tmscale, expect_wavelet):
+    # Confirm detail coeffs sum to original signal
+    csi_series = pd.Series(clear_sky_index, index=time)
+    wavelet, tmscale = scaling._compute_wavelet(csi_series)
+    assert_almost_equal(np.sum(wavelet, 0), csi_series)
+
+
 def test_wvm_series(clear_sky_index, time, positions, expect_cs_smooth):
     csi_series = pd.Series(clear_sky_index, index=time)
     cs_sm, _, _ = scaling.wvm(csi_series, positions, cloud_speed)