Pvlib-python: implement transient moving-average thermal model

[mikofski]

New Feature
Implement the transient moving-average model proposed by Matt Prillman

Describe the solution you'd like
a new function pvlib.temperature.prillman(times, ...)

Describe alternatives you've considered
see also #1080

Additional context
see https://github.com/NREL/ssc/issues/261

[kanderso-nrel]

@mjprilliman FYI

[kanderso-nrel]

I've taken a look at implementing this but am not satisfied with what I've come up with so far. Here are my notes in case anyone else is working on this.

Option 1: cell_temp_ss.rolling('1200s').apply(_prilliman_average)

My natural starting point was to use pandas to get the offset-based windows and apply a custom averaging function. The problem is that tseries.rolling(...).apply(...) is way too slow -- ~1 second computation time for three days of 1-minute data, even without a custom averaging function:

In [50]: N=60*24*3  # three days of 1-min data
    ...: s = pd.Series(range(N), pd.date_range('2019-01-01', freq='1T', periods=N))
    ...: %time _=s.rolling('1200s').apply(np.mean)
    ...: 
CPU times: user 1.08 s, sys: 2.41 ms, total: 1.08 s
Wall time: 1.02 s

It's also a little clumsy giving the averaging function access to the P time series. The best way I found involved lots of extra indexing and made it even slower.

Option 2: pandas EWM

