From 30181ca0ae73f5f7f1856ac289db8fcf849c9a84 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Sat, 20 May 2017 17:22:47 -0300 Subject: Electromechanical class and several methods implemented --- Project/ElectricCalculation.cpp | 241 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 218 insertions(+), 23 deletions(-) (limited to 'Project/ElectricCalculation.cpp') diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp index d024364..86653d5 100644 --- a/Project/ElectricCalculation.cpp +++ b/Project/ElectricCalculation.cpp @@ -7,6 +7,7 @@ ElectricCalculation::ElectricCalculation() {} ElectricCalculation::~ElectricCalculation() {} void ElectricCalculation::GetElementsFromList(std::vector elementList) { + m_powerElementList.clear(); m_busList.clear(); m_capacitorList.clear(); m_indMotorList.clear(); @@ -19,6 +20,8 @@ void ElectricCalculation::GetElementsFromList(std::vector elementList) // TODO: Bad design? for(auto it = elementList.begin(); it != elementList.end(); it++) { Element* element = *it; + m_powerElementList.push_back(static_cast(element)); + if(Bus* bus = dynamic_cast(element)) m_busList.push_back(bus); else if(Capacitor* capacitor = dynamic_cast(element)) @@ -41,9 +44,10 @@ void ElectricCalculation::GetElementsFromList(std::vector elementList) } bool ElectricCalculation::GetYBus(std::vector > >& yBus, - double systemPowerBase, - YBusSequence sequence, - bool includeSyncMachines) + double systemPowerBase, + YBusSequence sequence, + bool includeSyncMachines, + bool allLoadsAsImpedances) { if(m_busList.size() == 0) return false; @@ -73,7 +77,7 @@ bool ElectricCalculation::GetYBus(std::vector > if(load->IsOnline()) { int n = static_cast(load->GetParentList()[0])->GetEletricalData().number; LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); - if(data.loadType == CONST_IMPEDANCE) + if(data.loadType == CONST_IMPEDANCE || allLoadsAsImpedances) yBus[n][n] += std::complex(data.activePower, -data.reactivePower); } } @@ -160,8 +164,8 @@ bool ElectricCalculation::GetYBus(std::vector > else if(sequence != ZERO_SEQ) { // Complex turns ratio double radPhaseShift = wxDegToRad(data.phaseShift); - std::complex a = std::complex( - data.turnsRatio * std::cos(radPhaseShift), -data.turnsRatio * std::sin(radPhaseShift)); + 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); @@ -181,10 +185,11 @@ bool ElectricCalculation::GetYBus(std::vector > switch(data.connection) { case GWYE_GWYE: { std::complex y = - 1.0 / std::complex(data.zeroResistance + - 3.0 * (data.primaryGrndResistance + data.secondaryGrndResistance), - data.zeroIndReactance + - 3.0 * (data.primaryGrndReactance + data.secondaryGrndReactance)); + 1.0 / + std::complex( + data.zeroResistance + 3.0 * (data.primaryGrndResistance + data.secondaryGrndResistance), + data.zeroIndReactance + + 3.0 * (data.primaryGrndReactance + data.secondaryGrndReactance)); std::complex a = std::complex(data.turnsRatio, 0.0); yBus[n1][n1] += y / (a * a); @@ -195,14 +200,14 @@ bool ElectricCalculation::GetYBus(std::vector > case DELTA_GWYE: { std::complex y = 1.0 / std::complex(data.zeroResistance + 3.0 * (data.secondaryGrndResistance), - data.zeroIndReactance + 3.0 * (data.secondaryGrndReactance)); + data.zeroIndReactance + 3.0 * (data.secondaryGrndReactance)); yBus[n2][n2] += y; break; } case GWYE_DELTA: { std::complex y = 1.0 / std::complex(data.zeroResistance + 3.0 * (data.primaryGrndResistance), - data.zeroIndReactance + 3.0 * (data.primaryGrndReactance)); + data.zeroIndReactance + 3.0 * (data.primaryGrndReactance)); yBus[n1][n1] += y; break; } @@ -258,10 +263,10 @@ bool ElectricCalculation::GetYBus(std::vector > } void ElectricCalculation::UpdateElementsPowerFlow(std::vector > voltage, - std::vector > power, - std::vector busType, - std::vector reactiveLimit, - double systemPowerBase) + std::vector > power, + std::vector busType, + std::vector reactiveLimit, + double systemPowerBase) { for(int i = 0; i < (int)reactiveLimit.size(); ++i) { if(reactiveLimit[i].maxLimit > -1e-5 && reactiveLimit[i].maxLimit < 1e-5) reactiveLimit[i].maxLimit = 1e-5; @@ -287,9 +292,9 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector v2 = voltage[n2]; data.current[0] = (v1 - v2) / std::complex(data.resistance, data.indReactance) + - v1 * std::complex(0.0, data.capSusceptance / 2.0); + v1 * std::complex(0.0, data.capSusceptance / 2.0); data.current[1] = (v2 - v1) / std::complex(data.resistance, data.indReactance) + - v2 * std::complex(0.0, data.capSusceptance / 2.0); + v2 * std::complex(0.0, data.capSusceptance / 2.0); data.powerFlow[0] = v1 * std::conj(data.current[0]); data.powerFlow[1] = v2 * std::conj(data.current[1]); @@ -310,8 +315,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vectorGetElectricalData(); int n1 = static_cast(transformer->GetParentList()[0])->GetEletricalData().number; int n2 = static_cast(transformer->GetParentList()[1])->GetEletricalData().number; - std::complex v1 = voltage[n1]; // Primary voltage - std::complex v2 = voltage[n2]; // Secondary voltage + std::complex v1 = voltage[n1]; // Primary voltage + std::complex v2 = voltage[n2]; // Secondary voltage // Transformer admitance std::complex y = 1.0 / std::complex(data.resistance, data.indReactance); @@ -491,8 +496,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector > > matrix, - std::vector > >& inverse) + std::vector > >& inverse) { int order = static_cast(matrix.size()); @@ -599,3 +604,193 @@ bool ElectricCalculation::InvertMatrix(std::vector complexValue, double angle, double& dValue, double& qValue) +{ + dValue = -std::real(complexValue) * std::sin(angle) + std::imag(complexValue) * std::cos(angle); + qValue = std::real(complexValue) * std::cos(angle) + std::imag(complexValue) * std::sin(angle); +} + +void ElectricCalculation::DQ0toABC(double dValue, double qValue, double angle, std::complex& complexValue) +{ + double real = qValue * std::cos(angle) - dValue * std::sin(angle); + double imag = qValue * std::sin(angle) + dValue * std::cos(angle); + complexValue = std::complex(real, imag); +} + +std::vector > ElectricCalculation::GaussianElimination( + std::vector > > matrix, + std::vector > array) +{ + //[Ref] http://pt.wikipedia.org/wiki/Elimina%C3%A7%C3%A3o_de_Gauss + + std::vector > solution; + + std::vector > > triangMatrix; + triangMatrix.resize(matrix.size()); + for(unsigned int i = 0; i < matrix.size(); i++) { + triangMatrix[i].resize(matrix.size()); + } + + for(unsigned int i = 0; i < matrix.size(); i++) { + solution.push_back(array[i]); + } + + for(unsigned int i = 0; i < matrix.size(); i++) { + for(unsigned int j = 0; j < matrix.size(); j++) { + triangMatrix[i][j] = matrix[i][j]; + } + } + + for(unsigned int k = 0; k < matrix.size(); k++) { + unsigned int k1 = k + 1; + for(unsigned int i = k; i < matrix.size(); i++) { + if(triangMatrix[i][k] != std::complex(0.0, 0.0)) { + for(unsigned int j = k1; j < matrix.size(); j++) { + triangMatrix[i][j] = triangMatrix[i][j] / triangMatrix[i][k]; + } + solution[i] = solution[i] / triangMatrix[i][k]; + } + } + for(unsigned int i = k1; i < matrix.size(); i++) { + if(triangMatrix[i][k] != std::complex(0.0, 0.0)) { + for(unsigned int j = k1; j < matrix.size(); j++) { + triangMatrix[i][j] -= triangMatrix[k][j]; + } + solution[i] -= solution[k]; + } + } + } + for(unsigned int i = matrix.size() - 2; i >= 0; i--) { + for(unsigned int j = matrix.size() - 1; j >= i + 1; j--) { + solution[i] -= triangMatrix[i][j] * solution[j]; + } + } + + return solution; +} + +SyncMachineModel ElectricCalculation::GetMachineModel(SyncGenerator* generator) +{ + auto data = generator->GetElectricalData(); + if(data.transTd0 != 0.0) { + if(data.transTq0 != 0.0) { + if(data.subTd0 != 0.0) { + if(data.subTq0 != 0.0) + return SM_MODEL_5; + else + return SM_MODEL_4; + } else + return SM_MODEL_3; + } else + return SM_MODEL_2; + } + + return SM_MODEL_1; +} + +std::vector > ElectricCalculation::ComplexMatrixTimesVector( + std::vector > > matrix, + std::vector > vector) +{ + std::vector > solution; + for(unsigned int i = 0; i < matrix.size(); i++) { + solution.push_back(std::complex(0.0, 0.0)); + + for(unsigned int j = 0; j < matrix.size(); j++) { + solution[i] += matrix[i][j] * vector[j]; + } + } + + return solution; +} + +void ElectricCalculation::GetLUDecomposition(std::vector > > matrix, + std::vector > >& matrixL, + std::vector > >& matrixU) +{ + // Doolittle method + // [Ref] http://www3.nd.edu/~zxu2/acms40390F11/Alg-LU-Crout.pdf + // [Ref] http://www.engr.colostate.edu/~thompson/hPage/CourseMat/Tutorials/CompMethods/doolittle.pdf + + int size = static_cast(matrix.size()); // Decomposed matrix size. + + // Set upper and lower matrices sizes. + matrixL.resize(size); + matrixU.resize(size); + for(int i = 0; i < size; i++) { + matrixL[i].resize(size); + matrixU[i].resize(size); + } + + // First row of upper matrix and first column of lower matrix. + for(int i = 0; i < size; i++) { + matrixU[0][i] = matrix[0][i]; + matrixL[i][0] = matrix[i][0] / matrixU[0][0]; + } + + // Upper matrix main diagonal. + for(int i = 1; i < size; i++) { + matrixL[i][i] = std::complex(1.0, 0.0); + } + + for(int i = 1; i < size - 1; i++) { + // Upper matrix main diagonal. + matrixU[i][i] = matrix[i][i]; + for(int k = 0; k < i; k++) { + matrixU[i][i] -= matrixL[i][k] * matrixU[k][i]; + } + + // Others elements of upper matrix + for(int j = i + 1; j < size; j++) { + matrixU[i][j] = matrix[i][j]; + for(int k = 0; k < i; k++) { + matrixU[i][j] -= matrixL[i][k] * matrixU[k][j]; + } + } + + // Lower matrix elements + for(int j = i + 1; j < size; j++) { + matrixL[j][i] = matrix[j][i]; + for(int k = 0; k < i; k++) { + matrixL[j][i] -= matrixL[j][k] * matrixU[k][i]; + } + matrixL[j][i] = matrixL[j][i] / matrixU[i][i]; + } + } + + // Last element of upper matrix. + matrixU[size - 1][size - 1] = matrix[size - 1][size - 1]; + for(int k = 0; k < size - 1; k++) { + matrixU[size - 1][size - 1] -= matrixL[size - 1][k] * matrixU[k][size - 1]; + } +} + +std::vector > ElectricCalculation::LUEvaluate(std::vector > > u, + std::vector > > l, + std::vector > b) +{ + int size = static_cast(b.size()); + std::vector > x; + std::vector > y; + x.resize(size); + y.resize(size); + + // Forward + for(int i = 0; i < size; i++) { + y[i] = b[i]; + for(int j = 0; j < i; j++) { + y[i] -= l[i][j] * y[j]; + } + y[i] /= l[i][i]; + } + // Backward + for(int i = size - 1; i >= 0; i--) { + x[i] = y[i]; + for(int j = i + 1; j < size; j++) { + x[i] -= u[i][j] * x[j]; + } + x[i] /= u[i][i]; + } + return x; +} -- cgit From f995850b38916b38718b84f4b82948479a81855a Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Mon, 22 May 2017 20:17:41 -0300 Subject: Events implemented Effects on adimittance matrix and elements: ->Fault on bus; ->Generation/Load switching; ->Branch switching. --- Project/ElectricCalculation.cpp | 42 +++++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 18 deletions(-) (limited to 'Project/ElectricCalculation.cpp') diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp index 86653d5..847249d 100644 --- a/Project/ElectricCalculation.cpp +++ b/Project/ElectricCalculation.cpp @@ -65,7 +65,7 @@ bool ElectricCalculation::GetYBus(std::vector > int busNumber = 0; for(auto itb = m_busList.begin(); itb != m_busList.end(); itb++) { Bus* bus = *itb; - BusElectricalData data = bus->GetEletricalData(); + BusElectricalData data = bus->GetElectricalData(); data.number = busNumber; bus->SetElectricalData(data); busNumber++; @@ -75,10 +75,16 @@ bool ElectricCalculation::GetYBus(std::vector > 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])->GetEletricalData().number; + int n = static_cast(load->GetParentList()[0])->GetElectricalData().number; LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); - if(data.loadType == CONST_IMPEDANCE || allLoadsAsImpedances) - yBus[n][n] += std::complex(data.activePower, -data.reactivePower); + if(data.loadType == CONST_IMPEDANCE || allLoadsAsImpedances) { + std::complex yLoad = std::complex(data.activePower, -data.reactivePower); + if(allLoadsAsImpedances) { + std::complex v = static_cast(load->GetParentList()[0])->GetElectricalData().voltage; + yLoad /= (v * v); + } + yBus[n][n] += yLoad; + } } } @@ -86,7 +92,7 @@ bool ElectricCalculation::GetYBus(std::vector > 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])->GetEletricalData().number; + int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase); yBus[n][n] += std::complex(0.0, data.reactivePower); } @@ -96,7 +102,7 @@ bool ElectricCalculation::GetYBus(std::vector > 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])->GetEletricalData().number; + int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase); yBus[n][n] += std::complex(0.0, -data.reactivePower); } @@ -108,8 +114,8 @@ bool ElectricCalculation::GetYBus(std::vector > if(line->IsOnline()) { LineElectricalData data = line->GetElectricalData(); - int n1 = static_cast(line->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast(line->GetParentList()[1])->GetEletricalData().number; + int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(line->GetParentList()[1])->GetElectricalData().number; switch(sequence) { case POSITIVE_SEQ: @@ -144,8 +150,8 @@ bool ElectricCalculation::GetYBus(std::vector > if(transformer->IsOnline()) { TransformerElectricalData data = transformer->GetElectricalData(); - int n1 = static_cast(transformer->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast(transformer->GetParentList()[1])->GetEletricalData().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. @@ -223,7 +229,7 @@ bool ElectricCalculation::GetYBus(std::vector > 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])->GetEletricalData().number; + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); switch(sequence) { case POSITIVE_SEQ: { @@ -242,7 +248,7 @@ bool ElectricCalculation::GetYBus(std::vector > 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])->GetEletricalData().number; + int n = static_cast(syncMotor->GetParentList()[0])->GetElectricalData().number; SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); switch(sequence) { case POSITIVE_SEQ: { @@ -275,7 +281,7 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vectorGetEletricalData(); + BusElectricalData data = bus->GetElectricalData(); data.voltage = voltage[i]; bus->SetElectricalData(data); } @@ -284,8 +290,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vectorIsOnline()) { - int n1 = static_cast(line->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast(line->GetParentList()[1])->GetEletricalData().number; + int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(line->GetParentList()[1])->GetElectricalData().number; LineElectricalData data = line->GetElectricalData(); std::complex v1 = voltage[n1]; @@ -313,8 +319,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vectorIsOnline()) { TransformerElectricalData data = transformer->GetElectricalData(); - int n1 = static_cast(transformer->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast(transformer->GetParentList()[1])->GetEletricalData().number; + int n1 = static_cast(transformer->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(transformer->GetParentList()[1])->GetElectricalData().number; std::complex v1 = voltage[n1]; // Primary voltage std::complex v2 = voltage[n2]; // Secondary voltage @@ -348,7 +354,7 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vectorGetEletricalData(); + BusElectricalData data = bus->GetElectricalData(); // Get the synchronous machines connected and calculate the load power on the bus. std::vector syncGeneratorsOnBus; -- cgit From ec288f99d922ad81b560137c1e561ce02139690a Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Fri, 26 May 2017 17:05:54 -0300 Subject: Some bugs fixed --- Project/ElectricCalculation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Project/ElectricCalculation.cpp') diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp index 847249d..a0667d9 100644 --- a/Project/ElectricCalculation.cpp +++ b/Project/ElectricCalculation.cpp @@ -81,7 +81,7 @@ bool ElectricCalculation::GetYBus(std::vector > std::complex yLoad = std::complex(data.activePower, -data.reactivePower); if(allLoadsAsImpedances) { std::complex v = static_cast(load->GetParentList()[0])->GetElectricalData().voltage; - yLoad /= (v * v); + yLoad /= (std::abs(v) * std::abs(v)); } yBus[n][n] += yLoad; } @@ -735,7 +735,7 @@ void ElectricCalculation::GetLUDecomposition(std::vector(1.0, 0.0); } -- cgit From 4ca362b083d7b20adfd80f1ec4b46a52789cdeb7 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Mon, 17 Jul 2017 21:02:07 -0300 Subject: Some optimization, model 1 machine implementation start --- Project/ElectricCalculation.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'Project/ElectricCalculation.cpp') diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp index a0667d9..2fcab57 100644 --- a/Project/ElectricCalculation.cpp +++ b/Project/ElectricCalculation.cpp @@ -676,23 +676,23 @@ std::vector > ElectricCalculation::GaussianElimination( return solution; } -SyncMachineModel ElectricCalculation::GetMachineModel(SyncGenerator* generator) +Machines::SyncMachineModel ElectricCalculation::GetMachineModel(SyncGenerator* generator) { auto data = generator->GetElectricalData(); if(data.transTd0 != 0.0) { if(data.transTq0 != 0.0) { if(data.subTd0 != 0.0) { if(data.subTq0 != 0.0) - return SM_MODEL_5; + return Machines::SM_MODEL_5; else - return SM_MODEL_4; + return Machines::SM_MODEL_4; } else - return SM_MODEL_3; + return Machines::SM_MODEL_3; } else - return SM_MODEL_2; + return Machines::SM_MODEL_2; } - return SM_MODEL_1; + return Machines::SM_MODEL_1; } std::vector > ElectricCalculation::ComplexMatrixTimesVector( -- cgit From 49ac4c4fcedb03ef09f8faab176159f48896dc7e Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Wed, 26 Jul 2017 21:32:35 -0300 Subject: 5 machine models implemented --- Project/ElectricCalculation.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) (limited to 'Project/ElectricCalculation.cpp') diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp index 2fcab57..e67c9b4 100644 --- a/Project/ElectricCalculation.cpp +++ b/Project/ElectricCalculation.cpp @@ -681,15 +681,16 @@ Machines::SyncMachineModel ElectricCalculation::GetMachineModel(SyncGenerator* g auto data = generator->GetElectricalData(); if(data.transTd0 != 0.0) { if(data.transTq0 != 0.0) { - if(data.subTd0 != 0.0) { - if(data.subTq0 != 0.0) - return Machines::SM_MODEL_5; - else - return Machines::SM_MODEL_4; - } else - return Machines::SM_MODEL_3; - } else + if(data.subTd0 != 0.0 || data.subTq0 != 0.0) { + return Machines::SM_MODEL_5; + } + return Machines::SM_MODEL_3; + } else { + if(data.subTd0 != 0.0 || data.subTq0 != 0.0) { + return Machines::SM_MODEL_4; + } return Machines::SM_MODEL_2; + } } return Machines::SM_MODEL_1; -- cgit