Medición del factor de potencia en sinusoides

No estoy seguro de si este es el mejor foro para esto o si stackoverflow podría no ser mejor, pero preguntaré aquí y veré qué piensa la gente.

Tengo un sistema que mide voltaje de 50 Hz y corrientes a 1 kHz. Estoy usando demodulación en cuadratura para medir el cambio de fase y calcular cos(φ)a partir de eso, y también calcular la potencia activa y RMS e inferir cos(φ)de ellos. Estos cálculos se realizan en porciones de datos de un segundo.

Adjuntaré el código Python de ejemplo a continuación, pero básicamente el método se ve así:

• Multiplique las series temporales de tensión y corriente instantáneas para obtener la potencia instantánea, p_inst.
• Integre p_inst(usando una integración trapezoidal) y divida por la duración del segmento de datos para obtener la potencia activa, p_act.
• Cuadre cada muestra en p_inst, tome el promedio de p_inst**2las series temporales y luego saque la raíz cuadrada de ese valor para obtener la potencia RMS, p_rms.
• Calcular el factor de potencia, pf = p_act / p_rms.
• Usando un método de detección de cruce por cero, estime la frecuencia de la señal de voltaje.
• Genere señales de referencia en cuadratura a esa frecuencia (es decir cos(2πft), , sin(2πft)).
• Multiplique el voltaje por cada una de las señales de referencia e integre (nuevamente una aproximación trapezoidal) cada una de las series de tiempo resultantes, dando como resultado v_cosy v_sin. Entonces , el ángulo de fase de la señal de voltaje, v_phien relación con alguna referencia arbitraria, es atan2(v_sin, v_cos).
• Repita el paso anterior para el actual. Entonces, el desfase de la corriente en relación con el voltaje es φ = v_phi - i_phi.
• Mi segunda estimación del factor de potencia es entonces cos(φ).

Para una señal sinusoidal con un desfase de 30 grados, el valor correcto cos(φ)es aproximadamente 0,86602540378443871. El método de demodulación en cuadratura produce una muy buena aproximación, 0.86636025346085943. Pero el método de la razón de potencias produce una estimación muy incorrecta: 0,77542956418409648. Esto es equivalente a casi diez grados de error en el ángulo de fase.

Al principio supuse que me había equivocado en la demodulación en cuadratura (siendo el cálculo más complejo), pero produce la respuesta correcta. Luego asumí que la señal era muy no sinusoidal y que esto explicaría la diferencia, pero el código a continuación hace el mismo cálculo en sinusoides ideales.

¿Qué tengo mal aquí?

Código completo de Python que demuestra el problema:

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)))