Skip to content

Instantly share code, notes, and snippets.

@wholmgren

wholmgren/brentq.py

Last active March 6, 2018 21:42
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save wholmgren/6abe1fe9d4fe0db6985335dec750be05 to your computer and use it in GitHub Desktop.
Save wholmgren/6abe1fe9d4fe0db6985335dec750be05 to your computer and use it in GitHub Desktop.
Download ZIP
preliminary results for numba'ed maximum power point code. non-numba code was developed in https://github.com/pvlib/pvlib-python/pull/409
Raw
brentq.py
# based on scipy/optimize/Zeros/brentq.c
# modified to python by
# William Holmgren william.holmgren@gmail.com @wholmgren
# University of Arizona, 2018
# /* Written by Charles Harris charles.harris@sdl.usu.edu */
#
# #include <math.h>
# #include "zeros.h"
#
# #define MIN(a, b) ((a) < (b) ? (a) : (b))
#
# /*
# At the top of the loop the situation is the following:
#
# 1. the root is bracketed between xa and xb
# 2. xa is the most recent estimate
# 3. xp is the previous estimate
# 4. |fp| < |fb|
#
# The order of xa and xp doesn't matter, but assume xp < xb. Then xa lies to
# the right of xp and the assumption is that xa is increasing towards the root.
# In this situation we will attempt quadratic extrapolation as long as the
# condition
#
# * |fa| < |fp| < |fb|
#
# is satisfied. That is, the function value is decreasing as we go along.
# Note the 4 above implies that the right inequlity already holds.
#
# The first check is that xa is still to the left of the root. If not, xb is
# replaced by xp and the erval reverses, with xb < xa. In this situation
# we will try linear erpolation. That this has happened is signaled by the
# equality xb == xp;
#
# The second check is that |fa| < |fb|. If this is not the case, we swap
# xa and xb and resort to bisection.
#
# */
def brentq(f, xa, xb, xtol, rtol, iter, params):
xpre = xa
xcur = xb
xblk = 0.
# fpre
# fcur
fblk = 0.
spre = 0.
scur = 0.
# sbis
# the tolerance is 2*delta */
# delta
# stry, dpre, dblk
# i
fpre = f(xpre, *params)
fcur = f(xcur, *params)
# params->funcalls = 2
if fpre*fcur > 0:
# params->error_num = SIGNERR;
return 0.
if fpre == 0:
return xpre
if fcur == 0:
return xcur
iterations = 0
for i in range(iter):
iterations += 1
if fpre*fcur < 0:
xblk = xpre;
fblk = fpre;
spre = scur = xcur - xpre
if abs(fblk) < abs(fcur):
xpre = xcur
xcur = xblk
xblk = xpre
fpre = fcur
fcur = fblk
fblk = fpre
delta = (xtol + rtol*abs(xcur))/2
sbis = (xblk - xcur)/2
if (fcur == 0 | (abs(sbis) < delta)):
return xcur
if (abs(spre) > delta) & (abs(fcur) < abs(fpre)):
if (xpre == xblk):
# interpolate
stry = -fcur*(xcur - xpre)/(fcur - fpre)
else:
# extrapolate
dpre = (fpre - fcur)/(xpre - xcur)
dblk = (fblk - fcur)/(xblk - xcur)
stry = (-fcur*(fblk*dblk - fpre*dpre)
/(dblk*dpre*(fblk - fpre)))
if 2*abs(stry) < min(abs(spre), 3*abs(sbis) - delta):
#good short step
spre = scur
scur = stry
else:
# bisect
spre = sbis
scur = sbis
else:
# bisect
spre = sbis
scur = sbis
xpre = xcur
fpre = fcur
if abs(scur) > delta:
xcur += scur
else:
xcur += delta if sbis > 0 else -delta
fcur = f(xcur, *params)
i += 1
# params->funcalls++;
# params->error_num = CONVERR;
return xcur
Raw
brentq_bishop.py
# based on scipy/optimize/Zeros/brentq.c
# modified to python by
# William Holmgren william.holmgren@gmail.com @wholmgren
# University of Arizona, 2018
# /* Written by Charles Harris charles.harris@sdl.usu.edu */
#
# #include <math.h>
# #include "zeros.h"
#
# #define MIN(a, b) ((a) < (b) ? (a) : (b))
#
# /*
# At the top of the loop the situation is the following:
#
# 1. the root is bracketed between xa and xb
# 2. xa is the most recent estimate
# 3. xp is the previous estimate
# 4. |fp| < |fb|
#
# The order of xa and xp doesn't matter, but assume xp < xb. Then xa lies to
# the right of xp and the assumption is that xa is increasing towards the root.
# In this situation we will attempt quadratic extrapolation as long as the
# condition
#
# * |fa| < |fp| < |fb|
#
# is satisfied. That is, the function value is decreasing as we go along.
# Note the 4 above implies that the right inequlity already holds.
#
# The first check is that xa is still to the left of the root. If not, xb is
# replaced by xp and the erval reverses, with xb < xa. In this situation
# we will try linear erpolation. That this has happened is signaled by the
# equality xb == xp;
#
# The second check is that |fa| < |fb|. If this is not the case, we swap
# xa and xb and resort to bisection.
#
# */
import numpy as np
from numpy import abs
from numpy import minimum as min
from numba import njit, float64, int32
@njit([float64(float64, float64, float64, float64,
float64, float64)])
def bishop88_gradp(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
"""root finders only need dp/dv"""
a = np.exp(vd / nNsVth)
b = 1.0 / resistance_shunt
i = photocurrent - saturation_current * (a - 1.0) - vd * b
v = vd - i * resistance_series
c = saturation_current * a / nNsVth
grad_i = - c - b # di/dvd
grad_v = 1.0 - grad_i * resistance_series # dv/dvd
# dp/dv = d(iv)/dv = v * di/dv + i
grad = grad_i / grad_v # di/dv
grad_p = v * grad + i # dp/dv
return grad_p
@njit([float64(float64, float64, float64, float64, int32,
float64, float64, float64,
float64, float64)])
def brentq_bishop(xa, xb, xtol, rtol, iter,
photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
xpre = xa
xcur = xb
xblk = 0.
fblk = 0.
spre = 0.
scur = 0.
# the tolerance is 2*delta */
fpre = bishop88_gradp(xpre, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
fcur = bishop88_gradp(xcur, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
if fpre*fcur > 0:
return 0.
if fpre == 0:
return xpre
if fcur == 0:
return xcur
for i in range(iter):
if fpre*fcur < 0:
xblk = xpre
fblk = fpre
scur = xcur - xpre
spre = scur
if abs(fblk) < abs(fcur):
xpre = xcur
xcur = xblk
xblk = xpre
fpre = fcur
fcur = fblk
fblk = fpre
delta = (xtol + rtol*abs(xcur))/2
sbis = (xblk - xcur)/2
if (fcur == 0 | (abs(sbis) < delta)):
return xcur
if (abs(spre) > delta) & (abs(fcur) < abs(fpre)):
if xpre == xblk:
# interpolate
stry = -fcur*(xcur - xpre)/(fcur - fpre)
else:
# extrapolate
dpre = (fpre - fcur)/(xpre - xcur)
dblk = (fblk - fcur)/(xblk - xcur)
stry = (-fcur*(fblk*dblk - fpre*dpre)
/(dblk*dpre*(fblk - fpre)))
if 2*abs(stry) < min(abs(spre), 3*abs(sbis) - delta):
# good short step
spre = scur
scur = stry
else:
# bisect
spre = sbis
scur = sbis
else:
# bisect
spre = sbis
scur = sbis
xpre = xcur
fpre = fcur
if abs(scur) > delta:
xcur += scur
else:
if sbis > 0:
xcur += delta
else:
xcur -= delta
fcur = bishop88_gradp(xcur, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
return xcur
Raw
numba_mpp_testing.py
from time import clock
import numpy as np
from numba import jit, njit, vectorize, float64, typeof
import pandas as pd
import scipy.optimize
import pvlib
from pvlib import pvsystem
# comment out @jit, @vectorize, and uncomment slow_vd_jit_vec if no numba
import brentq
# hangs if nopython=True
brentq_jit = jit(nopython=False)(brentq.brentq)
from brentq_bishop import brentq_bishop
#
# def jit_in_context(fun, d, *args, **kwargs):
# # to get an AST of fun, without reference to initial module
# import ast
# import inspect
# source = inspect.getsource(fun)
# tree = ast.parse(source)
# name = tree.body[0].name
#
# # eval the AST in context d
# code = compile(tree, "<string>", 'exec')
# eval(code, d)
#
# # jit the function in context d
# fun = d[name]
# jitted_fun = jit(fun, *args, **kwargs)
# d[name] = jitted_fun
#
# return jitted_fun
#
#
# def factory(f, g, *args, **kwargs):
# d = {}
# jit_in_context(f, d, *args, **kwargs)
# myfun = jit_in_context(g, d, *args, **kwargs)
# return myfun
@njit
def bishop88_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
"""
Explicit calculation single-diode-model (SDM) currents and voltages using
diode junction voltages [1].
[1] "Computer simulation of the effects of electrical mismatches in
photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
https://doi.org/10.1016/0379-6787(88)90059-2
:param numeric vd: diode voltages [V]
:param numeric photocurrent: photo-generated current [A]
:param numeric saturation_current: diode one reverse saturation current [A]
:param numeric resistance_series: series resitance [ohms]
:param numeric resistance_shunt: shunt resitance [ohms]
:param numeric nNsVth: product of thermal voltage ``Vth`` [V], diode
ideality factor ``n``, and number of series cells ``Ns``
:returns: tuple containing currents [A], voltages [V], power [W],
"""
a = np.exp(vd / nNsVth)
b = 1.0 / resistance_shunt
i = photocurrent - saturation_current * (a - 1.0) - vd * b
v = vd - i * resistance_series
retval = i, v, i*v
return retval
@njit([float64(float64, float64, float64, float64, float64, float64)])
def bishop88_gradp_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
"""root finders only need dp/dv"""
a = np.exp(vd / nNsVth)
b = 1.0 / resistance_shunt
i = photocurrent - saturation_current * (a - 1.0) - vd * b
v = vd - i * resistance_series
c = saturation_current * a / nNsVth
grad_i = - c - b # di/dvd
grad_v = 1.0 - grad_i * resistance_series # dv/dvd
# dp/dv = d(iv)/dv = v * di/dv + i
grad = grad_i / grad_v # di/dv
grad_p = v * grad + i # dp/dv
return grad_p
@njit
def est_voc_jit(photocurrent, saturation_current, nNsVth):
"""
Rough estimate of open circuit voltage useful for bounding searches for
``i`` of ``v`` when using :func:`~pvlib.pvsystem.singlediode`.
:param numeric photocurrent: photo-generated current [A]
:param numeric saturation_current: diode one reverse saturation current [A]
:param numeric nNsVth: product of thermal voltage ``Vth`` [V], diode
ideality factor ``n``, and number of series cells ``Ns``
:returns: rough estimate of open circuit voltage [V]
Calculating the open circuit voltage, :math:`V_{oc}`, of an ideal device
with infinite shunt resistance, :math:`R_{sh} \\to \\infty`, and zero series
resistance, :math:`R_s = 0`, yields the following equation [1]. As an
estimate of :math:`V_{oc}` it is useful as an upper bound for the bisection
method.
.. math::
V_{oc, est}=n Ns V_{th} \\log \\left( \\frac{I_L}{I_0} + 1 \\right)
[1] http://www.pveducation.org/pvcdrom/open-circuit-voltage
"""
return nNsVth * np.log(photocurrent / saturation_current + 1.0)
@vectorize([float64(float64, float64, float64, float64, float64, float64)], target='cpu')
def slow_vd_jit_vec(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth, voc_est):
"""
This is a slow but reliable way to find mpp.
"""
# collect args
args = (photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
# first bound the search using voc
vd = scipy.optimize.brentq(bishop88_gradp_jit, 0.0, voc_est, args)
return vd
@jit([float64(float64, float64, float64, float64, float64)])
def slow_mpp_jit(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
voc_est = est_voc_jit(photocurrent, saturation_current, nNsVth)
# root finder fails if bounds are both 0
nonzeros = voc_est != 0
IL_pos = photocurrent[nonzeros]
RSH_pos = resistance_shunt[nonzeros]
voc_est_pos = voc_est[nonzeros]
vd_pos = slow_vd_jit_vec(IL_pos, saturation_current, resistance_series,
RSH_pos, nNsVth, voc_est_pos)
vd = np.zeros_like(photocurrent)
vd[nonzeros] = vd_pos
mpp = bishop88_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
# guessing that some code is needed here to handle nans and/or
# differences between pandas/numpy nonzeros indexing
return mpp
#brentq_jit = factory(brentq.brentq, bishop88_gradp_jit)
@vectorize([float64(float64, float64, float64, float64, float64, float64)], target='cpu')
def slow_vd_jit_vec_brentq_jit(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth, voc_est):
"""
This is a slow but reliable way to find mpp.
"""
# collect args
args = (photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
# first bound the search using voc
vd = brentq_jit(bishop88_gradp_jit, 0.0, voc_est, 2e-12, 1e-15, 100, args)
return vd
@jit([float64(float64, float64, float64, float64, float64)])
def slow_mpp_jit_brentq_jit(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
voc_est = est_voc_jit(photocurrent, saturation_current, nNsVth)
# root finder fails if bounds are both 0
nonzeros = voc_est != 0
IL_pos = photocurrent[nonzeros]
RSH_pos = resistance_shunt[nonzeros]
voc_est_pos = voc_est[nonzeros]
vd_pos = slow_vd_jit_vec_brentq_jit(IL_pos, saturation_current, resistance_series,
RSH_pos, nNsVth, voc_est_pos)
vd = np.zeros_like(photocurrent)
vd[nonzeros] = vd_pos
mpp = bishop88_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
# guessing that some code is needed here to handle nans and/or
# differences between pandas/numpy nonzeros indexing
return mpp
@jit([float64(float64, float64, float64, float64, float64)])
def slow_mpp_jit_brentq_bishop(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
voc_est = est_voc_jit(photocurrent, saturation_current, nNsVth)
# root finder fails if bounds are both 0
nonzeros = voc_est != 0
IL_pos = photocurrent[nonzeros]
RSH_pos = resistance_shunt[nonzeros]
voc_est_pos = voc_est[nonzeros]
vd_pos = slow_vd_jit_vec_brentq_bishop(IL_pos, saturation_current, resistance_series,
RSH_pos, nNsVth, voc_est_pos)
vd = np.zeros_like(photocurrent)
vd[nonzeros] = vd_pos
mpp = bishop88_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
# guessing that some code is needed here to handle nans and/or
# differences between pandas/numpy nonzeros indexing
return mpp
@vectorize([float64(float64, float64, float64, float64, float64, float64)], nopython=True, target='cpu')
def slow_vd_jit_vec_brentq_bishop(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth, voc_est):
"""
This is a slow but reliable way to find mpp.
"""
# first bound the search using voc
vd = brentq_bishop(0.0, voc_est, 2e-12, 1e-15, 100,
photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
return vd
# comment out if using numba
#slow_vd_jit_vec = np.vectorize(slow_vd_jit_vec)
def prepare_data():
print('preparing single diode data from clear sky ghi...')
# adjust values to change length of test data
times = pd.DatetimeIndex(start='20180101', end='20190101', freq='15min', tz='America/Phoenix')
location = pvlib.location.Location(32.2, -110.9, altitude=710)
cs = location.get_clearsky(times)
poa_data = cs['ghi']
cec_modules = pvsystem.retrieve_sam('cecmod')
cec_module_params = cec_modules['Example_Module']
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)
return IL, I0, Rs, Rsh, nNsVth
# typeof_params = typeof((float64, float64, float64, float64, float64))
# typeof_f = typeof(bishop88_gradp_jit)
# brentq_jit = jit(
# brentq.brentq,
# [float64(typeof_f, float64, float64, float64, float64, float64, typeof_params)],
# nopython=True)
if __name__ == '__main__':
IL, I0, Rs, Rsh, nNsVth = prepare_data()
print('number of points = %s' % len(IL))
for n in range(4):
tstart = clock()
singlediode_out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth)
tstop = clock()
dt_slow = tstop - tstart
print('%s singlediode elapsed time = %g[s]' % (n, dt_slow))
for n in range(4):
tstart = clock()
i_mp, v_mp, p_mp = slow_mpp_jit(IL.values, I0, Rs, Rsh.values, nNsVth)
i_mp = pd.Series(i_mp, index=IL.index)
v_mp = pd.Series(v_mp, index=IL.index)
p_mp = pd.Series(p_mp, index=IL.index)
tstop = clock()
dt_slow = tstop - tstart
print('%s slow_mpp_jit elapsed time = %g[s]' % (n, dt_slow))
print('(singlediode - slow_mpp_jit).abs()\n',
(singlediode_out['p_mp'] - p_mp.fillna(0)).describe())
for n in range(4):
tstart = clock()
i_mp, v_mp, p_mp = slow_mpp_jit_brentq_bishop(IL.values, I0, Rs, Rsh.values, nNsVth)
i_mp = pd.Series(i_mp, index=IL.index)
v_mp = pd.Series(v_mp, index=IL.index)
p_mp = pd.Series(p_mp, index=IL.index)
tstop = clock()
dt_slow = tstop - tstart
print('%s slow_mpp_jit_brentq_bishop elapsed time = %g[s]' % (n, dt_slow))
print('(singlediode - slow_mpp_jit_brentq_bishop).abs()\n',
(singlediode_out['p_mp'] - p_mp.fillna(0)).describe())
for n in range(4):
tstart = clock()
i_mp, v_mp, p_mp = slow_mpp_jit_brentq_jit(IL.values, I0, Rs, Rsh.values, nNsVth)
i_mp = pd.Series(i_mp, index=IL.index)
v_mp = pd.Series(v_mp, index=IL.index)
p_mp = pd.Series(p_mp, index=IL.index)
tstop = clock()
dt_slow = tstop - tstart
print('%s slow_mpp_jit_brentq_jit elapsed time = %g[s]' % (n, dt_slow))
print('(singlediode - slow_mpp_jit_brentq_jit).abs()\n',
(singlediode_out['p_mp'] - p_mp.fillna(0)).describe())
# print(brentq_bishop.inspect_types())
Raw
numba_mpp_testing_clean.py
from time import clock
import numpy as np
from numba import jit, njit, vectorize, float64
import pandas as pd
import scipy.optimize
import pvlib
from pvlib import pvsystem
@njit
def bishop88_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
"""
Explicit calculation single-diode-model (SDM) currents and voltages using
diode junction voltages [1].
[1] "Computer simulation of the effects of electrical mismatches in
photovoltaic cell interconnection circuits" JW Bishop, Solar Cell (1988)
https://doi.org/10.1016/0379-6787(88)90059-2
:param numeric vd: diode voltages [V]
:param numeric photocurrent: photo-generated current [A]
:param numeric saturation_current: diode one reverse saturation current [A]
:param numeric resistance_series: series resitance [ohms]
:param numeric resistance_shunt: shunt resitance [ohms]
:param numeric nNsVth: product of thermal voltage ``Vth`` [V], diode
ideality factor ``n``, and number of series cells ``Ns``
:returns: tuple containing currents [A], voltages [V], power [W],
"""
a = np.exp(vd / nNsVth)
b = 1.0 / resistance_shunt
i = photocurrent - saturation_current * (a - 1.0) - vd * b
v = vd - i * resistance_series
retval = i, v, i*v
return retval
@njit([float64(float64, float64, float64, float64, float64, float64)])
def bishop88_gradp_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
"""root finders only need dp/dv"""
a = np.exp(vd / nNsVth)
b = 1.0 / resistance_shunt
i = photocurrent - saturation_current * (a - 1.0) - vd * b
v = vd - i * resistance_series
c = saturation_current * a / nNsVth
grad_i = - c - b # di/dvd
grad_v = 1.0 - grad_i * resistance_series # dv/dvd
# dp/dv = d(iv)/dv = v * di/dv + i
grad = grad_i / grad_v # di/dv
grad_p = v * grad + i # dp/dv
return grad_p
@njit
def est_voc_jit(photocurrent, saturation_current, nNsVth):
"""
Rough estimate of open circuit voltage useful for bounding searches for
``i`` of ``v`` when using :func:`~pvlib.pvsystem.singlediode`.
:param numeric photocurrent: photo-generated current [A]
:param numeric saturation_current: diode one reverse saturation current [A]
:param numeric nNsVth: product of thermal voltage ``Vth`` [V], diode
ideality factor ``n``, and number of series cells ``Ns``
:returns: rough estimate of open circuit voltage [V]
Calculating the open circuit voltage, :math:`V_{oc}`, of an ideal device
with infinite shunt resistance, :math:`R_{sh} \\to \\infty`, and zero series
resistance, :math:`R_s = 0`, yields the following equation [1]. As an
estimate of :math:`V_{oc}` it is useful as an upper bound for the bisection
method.
.. math::
V_{oc, est}=n Ns V_{th} \\log \\left( \\frac{I_L}{I_0} + 1 \\right)
[1] http://www.pveducation.org/pvcdrom/open-circuit-voltage
"""
return nNsVth * np.log(photocurrent / saturation_current + 1.0)
@vectorize([float64(float64, float64, float64, float64, float64, float64)], target='cpu')
def slow_vd_jit_vec(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth, voc_est):
"""
This is a slow but reliable way to find mpp.
"""
# collect args
args = (photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
# first bound the search using voc
vd = scipy.optimize.brentq(bishop88_gradp_jit, 0.0, voc_est, args)
return vd
@jit([float64(float64, float64, float64, float64, float64)])
def slow_mpp_jit(photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth):
voc_est = est_voc_jit(photocurrent, saturation_current, nNsVth)
# root finder fails if bounds are both 0
nonzeros = voc_est != 0
IL_pos = photocurrent[nonzeros]
RSH_pos = resistance_shunt[nonzeros]
voc_est_pos = voc_est[nonzeros]
vd_pos = slow_vd_jit_vec(IL_pos, saturation_current, resistance_series,
RSH_pos, nNsVth, voc_est_pos)
vd = np.zeros_like(photocurrent)
vd[nonzeros] = vd_pos
mpp = bishop88_jit(vd, photocurrent, saturation_current, resistance_series,
resistance_shunt, nNsVth)
# guessing that some code is needed here to handle nans and/or
# differences between pandas/numpy nonzeros indexing
return mpp
def prepare_data():
print('preparing single diode data from clear sky ghi...')
# adjust values to change length of test data
times = pd.DatetimeIndex(start='20180101', end='20190101', freq='1min', tz='America/Phoenix')
location = pvlib.location.Location(32.2, -110.9, altitude=710)
cs = location.get_clearsky(times)
poa_data = cs['ghi']
cec_modules = pvsystem.retrieve_sam('cecmod')
cec_module_params = cec_modules['Example_Module']
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)
return IL, I0, Rs, Rsh, nNsVth
if __name__ == '__main__':
IL, I0, Rs, Rsh, nNsVth = prepare_data()
print('number of points = %s' % len(IL))
for n in range(4):
tstart = clock()
singlediode_out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth)
tstop = clock()
dt_slow = tstop - tstart
print('%s singlediode elapsed time = %g[s]' % (n, dt_slow))
for n in range(4):
tstart = clock()
i_mp, v_mp, p_mp = slow_mpp_jit(IL.values, I0, Rs, Rsh.values, nNsVth)
i_mp = pd.Series(i_mp, index=IL.index)
v_mp = pd.Series(v_mp, index=IL.index)
p_mp = pd.Series(p_mp, index=IL.index)
tstop = clock()
dt_slow = tstop - tstart
print('%s slow_mpp_jit elapsed time = %g[s]' % (n, dt_slow))
print('(singlediode - slow_mpp_jit).abs()\n',
(singlediode_out['p_mp'] - p_mp.fillna(0)).describe())
Raw
numba_mpp_testing_clean_results.txt
(numba) 10:51:43 holmgren@holmgren-mbp.local numba_mmp_testing master ? python numba_mpp_testing_clean.py
preparing single diode data from clear sky ghi...
number of points = 525601
0 singlediode elapsed time = 3.32074[s]
1 singlediode elapsed time = 3.45139[s]
2 singlediode elapsed time = 3.30214[s]
3 singlediode elapsed time = 3.2279[s]
0 slow_mpp_jit elapsed time = 2.82168[s]
1 slow_mpp_jit elapsed time = 2.30675[s]
2 slow_mpp_jit elapsed time = 2.20894[s]
3 slow_mpp_jit elapsed time = 2.31824[s]
(singlediode - slow_mpp_jit).abs()
count 5.256010e+05
mean -9.334509e-05
std 2.109775e-04
min -1.557723e-03
25% -6.835588e-05
50% -2.328990e-08
75% 0.000000e+00
max 0.000000e+00
dtype: float64
(numba) 11:01:05 holmgren@holmgren-mbp.local numba_mmp_testing master ? NUMBA_DISABLE_JIT=1 python numba_mpp_testing_clean.py
preparing single diode data from clear sky ghi...
number of points = 525601
0 singlediode elapsed time = 3.12153[s]
1 singlediode elapsed time = 3.20733[s]
2 singlediode elapsed time = 3.1084[s]
3 singlediode elapsed time = 3.14283[s]
0 slow_mpp_jit elapsed time = 10.6972[s]
1 slow_mpp_jit elapsed time = 10.8363[s]
2 slow_mpp_jit elapsed time = 10.8462[s]
3 slow_mpp_jit elapsed time = 10.7065[s]
(singlediode - slow_mpp_jit).abs()
count 5.256010e+05
mean -9.334509e-05
std 2.109775e-04
min -1.557723e-03
25% -6.835588e-05
50% -2.328990e-08
75% 0.000000e+00
max 0.000000e+00
dtype: float64
Raw
numba_mpp_testing_results.txt
number of points = 731
0 slow_mpp_jit elapsed time = 0.267107[s]
1 slow_mpp_jit elapsed time = 0.010541[s]
2 slow_mpp_jit elapsed time = 0.011016[s]
3 slow_mpp_jit elapsed time = 0.011841[s]
0 singlediode elapsed time = 0.013187[s]
1 singlediode elapsed time = 0.011986[s]
2 singlediode elapsed time = 0.012876[s]
3 singlediode elapsed time = 0.013185[s]
(singlediode - slow_mpp_jit).abs()
count 731.000000
mean -0.000140
std 0.000268
min -0.001539
25% -0.000163
50% 0.000000
75% 0.000000
max 0.000000
dtype: float64
number of points = 8761
0 slow_mpp_jit elapsed time = 0.304554[s]
1 slow_mpp_jit elapsed time = 0.052078[s]
2 slow_mpp_jit elapsed time = 0.050663[s]
3 slow_mpp_jit elapsed time = 0.050045[s]
0 singlediode elapsed time = 0.060079[s]
1 singlediode elapsed time = 0.052729[s]
2 singlediode elapsed time = 0.053135[s]
3 singlediode elapsed time = 0.055275[s]
(singlediode - slow_mpp_jit).abs()
count 8.761000e+03
mean -9.067401e-05
std 2.059644e-04
min -1.549522e-03
25% -6.873591e-05
50% -2.308914e-08
75% 0.000000e+00
max 0.000000e+00
dtype: float64
number of points = 525601
0 slow_mpp_jit elapsed time = 2.46883[s]
1 slow_mpp_jit elapsed time = 2.3082[s]
2 slow_mpp_jit elapsed time = 2.19921[s]
3 slow_mpp_jit elapsed time = 2.33642[s]
0 singlediode elapsed time = 3.41547[s]
1 singlediode elapsed time = 3.64378[s]
2 singlediode elapsed time = 3.54097[s]
3 singlediode elapsed time = 3.193[s]
count 5.256010e+05
mean 9.334509e-05
std 2.109775e-04
min 0.000000e+00
25% 0.000000e+00
50% 2.328990e-08
75% 6.835588e-05
max 1.557723e-03
dtype: float64
number of points = 1051201
0 slow_mpp_jit elapsed time = 25.1237[s]
1 slow_mpp_jit elapsed time = 19.005[s]
2 slow_mpp_jit elapsed time = 18.9404[s]
3 slow_mpp_jit elapsed time = 18.9563[s]
0 singlediode elapsed time = 6.71949[s]
1 singlediode elapsed time = 6.79304[s]
2 singlediode elapsed time = 7.28941[s]
3 singlediode elapsed time = 7.42004[s]
(singlediode - slow_mpp_jit).abs()
count 1.051201e+06
mean -9.248689e-05
std 2.105602e-04
min -1.557463e-03
25% -6.626797e-05
50% 0.000000e+00
75% 0.000000e+00
max 6.313542e-06
dtype: float64
####################################################################
####################################################################
####################################################################
####################################################################
USE_NUMBA = False
number of points = 8761
0 slow_mpp_jit elapsed time = 0.163206[s]
1 slow_mpp_jit elapsed time = 0.167645[s]
2 slow_mpp_jit elapsed time = 0.159052[s]
3 slow_mpp_jit elapsed time = 0.166102[s]
0 singlediode elapsed time = 0.058287[s]
1 singlediode elapsed time = 0.054926[s]
2 singlediode elapsed time = 0.057167[s]
3 singlediode elapsed time = 0.054014[s]
(singlediode - slow_mpp_jit).abs()
count 8.761000e+03
mean -9.067401e-05
std 2.059644e-04
min -1.549522e-03
25% -6.873591e-05
50% -2.308914e-08
75% 0.000000e+00
max 0.000000e+00
dtype: float64
number of points = 525601
0 slow_mpp_jit elapsed time = 9.03347[s]
1 slow_mpp_jit elapsed time = 9.13156[s]
2 slow_mpp_jit elapsed time = 9.04388[s]
3 slow_mpp_jit elapsed time = 9.06983[s]
0 singlediode elapsed time = 3.38312[s]
1 singlediode elapsed time = 3.04524[s]
2 singlediode elapsed time = 3.03787[s]
3 singlediode elapsed time = 3.08195[s]
(singlediode - slow_mpp_jit).abs()
count 5.256010e+05
mean -9.334509e-05
std 2.109775e-04
min -1.557723e-03
25% -6.835588e-05
50% -2.328990e-08
75% 0.000000e+00
max 0.000000e+00
dtype: float64
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment