diff options
Diffstat (limited to 'Project/PowerQuality.cpp')
| -rw-r--r-- | Project/PowerQuality.cpp | 266 |
1 files changed, 158 insertions, 108 deletions
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<Bus*>(load->GetParentList()[0])->GetElectricalData().number; - LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); - wxMessageBox(wxString::Format("%f %f", data.activePower, data.reactivePower)); - std::complex<double> yLoad = - std::complex<double>(data.activePower, -data.reactivePower / harmYBus.order); - std::complex<double> v = static_cast<Bus*>(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<std::vector<std::complex<double> > >& 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<Bus*>(load->GetParentList()[0])->GetElectricalData().number; + LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); + std::complex<double> yLoad = std::complex<double>(data.activePower, -data.reactivePower / order); + std::complex<double> v = static_cast<Bus*>(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<Bus*>(capacitor->GetParentList()[0])->GetElectricalData().number; - CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += std::complex<double>(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<Bus*>(capacitor->GetParentList()[0])->GetElectricalData().number; + CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase); + yBus[n][n] += std::complex<double>(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<Bus*>(inductor->GetParentList()[0])->GetElectricalData().number; - InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += std::complex<double>(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<Bus*>(inductor->GetParentList()[0])->GetElectricalData().number; + InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase); + yBus[n][n] += std::complex<double>(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<Bus*>(line->GetParentList()[0])->GetElectricalData().number; - int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number; + int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number; - harmYBus.yBus[n1][n2] -= - 1.0 / std::complex<double>(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n1] -= - 1.0 / std::complex<double>(data.resistance, data.indReactance * harmYBus.order); + yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance * order); + yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance * order); - harmYBus.yBus[n1][n1] += - 1.0 / std::complex<double>(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n2] += - 1.0 / std::complex<double>(data.resistance, data.indReactance * harmYBus.order); + yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance * order); + yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance * order); - harmYBus.yBus[n1][n1] += std::complex<double>(0.0, (data.capSusceptance * harmYBus.order) / 2.0); - harmYBus.yBus[n2][n2] += std::complex<double>(0.0, (data.capSusceptance * harmYBus.order) / 2.0); - } + yBus[n1][n1] += std::complex<double>(0.0, (data.capSusceptance * order) / 2.0); + yBus[n2][n2] += std::complex<double>(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<Bus*>(transformer->GetParentList()[0])->GetElectricalData().number; - int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData().number; + int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast<Bus*>(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<double>(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n1] += - -1.0 / std::complex<double>(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<double>(data.resistance, data.indReactance * order); + yBus[n2][n1] += -1.0 / std::complex<double>(data.resistance, data.indReactance * order); - harmYBus.yBus[n1][n1] += - 1.0 / std::complex<double>(data.resistance, data.indReactance * harmYBus.order); - harmYBus.yBus[n2][n2] += - 1.0 / std::complex<double>(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<double> a = std::complex<double>(data.turnsRatio * std::cos(radPhaseShift), - -data.turnsRatio * std::sin(radPhaseShift)); - - // Transformer admitance - std::complex<double> y = - 1.0 / std::complex<double>(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<double>(data.resistance, data.indReactance * order); + yBus[n2][n2] += 1.0 / std::complex<double>(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<double> a = std::complex<double>(data.turnsRatio * std::cos(radPhaseShift), + -data.turnsRatio * std::sin(radPhaseShift)); + + // Transformer admitance + std::complex<double> y = 1.0 / std::complex<double>(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<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number; - SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += - 1.0 / std::complex<double>(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<Bus*>(syncMotor->GetParentList()[0])->GetElectricalData().number; - SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); - harmYBus.yBus[n][n] += - 1.0 / std::complex<double>(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<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number; + SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); + yBus[n][n] += 1.0 / std::complex<double>(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<Bus*>(syncMotor->GetParentList()[0])->GetElectricalData().number; + SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); + yBus[n][n] += 1.0 / std::complex<double>(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<int>(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<double> 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<std::vector<std::complex<double> > > yBus; + for(unsigned int i = 0; i < m_busList.size(); i++) { + std::vector<std::complex<double> > line; + for(unsigned int j = 0; j < m_busList.size(); j++) { line.push_back(std::complex<double>(0.0, 0.0)); } + yBus.push_back(line); + } + // Create and fill with zoros the injected current vector + std::vector<std::complex<double> > iInj; + for(unsigned int i = 0; i < m_busList.size(); i++) { iInj.push_back(std::complex<double>(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<double>(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<double>(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<double>(0.0, 0.0); + } + + currentFreq += stepFreq; + } + return false; +} |