Сообщество Engee
Engee Community
/
Примеры и проекты
/
QR разложение при оценке первой гармоники тока КЗ

QR разложение при оценке первой гармоники тока КЗ

Автор
avatar-andreyilinskiyandreyilinskiy
Notebook

QR разложение при оценке первой гармоники тока КЗ

Одна из основных задач релейной защиты - производить классификацию наблюдаемых в энергосистеме (ЭС) режимов на аварийные и не аварийные. В простейшем случае эта задача решается по величине протекающего тока через защищаемый элемент сети. Под величиной тока можно понимать действующее значение основной (первой) гармоники тока (хотя сам сигнал часто более сложный и может содержать как кратные высшие гармоники, так и затухающую экспоненциальную составляющую).

BASELINE

Базовым подходом для выделения первой гармонической составляющей является использование формул дискретного преобразования Фурье. Используется вычисление лишь одного спектрального отсчета. ДПФ применяется к данным, длиной один период на каждом такте расчета.
Предположим, для простоты, что входной сигнал оцифрован с частотой 400 Гц (8 точек на период 20 мс). Тогда входной сигнал для каждого момента времени является вектор-столбцом с восьмью элементами.

Также введем вектор-столбцы кратных гармоник, то есть, ортогональный базис, по которому будет раскладываться вектор входного сигнала:

Столбец неизвестных составляющих гармоник (коэффициенты при базисных векторах):

Матрицу , столбцами которой будут базисные вектора:

In [ ]: 
using Printf
In [ ]: 
# Сформируем основную матрицу системы H
n = 8
dpsi = 2pi/n
psi = [0: dpsi: (n-1)*dpsi;]
H = zeros(n, n)
for k in 1:n-1
    if isodd(k)
        H[:, k] = cos.((div(k, 2)+1)*psi)
    else
        H[:, k] = sin.(div(k, 2)*psi)
    end
end
H[:, n] = ones(n, 1);

Я, к сожалению, быстро не нашел функцию для отображения матриц без значений типа 6.12323e-17, поэтому пришлось написать собственную функцию.

In [ ]: 
function showm(matrix)
    n, m = size(matrix)
    for i in 1:n
        for j in 1:m
            @printf("%9.2f  ", matrix[i, j])
        end
        print("\n")
    end
    print("\n")
    return nothing
end;

Тогда задача нахождения составляющих гармоник в сигнале сводится к СЛАУ вида:

Принципиальным моментом является ортогональность введенного базиса:

где - диагональная матрица.