pandas does provide exponential-weighted functions via s.ewm, but as far as I can tell it considers the entire history of the time series; there is no way to do a moving window (https://github.com/pandas-dev/pandas/issues/11288). That's maybe an acceptable error because the exponential weighting would minimize the contribution of points outside the window, but according to its docstring, the weight must be a float, i.e. it is a constant, which is not the case for the Prilliman et al model (Eq. 11; P depends on the wind speed time series). So I don't think the pandas EWM functions are flexible enough for this.

Option 3: constrain to regular time series

I think the inefficiency of tseries.rolling(...).apply(...) can be avoided if the input time series index is evenly spaced -- something like the Hankel matrix indexing approach of the pre-#1074 detect_clearsky. But it's a shame to require even spacing, especially when the paper specifically highlights the algorithm's ability to handle data with nonuniform sampling.

Option 4: accelerate it using numba, or maybe cython

Haven't tried it but I'm guessing this is the answer, in coordination with #1060. Like in #1037, I can't find the trick to get these transient models running efficiently with base pandas/numpy.

[mjprilliman]

Hi Kevin,

I'd like to help with this implementation. I had previously written a standalone Python script that was running the model accurately but I will have to revisit to see how or if it can tie in with the existing framework. I am working on SAM features this week but I can start working on this next week (possibly late this week) if that works.

[mjprilliman]

def pvl_MAmodel_2(Tmstamp, SS, WS, m_u, a=np.array([.0046,4.5537e-4,-2.2586e-4,-1.5661e-5])):
    #Numpy array MA model
    #convert datetime to datenum
    Tmstamp['Datenum'] = Tmstamp[['Times']].apply(pd.to_numeric)/(1e9)
    Tmstamp2 = np.array([Tmstamp.Datenum])

    #Separate bilinear interpolation coefficients into 4 variables
    a0 = a[0]
    a1 = a[1]
    a2 = a[2]
    a3 = a[3]
    #P = pd.DataFrame(a0+a1*WS+a2*m_u+a3*m_u*WS,dtype=float) # Power parameter for exponential weighting function
    P = np.array(a0+a1*WS+a2*m_u+a3*m_u*WS)
    #Initialize Results variable T_MA
    #T_MA = pd.DataFrame(SS.values[0]) #Initialize Moving Average with 
    #T_MA = np.arange(len(Tmstamp))
    #T_MA = np.full_like(T_MA,np.nan,dtype=np.double)
    #T_MA = np.array(SS[0,0])
    T_MA = np.array(SS[0])
    #T_MA[0] = SS[0,0]

    #Set constant window length
    WindowLength = 20*60
    #Initialize Back Counter
    I_B = 0

    cntr1 = 1
    for i in np.linspace(1,len(Tmstamp)-1,len(Tmstamp)-1):

        #Set cntr1 to iterative i to match Matlab code
        I_F = cntr1 - 1 #Front indice always set one timestep behind current value
        deltaT_I_F = Tmstamp2[0,cntr1] - Tmstamp2[0,I_F] #Time difference between current time and front indice (in case time series isn't uniform)
        if deltaT_I_F > WindowLength:
            #T_MA[cntr1] = SS[0,cntr1] #if the front indice is more than 20 minutes behind the current value, just use the steady-state approximation
            #T_MA = np.append(T_MA,SS[0,cntr1])
            T_MA = np.append(T_MA,SS[cntr1])

        else:
            while((Tmstamp2[0,cntr1] - Tmstamp2[0,I_B]) > WindowLength) & (I_B<I_F):
                     I_B = I_B+1
                     #If the front indice is within 20 minutes and the back indice is not going to bump into the front, bump up the back indice by one until it is within 20 minutes (window length)
            TimeBack = (Tmstamp2[0,cntr1] - Tmstamp2[0,I_B:I_F+1]) #Calculate the time in seconds back from the current value for each value between the front and back indices
            #TempsInWindow = SS[0,I_B:I_F+1] #find the steady-state temperature approximations for the times within the indices (include front indice)
            TempsInWindow = SS[I_B:I_F+1]
            #Weight = np.exp(-P[0,cntr1]*TimeBack.astype(float)) #Calculate the weight of each element in indice as function of power parameter P = fcn(WS, unit mass) and time back          
            Weight = np.exp(-P[cntr1]*TimeBack) #Calculate the weight of each element in indice as function of power parameter P = fcn(WS, unit mass) and time back          

            Rel_weight = Weight/np.sum(Weight) #Calculate weights relative to the total (weighted average)
            TempsWeighted = (TempsInWindow)*Rel_weight;
            Temp = np.sum(TempsWeighted);
            #Temp = pd.DataFrame(TempsInWindow.dot(np.transpose(Rel_weight))) #Dot product to find the temperature prediction at the current timestep
            #Temp = (TempsInWindow.dot(np.transpose(Rel_weight))) #Dot product to find the temperature prediction at the current timestep
            #Temp = (np.transpose(TempsInWindow).dot(Rel_weight))

            T_MA = np.append(T_MA,Temp) #append Temp to the ongoing T_MA dataframe (contains values for the whole year)
            #T_MA[cntr1] = Temp

        cntr1 += 1 


    return(T_MA)

For reference this is the Python function I was using in my grad studies, although most of the model development and validation was done in Matlab.

[cwhanse]

@kanderso-nrel I didn't realize that rolling.apply() was so slow (compared to built-in methods). I'll take a second look at it's use in #1074.

If performance is a roadblock, for irregular data on short intervals, resampling seems appropriate.

[kanderso-nrel]

Thanks @mjprilliman, it'll be interesting to see how your implementation compares. It's good to have a reference implementation too. I don't think this function is tied to a deadline so of course whenever you can find the time for it is fine.

I didn't realize that rolling.apply() was so slow (compared to built-in methods)

@cwhanse I forgot to mention this earlier but you can specify raw=True in apply to get an order of magnitude or two speedup. The difference is whether the applied function is passed Series objects (raw=False) or numpy arrays (raw=True). In this application I needed to keep it as a Series for the timedelta calculations so I used the slow route. Maybe the fast path could work in #1074, not sure.

[mikofski]

My natural starting point was to use pandas to get the offset-based windows and apply a custom averaging function. The problem is that tseries.rolling(...).apply(...) is way too slow

• It maybe b/c you are using apply. Have you tried: tseries.rolling('1200s').mean()? it may be faster

• A few times I have also found Pandas operations to be convenient but slow compared to NumPy, but AFAIK NumPy doesn't have rolling(), but I think I found this online somewhere:

def moving_average(x, window=5):
    """
    Moving average of numpy array

    Parameters
    ----------
    x : numeric
        a numpy array to average
    window : int
        the window over which to average
    Returns
    -------
    an array of the same size with index at beginning of window

    If ``window <= 1`` then do nothing and return ``x``.
    """
    if window <= 1: return x
    m = window - 1
    y = np.pad(x, (0, m), 'edge')
    z = 0
    for n in range(window):
        z += np.roll(y, -n)
    return z[:-m] / window

It's right bounded, so not sure how to backfill or center. I using numpy.roll and numpy.pad

[kanderso-nrel]

It maybe b/c you are using apply. Have you tried: tseries.rolling('1200s').mean()? it may be faster

The built-in functions are definitely way faster. The problem is that I didn't see a way to implement the exponential weighting using only the built-in functions, so I went with apply and a custom function. It's tricky not only because it's a weighted average of all points in the last 1200s but also because the weights vary for each window based on wind speed (and also the index if the series is irregularly-sampled).

I think I found this online somewhere

The issue there, assuming we want to accommodate irregular sampling, is that the window size can vary based on how many timestamps fall into the last 1200s -- it's an offset-based window, not a count-based window. I couldn't see an elegant/efficient way to recreate pandas's offset-based window with numpy, but would love to see one!

[jforbess]

I recently became aware of np.nditer, which sped up my rolling window apply greatly. I think you should be able to use it for this exponential weighted function, basically applied as a for loop with an efficiently stored/accessed array. I was doing a simple calc on 20 years of 525600 long series, and it was quite snappy, relative to the rolling window apply which was a long enough process that the full 20 year calc was at least tens of minutes if not an hour. Sorry I don't have more quantified values.