From 2771fff79ac9c3c09b70f4668e7142b2e944d1f2 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 25 Apr 2019 01:25:41 -0300 Subject: Matpower Importer and power quality calculation Power quality in implementation --- Project/PowerQuality.cpp | 250 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 250 insertions(+) create mode 100644 Project/PowerQuality.cpp (limited to 'Project/PowerQuality.cpp') diff --git a/Project/PowerQuality.cpp b/Project/PowerQuality.cpp new file mode 100644 index 0000000..33acbe6 --- /dev/null +++ b/Project/PowerQuality.cpp @@ -0,0 +1,250 @@ +#include "PowerQuality.h" + +PowerQuality::PowerQuality() {} + +PowerQuality::~PowerQuality() {} + +PowerQuality::PowerQuality(std::vector elementList) { GetElementsFromList(elementList); } + +void PowerQuality::CalculateHarmonicYbusList(double systemPowerBase) +{ + // Clear and fill with zeros all the harmonic Ybuses + for(auto it = m_harmYbusList.begin(), itEnd = m_harmYbusList.end(); it != itEnd; ++it) { + HarmonicYbus harmYBus = *it; + harmYBus.yBus.clear(); + for(unsigned int i = 0; i < m_busList.size(); i++) { + std::vector > line; + for(unsigned int j = 0; j < m_busList.size(); j++) { line.push_back(std::complex(0.0, 0.0)); } + harmYBus.yBus.push_back(line); + } + *it = harmYBus; + } + + // Fill all Ybuses + for(auto itYbus = m_harmYbusList.begin(), itYbusEnd = m_harmYbusList.end(); itYbus != itYbusEnd; ++itYbus) { + HarmonicYbus harmYBus = *itYbus; + // Load + for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) { + Load* load = *it; + if(load->IsOnline()) { + int n = static_cast(load->GetParentList()[0])->GetElectricalData().number; + LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); + if(data.loadType == CONST_IMPEDANCE) { + std::complex yLoad = + std::complex(data.activePower, -data.reactivePower / harmYBus.order); + harmYBus.yBus[n][n] += yLoad; + } + } + } + + // Capacitor + for(auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) { + Capacitor* capacitor = *it; + if(capacitor->IsOnline()) { + int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; + CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase); + harmYBus.yBus[n][n] += std::complex(0.0, data.reactivePower) * harmYBus.order; + } + } + + // Inductor + for(auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) { + Inductor* inductor = *it; + if(inductor->IsOnline()) { + int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; + InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase); + harmYBus.yBus[n][n] += std::complex(0.0, -data.reactivePower) / harmYBus.order; + } + } + + // Power line + for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) { + Line* line = *it; + if(line->IsOnline()) { + LineElectricalData data = line->GetPUElectricalData(systemPowerBase); + + int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(line->GetParentList()[1])->GetElectricalData().number; + + harmYBus.yBus[n1][n2] -= + 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + harmYBus.yBus[n2][n1] -= + 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + + harmYBus.yBus[n1][n1] += + 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + harmYBus.yBus[n2][n2] += + 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + + harmYBus.yBus[n1][n1] += std::complex(0.0, (data.capSusceptance * harmYBus.order) / 2.0); + harmYBus.yBus[n2][n2] += std::complex(0.0, (data.capSusceptance * harmYBus.order) / 2.0); + } + } + + // Transformer + for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) { + Transformer* transformer = *it; + + if(transformer->IsOnline()) { + TransformerElectricalData data = transformer->GetPUElectricalData(systemPowerBase); + + int n1 = static_cast(transformer->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(transformer->GetParentList()[1])->GetElectricalData().number; + + // If the transformer have nominal turns ratio (1.0) and no phase shifting, it will be modelled as + // series impedance. + if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) { + harmYBus.yBus[n1][n2] += + -1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + harmYBus.yBus[n2][n1] += + -1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + + harmYBus.yBus[n1][n1] += + 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + harmYBus.yBus[n2][n2] += + 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + } + // If the transformer have no-nominal turn ratio and/or phase shifting, it will be modelled in a + // different way (see references). + //[Ref. 1: Elementos de analise de sistemas de potencia - Stevenson - pg. 232] + //[Ref. 2: + // http://www.ee.washington.edu/research/real/Library/Reports/Tap_Adjustments_in_AC_Load_Flows.pdf] + // [Ref. 3: http://www.columbia.edu/~dano/courses/power/notes/power/andersson1.pdf] + else { + // Complex turns ratio + double radPhaseShift = wxDegToRad(data.phaseShift); + std::complex a = std::complex(data.turnsRatio * std::cos(radPhaseShift), + -data.turnsRatio * std::sin(radPhaseShift)); + + // Transformer admitance + std::complex y = 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + + harmYBus.yBus[n1][n1] += y / (std::pow(std::abs(a), 2.0)); + harmYBus.yBus[n1][n2] += -(y / std::conj(a)); + harmYBus.yBus[n2][n1] += -(y / a); + harmYBus.yBus[n2][n2] += y; + } + } + } + + // Synchronous generator + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + if(syncGenerator->IsOnline()) { + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); + harmYBus.yBus[n][n] += + 1.0 / std::complex(data.positiveResistance, data.positiveReactance * harmYBus.order); + } + } + // Synchronous motor + for(auto it = m_syncMotorList.begin(), itEnd = m_syncMotorList.end(); it != itEnd; ++it) { + SyncMotor* syncMotor = *it; + if(syncMotor->IsOnline()) { + int n = static_cast(syncMotor->GetParentList()[0])->GetElectricalData().number; + SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); + harmYBus.yBus[n][n] += + 1.0 / std::complex(data.positiveResistance, data.positiveReactance * harmYBus.order); + } + } + *itYbus = harmYBus; + } +} + +bool PowerQuality::CalculateDistortions(double systemPowerBase) +{ + // Get harmonic orders + m_harmYbusList.clear(); + std::vector harmOrders = GetHarmonicOrdersList(); + + // Fill the Ybuses + for(unsigned int i = 0; i < harmOrders.size(); ++i) { + HarmonicYbus newYbus; + newYbus.order = harmOrders[i]; + m_harmYbusList.push_back(newYbus); + } + CalculateHarmonicYbusList(systemPowerBase); + + // Initialize current arrays with zeros + std::vector > > iHarmInjList; + for(unsigned int i = 0; i < harmOrders.size(); i++) { + std::vector > line; + for(unsigned int j = 0; j < m_busList.size(); j++) { line.push_back(std::complex(0.0, 0.0)); } + iHarmInjList.push_back(line); + } + + // Fill the current array + for(unsigned int i = 0; i < harmOrders.size(); ++i) { + for(auto it = m_harmCurrentList.begin(), itEnd = m_harmCurrentList.end(); it != itEnd; ++it) { + HarmCurrent* harmCurrent = *it; + if(harmCurrent->IsOnline()) { + // Get only the current order in analysis + for(unsigned int k = 0; k < harmCurrent->GetElectricalData().harmonicOrder.size(); ++k) { + if(harmCurrent->GetElectricalData().harmonicOrder[k] == static_cast(harmOrders[i])) { + Bus* parentBus = static_cast(harmCurrent->GetParentList()[0]); + auto busData = parentBus->GetElectricalData(); + int j = busData.number; + + // Bus voltage + double voltage = busData.nominalVoltage; + if(busData.nominalVoltageUnit == UNIT_kV) voltage *= 1e3; + + auto puData = harmCurrent->GetPUElectricalData(systemPowerBase, voltage); + + iHarmInjList[i][j] += std::complex( + puData.injHarmCurrent[k] * std::cos(wxDegToRad(puData.injHarmAngle[k])), + puData.injHarmCurrent[k] * std::sin(wxDegToRad(puData.injHarmAngle[k]))); + } + } + } + } + } + + // Calculate harmonic voltages + std::vector > > vHarmList; + for(unsigned int i = 0; i < m_harmYbusList.size(); ++i) { + vHarmList.push_back(GaussianElimination(m_harmYbusList[i].yBus, iHarmInjList[i])); + } + + wxString volt = ""; + for(unsigned int i = 0; i < vHarmList.size(); ++i) { + wxString orderStr = ""; + for(unsigned int j = 0; j < vHarmList[i].size(); ++j) { + orderStr += wxString::Format("%f; ", std::abs(vHarmList[i][j])); + } + volt += orderStr + "\n"; + } + wxMessageBox(volt); + + return true; +} + +std::vector PowerQuality::GetHarmonicOrdersList() +{ + std::vector harmOrders; + auto harmCurrentList = GetHarmCurrentList(); + // Check all harmonic sources and get all harmonic orders in the system + for(auto it = harmCurrentList.begin(), itEnd = harmCurrentList.end(); it != itEnd; ++it) { + HarmCurrent* harmCurrent = *it; + if(harmCurrent->IsOnline()) { + auto data = harmCurrent->GetElectricalData(); + for(unsigned int i = 0; i < data.harmonicOrder.size(); ++i) { + int order = data.harmonicOrder[i]; + // Check if this harmonic order have been added already + bool newOrder = true; + for(unsigned int j = 0; j < harmOrders.size(); ++j) { + if(order == harmOrders[j]) { + newOrder = false; + break; + } + } + if(newOrder) harmOrders.push_back(order); + } + } + } + std::vector doubleHarmOrder; + for(unsigned int i = 0; i < harmOrders.size(); ++i) { + doubleHarmOrder.push_back(static_cast(harmOrders[i])); + } + return doubleHarmOrder; +} -- cgit From 4dabf27f998db83e20bc0eca7e18672777f0bf5b Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 25 Apr 2019 19:26:15 -0300 Subject: THD implemented --- Project/PowerQuality.cpp | 36 ++++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) (limited to 'Project/PowerQuality.cpp') diff --git a/Project/PowerQuality.cpp b/Project/PowerQuality.cpp index 33acbe6..4efdc49 100644 --- a/Project/PowerQuality.cpp +++ b/Project/PowerQuality.cpp @@ -29,11 +29,12 @@ void PowerQuality::CalculateHarmonicYbusList(double systemPowerBase) if(load->IsOnline()) { int n = static_cast(load->GetParentList()[0])->GetElectricalData().number; LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); - if(data.loadType == CONST_IMPEDANCE) { - std::complex yLoad = - std::complex(data.activePower, -data.reactivePower / harmYBus.order); - harmYBus.yBus[n][n] += yLoad; - } + wxMessageBox(wxString::Format("%f %f", data.activePower, data.reactivePower)); + std::complex yLoad = + std::complex(data.activePower, -data.reactivePower / harmYBus.order); + std::complex v = static_cast(load->GetParentList()[0])->GetElectricalData().voltage; + yLoad /= (std::abs(v) * std::abs(v)); + harmYBus.yBus[n][n] += yLoad; } } @@ -117,7 +118,8 @@ void PowerQuality::CalculateHarmonicYbusList(double systemPowerBase) -data.turnsRatio * std::sin(radPhaseShift)); // Transformer admitance - std::complex y = 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + std::complex y = + 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); harmYBus.yBus[n1][n1] += y / (std::pow(std::abs(a), 2.0)); harmYBus.yBus[n1][n2] += -(y / std::conj(a)); @@ -205,16 +207,22 @@ bool PowerQuality::CalculateDistortions(double systemPowerBase) for(unsigned int i = 0; i < m_harmYbusList.size(); ++i) { vHarmList.push_back(GaussianElimination(m_harmYbusList[i].yBus, iHarmInjList[i])); } - - wxString volt = ""; - for(unsigned int i = 0; i < vHarmList.size(); ++i) { - wxString orderStr = ""; - for(unsigned int j = 0; j < vHarmList[i].size(); ++j) { - orderStr += wxString::Format("%f; ", std::abs(vHarmList[i][j])); + + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + auto data = bus->GetElectricalData(); + data.harmonicOrder.clear(); + data.harmonicVoltage.clear(); + double thd = 0.0; + for(unsigned int i = 0; i < vHarmList.size(); ++i) { + thd += std::abs(vHarmList[i][data.number]) * std::abs(vHarmList[i][data.number]); + data.harmonicVoltage.push_back(vHarmList[i][data.number]); + data.harmonicOrder.push_back(static_cast(harmOrders[i])); } - volt += orderStr + "\n"; + //distortion = std::sqrt(distortion) / std::abs(data.voltage); + thd = std::sqrt(thd) * 100.0; + data.thd = thd; } - wxMessageBox(volt); return true; } -- cgit From a40d5a405d60b4e429f6f578dcfe3c33ab5ad81a Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Fri, 26 Apr 2019 02:29:47 -0300 Subject: Frequency response implemented Need a form --- Project/PowerQuality.cpp | 266 ++++++++++++++++++++++++++++------------------- 1 file changed, 158 insertions(+), 108 deletions(-) (limited to 'Project/PowerQuality.cpp') diff --git a/Project/PowerQuality.cpp b/Project/PowerQuality.cpp index 4efdc49..30969d9 100644 --- a/Project/PowerQuality.cpp +++ b/Project/PowerQuality.cpp @@ -23,133 +23,127 @@ void PowerQuality::CalculateHarmonicYbusList(double systemPowerBase) // Fill all Ybuses for(auto itYbus = m_harmYbusList.begin(), itYbusEnd = m_harmYbusList.end(); itYbus != itYbusEnd; ++itYbus) { HarmonicYbus harmYBus = *itYbus; - // Load - for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) { - Load* load = *it; - if(load->IsOnline()) { - int n = static_cast(load->GetParentList()[0])->GetElectricalData().number; - LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); - wxMessageBox(wxString::Format("%f %f", data.activePower, data.reactivePower)); - std::complex yLoad = - std::complex(data.activePower, -data.reactivePower / harmYBus.order); - std::complex v = static_cast(load->GetParentList()[0])->GetElectricalData().voltage; - yLoad /= (std::abs(v) * std::abs(v)); - harmYBus.yBus[n][n] += yLoad; - } + CalculateHarmonicYbus(harmYBus.yBus, systemPowerBase, harmYBus.order); + *itYbus = harmYBus; + } +} + +void PowerQuality::CalculateHarmonicYbus(std::vector > >& yBus, + double systemPowerBase, + double order) +{ + // Load + for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) { + Load* load = *it; + if(load->IsOnline()) { + int n = static_cast(load->GetParentList()[0])->GetElectricalData().number; + LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); + std::complex yLoad = std::complex(data.activePower, -data.reactivePower / order); + std::complex v = static_cast(load->GetParentList()[0])->GetElectricalData().voltage; + yLoad /= (std::abs(v) * std::abs(v)); + yBus[n][n] += yLoad; } + } - // Capacitor - for(auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) { - Capacitor* capacitor = *it; - if(capacitor->IsOnline()) { - int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; - CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += std::complex(0.0, data.reactivePower) * harmYBus.order; - } + // Capacitor + for(auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) { + Capacitor* capacitor = *it; + if(capacitor->IsOnline()) { + int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; + CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase); + yBus[n][n] += std::complex(0.0, data.reactivePower) * order; } + } - // Inductor - for(auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) { - Inductor* inductor = *it; - if(inductor->IsOnline()) { - int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; - InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += std::complex(0.0, -data.reactivePower) / harmYBus.order; - } + // Inductor + for(auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) { + Inductor* inductor = *it; + if(inductor->IsOnline()) { + int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; + InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase); + yBus[n][n] += std::complex(0.0, -data.reactivePower) / order; } + } - // Power line - for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) { - Line* line = *it; - if(line->IsOnline()) { - LineElectricalData data = line->GetPUElectricalData(systemPowerBase); + // Power line + for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) { + Line* line = *it; + if(line->IsOnline()) { + LineElectricalData data = line->GetPUElectricalData(systemPowerBase); - int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; - int n2 = static_cast(line->GetParentList()[1])->GetElectricalData().number; + int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(line->GetParentList()[1])->GetElectricalData().number; - harmYBus.yBus[n1][n2] -= - 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n1] -= - 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + yBus[n1][n2] -= 1.0 / std::complex(data.resistance, data.indReactance * order); + yBus[n2][n1] -= 1.0 / std::complex(data.resistance, data.indReactance * order); - harmYBus.yBus[n1][n1] += - 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n2] += - 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + yBus[n1][n1] += 1.0 / std::complex(data.resistance, data.indReactance * order); + yBus[n2][n2] += 1.0 / std::complex(data.resistance, data.indReactance * order); - harmYBus.yBus[n1][n1] += std::complex(0.0, (data.capSusceptance * harmYBus.order) / 2.0); - harmYBus.yBus[n2][n2] += std::complex(0.0, (data.capSusceptance * harmYBus.order) / 2.0); - } + yBus[n1][n1] += std::complex(0.0, (data.capSusceptance * order) / 2.0); + yBus[n2][n2] += std::complex(0.0, (data.capSusceptance * order) / 2.0); } + } - // Transformer - for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) { - Transformer* transformer = *it; + // Transformer + for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) { + Transformer* transformer = *it; - if(transformer->IsOnline()) { - TransformerElectricalData data = transformer->GetPUElectricalData(systemPowerBase); + if(transformer->IsOnline()) { + TransformerElectricalData data = transformer->GetPUElectricalData(systemPowerBase); - int n1 = static_cast(transformer->GetParentList()[0])->GetElectricalData().number; - int n2 = static_cast(transformer->GetParentList()[1])->GetElectricalData().number; + int n1 = static_cast(transformer->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(transformer->GetParentList()[1])->GetElectricalData().number; - // If the transformer have nominal turns ratio (1.0) and no phase shifting, it will be modelled as - // series impedance. - if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) { - harmYBus.yBus[n1][n2] += - -1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n1] += - -1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); + // If the transformer have nominal turns ratio (1.0) and no phase shifting, it will be modelled as + // series impedance. + if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) { + yBus[n1][n2] += -1.0 / std::complex(data.resistance, data.indReactance * order); + yBus[n2][n1] += -1.0 / std::complex(data.resistance, data.indReactance * order); - harmYBus.yBus[n1][n1] += - 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n2] += - 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); - } - // If the transformer have no-nominal turn ratio and/or phase shifting, it will be modelled in a - // different way (see references). - //[Ref. 1: Elementos de analise de sistemas de potencia - Stevenson - pg. 232] - //[Ref. 2: - // http://www.ee.washington.edu/research/real/Library/Reports/Tap_Adjustments_in_AC_Load_Flows.pdf] - // [Ref. 3: http://www.columbia.edu/~dano/courses/power/notes/power/andersson1.pdf] - else { - // Complex turns ratio - double radPhaseShift = wxDegToRad(data.phaseShift); - std::complex a = std::complex(data.turnsRatio * std::cos(radPhaseShift), - -data.turnsRatio * std::sin(radPhaseShift)); - - // Transformer admitance - std::complex y = - 1.0 / std::complex(data.resistance, data.indReactance * harmYBus.order); - - harmYBus.yBus[n1][n1] += y / (std::pow(std::abs(a), 2.0)); - harmYBus.yBus[n1][n2] += -(y / std::conj(a)); - harmYBus.yBus[n2][n1] += -(y / a); - harmYBus.yBus[n2][n2] += y; - } + yBus[n1][n1] += 1.0 / std::complex(data.resistance, data.indReactance * order); + yBus[n2][n2] += 1.0 / std::complex(data.resistance, data.indReactance * order); } - } + // If the transformer have no-nominal turn ratio and/or phase shifting, it will be modelled in a + // different way (see references). + //[Ref. 1: Elementos de analise de sistemas de potencia - Stevenson - pg. 232] + //[Ref. 2: + // http://www.ee.washington.edu/research/real/Library/Reports/Tap_Adjustments_in_AC_Load_Flows.pdf] + // [Ref. 3: http://www.columbia.edu/~dano/courses/power/notes/power/andersson1.pdf] + else { + // Complex turns ratio + double radPhaseShift = wxDegToRad(data.phaseShift); + std::complex a = std::complex(data.turnsRatio * std::cos(radPhaseShift), + -data.turnsRatio * std::sin(radPhaseShift)); + + // Transformer admitance + std::complex y = 1.0 / std::complex(data.resistance, data.indReactance * order); - // Synchronous generator - for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { - SyncGenerator* syncGenerator = *it; - if(syncGenerator->IsOnline()) { - int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; - SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += - 1.0 / std::complex(data.positiveResistance, data.positiveReactance * harmYBus.order); + yBus[n1][n1] += y / (std::pow(std::abs(a), 2.0)); + yBus[n1][n2] += -(y / std::conj(a)); + yBus[n2][n1] += -(y / a); + yBus[n2][n2] += y; } } - // Synchronous motor - for(auto it = m_syncMotorList.begin(), itEnd = m_syncMotorList.end(); it != itEnd; ++it) { - SyncMotor* syncMotor = *it; - if(syncMotor->IsOnline()) { - int n = static_cast(syncMotor->GetParentList()[0])->GetElectricalData().number; - SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += - 1.0 / std::complex(data.positiveResistance, data.positiveReactance * harmYBus.order); - } + } + + // Synchronous generator + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + if(syncGenerator->IsOnline()) { + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); + yBus[n][n] += 1.0 / std::complex(data.positiveResistance, data.positiveReactance * order); + } + } + // Synchronous motor + for(auto it = m_syncMotorList.begin(), itEnd = m_syncMotorList.end(); it != itEnd; ++it) { + SyncMotor* syncMotor = *it; + if(syncMotor->IsOnline()) { + int n = static_cast(syncMotor->GetParentList()[0])->GetElectricalData().number; + SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); + yBus[n][n] += 1.0 / std::complex(data.positiveResistance, data.positiveReactance * order); } - *itYbus = harmYBus; } } @@ -207,7 +201,7 @@ bool PowerQuality::CalculateDistortions(double systemPowerBase) for(unsigned int i = 0; i < m_harmYbusList.size(); ++i) { vHarmList.push_back(GaussianElimination(m_harmYbusList[i].yBus, iHarmInjList[i])); } - + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { Bus* bus = *it; auto data = bus->GetElectricalData(); @@ -219,9 +213,10 @@ bool PowerQuality::CalculateDistortions(double systemPowerBase) data.harmonicVoltage.push_back(vHarmList[i][data.number]); data.harmonicOrder.push_back(static_cast(harmOrders[i])); } - //distortion = std::sqrt(distortion) / std::abs(data.voltage); + // distortion = std::sqrt(distortion) / std::abs(data.voltage); thd = std::sqrt(thd) * 100.0; data.thd = thd; + bus->SetElectricalData(data); } return true; @@ -256,3 +251,58 @@ std::vector PowerQuality::GetHarmonicOrdersList() } return doubleHarmOrder; } + +bool PowerQuality::CalculateFrequencyResponse(double systemFreq, + double initFreq, + double endFreq, + double stepFreq, + double systemPowerBase) +{ + // Clear all previous data + for(unsigned int i = 0; i < m_busList.size(); i++) { + auto data = m_busList[i]->GetElectricalData(); + data.absImpedanceVector.clear(); + data.absImpedanceVector.shrink_to_fit(); + m_busList[i]->SetElectricalData(data); + } + + // Create and fill with zeros the YBus + std::vector > > yBus; + for(unsigned int i = 0; i < m_busList.size(); i++) { + std::vector > line; + for(unsigned int j = 0; j < m_busList.size(); j++) { line.push_back(std::complex(0.0, 0.0)); } + yBus.push_back(line); + } + // Create and fill with zoros the injected current vector + std::vector > iInj; + for(unsigned int i = 0; i < m_busList.size(); i++) { iInj.push_back(std::complex(10.0, 0.0)); } + + if(initFreq < 1e-6) initFreq = stepFreq; + double currentFreq = initFreq; + while(currentFreq <= endFreq) { + m_frequencyList.push_back(currentFreq); + + double order = currentFreq / systemFreq; + + // Fill YBus with zeros + for(unsigned int i = 0; i < m_busList.size(); i++) { + for(unsigned int j = 0; j < m_busList.size(); j++) { yBus[i][j] = std::complex(0.0, 0.0); } + } + + CalculateHarmonicYbus(yBus, systemPowerBase, order); + + for(unsigned int i = 0; i < m_busList.size(); i++) { + auto data = m_busList[i]->GetElectricalData(); + iInj[data.number] = std::complex(1.0, 0.0); + + auto zh = GaussianElimination(yBus, iInj); + + data.absImpedanceVector.push_back(std::abs(zh[data.number])); + m_busList[i]->SetElectricalData(data); + iInj[data.number] = std::complex(0.0, 0.0); + } + + currentFreq += stepFreq; + } + return false; +} -- cgit From 2b02ef22cc5f2025b09b700f1cb6e1cec94d80f6 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Fri, 26 Apr 2019 17:52:32 -0300 Subject: Power quality fully implemented A formal filter element must be implemented in future --- Project/PowerQuality.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) (limited to 'Project/PowerQuality.cpp') diff --git a/Project/PowerQuality.cpp b/Project/PowerQuality.cpp index 30969d9..c8fcd91 100644 --- a/Project/PowerQuality.cpp +++ b/Project/PowerQuality.cpp @@ -256,6 +256,7 @@ bool PowerQuality::CalculateFrequencyResponse(double systemFreq, double initFreq, double endFreq, double stepFreq, + int injBusNumber, double systemPowerBase) { // Clear all previous data @@ -273,9 +274,10 @@ bool PowerQuality::CalculateFrequencyResponse(double systemFreq, for(unsigned int j = 0; j < m_busList.size(); j++) { line.push_back(std::complex(0.0, 0.0)); } yBus.push_back(line); } - // Create and fill with zoros the injected current vector + // Create and fill with zeros the injected current vector std::vector > iInj; - for(unsigned int i = 0; i < m_busList.size(); i++) { iInj.push_back(std::complex(10.0, 0.0)); } + for(unsigned int i = 0; i < m_busList.size(); i++) { iInj.push_back(std::complex(0.0, 0.0)); } + iInj[injBusNumber] = std::complex(1.0, 0.0); if(initFreq < 1e-6) initFreq = stepFreq; double currentFreq = initFreq; @@ -293,13 +295,12 @@ bool PowerQuality::CalculateFrequencyResponse(double systemFreq, for(unsigned int i = 0; i < m_busList.size(); i++) { auto data = m_busList[i]->GetElectricalData(); - iInj[data.number] = std::complex(1.0, 0.0); + if(data.plotPQData) { + auto zh = GaussianElimination(yBus, iInj); - auto zh = GaussianElimination(yBus, iInj); - - data.absImpedanceVector.push_back(std::abs(zh[data.number])); - m_busList[i]->SetElectricalData(data); - iInj[data.number] = std::complex(0.0, 0.0); + data.absImpedanceVector.push_back(std::abs(zh[data.number])); + m_busList[i]->SetElectricalData(data); + } } currentFreq += stepFreq; -- cgit