"""
The ``pvlib.ivtools.utils.py`` module contains utility functions related to
working with IV curves, or fitting equations to IV curve data.
"""
import numpy as np
import pandas as pd
# A small number used to decide when a slope is equivalent to zero
EPS = np.finfo('float').eps**(1/3)
def _numdiff(x, f):
"""
Compute first and second order derivative using possibly unequally
spaced data.
Parameters
----------
x : numeric
a numpy array of values of x
f : numeric
a numpy array of values of the function f for which derivatives are to
be computed. Must be the same length as x.
Returns
-------
df : numeric
a numpy array of len(x) containing the first derivative of f at each
point x except at the first 2 and last 2 points
df2 : numeric
a numpy array of len(x) containing the second derivative of f at each
point x except at the first 2 and last 2 points.
Notes
-----
``numdiff`` computes first and second order derivatives using a 5th order
formula that accounts for possibly unequally spaced data [1]_. Because a
5th order centered difference formula is used, ``numdiff`` returns NaNs
for the first 2 and last 2 points in the input vector for x. Ported from
PVLib Matlab [2]_.
References
----------
.. [1] M. K. Bowen, R. Smith, "Derivative formulae and errors for
non-uniformly spaced points", Proceedings of the Royal Society A, vol.
461 pp 1975 - 1997, July 2005. DOI: 10.1098/rpsa.2004.1430
.. [2] PVLib MATLAB https://github.com/sandialabs/MATLAB_PV_LIB
"""
n = len(f)
df = np.zeros(n)
df2 = np.zeros(n)
# first two points are special
df[:2] = float("Nan")
df2[:2] = float("Nan")
# Last two points are special
df[-2:] = float("Nan")
df2[-2:] = float("Nan")
# Rest of points. Take reference point to be the middle of each group of 5
# points. Calculate displacements
ff = np.vstack((f[:-4], f[1:-3], f[2:-2], f[3:-1], f[4:])).T
a0 = (np.vstack((x[:-4], x[1:-3], x[2:-2], x[3:-1], x[4:])).T
- np.tile(x[2:-2], [5, 1]).T)
u1 = np.zeros(a0.shape)
left = np.zeros(a0.shape)
u2 = np.zeros(a0.shape)
u1[:, 0] = (
a0[:, 1] * a0[:, 2] * a0[:, 3] + a0[:, 1] * a0[:, 2] * a0[:, 4]
+ a0[:, 1] * a0[:, 3] * a0[:, 4] + a0[:, 2] * a0[:, 3] * a0[:, 4])
u1[:, 1] = (
a0[:, 0] * a0[:, 2] * a0[:, 3] + a0[:, 0] * a0[:, 2] * a0[:, 4]
+ a0[:, 0] * a0[:, 3] * a0[:, 4] + a0[:, 2] * a0[:, 3] * a0[:, 4])
u1[:, 2] = (
a0[:, 0] * a0[:, 1] * a0[:, 3] + a0[:, 0] * a0[:, 1] * a0[:, 4]
+ a0[:, 0] * a0[:, 3] * a0[:, 4] + a0[:, 1] * a0[:, 3] * a0[:, 4])
u1[:, 3] = (
a0[:, 0] * a0[:, 1] * a0[:, 2] + a0[:, 0] * a0[:, 1] * a0[:, 4]
+ a0[:, 0] * a0[:, 2] * a0[:, 4] + a0[:, 1] * a0[:, 2] * a0[:, 4])
u1[:, 4] = (
a0[:, 0] * a0[:, 1] * a0[:, 2] + a0[:, 0] * a0[:, 1] * a0[:, 3]
+ a0[:, 0] * a0[:, 2] * a0[:, 3] + a0[:, 1] * a0[:, 2] * a0[:, 3])
left[:, 0] = (a0[:, 0] - a0[:, 1]) * (a0[:, 0] - a0[:, 2]) * \
(a0[:, 0] - a0[:, 3]) * (a0[:, 0] - a0[:, 4])
left[:, 1] = (a0[:, 1] - a0[:, 0]) * (a0[:, 1] - a0[:, 2]) * \
(a0[:, 1] - a0[:, 3]) * (a0[:, 1] - a0[:, 4])
left[:, 2] = (a0[:, 2] - a0[:, 0]) * (a0[:, 2] - a0[:, 1]) * \
(a0[:, 2] - a0[:, 3]) * (a0[:, 2] - a0[:, 4])
left[:, 3] = (a0[:, 3] - a0[:, 0]) * (a0[:, 3] - a0[:, 1]) * \
(a0[:, 3] - a0[:, 2]) * (a0[:, 3] - a0[:, 4])
left[:, 4] = (a0[:, 4] - a0[:, 0]) * (a0[:, 4] - a0[:, 1]) * \
(a0[:, 4] - a0[:, 2]) * (a0[:, 4] - a0[:, 3])
df[2:-2] = np.sum(-(u1 / left) * ff, axis=1)
# second derivative
u2[:, 0] = (
a0[:, 1] * a0[:, 2] + a0[:, 1] * a0[:, 3] + a0[:, 1] * a0[:, 4]
+ a0[:, 2] * a0[:, 3] + a0[:, 2] * a0[:, 4] + a0[:, 3] * a0[:, 4])
u2[:, 1] = (
a0[:, 0] * a0[:, 2] + a0[:, 0] * a0[:, 3] + a0[:, 0] * a0[:, 4]
+ a0[:, 2] * a0[:, 3] + a0[:, 2] * a0[:, 4] + a0[:, 3] * a0[:, 4])
u2[:, 2] = (
a0[:, 0] * a0[:, 1] + a0[:, 0] * a0[:, 3] + a0[:, 0] * a0[:, 4]
+ a0[:, 1] * a0[:, 3] + a0[:, 1] * a0[:, 3] + a0[:, 3] * a0[:, 4])
u2[:, 3] = (
a0[:, 0] * a0[:, 1] + a0[:, 0] * a0[:, 2] + a0[:, 0] * a0[:, 4]
+ a0[:, 1] * a0[:, 2] + a0[:, 1] * a0[:, 4] + a0[:, 2] * a0[:, 4])
u2[:, 4] = (
a0[:, 0] * a0[:, 1] + a0[:, 0] * a0[:, 2] + a0[:, 0] * a0[:, 3]
+ a0[:, 1] * a0[:, 2] + a0[:, 1] * a0[:, 4] + a0[:, 2] * a0[:, 3])
df2[2:-2] = 2. * np.sum(u2 * ff, axis=1)
return df, df2
[docs]def rectify_iv_curve(voltage, current, decimals=None):
"""
Sort the IV curve data, remove NaNs and negative
values, and combine points with duplicate voltage.
Parameters
----------
voltage : numeric [V]
current : numeric [A]
decimals : int or None, default None
number of decimal places to which voltage is rounded to remove
duplicated points. If None, no rounding is done.
Returns
-------
voltage : numeric [V]
current : numeric [A]
Notes
-----
``rectify_iv_curve`` ensures that the IV curve lies in the first quadrant
of the (voltage, current) plane. The returned IV curve:
* increases in voltage
* contains no negative current or voltage values
* contains no NaNs
* contains no points with duplicate voltage values. Where voltage
values are repeated, a single data point is substituted with current
equal to the average of current at duplicated voltages.
"""
df = pd.DataFrame(data=np.vstack((voltage, current)).T, columns=['v', 'i'])
# restrict to first quadrant
df.dropna(inplace=True)
df = df[(df['v'] >= 0) & (df['i'] >= 0)]
# sort pairs on voltage, then current
df = df.sort_values(by=['v', 'i'], ascending=[True, False])
# eliminate duplicate voltage points
if decimals is not None:
df['v'] = np.round(df['v'], decimals=decimals)
_, inv = np.unique(df['v'], return_inverse=True)
df.index = inv
# average current at each common voltage
df = df.groupby(by=inv).mean()
tmp = np.array(df).T
return tmp[0, ], tmp[1, ]
def _schumaker_qspline(x, y):
"""
Fit a quadratic spline which preserves monotonicity and
convexity in the data.
Parameters
----------
x : numeric
independent points between which the spline will interpolate.
y : numeric
dependent points between which the spline will interpolate.
Returns
-------
t : array
an ordered vector of knots, i.e., X values where the spline
changes coefficients. All values in x are used as knots.
The algorithm may insert additional knots between data points in x
where changes in convexity are indicated by the (numerical)
derivative. Consequently len(t) >= len(x).
c : array
a Nx3 matrix of coefficients where the kth row defines the quadratic
interpolant between t_k and t_(k+1), i.e., y = c[i, 0] *
(x - t_k)^2 + c[i, 1] * (x - t_k) + c[i, 2]
yhat : array
y values corresponding to the knots in t. Contains the original
data points, y, and also y values estimated from the spline at the
inserted knots.
kflag : array
a vector of len(t) of logicals, which are set to true for
elements of t that are knots inserted by the algorithm.
Notes
-----
Algorithm is taken from [1]_, which relies on prior work described in [2]_.
Ported from PVLib Matlab [3]_.
References
----------
.. [1] L. L. Schumaker, "On Shape Preserving Quadratic Spline
Interpolation", SIAM Journal on Numerical Analysis 20(4), August 1983,
pp 854 - 864
.. [2] M. H. Lam, "Monotone and Convex Quadratic Spline Interpolation",
Virginia Journal of Science 41(1), Spring 1990
.. [3] PVLib MATLAB https://github.com/sandialabs/MATLAB_PV_LIB
"""
# Make sure vectors are 1D arrays
x = x.flatten()
y = y.flatten()
n = x.size
# compute various values used by the algorithm: differences, length of line
# segments between data points, and ratios of differences.
delx = np.diff(x) # delx[i] = x[i + 1] - x[i]
dely = np.diff(y)
delta = dely / delx
# Calculate first derivative at each x value per [3]
s = np.zeros_like(x)
left = np.append(0.0, delta)
right = np.append(delta, 0.0)
pdelta = left * right
u = pdelta > 0
# [3], Eq. 9 for interior points
# fix tuning parameters in [2], Eq 9 at chi = .5 and eta = .5
s[u] = pdelta[u] / (0.5*left[u] + 0.5*right[u])
# [3], Eq. 7 for left endpoint
left_end = 2.0 * delta[0] - s[1]
if delta[0] * left_end > 0:
s[0] = left_end
# [3], Eq. 8 for right endpoint
right_end = 2.0 * delta[-1] - s[-2]
if delta[-1] * right_end > 0:
s[-1] = right_end
# determine knots. Start with initial points x
# [2], Algorithm 4.1 first 'if' condition of step 5 defines intervals
# which won't get internal knots
tests = s[:-1] + s[1:]
u = np.isclose(tests, 2.0 * delta, atol=EPS)
# u = true for an interval which will not get an internal knot
k = n + sum(~u) # total number of knots = original data + inserted knots
# set up output arrays
# knot locations, first n - 1 and very last (n + k) are original data
xk = np.zeros(k)
yk = np.zeros(k) # function values at knot locations
# logicals that will indicate where additional knots are inserted
flag = np.zeros(k, dtype=bool)
a = np.zeros((k, 3))
# structures needed to compute coefficients, have to be maintained in
# association with each knot
tmpx = x[:-1]
tmpy = y[:-1]
tmpx2 = x[1:]
tmps = s[:-1]
tmps2 = s[1:]
diffs = np.diff(s)
# structure to contain information associated with each knot, used to
# calculate coefficients
uu = np.zeros((k, 6))
uu[:(n - 1), :] = np.array([tmpx, tmpx2, tmpy, tmps, tmps2, delta]).T
# [2], Algorithm 4.1 subpart 1 of Step 5
# original x values that are left points of intervals without internal
# knots
# MATLAB differs from NumPy, boolean indices must be same size as
# array
xk[:(n - 1)][u] = tmpx[u]
yk[:(n - 1)][u] = tmpy[u]
# constant term for each polynomial for intervals without knots
a[:(n - 1), 2][u] = tmpy[u]
a[:(n - 1), 1][u] = s[:-1][u]
a[:(n - 1), 0][u] = 0.5 * diffs[u] / delx[u] # leading coefficients
# [2], Algorithm 4.1 subpart 2 of Step 5
# original x values that are left points of intervals with internal knots
xk[:(n-1)][~u] = tmpx[~u]
yk[:(n-1)][~u] = tmpy[~u]
aa = s[:-1] - delta
b = s[1:] - delta
sbar = np.zeros(k)
eta = np.zeros(k)
# will contain mapping from the left points of intervals containing an
# added knot to each inverval's internal knot value
xi = np.zeros(k)
t0 = aa * b >= 0
# first 'else' in Algorithm 4.1 Step 5
v = np.logical_and(~u, t0) # len(u) == (n - 1) always
q = np.sum(v) # number of this type of knot to add
if q > 0.:
xk[(n - 1):(n + q - 1)] = .5 * (tmpx[v] + tmpx2[v]) # knot location
uu[(n - 1):(n + q - 1), :] = np.array([tmpx[v], tmpx2[v], tmpy[v],
tmps[v], tmps2[v], delta[v]]).T
xi[:(n-1)][v] = xk[(n - 1):(n + q - 1)]
t1 = np.abs(aa) > np.abs(b)
w = np.logical_and(~u, ~v) # second 'else' in Algorithm 4.1 Step 5
w = np.logical_and(w, t1)
r = np.sum(w)
if r > 0.:
xk[(n + q - 1):(n + q + r - 1)] = tmpx2[w] + aa[w] * delx[w] / diffs[w]
uu[(n + q - 1):(n + q + r - 1), :] = np.array([tmpx[w], tmpx2[w],
tmpy[w], tmps[w],
tmps2[w], delta[w]]).T
xi[:(n - 1)][w] = xk[(n + q - 1):(n + q + r - 1)]
z = np.logical_and(~u, ~v) # last 'else' in Algorithm 4.1 Step 5
z = np.logical_and(z, ~w)
ss = np.sum(z)
if ss > 0.:
xk[(n + q + r - 1):(n + q + r + ss - 1)] = \
tmpx[z] + b[z] * delx[z] / diffs[z]
uu[(n + q + r - 1):(n + q + r + ss - 1), :] = \
np.array([tmpx[z], tmpx2[z], tmpy[z], tmps[z], tmps2[z],
delta[z]]).T
xi[:(n-1)][z] = xk[(n + q + r - 1):(n + q + r + ss - 1)]
# define polynomial coefficients for intervals with added knots
ff = ~u
sbar[:(n-1)][ff] = (
(2 * uu[:(n - 1), 5][ff] - uu[:(n-1), 4][ff])
+ (uu[:(n - 1), 4][ff] - uu[:(n-1), 3][ff])
* (xi[:(n - 1)][ff] - uu[:(n-1), 0][ff])
/ (uu[:(n - 1), 1][ff] - uu[:(n-1), 0][ff]))
eta[:(n-1)][ff] = (
(sbar[:(n - 1)][ff] - uu[:(n-1), 3][ff])
/ (xi[:(n - 1)][ff] - uu[:(n-1), 0][ff]))
sbar[(n - 1):(n + q + r + ss - 1)] = \
(2 * uu[(n - 1):(n + q + r + ss - 1), 5] -
uu[(n - 1):(n + q + r + ss - 1), 4]) + \
(uu[(n - 1):(n + q + r + ss - 1), 4] -
uu[(n - 1):(n + q + r + ss - 1), 3]) * \
(xk[(n - 1):(n + q + r + ss - 1)] -
uu[(n - 1):(n + q + r + ss - 1), 0]) / \
(uu[(n - 1):(n + q + r + ss - 1), 1] -
uu[(n - 1):(n + q + r + ss - 1), 0])
eta[(n - 1):(n + q + r + ss - 1)] = \
(sbar[(n - 1):(n + q + r + ss - 1)] -
uu[(n - 1):(n + q + r + ss - 1), 3]) / \
(xk[(n - 1):(n + q + r + ss - 1)] -
uu[(n - 1):(n + q + r + ss - 1), 0])
# constant term for polynomial for intervals with internal knots
a[:(n - 1), 2][~u] = uu[:(n - 1), 2][~u]
a[:(n - 1), 1][~u] = uu[:(n - 1), 3][~u]
a[:(n - 1), 0][~u] = 0.5 * eta[:(n - 1)][~u] # leading coefficient
a[(n - 1):(n + q + r + ss - 1), 2] = \
uu[(n - 1):(n + q + r + ss - 1), 2] + \
uu[(n - 1):(n + q + r + ss - 1), 3] * \
(xk[(n - 1):(n + q + r + ss - 1)] -
uu[(n - 1):(n + q + r + ss - 1), 0]) + \
.5 * eta[(n - 1):(n + q + r + ss - 1)] * \
(xk[(n - 1):(n + q + r + ss - 1)] -
uu[(n - 1):(n + q + r + ss - 1), 0]) ** 2.
a[(n - 1):(n + q + r + ss - 1), 1] = sbar[(n - 1):(n + q + r + ss - 1)]
a[(n - 1):(n + q + r + ss - 1), 0] = \
.5 * (uu[(n - 1):(n + q + r + ss - 1), 4] -
sbar[(n - 1):(n + q + r + ss - 1)]) / \
(uu[(n - 1):(n + q + r + ss - 1), 1] -
uu[(n - 1):(n + q + r + ss - 1), 0])
yk[(n - 1):(n + q + r + ss - 1)] = a[(n - 1):(n + q + r + ss - 1), 2]
xk[n + q + r + ss - 1] = x[n - 1]
yk[n + q + r + ss - 1] = y[n - 1]
flag[(n - 1):(n + q + r + ss - 1)] = True # these are all inserted knots
tmp = np.vstack((xk, a.T, yk, flag)).T
# sort output in terms of increasing x (original plus added knots)
tmp2 = tmp[tmp[:, 0].argsort(kind='mergesort')]
t = tmp2[:, 0]
outn = len(t)
c = tmp2[0:(outn - 1), 1:4]
yhat = tmp2[:, 4]
kflag = tmp2[:, 5]
return t, c, yhat, kflag