Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 7 commits into from
May 7, 2021
Merged

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

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
b9e57f0
Add theoretical DWT test for scaling.py
jranalli Apr 22, 2021
73e7d08
Adjust test fixtures relative to outputs from corrected MATLAB version.
jranalli Apr 23, 2021
a6efc70
Implement fix to Python version of code. Some additional work needed …
jranalli Apr 23, 2021
979cbdb
Update docs
jranalli Apr 23, 2021
54928eb
Update reference to WVM in scaling.py and formatting thereof.
jranalli Apr 26, 2021
56a8928
Suggested update to the documentation language.
jranalli May 6, 2021
0053e59
docstring and whatsnew adjustments
cwhanse May 7, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions docs/sphinx/source/whatsnew/v0.9.0.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 scaling.py WVM. Tracks with fix in MATLAB ver.
(:issue:`1206`, :pull:`1213`)

Testing
~~~~~~~
Expand Down Expand Up @@ -148,3 +150,4 @@ Contributors
* Joshua Stein (:ghuser:`jsstein`)
* Tony Lorenzo (:ghuser:`alorenzo175`)
* Damjan Postolovski (:ghuser:`dpostolovski`)
* Joe Ranalli (:ghuser:`jranalli`)
68 changes: 40 additions & 28 deletions pvlib/scaling.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
----------
Expand All @@ -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
Expand All @@ -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
27 changes: 19 additions & 8 deletions pvlib/tests/test_scaling.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand Down Expand Up @@ -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,
Expand All @@ -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)
Expand Down