Medición del factor de potencia en sinusoides

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

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

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

• Multiplique el voltaje instantáneo y la serie temporal actual para obtener potencia instantánea, p_inst .
• Integre p_inst (utilizando una integración trapezoidal) y divida por la longitud de tiempo del segmento de datos para obtener la potencia activa, p_act .
• Escuadrar cada muestra en p_inst , tomar el promedio de p_inst**2 en la serie temporal y luego tomar la raíz cuadrada de ese valor para obtener la potencia RMS, p_rms .
• Calcule 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 (de nuevo una aproximación trapezoidal) cada una de las series de tiempo resultantes, dando como resultado v_cos y v_sin . El ángulo de fase de la señal de voltaje, v_phi , relativo a alguna referencia arbitraria, es entonces atan2(v_sin, v_cos) .
• Repita el paso anterior para el actual. El retraso de fase de la corriente en relación con el voltaje es entonces φ = 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 cos(φ) correcto 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 relación de potencias produce una estimación muy errónea: 0.77542956418409648. Esto equivale a un error de casi diez grados en el ángulo de fase.

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

¿En qué me he equivocado aquí?

Código Python completo 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)))