Electronic – Measuring power factor on sinusoids

I'm not sure if this is the best forum for this or whether stackoverflow mightn't be better, but I'll ask here and see what people think.

I have a system measuring 50Hz voltage and currents at 1kHz. I'm using quadrature demodulation to measure the phase shift and calculating cos(φ) from that, and also calculating the active and RMS power and inferring cos(φ) from them. These calculations are done over one-second slices of the data.

I'll attach example Python code below but basically the method looks like this:

• Multiply the instantaneous voltage and current timeseries together to get instantaneous power, p_inst.
• Integrate p_inst (using a trapezoidal integration) and divide by the time length of the data slice to get the active power, p_act.
• Square each sample in p_inst, take the average of p_inst**2 across the timeseries and then take the square root of that value to get the RMS power, p_rms.
• Calculate the power factor, pf = p_act / p_rms.
• Using a zero-crossing-detection method, estimate the frequency of the voltage signal.
• Generate quadrature reference signals at that frequency (ie cos(2πft), sin(2πft)).
• Multiply the voltage by each of the reference signals and integrate (again a trapezoidal approximation) each of the resulting timeseries, resulting in v_cos and v_sin. The phase angle of the voltage signal, v_phi, relative to some arbitrary reference, is then atan2(v_sin, v_cos).
• Repeat the previous step for the current. The phase lag of the current relative to the voltage is then φ = v_phi - i_phi.
• My second estimate of the power factor is then cos(φ).

For a sinusoidal signal with 30 degrees phase lag, the correct cos(φ) is approximately 0.86602540378443871. The quadrature demodulation method produces a very good approximation, 0.86636025346085943. But the ratio-of-powers method produces a very wrong estimate – 0.77542956418409648. This is equivalent to nearly ten degrees error in the phase angle.

I assumed at first that I'd got the quadrature demodulation wrong (being the more complex calculation) but it produces the right answer. I then assumed that the signal was badly non-sinusoidal and that this would explain the difference – but the code below does the same calculation on ideal sinusoids.

What have I got wrong here?

Full Python code demonstrating the problem:

import numpy as np, pandas as pd
from matplotlib.pyplot import *

t = np.arange(0, 600, 0.001)
t_rad = 2 * np.pi * 50 * t
phi = 30 * np.pi / 180

v = 230 * np.sin(t_rad)
i = 200 * np.sin(t_rad + phi)

p_inst = v * i

data = pd.DataFrame(np.array([t, v, i, p_inst]).transpose(), columns=['t', 'v', 'i', 'p_inst'], index=t)

def gen_act(x):
    return np.trapz(x.p_inst, x=x.index) / (x.index[-1] - x.index[0])

def gen_rms(x):
    return np.sqrt(np.mean(x.p_inst**2))

second_bins = data.groupby(lambda x: int(x))
p_act = second_bins.apply(gen_act)
p_rms = second_bins.apply(gen_rms)

def estimate_frequency(t, v):
    t = t.values
    v = v.values
    try:
        zero_crossings = np.where(np.diff(np.sign(v)) > 0.5)[0]
        diffs = np.diff(t[zero_crossings])
        accum_intervals = np.copy(diffs)
        accum_interval = 0
        ii = 0
        for ii in range(len(accum_intervals)):
            accum_interval += diffs[ii]
            accum_intervals[ii] = accum_interval
            if accum_interval >= 0.01:
                accum_interval = 0

        zero_crossings = zero_crossings[np.hstack([np.where(accum_intervals > 0.01)[0], -1])]
        ends = zero_crossings[[0, -1]]
        frequency = (len(zero_crossings) - 1) / (t[ends[1]] - t[ends[0]])
        return frequency
    except:
        return float('NaN')

def generate_reference_signals(t, f):
    freq = np.mean(f)
    if np.isnan(freq):
        return None
    x = t * 2 * np.pi * freq
    return np.sin(x), np.cos(x)

def estimate_phi(t, v, i, refs):
    def angle_from_refs(t, x, refs):
        s_i = np.trapz(x * refs[0], t)
        c_i = np.trapz(x * refs[1], t)
        return np.arctan2(c_i, s_i)

    angle = angle_from_refs(t, i, refs) - angle_from_refs(t, v, refs)
    while angle > np.pi:
        angle -= 2*np.pi
    return angle

def gen_phi(x):
    f = estimate_frequency(x.t, x.v)
    refs = generate_reference_signals(x.t, f)
    return estimate_phi(x.t, x.v, x.i, refs)

phi = second_bins.apply(gen_phi)
cos_phi = np.cos(phi)
pf = p_act / p_rms

print('Power factor estimated from quadrature demodulation: {}'.format(np.mean(cos_phi)))
print('Power factor estimated from measured power: {}'.format(np.mean(pf)))
print('True power factor: {}'.format(np.cos(30 * np.pi / 180)))