Source code for uncertainty_wrapper.tests.test_uncertainty_wrapper

"""
Tests for :func:`~uncertainty_wrapper.unc_wrapper` and
:func:`~uncertainty_wrapper.unc_wrapper_args`
"""

# import ok_, np, pd,unc_wrapper, unc_wrapper_args, KB, QE,
# pvlib, plt, UREG, PST and LOGGER from .tests
from uncertainty_wrapper.tests import *
from uncertainty_wrapper.tests.test_algopy import IV_algopy_jac, solpos_nd_jac

LOGGER = logging.getLogger(__name__)
LOGGER.setLevel(logging.DEBUG)


[docs]def test_unc_wrapper(): """ Test uncertainty wrapper with grouped arguments. """ x, cov = np.array([[1.0]]), np.array([[0.1]]) @unc_wrapper def f(y): return np.exp(y) avg, var, jac = f(x, __covariance__=cov, __method__='dense') LOGGER.debug("average = %g", avg) LOGGER.debug("variance = %g", var) ok_(np.isclose(avg, np.exp(x))) ok_(np.isclose(var, cov * np.exp(x) ** 2)) ok_(np.isclose(jac, np.exp(x))) return avg, var, jac
def test_unc_wrapper_args(): """ Test uncertainty wrapper with ungrouped arguments. """ x, cov = 1.0, np.array([[0.1]]) @unc_wrapper_args(None) def f(y): return np.exp(y) avg, var, jac = f(x, __covariance__=cov, __method__='dense') LOGGER.debug("average = %g", avg) LOGGER.debug("variance = %g", var) ok_(np.isclose(avg, np.exp(x))) ok_(np.isclose(var, cov * np.exp(x) ** 2)) ok_(np.isclose(jac, np.exp(x))) return avg, var, jac def test_multiple_observations(): """ Test multiple observations. """ x, cov = [1.0, 1.0], np.array([[[0.1]], [[0.1]]]) @unc_wrapper_args(None) def f(y): return np.exp(y).reshape(1, -1) avg, var, jac = f(x, __covariance__=cov, __method__='dense') LOGGER.debug("average = %r", avg) LOGGER.debug("variance = %r", var) ok_(np.allclose(avg, np.exp(x))) ok_(np.allclose(var, cov * np.exp(x) ** 2)) ok_(np.allclose(jac, np.exp(x))) return avg, var, jac def test_jagged(): """ Test jagged array. """ w, x = 1.0, [1.0, 1.0] cov = np.array([[[0.1, 0.], [0., 0.1]], [[0.1, 0.], [0., 0.1]]]) @unc_wrapper_args(None) def f(y, z): return (y * np.exp(z)).reshape(1, -1) avg, var, jac = f(w, x, __covariance__=cov, __method__='dense') LOGGER.debug("average = %r", avg) LOGGER.debug("jacobian = %r", jac) LOGGER.debug("variance = %r", var) ok_(np.allclose(avg, w * np.exp(x))) ok_(np.allclose(jac, [np.exp(x), np.exp(x)])) var_calc = np.concatenate( [[np.exp(x) * cov[:, 0, 0] + np.exp(x) * cov[:, 1, 0]], [np.exp(x) * cov[:, 0, 1] + np.exp(x) * cov[:, 1, 1]]], 0 ).reshape(2, 1, 2) var_calc = var_calc[:, 0, 0] * np.exp(x) + var_calc[:, 0, 1] * np.exp(x) ok_(np.allclose(var, var_calc)) return avg, var, jac def IV(x, Vd): """ Calculate IV curve using 2-diode model. :param x: independent variables: :type x: sequence :param Vd: diode voltages :type Vd: :class:`numpy.ndarray` :returns: current [A], voltage [V] and power [W] :rtype: :class:`numpy.ndarray` The sequence of independent variables must contain the following in the specified order:: [Ee, Tc, Rs, Rsh, Isat1_0, Isat2, Isc0, alpha_Isc, Eg] This function is an example of grouping the independent variables together so that :class:~`uncertianty_wrapper.core.unc_wrapper` can be used. """ Ee, Tc, Rs, Rsh, Isat1_0, Isat2, Isc0, alpha_Isc, Eg = x Vt = Tc * KB / QE Isc = Ee * Isc0 * (1.0 + (Tc - T0) * alpha_Isc) Isat1 = ( Isat1_0 * (Tc ** 3.0 / T0 ** 3.0) * np.exp(Eg * QE / KB * (1.0 / T0 - 1.0 / Tc)) ) Vd_sc = Isc * Rs # at short circuit Vc = 0 Id1_sc = Isat1 * (np.exp(Vd_sc / Vt) - 1.0) Id2_sc = Isat2 * (np.exp(Vd_sc / 2.0 / Vt) - 1.0) Ish_sc = Vd_sc / Rsh Iph = Isc + Id1_sc + Id2_sc + Ish_sc Id1 = Isat1 * (np.exp(Vd / Vt) - 1.0) Id2 = Isat2 * (np.exp(Vd / 2.0 / Vt) - 1.0) Ish = Vd / Rsh Ic = Iph - Id1 - Id2 - Ish Vc = Vd - Ic * Rs return np.array([Ic, Vc, Ic * Vc]) def Voc(x): """ Estimate open circuit voltage (Voc). """ Ee, Tc, Rs, Rsh, Isat1_0, Isat2, Isc0, alpha_Isc, Eg = x msg = ['Ee=%g[suns]','Tc=%g[K]','Rs=%g[ohms]','Rsh=%g[ohms]', 'Isat1_0=%g[A]','Isat2=%g[A]','Isc0=%g[A]','alpha_Isc=%g[]', 'Eg=%g[eV]'] LOGGER.debug('\n' + '\n'.join(msg) + '\n', *x) Vt = Tc * KB / QE LOGGER.debug('Vt=%g[V]', Vt) Isc = Ee * Isc0 * (1.0 + (Tc - T0) * alpha_Isc) LOGGER.debug('Isc=%g[A]', Isc) Isat1 = ( Isat1_0 * (Tc ** 3.0 / T0 ** 3.0) * np.exp(Eg * QE / KB * (1.0 / T0 - 1.0 / Tc)) ) LOGGER.debug('Isat1=%g[A]', Isat1) Vd_sc = Isc * Rs # at short circuit Vc = 0 Id1_sc = Isat1 * (np.exp(Vd_sc / Vt) - 1.0) Id2_sc = Isat2 * (np.exp(Vd_sc / 2.0 / Vt) - 1.0) Ish_sc = Vd_sc / Rsh Iph = Isc + Id1_sc + Id2_sc + Ish_sc # estimate Voc delta = Isat2 ** 2.0 + 4.0 * Isat1 * (Iph + Isat1 + Isat2) return Vt * np.log(((-Isat2 + np.sqrt(delta)) / 2.0 / Isat1) ** 2.0) # constants for IV test RS = 0.004267236774264931 # [ohm] series resistance RSH = 10.01226369025448 # [ohm] shunt resistance ISAT1_0 = 2.286188161253440E-11 # [A] diode one saturation current ISAT2 = 1.117455042372326E-6 # [A] diode two saturation current ISC0 = 6.3056 # [A] reference short circuit current EE = 0.8 # [suns] effective irradiance TC = 323.15 # [K] cell temperature EG = 1.1 # [eV] c-Si band gap ALPHA_ISC = 0.0003551 # [1/degC] short circuit current temp co # [V] open circuit voltage VOC = Voc((EE, TC, RS, RSH, ISAT1_0, ISAT2, ISC0, ALPHA_ISC, EG)) assert np.isclose(VOC, 0.62816490891656673) LOGGER.debug('Voc = %g[V]', VOC) VD = np.arange(0, VOC, 0.005) # [V] diode voltages X = np.array([EE, TC, RS, RSH, ISAT1_0, ISAT2, ISC0, ALPHA_ISC, EG]) X = X.reshape(-1, 1) # covariance equivalent to standard deviation of 1.0 [%] COV = np.diag([1e-4] * X.size) X_algopy = X.repeat(VD.size, axis=1)
[docs]def test_IV(method='sparse'): """ Test calculation of photovoltaic cell IV curve using 2-diode model and and compare Jacobian estimated by finite central difference to AlgoPy automatic differentiation. """ f = unc_wrapper(IV) pv, pv_cov, pv_jac = f(X, VD, __covariance__=COV, __method__=method) pv_cov = jflatten(pv_cov) pv_jac = jflatten(pv_jac) pv_jac_algopy = IV_algopy_jac(*X_algopy, Vd=VD) nVd = pv_jac_algopy.shape[1] for n in xrange(nVd // 2, nVd): irow, icol = 3 * n, 9 * n jrow, jcol = 3 + irow, 9 +icol pv_jac_n = pv_jac[irow:jrow, icol:jcol] pv_jac_algopy_n = pv_jac_algopy[:, n, n::VD.size] LOGGER.debug('pv jac at Vd = %g[V]:\n%r', VD[n], pv_jac_n) LOGGER.debug('pv jac AlgoPy at Vd = %g[V]:\n%r', VD[n], pv_jac_algopy_n) reldiff = pv_jac_n / pv_jac_algopy_n - 1.0 LOGGER.debug('reldiff at Vd = %g[V]:\n%r', VD[n], reldiff) resnorm = np.linalg.norm(reldiff) LOGGER.debug('resnorm at Vd = %g[V]: %r', VD[n], resnorm) rms = np.sqrt(np.sum(reldiff ** 2.0) / 9.0/ 3.0) LOGGER.debug('rms at Vd = %g[V]: %r', VD[n], rms) ok_(np.allclose(pv_jac_n, pv_jac_algopy_n, rtol=1e-3, atol=1e-3)) return pv, pv_cov, pv_jac, pv_jac_algopy
def plot_pv(pv, pv_cov): """ IV and PV 2-axis plot with errorbars """ i_pv, v_pv, p_pv = pv i_stdev = np.sqrt(pv_cov.diagonal()[::3]) v_stdev = np.sqrt(pv_cov.diagonal()[1::3]) p_stdev = np.sqrt(pv_cov.diagonal()[2::3]) fig, ax1 = plt.subplots() ax1.errorbar(v_pv, i_pv, i_stdev, v_stdev) ax1.grid() ax1.set_xlabel('voltage [V]') ax1.set_ylabel('current [A]', color='b') ax1.set_ylim([0, 6.0]) ax2 = ax1.twinx() ax2.errorbar(v_pv, p_pv, p_stdev, v_stdev, fmt='r') ax2.grid() ax2.set_ylabel('power [W]', color='r') ax2.set_ylim([0, 3.0]) ax1.set_title('IV and PV curves') return fig def plot_pv_jac(pv_jac, pv_jac_algopy, Vd=VD): """ Log plot of relative difference between AlgoPy and central finite difference approximations :param pv_jac: central finite approximations :param pv_jac_algopy: automatic differentiation :param Vd: voltages :return: fig """ fn = ['Cell Current, Ic [A]', 'Cell Voltage, Vc [V]', 'Cell Power, Pc [W]'] fig, ax = plt.subplots(3, 1, **{'figsize': (8.0, 18.0)}) colorcycle = [ 'firebrick', 'goldenrod', 'sage', 'lime', 'seagreen', 'turquoise', 'royalblue', 'indigo', 'fuchsia' ] for m in xrange(3): for n in xrange(9): pv_jac_n = pv_jac[m::3, n::9].diagonal() pv_jac_algopy_n = pv_jac_algopy[ m, :, n * 126:(n + 1) * 126 ].diagonal() reldiff = np.abs(pv_jac_n / pv_jac_algopy_n - 1.0) ax[m].semilogy(Vd, reldiff, colorcycle[n]) ax[m].grid() ax[m].legend( ['Ee', 'Tc', 'Rs', 'Rsh', 'Isat1_0', 'Isat2', 'Isc0', 'alpha_Isc', 'Eg'], fancybox=True, framealpha=0.5 ) ax[m].set_xlabel('Diode Voltage, Vd [V]') ax[m].set_ylabel('Relative Difference') ax[m].set_title(fn[m]) plt.tight_layout() return fig @UREG.wraps(('deg', 'deg', None, None), (None, 'deg', 'deg', 'Pa', 'm', 'degC')) @unc_wrapper_args(1, 2, 3, 4, 5) # indices specify positions of independent variables: # 1: latitude, 2: longitude, 3: pressure, 4: altitude, 5: temperature def spa(times, latitude, longitude, pressure, altitude, temperature): """ Calculate solar position using PVLIB Cython wrapper around NREL SPA. :param times: list of times, must be localized as UTC :type times: :class:`pandas.DatetimeIndex` :param latitude: latitude [deg] :param latitude: longitude [deg] :param pressure: pressure [Pa] :param latitude: altitude [m] :param temperature: temperature [degC] :returns: zenith, azimuth """ dataframe = pvlib.solarposition.spa_c(times, latitude, longitude, pressure, temperature) retvals = dataframe.to_records() zenith = retvals['apparent_zenith'] zenith = np.where(zenith < 90, zenith, np.nan) azimuth = retvals['azimuth'] return zenith, azimuth
[docs]def test_solpos(method='loop'): """ Test solar position calculation using NREL's SOLPOS. """ times = pd.DatetimeIndex(start='2015/1/1', end='2015/1/2', freq='1h', tz=PST).tz_convert(UTC) latitude, longitude = 37.0 * UREG.deg, -122.0 * UREG.deg pressure, temperature = 101325.0 * UREG.Pa, UREG.Quantity(22.0, UREG.degC) altitude = 0.0 * UREG.m # standard deviation of 1% assuming normal distribution covariance = np.diag([0.0001] * 5) ze, az, cov, jac = spa(times, latitude, longitude, pressure, altitude, temperature, __covariance__=covariance, __method__=method) cov = jflatten(cov) jac = jflatten(jac) jac_nd = solpos_nd_jac(times, latitude, longitude, pressure, altitude, temperature) for n in xrange(times.size): r, c = 2 * n, 5 * n # some rows which numdifftools returned nan if n in [0, 8, 17, 24]: continue ok_(np.allclose(jac[r:(r + 2), c:(c + 5)], jac_nd[n], equal_nan=True)) return ze, az, cov, jac, jac_nd
if __name__ == '__main__': test_unc_wrapper() pv, pv_cov, pv_jac, pv_jac_algopy = test_IV() test_solpos() fig1 = plot_pv(pv, pv_cov) fig1.show() fig2 = plot_pv_jac(pv_jac, pv_jac_algopy) fig2.savefig('IV-PV-jac-errors.png') fig2.show()