From addef9a1e5bdc3f7951c35be41437e1ae13f1715 Mon Sep 17 00:00:00 2001 From: Manuel Pitz Date: Wed, 1 Jun 2022 18:15:29 +0200 Subject: [PATCH] Hook pmu dft classes --- common | 2 +- etc/examples/hooks/ip-dft-pmu.conf | 27 +++ etc/examples/hooks/pmu.conf | 26 +++ include/villas/hooks/pmu.hpp | 103 +++++++++++ lib/hooks/CMakeLists.txt | 2 + lib/hooks/dp.cpp | 2 +- lib/hooks/pmu.cpp | 270 +++++++++++++++++++++++++++++ lib/hooks/pmu_ipdft.cpp | 173 ++++++++++++++++++ 8 files changed, 603 insertions(+), 2 deletions(-) create mode 100644 etc/examples/hooks/ip-dft-pmu.conf create mode 100644 etc/examples/hooks/pmu.conf create mode 100644 include/villas/hooks/pmu.hpp create mode 100644 lib/hooks/pmu.cpp create mode 100644 lib/hooks/pmu_ipdft.cpp diff --git a/common b/common index 0368040dd..25cd53ee6 160000 --- a/common +++ b/common @@ -1 +1 @@ -Subproject commit 0368040ddaef8cfeeb64f47fb16b6b1333b18a52 +Subproject commit 25cd53ee6882c3f66746d6d8c27790ef22d18322 diff --git a/etc/examples/hooks/ip-dft-pmu.conf b/etc/examples/hooks/ip-dft-pmu.conf new file mode 100644 index 000000000..e7c4bb0d5 --- /dev/null +++ b/etc/examples/hooks/ip-dft-pmu.conf @@ -0,0 +1,27 @@ +@include "hook-nodes.conf" + +paths = ( + { + in = "signal_node" + out = "file_node" + + hooks = ( + { + type = "pmu", + + signals = ( + "sine" + ) + + sample_rate = 1000, # sample rate of the input signal + dft_rate = 10, # number of phasors calculated per second + + estimation_range = 10, # the range around the nominal_freq in with the estimation is done + nominal_freq = 50, # the nominal grid frequnecy + number_plc = 10., # the number of power line cylces stored in the buffer + + angle_unit = "rad" # one of: rad, degree + } + ) + } +) diff --git a/etc/examples/hooks/pmu.conf b/etc/examples/hooks/pmu.conf new file mode 100644 index 000000000..ab36e9d2e --- /dev/null +++ b/etc/examples/hooks/pmu.conf @@ -0,0 +1,26 @@ +@include "hook-nodes.conf" + +paths = ( + { + in = "signal_node" + out = "file_node" + + hooks = ( + { + type = "pmu", + + signals = ( + "sine" + ) + + sample_rate = 1000, # sample rate of the input signal + dft_rate = 10, # number of phasors calculated per second + + nominal_freq = 50, # the nominal grid frequnecy + number_plc = 10., # the number of power line cylces stored in the buffer + + angle_unit = "rad" # one of: rad, degree + } + ) + } +) diff --git a/include/villas/hooks/pmu.hpp b/include/villas/hooks/pmu.hpp new file mode 100644 index 000000000..a537b0ca8 --- /dev/null +++ b/include/villas/hooks/pmu.hpp @@ -0,0 +1,103 @@ +/** PMU hook. + * + * @author Manuel Pitz + * @copyright 2014-2022, Institute for Automation of Complex Power Systems, EONERC + * @license GNU General Public License (version 3) + * + * VILLASnode + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + *********************************************************************************/ + +#include +#include +#include + +namespace villas { +namespace node { + +class PmuHook : public MultiSignalHook { + + +protected: + enum class Status { + INVALID, + VALID + }; + + struct Phasor { + double frequency; + double amplitude; + double phase; + double rocof; /* Rate of change of frequency. */ + Status valid; + }; + + enum class WindowType { + NONE, + FLATTOP, + HANN, + HAMMING, + NUTTAL, + BLACKMAN + }; + + enum class TimeAlign { + LEFT, + CENTER, + RIGHT, + }; + + std::vector*> windows; + dsp::Window *windowsTs; + std::vector lastPhasors; + + enum TimeAlign timeAlignType; + enum WindowType windowType; + + unsigned sampleRate; + double phasorRate; + double nominalFreq; + double numberPlc; + unsigned windowSize; + bool channelNameEnable; + double angleUnitFactor; + uint64_t lastSequence; + timespec nextRun; + bool init; + unsigned initSampleCount; + /* Correction factors. */ + double phaseOffset; + double amplitudeOffset; + double frequencyOffset; + double rocofOffset; + +public: + PmuHook(Path *p, Node *n, int fl, int prio, bool en = true); + + virtual + void prepare(); + + virtual + void parse(json_t *json); + + virtual + Hook::Reason process(struct Sample *smp); + + virtual + Phasor estimatePhasor(dsp::CosineWindow *window, Phasor lastPhasor); +}; + +} /* namespace node */ +} /* namespace villas */ diff --git a/lib/hooks/CMakeLists.txt b/lib/hooks/CMakeLists.txt index fefb2ac78..e0e634833 100644 --- a/lib/hooks/CMakeLists.txt +++ b/lib/hooks/CMakeLists.txt @@ -35,7 +35,9 @@ set(HOOK_SRC limit_value.cpp ma.cpp power.cpp + pmu.cpp pmu_dft.cpp + pmu_ipdft.cpp pps_ts.cpp print.cpp restart.cpp diff --git a/lib/hooks/dp.cpp b/lib/hooks/dp.cpp index 125a561c5..919be7cf0 100644 --- a/lib/hooks/dp.cpp +++ b/lib/hooks/dp.cpp @@ -58,7 +58,7 @@ protected: void step(double *in, std::complex *out) { - int N = window.getLength(); + int N = window.size(); __attribute__((unused)) std::complex om_k, corr; double newest = *in; __attribute__((unused)) double oldest = window.update(newest); diff --git a/lib/hooks/pmu.cpp b/lib/hooks/pmu.cpp new file mode 100644 index 000000000..629f3d975 --- /dev/null +++ b/lib/hooks/pmu.cpp @@ -0,0 +1,270 @@ +/** PMU hook. + * + * @author Manuel Pitz + * @copyright 2014-2022, Institute for Automation of Complex Power Systems, EONERC + * @license GNU General Public License (version 3) + * + * VILLASnode + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + *********************************************************************************/ + +#include +#include + +namespace villas { +namespace node { + +PmuHook::PmuHook(Path *p, Node *n, int fl, int prio, bool en): + MultiSignalHook(p, n, fl, prio, en), + windows(), + windowsTs(), + timeAlignType(TimeAlign::CENTER), + windowType(WindowType::NONE), + sampleRate(1), + phasorRate(1.0), + nominalFreq(1.0), + numberPlc(1.), + windowSize(1), + channelNameEnable(true), + angleUnitFactor(1.0), + lastSequence(0), + nextRun({0}), + init(false), + initSampleCount(0), + phaseOffset(0.0), + amplitudeOffset(0.0), + frequencyOffset(0.0), + rocofOffset(0.0) +{ } + +void PmuHook::prepare() +{ + MultiSignalHook::prepare(); + + signals->clear(); + for (unsigned i = 0; i < signalIndices.size(); i++) { + /* Add signals */ + auto freqSig = std::make_shared("frequency", "Hz", SignalType::FLOAT); + auto amplSig = std::make_shared("amplitude", "V", SignalType::FLOAT); + auto phaseSig = std::make_shared("phase", (angleUnitFactor)?"rad":"deg", SignalType::FLOAT);//angleUnitFactor==1 means rad + auto rocofSig = std::make_shared("rocof", "Hz/s", SignalType::FLOAT); + + if (!freqSig || !amplSig || !phaseSig || !rocofSig) + throw RuntimeError("Failed to create new signals"); + + if (channelNameEnable) { + auto suffix = fmt::format("_{}", signalNames[i]); + + freqSig->name += suffix; + amplSig->name += suffix; + phaseSig->name += suffix; + rocofSig->name += suffix; + } + + signals->push_back(freqSig); + signals->push_back(amplSig); + signals->push_back(phaseSig); + signals->push_back(rocofSig); + + lastPhasors.push_back({0.,0.,0.,0.}); + } + + windowSize = ceil(sampleRate * numberPlc / nominalFreq); + + for (unsigned i = 0; i < signalIndices.size(); i++) { + if (windowType == WindowType::NONE) + windows.push_back(new dsp::RectangularWindow(windowSize, 0.0)); + else if (windowType == WindowType::FLATTOP) + windows.push_back(new dsp::FlattopWindow(windowSize, 0.0)); + else if (windowType == WindowType::HAMMING) + windows.push_back(new dsp::HammingWindow(windowSize, 0.0)); + else if (windowType == WindowType::HANN) + windows.push_back(new dsp::HannWindow(windowSize, 0.0)); + else if (windowType == WindowType::NUTTAL) + windows.push_back(new dsp::NuttallWindow(windowSize, 0.0)); + else if (windowType == WindowType::BLACKMAN) + windows.push_back(new dsp::BlackmanWindow(windowSize, 0.0)); + } + + windowsTs = new dsp::Window(windowSize, timespec{0}); +} + +void PmuHook::parse(json_t *json) +{ + MultiSignalHook::parse(json); + int ret; + + const char *windowTypeC = nullptr; + const char *angleUnitC = nullptr; + const char *timeAlignC = nullptr; + + json_error_t err; + + assert(state != State::STARTED); + + Hook::parse(json); + + ret = json_unpack_ex(json, &err, 0, "{ s?: i, s?: F, s?: F, s?: F, s?: s, s?: s, s?: b, s?: s, s?: F, s?: F, s?: F, s?: F}", + "sample_rate", &sampleRate, + "dft_rate", &phasorRate, + "nominal_freq", &nominalFreq, + "number_plc", &numberPlc, + "window_type", &windowTypeC, + "angle_unit", &angleUnitC, + "add_channel_name", &channelNameEnable, + "timestamp_align", &timeAlignC, + "phase_offset", &phaseOffset, + "amplitude_offset", &litudeOffset, + "frequency_offset", &frequencyOffset, + "rocof_offset", &rocofOffset + ); + + if (ret) + throw ConfigError(json, err, "node-config-hook-pmu"); + + if (sampleRate <= 0) + throw ConfigError(json, "node-config-hook-pmu-sample_rate", "Sample rate cannot be less than 0 tried to set {}", sampleRate); + + if (phasorRate <= 0) + throw ConfigError(json, "node-config-hook-pmu-phasor_rate", "Phasor rate cannot be less than 0 tried to set {}", phasorRate); + + if (nominalFreq <= 0) + throw ConfigError(json, "node-config-hook-pmu-nominal_freq", "Nominal frequency cannot be less than 0 tried to set {}", nominalFreq); + + if (numberPlc <= 0) + throw ConfigError(json, "node-config-hook-pmu-number_plc", "Number of power line cycles cannot be less than 0 tried to set {}", numberPlc); + + if (!windowTypeC) + logger->info("No Window type given, assume no windowing"); + else if (strcmp(windowTypeC, "flattop") == 0) + windowType = WindowType::FLATTOP; + else if (strcmp(windowTypeC, "hamming") == 0) + windowType = WindowType::HAMMING; + else if (strcmp(windowTypeC, "hann") == 0) + windowType = WindowType::HANN; + else if (strcmp(windowTypeC, "nuttal") == 0) + windowType = WindowType::NUTTAL; + else if (strcmp(windowTypeC, "blackman") == 0) + windowType = WindowType::BLACKMAN; + else + throw ConfigError(json, "node-config-hook-pmu-window-type", "Invalid window type: {}", windowTypeC); + + if (!angleUnitC) + logger->info("No angle type given, assume rad"); + else if (strcmp(angleUnitC, "rad") == 0) + angleUnitFactor = 1; + else if (strcmp(angleUnitC, "degree") == 0) + angleUnitFactor = 180 / M_PI; + else + throw ConfigError(json, "node-config-hook-pmu-angle-unit", "Angle unit {} not recognized", angleUnitC); + + if (!timeAlignC) + logger->info("No timestamp alignment defined. Assume alignment center"); + else if (strcmp(timeAlignC, "left") == 0) + timeAlignType = TimeAlign::LEFT; + else if (strcmp(timeAlignC, "center") == 0) + timeAlignType = TimeAlign::CENTER; + else if (strcmp(timeAlignC, "right") == 0) + timeAlignType = TimeAlign::RIGHT; + else + throw ConfigError(json, "node-config-hook-dft-timestamp-alignment", "Timestamp alignment {} not recognized", timeAlignC); + + +} + +Hook::Reason PmuHook::process(struct Sample *smp) +{ + assert(state == State::STARTED); + + initSampleCount++; + + if (smp->sequence - lastSequence > 1) + logger->warn("Calculation is not Realtime. {} sampled missed", smp->sequence - lastSequence); + lastSequence = smp->sequence; + + if (!init && initSampleCount > windowSize) + init = true; + + timespec timeDiff = time_diff(&nextRun, &smp->ts.origin); + double tmpTimeDiff = time_to_double(&timeDiff); + bool run = false; + if (tmpTimeDiff > 0. && init) + run = true; + + Status phasorStatus = Status::VALID; + timespec phasorTimestamp = {0}; + if (run) { + for (unsigned i = 0; i < signalIndices.size(); i++) { + lastPhasors[i] = estimatePhasor(windows[i], lastPhasors[i]); + if (lastPhasors[i].valid != Status::VALID) + phasorStatus = Status::INVALID; + } + + double m = pow(10, floor(log10(phasorRate) + 1)); + if (phasorRate <= 1) // For less then 1 phasor per second + m = pow(10, floor(log10(phasorRate))); + double nextRunDouble = (floor(time_to_double(&smp->ts.origin)*m)/m) + 1.0 / phasorRate; + double r = fmod(nextRunDouble, 1 / phasorRate); + if( (r > 1 / (4 * phasorRate)) && (r < 3 / (4 * phasorRate))) + nextRunDouble += (1.0 / phasorRate) - fmod(nextRunDouble, 1 / phasorRate); + nextRun = time_from_double(nextRunDouble); + + unsigned tsPos = 0; + if (timeAlignType == TimeAlign::RIGHT) + tsPos = windowSize; + else if (timeAlignType == TimeAlign::LEFT) + tsPos = 0; + else if (timeAlignType == TimeAlign::CENTER) + tsPos = windowSize / 2; + phasorTimestamp = (*windowsTs)[tsPos]; + } + + /* Update sample memory */ + unsigned i = 0; + for (auto index : signalIndices) + windows[i++]->update(smp->data[index].f); + windowsTs->update(smp->ts.origin); + + /* Make sure to update phasors after window update but estimate them before */ + if(run) { + for (unsigned i = 0; i < signalIndices.size(); i++) { + smp->data[i * 4 + 0].f = lastPhasors[i].frequency + frequencyOffset; /* Frequency */ + smp->data[i * 4 + 1].f = (lastPhasors[i].amplitude / pow(2, 0.5)) + amplitudeOffset; /* Amplitude */ + smp->data[i * 4 + 2].f = (lastPhasors[i].phase * 180 / M_PI) + phaseOffset; /* Phase */ + smp->data[i * 4 + 3].f = lastPhasors[i].rocof + rocofOffset; /* ROCOF */; + } + smp->ts.origin = phasorTimestamp; + + smp->length = signalIndices.size() * 4; + } + + if (!run || phasorStatus != Status::VALID) + return Reason::SKIP_SAMPLE; + + return Reason::OK; +} + +PmuHook::Phasor PmuHook::estimatePhasor(dsp::CosineWindow *window, Phasor lastPhasor) +{ + return {0., 0., 0., 0., Status::INVALID}; +} + +/* Register hook */ +static char n[] = "pmu"; +static char d[] = "This hook estimates a phsor"; +static HookPlugin p; + +} /* namespace node */ +} /* namespace villas */ diff --git a/lib/hooks/pmu_ipdft.cpp b/lib/hooks/pmu_ipdft.cpp new file mode 100644 index 000000000..2c1e22b26 --- /dev/null +++ b/lib/hooks/pmu_ipdft.cpp @@ -0,0 +1,173 @@ +/** ipDFT PMU hook. + * + * @author Manuel Pitz + * @copyright 2014-2022, Institute for Automation of Complex Power Systems, EONERC + * @license GNU General Public License (version 3) + * + * VILLASnode + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + *********************************************************************************/ + +#include + +namespace villas { +namespace node { + +class IpDftPmuHook : public PmuHook { + +protected: + std::complex omega; + std::vector>> dftMatrix; + std::vector> dftResult; + + unsigned frequencyCount; /* Number of requency bins that are calculated */ + double estimationRange; /* The range around nominalFreq used for estimation */ + +public: + IpDftPmuHook(Path *p, Node *n, int fl, int prio, bool en = true): + PmuHook(p, n, fl, prio, en), + frequencyCount(0), + estimationRange(0) + + { } + + void prepare() + { + PmuHook::prepare(); + + const double startFrequency = nominalFreq - estimationRange; + const double endFrequency = nominalFreq + estimationRange; + const double frequencyResolution = (double)sampleRate / windowSize; + + frequencyCount = ceil((endFrequency - startFrequency) / frequencyResolution); + + /* Initialize matrix of dft coeffients */ + dftMatrix.clear(); + for (unsigned i = 0; i < frequencyCount; i++) + dftMatrix.emplace_back(windowSize, 0.0); + + using namespace std::complex_literals; + + omega = exp((-2i * M_PI) / (double)(windowSize)); + unsigned startBin = floor(startFrequency / frequencyResolution); + + for (unsigned i = 0; i < frequencyCount ; i++) { + for (unsigned j = 0 ; j < windowSize; j++) + dftMatrix[i][j] = pow(omega, (i + startBin) * j); + dftResult.push_back(0.0); + } + } + + void parse(json_t *json) + { + PmuHook::parse(json); + int ret; + + json_error_t err; + + assert(state != State::STARTED); + + Hook::parse(json); + + ret = json_unpack_ex(json, &err, 0, "{ s?: F}", + "estimation_range", &estimationRange + ); + + if (ret) + throw ConfigError(json, err, "node-config-hook-ip-dft-pmu"); + + if (estimationRange <= 0) + throw ConfigError(json, "node-config-hook-ip-dft-pmu-estimation_range", "Estimation range cannot be less or equal than 0 tried to set {}", estimationRange); + + } + + PmuHook::Phasor estimatePhasor(dsp::CosineWindow *window, PmuHook::Phasor lastPhasor) + { + PmuHook::Phasor phasor = {0}; + + /* Calculate DFT */ + for (unsigned i = 0; i < frequencyCount; i++) { + dftResult[i] = 0; + + const unsigned size = (*window).size(); + for (unsigned j = 0; j < size; j++) + dftResult[i] += (*window)[j] * dftMatrix[i][j]; + } + /* End calculate DFT */ + + /* Find max bin */ + unsigned maxBin = 0; + double absAmplitude = 0; + + for(unsigned j = 0; j < frequencyCount; j++) { + if (absAmplitude < abs(dftResult[j])) { + absAmplitude = abs(dftResult[j]); + maxBin = j; + } + } + /* End find max bin */ + + if (maxBin == 0 || maxBin == (frequencyCount - 1)) { + logger->warn("Maximum frequency bin lies on window boundary. Using non-estimated results!"); + //@todo add handling to not forward this phasor!! + } else { + const double startFrequency = nominalFreq - estimationRange; + const double frequencyResolution = (double)sampleRate / windowSize; + + double a = abs(dftResult[ maxBin - 1 ]); + double b = abs(dftResult[ maxBin + 0 ]); + double c = abs(dftResult[ maxBin + 1 ]); + double bPhase = atan2(dftResult[maxBin].imag(), dftResult[maxBin].real()); + + /* Estimate phasor */ + /* Based on https://ieeexplore.ieee.org/document/7980868 */ + double delta = 0; + + /* Paper eq 8 */ + if (c > a) { + delta = 1. * (2. * c - b) / (b + c); + } else { + delta = -1. * (2. * a - b) / (b + a); + } + + /* Frequency estimation (eq 4) */ + phasor.frequency = startFrequency + ( (double) maxBin + delta) * frequencyResolution; + + /* Amplitude estimation (eq 9) */ + phasor.amplitude = b * abs( (M_PI * delta) / sin( M_PI * delta) ) * abs( pow(delta, 2) - 1); + phasor.amplitude *= 2 / (windowSize * window->getCorrectionFactor()); + + /* Phase estimation (eq 10) */ + phasor.phase = bPhase - M_PI * delta; + + /* ROCOF estimation */ + phasor.rocof = ((phasor.frequency - lastPhasor.frequency) * (double)phasorRate); + /* End estimate phasor */ + } + + if (lastPhasor.frequency != 0) /* Check if we already calculated a phasor before */ + phasor.valid = Status::VALID; + + return phasor; + } +}; + +/* Register hook */ +static char n[] = "ip-dft-pmu"; +static char d[] = "This hook calculates a phasor based on ipDFT"; +static HookPlugin p; + +} /* namespace node */ +} /* namespace villas */