In [ ]: 
# убедимся в диагональности матрицы
showm(H'*H)

     4.00       0.00      -0.00       0.00       0.00       0.00       0.00      -0.00  
     0.00       4.00       0.00      -0.00       0.00       0.00       0.00      -0.00  
    -0.00       0.00       4.00       0.00      -0.00       0.00       0.00      -0.00  
     0.00      -0.00       0.00       4.00      -0.00      -0.00       0.00       0.00  
     0.00       0.00      -0.00      -0.00       4.00       0.00       0.00      -0.00  
     0.00       0.00       0.00      -0.00       0.00       4.00       0.00      -0.00  
     0.00       0.00       0.00       0.00       0.00       0.00       8.00       0.00  
    -0.00      -0.00      -0.00       0.00      -0.00      -0.00       0.00       8.00

Это позволяет находить проекции входного вектора независимо друг от друга, и формула для вычисления одного отсчета ДПФ становится просто переписанной формулой нахождения проекции вектора на вектор из курса линейной алгебры.

Например, для :

Затруднения

При КЗ в сигналах от энергосистемы появляются свободные составляющие. В простейшем случае это одна экспоненциально затухающая составляющая. В этом случае базис константа + кратные гармоники уже подходит не так хорошо.

Зная, что в сигнале впредь есть экспонента, можно попробовать представить вектор сигнала как

где, например,

Это сработает, но при добавлении в эту линейную регрессию еще членов (вторую, третью гармонику) придется пересчитывать и те коэффициенты, которые были рассчитаны ранее. Потеряна возможность находить ответы независимо друг от друга. Это связано с тем, что базисные вектора из-за замены на более в совокупности не ортогональны.

In [ ]: 
# Создадим матрицу He, где заменим последний столбец
q = 0.8 # при частоте дискретизации 400 Гц данное значение соответствует экспоненте с постоянной времени 11 мс
expv = [q^k for k in 0:n-1]
He = hcat(H[:, 1:n-1], expv)

# Убедимся что матрица He'*He не диагональна (но определенная структура всё же есть)
showm(He'*He)

     4.00       0.00      -0.00       0.00       0.00       0.00       0.00       0.71  
     0.00       4.00       0.00      -0.00       0.00       0.00       0.00       0.93  
    -0.00       0.00       4.00       0.00      -0.00       0.00       0.00       0.51  
     0.00      -0.00       0.00       4.00      -0.00      -0.00       0.00       0.41  
     0.00       0.00      -0.00      -0.00       4.00       0.00       0.00       0.47  
     0.00       0.00       0.00      -0.00       0.00       4.00       0.00       0.17  
     0.00       0.00       0.00       0.00       0.00       0.00       8.00       0.46  
     0.71       0.93       0.51       0.41       0.47       0.17       0.46       2.70

Решение

В статье [1] предложено перейти к взвешенному МНК и решать задачу вида

где - матрица весов.
Согласно [1] она вводится формально без физического смысла для целей приведения к диаогональному виду (это приводит к независимости уравнений СЛАУ и позволяет находить каждое неизвестное из одного уравнения, а не из полного решения СЛАУ). Это делается за счет того что строки матрицы равны столбцам матрицы смещенным на определенную константу.

Обоснование

Дадим обоснование вышеописанному способу через QR разложение. Рассмотрим матрицу . Её столбцы ортогональны между собой, за исключением последнего, который экспонента. Давайте ортогонализируем весь базис. Применим к последнему столбцу алгоритм Грамма-Шмидта.

In [ ]: 
v = expv
for k in 1:n-1    
    v -= (v'*H[:, k]) / (H[:, k]'*H[:, k]) * H[:, k]
end
display(v)

8-element Vector{Float64}:
 0.5201424000000001
 0.5201424000000002
 0.5201423999999999
 0.5201424000000001
 0.5201424000000001
 0.5201424000000001
 0.5201424000000003
 0.5201424000000001

Результат, по сути, вектор (пропорционален ему с коэффициентом 0,520142). Можно было предугадать этот результат, так как это последний вектор ортогональнгого базиса которого нам не хватало. То есть мы представили матрицу как , где - верхне-треугольная матрица. Мы произвели QR разложение с оговоркой что столбцы матрицы ортогональны, но не нормированны.

где - коэффициенты разложения вектора по ортогональному базису кратных гармоник.
Тогда
\begin{aligned}
\mathbf{He}\cdot \mathbf{x}&=\mathbf{y} \
\mathbf{HR}\cdot \mathbf{x}&=\mathbf{y} \
\mathbf{H}^T\mathbf{HR}\cdot \mathbf{x}&=\mathbf{H}^T \mathbf{y}\
\mathbf{DR}\cdot \mathbf{x}&=\mathbf{H}^T \mathbf{y}\
\mathbf{R}\cdot \mathbf{x}&=\mathbf{D}{-1}\mathbf{H}T \mathbf{y}\
\end{aligned}
Из последнего уравнения можно сразу найти коэффициент при экспоненте и далее уже независимо друг от друга можно находить коэффициенты при любой гармонике. Так как в каждое уравнение входит только две неизвестных (искомый коэффициент при синусе/косинусе и величина экспоненты).
Для коэффициента теперь будем иметь

В скобках как раз получили значения косинуса сдвинутые на постоянную величину. Тем самым пришли к описанному в статье алгоритму.

Частотные характеристики

Сравним частотные характеристки и характеристики подавления вещественных компонент. Используем библиотеку DSP.jl

In [ ]: 
using DSP

c1 = H[:, 1]
s1 = H[:, 2]
one = H[:, n]

fs = 400
filterCos = PolynomialRatio(2c1/n, [1])
filterSin = PolynomialRatio(2s1/n, [1])
c1FreqResp, w = freqresp(filterCos)
s1FreqResp, w = freqresp(filterSin)

filterRe = PolynomialRatio(2/n*(c1 - ones(n)*(expv'*c1) / (expv'*one)), [1])
filterIm = PolynomialRatio(2/n*(s1 - ones(n)*(expv'*s1) / (expv'*one)), [1])
c1cFreqResp, w = freqresp(filterRe)
s1cFreqResp, w = freqresp(filterIm)

freq = collect(w)/pi*fs/2

plot(freq, abs.(c1FreqResp), label="cos", lw=3, lc="purple")
plot!(freq, abs.(s1FreqResp), label="sin", lw=3, lc="purple", ls=:dash)
plot!(freq, abs.(c1cFreqResp), label="cos-const", lw=3, lc="pink")
plot!(freq, abs.(s1cFreqResp), label="sin-const", lw=3, lc="pink", ls=:dash)
title!("АЧХ")
xlabel!("F, Гц")
ylabel!("Коэф. передачи")
Out[0]:

На графиках видно, что все фильтры без изменения передают первую гармонику (50 Гц), давят кратные частоты. Отличие возникает при передаче постоянной составляющей. Это связано с тем, что модифицированные фильтры реализованы в предположении, что в сигнале есть экспонента. И эта модификация сделана как раз ценой отказа от подавления постоянной составляющей.

Характеристики подавления апериодической составляющей

Сравним те же фильтры по подавлению конкретной апериодической составляющей.

In [ ]: 
# Напишем функцию для построения характеристики
function presp(coefFIR, p, fs)
    Hp = []
    for x in p
        Hp = push!(Hp, coefFIR' * exp.([-x/fs*k for k in 0:length(coefFIR)-1]))
    end

    return Hp
end

p = 0: .1: 150
Hpc1 = presp(2c1/n, p, fs);
Hps1 = presp(2s1/n, p, fs);

Hpc1c = presp(2/n*(c1 - ones(n)*(expv'*c1) / (expv'*one)), p, fs);
Hps1c = presp(2/n*(s1 - ones(n)*(expv'*s1) / (expv'*one)), p, fs);

plot(p, abs.(Hpc1), label="cos", lw=3, lc="purple")
plot!(p, abs.(Hps1), label="sin", lw=3, lc="purple", ls=:dash)
plot!(p, abs.(Hpc1c), label="cos-const", lw=3, lc="pink")
plot!(p, abs.(Hps1c), label="sin-const", lw=3, lc="pink", ls=:dash)
title!("Амплитудно \"затухательная\" характеристика")
xlabel!("Коэффициент затухания, 1/с")
ylabel!("Коэф. передачи")
Out[0]:

По графиками видно, что классические фильтры (sin, cos) подавляют экпоненту с нулевым затуханием (то есть константу). Это было видно и из АЧХ. А модифицированные зануляют экспоненту с коэффициентом затухания 89.26 (по модулю), что соответствует введенному вектору и принятой частоте дискретизации.

Иллюстрация работы алгоритма

Тестовый сигнал

Напишем функцию, реализующую вышеописанный алгоритм, и проверим ее работу на тестовом сигнале и на токе, записанном устройством РЗА в реальной ЭС.

In [ ]: 
function calcAmp(signal, samplesInPer, useModAlg)
    N = length(signal)

    dphi = 2pi/samplesInPer
    phi = (0: (samplesInPer-1))*dphi
    c1Base = 2/samplesInPer*cos.(phi)
    s1Base = 2/samplesInPer*sin.(phi)

    mesData = zeros(samplesInPer)
    prevSum = 1
    levelOfConst = 2
    minq = exp(-1/50/samplesInPer/0.005)
    vector = zeros(N)+im*zeros(N)

    for k in 1:N
        mesData = vcat(mesData[2:samplesInPer], current[k]) #скользящее окно по данным длиной 1 период
        
        if useModAlg
            s = sum(mesData) # В этом сигнале останется только экспонента и константа
        else
            s = 0
        end

        if abs(s) > levelOfConst # Запускаем удаление апериодики только если есть постоянная составялющая сигнала достаточной величины
            q = s / prevSum # Оцениваем exp(pT), где Т - период дискретризации, р - коэффициент затухания
            prevSum = s

            if q > 1 # Проверка что экспонента затухающая
                q = 1
            elseif q < minq # Слишком быстро затухающие экспоненты тоже странно для реальных сигналов
                q = minq
            end

            expv = [q^j for j in 0:samplesInPer-1]
            c1Mod = c1Base - ones(samplesInPer) * (c1Base'*expv/sum(expv)) # Формируем сдвинутые коэффициенты фильтров для подавления оцененной выше апериодики
            s1Mod = s1Base - ones(samplesInPer) * (s1Base'*expv/sum(expv))

            vector[k] = mesData'*(c1Mod - im*s1Mod)
        else
            prevSum = 1

            vector[k] = mesData'*(c1Base - im*s1Base)
        end
    end

    return vector
end;
In [ ]: 
#Для начала посмотрим работу на "стерильных" сигналах
n = 48
N = 5*n
dt = 1/50/n

time_ = (0: (N-1))*dt
current = (exp.(-time_/0.015)-cos.(100pi*time_))*10

vectorMod = calcAmp(current, n, true)

# Для сравнения расчитаем действующее значение тока также и базовым алгоритмом
vectorBase = calcAmp(current, n, false)

plot(time_, current, lw = 2, label="Ток")
plot!(time_, abs.(vectorBase), lw=2, label="Действующее значение по классич. алгоритму")
plot!(time_, abs.(vectorMod), lw=2, label="Действующее значение по модифицир. алгоритму")
xlabel!("Время, с")
ylabel!("Ток, А")
title!("Работа алгоритма на чистых сигналах")
Out[0]:

По графикам хорошо видно, что на зеленую кривую наличие затухающей апериодической составляющей в сигнале не влияет. Алгоритм сохранил частотные свойства базисного алгоритма (с заменой подавления константы на экспоненту). При этом количество точек, необходимых для расчета не увеличивается или увеличиается незначительно (для целей оценки коэффициента затухания апериодической состалвяющей)

Реальный сигнал

Я специально выбрал сигнал "покривее" чтобы было видно, что реальные сигналы имееют более сложную структуру. Хотя и чистые, как в книжке, осциллограммы также встречаются.

In [ ]: 
# Реальный ток из осциллограммы аварии
current = [0.32, 0.32, 0.32, 0.32, 0.32, 0.32, 0.32, 0.32, 0.32, 0.9, 1.3, 1.9000000000000001, 2.62, 3.42, 4.26, 5.1000000000000005, 5.9, 6.72, 7.640000000000001, 8.9, 10.200000000000001, 11.44, 12.48, 13.82, 14.8, 15.72, 16.32, 16.52, 
16.34, 15.92, 15.22, 14.120000000000001, 12.64, 10.92, 9.02, 7.38, 5.76, 4.04, 2.38, 0.92, -0.38, -1.58, -2.58, -3.3000000000000003, -2.14, -3.84, -3.96, -3.98, -3.96, -3.88, -3.84, -3.54, -3.2600000000000002, -2.94, -2.58, -2.04, -1.5, -0.86, -0.16, 0.48, 1.26, 2.06, 2.92, 3.68, 4.5200000000000005, 5.28, 5.96, 6.62, 7.36, 8.36, 9.28, 9.92, 10.28, 10.76, 11.26, 11.52, 11.52, 11.36, 10.68, 9.78, 8.540000000000001, 7.28, 6.140000000000001, 5.0600000000000005, 3.88, 2.42, 0.96, -0.4, -1.6400000000000001, -2.64, -3.5, -4.1, -4.44, -4.5600000000000005, -4.8, -4.8, -4.84, -4.8, -4.8, -4.5, -4.18, -3.88, -3.38, -2.9, -2.36, -1.6600000000000001, -0.92, -0.14, 0.68, 1.54, 2.36, 
3.2, 3.96, 4.7, 5.4, 6.0600000000000005, 6.62, 7.22, 7.84, 8.52, 8.9, 9.1, 9.16, 9.28, 9.3, 9.08, 8.66, 7.84, 6.94, 6.08, 5.28, 4.28, 3.18, 1.92, 0.5, -0.76, -1.86, -2.8000000000000003, -3.56, -4.08, -4.42, -4.68, -4.8, -5.0200000000000005, -5.12, -5.1000000000000005, -4.84, -4.76, -4.46, -4.12, -3.6, -3.12, -2.54, -1.8800000000000001, -1.12, -0.36, 0.46, 1.26, 2.1, 2.9, 3.88, 5.32, 5.42, 5.76, 6.34, 6.8, 7.26, 7.62, 8.0, 8.14, 8.1, 7.98, 7.76, 7.640000000000001, 7.2, 6.72, 6.04, 5.3, 4.58, 3.6, 2.5, 1.3800000000000001, 0.08, -1.2, -2.32, -3.1, -3.7, -4.2, -4.62, -4.88, -5.14, -5.38, -5.44, -5.44, -5.24, -5.08, -4.8, -4.36, -3.92, -3.4, -2.82, -2.18, -1.44, -0.68, 0.1, 
0.9, 1.72, 2.48, 3.3000000000000003, 4.0200000000000005, 4.66, 5.26, 5.76, 6.26, 6.66, 6.96, 7.2, 7.36, 7.32, 7.04, 6.68, 6.4, 6.08, 5.68, 5.1000000000000005, 4.48, 3.66, 2.7, 1.6400000000000001, 0.58, -0.54, -1.76, -2.68, -3.48, -4.2, -4.78, -5.12, -5.46, -5.76, -6.0200000000000005, -6.08, -6.08, -5.94, -5.74, -5.44, -4.94, -4.54, -4.0, -3.36, -2.7, -2.04, -1.26, -0.38, 0.38, 1.18, 1.98, 2.66, 3.36, 4.04, 4.58, 5.0600000000000005, 5.54, 5.9, 6.18, 6.4, 6.5, 6.42, 6.2, 5.8, 5.46, 5.12]

N = length(current)

time_ = (0: (N-1))*dt

vectorMod = calcAmp(current, 48, true)

# Для сравнения расчитаем действующее значение тока также и базовым алгоритмом
vectorBase = calcAmp(current, 48, false)

plot(time_, current, lw = 2, label="Ток")
plot!(time_, abs.(vectorBase), lw=2, label="Действующее значение по классич. алгоритму")
plot!(time_, abs.(vectorMod), lw=2, label="Действующее значение по модифицир. алгоритму")
xlabel!("Время, с")
ylabel!("Ток, А")
title!("Работа алгоритма на реальном сигнале")
Out[0]:

Результат не столь впечатляющий, хотя и тут заметно что зеленая линия меньше колеблется относительно некоторого умозрительного среднего.

Заключение

Мы посмотрели как можно объяснить предложение статьи [1] по смещению коэффициентов классического фильтра Фурье на константу для удаления апериодической составляющей. Я потренировался в базовом синтаксисе Julia.

Литература

[1] Eugeniusz Rosołowski, Jan Izykowski, Bogdan Kasztenny. Adaptive measuring algorithm suppressing a decaying DC component for digital protective relays

UPD

Тот же результат в плане математических формул можно получить через понятие взаимного базиса. В нем формула для проекции вектора на вектор работает в своем простейшем виде с заменой координат одного из сомножителей на координаты в взаимном базисе.

Группа видимости
ПубличноПросматривать контент могут все пользователи сообщества

Тип

Интерактивные скрипты

Краткое описание

Дано пояснение смещения коэффициентов фильтров Фурье на константу для удаления апериодической составляющей

Связанные материалы

Пусто

Теги

Категории

Отраслевые примененияОбработка сигналов
МероприятияКонкурс работ Engee 2024

Языки

    Julia

Форматы

    ngscript

Уровень

Пусто

Статус

Опубликовано

Обновлено

Источник

Сообщество