diff options
Diffstat (limited to 'Project/ElectricCalculation.cpp')
| -rw-r--r-- | Project/ElectricCalculation.cpp | 282 |
1 files changed, 242 insertions, 40 deletions
diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp index d024364..e67c9b4 100644 --- a/Project/ElectricCalculation.cpp +++ b/Project/ElectricCalculation.cpp @@ -7,6 +7,7 @@ ElectricCalculation::ElectricCalculation() {} ElectricCalculation::~ElectricCalculation() {} void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList) { + m_powerElementList.clear(); m_busList.clear(); m_capacitorList.clear(); m_indMotorList.clear(); @@ -19,6 +20,8 @@ void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList) // TODO: Bad design? for(auto it = elementList.begin(); it != elementList.end(); it++) { Element* element = *it; + m_powerElementList.push_back(static_cast<PowerElement*>(element)); + if(Bus* bus = dynamic_cast<Bus*>(element)) m_busList.push_back(bus); else if(Capacitor* capacitor = dynamic_cast<Capacitor*>(element)) @@ -41,9 +44,10 @@ void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList) } bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, - double systemPowerBase, - YBusSequence sequence, - bool includeSyncMachines) + double systemPowerBase, + YBusSequence sequence, + bool includeSyncMachines, + bool allLoadsAsImpedances) { if(m_busList.size() == 0) return false; @@ -61,7 +65,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > 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++; @@ -71,10 +75,16 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > 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])->GetEletricalData().number; + int n = static_cast<Bus*>(load->GetParentList()[0])->GetElectricalData().number; LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); - if(data.loadType == CONST_IMPEDANCE) - yBus[n][n] += std::complex<double>(data.activePower, -data.reactivePower); + if(data.loadType == CONST_IMPEDANCE || allLoadsAsImpedances) { + std::complex<double> yLoad = std::complex<double>(data.activePower, -data.reactivePower); + if(allLoadsAsImpedances) { + std::complex<double> v = static_cast<Bus*>(load->GetParentList()[0])->GetElectricalData().voltage; + yLoad /= (std::abs(v) * std::abs(v)); + } + yBus[n][n] += yLoad; + } } } @@ -82,7 +92,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > 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])->GetEletricalData().number; + 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); } @@ -92,7 +102,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > 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])->GetEletricalData().number; + 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); } @@ -104,8 +114,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > if(line->IsOnline()) { LineElectricalData data = line->GetElectricalData(); - int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetEletricalData().number; + int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number; switch(sequence) { case POSITIVE_SEQ: @@ -140,8 +150,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > if(transformer->IsOnline()) { TransformerElectricalData data = transformer->GetElectricalData(); - int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetEletricalData().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. @@ -160,8 +170,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > else if(sequence != ZERO_SEQ) { // 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)); + 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); @@ -181,10 +191,11 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > switch(data.connection) { case GWYE_GWYE: { std::complex<double> y = - 1.0 / std::complex<double>(data.zeroResistance + - 3.0 * (data.primaryGrndResistance + data.secondaryGrndResistance), - data.zeroIndReactance + - 3.0 * (data.primaryGrndReactance + data.secondaryGrndReactance)); + 1.0 / + std::complex<double>( + data.zeroResistance + 3.0 * (data.primaryGrndResistance + data.secondaryGrndResistance), + data.zeroIndReactance + + 3.0 * (data.primaryGrndReactance + data.secondaryGrndReactance)); std::complex<double> a = std::complex<double>(data.turnsRatio, 0.0); yBus[n1][n1] += y / (a * a); @@ -195,14 +206,14 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > case DELTA_GWYE: { std::complex<double> y = 1.0 / std::complex<double>(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<double> y = 1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.primaryGrndResistance), - data.zeroIndReactance + 3.0 * (data.primaryGrndReactance)); + data.zeroIndReactance + 3.0 * (data.primaryGrndReactance)); yBus[n1][n1] += y; break; } @@ -218,7 +229,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > 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])->GetEletricalData().number; + int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetElectricalData().number; SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); switch(sequence) { case POSITIVE_SEQ: { @@ -237,7 +248,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > 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])->GetEletricalData().number; + int n = static_cast<Bus*>(syncMotor->GetParentList()[0])->GetElectricalData().number; SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); switch(sequence) { case POSITIVE_SEQ: { @@ -258,10 +269,10 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > } void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<double> > voltage, - std::vector<std::complex<double> > power, - std::vector<BusType> busType, - std::vector<ReactiveLimits> reactiveLimit, - double systemPowerBase) + std::vector<std::complex<double> > power, + std::vector<BusType> busType, + std::vector<ReactiveLimits> 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; @@ -270,7 +281,7 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl // Buses voltages for(int i = 0; i < (int)m_busList.size(); i++) { Bus* bus = m_busList[i]; - BusElectricalData data = bus->GetEletricalData(); + BusElectricalData data = bus->GetElectricalData(); data.voltage = voltage[i]; bus->SetElectricalData(data); } @@ -279,17 +290,17 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl for(int i = 0; i < (int)m_lineList.size(); i++) { Line* line = m_lineList[i]; if(line->IsOnline()) { - int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetEletricalData().number; + int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetElectricalData().number; LineElectricalData data = line->GetElectricalData(); std::complex<double> v1 = voltage[n1]; std::complex<double> v2 = voltage[n2]; data.current[0] = (v1 - v2) / std::complex<double>(data.resistance, data.indReactance) + - v1 * std::complex<double>(0.0, data.capSusceptance / 2.0); + v1 * std::complex<double>(0.0, data.capSusceptance / 2.0); data.current[1] = (v2 - v1) / std::complex<double>(data.resistance, data.indReactance) + - v2 * std::complex<double>(0.0, data.capSusceptance / 2.0); + v2 * std::complex<double>(0.0, data.capSusceptance / 2.0); data.powerFlow[0] = v1 * std::conj(data.current[0]); data.powerFlow[1] = v2 * std::conj(data.current[1]); @@ -308,10 +319,10 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl Transformer* transformer = m_transformerList[i]; if(transformer->IsOnline()) { TransformerElectricalData data = transformer->GetElectricalData(); - int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetEletricalData().number; - int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetEletricalData().number; - std::complex<double> v1 = voltage[n1]; // Primary voltage - std::complex<double> v2 = voltage[n2]; // Secondary voltage + int n1 = static_cast<Bus*>(transformer->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast<Bus*>(transformer->GetParentList()[1])->GetElectricalData().number; + std::complex<double> v1 = voltage[n1]; // Primary voltage + std::complex<double> v2 = voltage[n2]; // Secondary voltage // Transformer admitance std::complex<double> y = 1.0 / std::complex<double>(data.resistance, data.indReactance); @@ -343,7 +354,7 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl // Synchronous machines for(int i = 0; i < (int)m_busList.size(); i++) { Bus* bus = m_busList[i]; - BusElectricalData data = bus->GetEletricalData(); + BusElectricalData data = bus->GetElectricalData(); // Get the synchronous machines connected and calculate the load power on the bus. std::vector<SyncGenerator*> syncGeneratorsOnBus; @@ -491,8 +502,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl reactivePower = childData_PU.minReactive; reachedMachineLimit = true; } else if((!childData_PU.haveMaxReactive && reactiveLimit[i].limitReached == RL_MAX_REACHED) || - (!childData_PU.haveMinReactive && reactiveLimit[i].limitReached == RL_MIN_REACHED) || - (!childData_PU.haveMaxReactive && !childData_PU.haveMaxReactive)) { + (!childData_PU.haveMinReactive && reactiveLimit[i].limitReached == RL_MIN_REACHED) || + (!childData_PU.haveMaxReactive && !childData_PU.haveMaxReactive)) { reactivePower += exceededReactive; exceededReactive = 0.0; } @@ -532,7 +543,7 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl } bool ElectricCalculation::InvertMatrix(std::vector<std::vector<std::complex<double> > > matrix, - std::vector<std::vector<std::complex<double> > >& inverse) + std::vector<std::vector<std::complex<double> > >& inverse) { int order = static_cast<int>(matrix.size()); @@ -599,3 +610,194 @@ bool ElectricCalculation::InvertMatrix(std::vector<std::vector<std::complex<doub return true; } + +void ElectricCalculation::ABCtoDQ0(std::complex<double> 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<double>& complexValue) +{ + double real = qValue * std::cos(angle) - dValue * std::sin(angle); + double imag = qValue * std::sin(angle) + dValue * std::cos(angle); + complexValue = std::complex<double>(real, imag); +} + +std::vector<std::complex<double> > ElectricCalculation::GaussianElimination( + std::vector<std::vector<std::complex<double> > > matrix, + std::vector<std::complex<double> > array) +{ + //[Ref] http://pt.wikipedia.org/wiki/Elimina%C3%A7%C3%A3o_de_Gauss + + std::vector<std::complex<double> > solution; + + std::vector<std::vector<std::complex<double> > > 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<double>(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<double>(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; +} + +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 || 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; +} + +std::vector<std::complex<double> > ElectricCalculation::ComplexMatrixTimesVector( + std::vector<std::vector<std::complex<double> > > matrix, + std::vector<std::complex<double> > vector) +{ + std::vector<std::complex<double> > solution; + for(unsigned int i = 0; i < matrix.size(); i++) { + solution.push_back(std::complex<double>(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<std::vector<std::complex<double> > > matrix, + std::vector<std::vector<std::complex<double> > >& matrixL, + std::vector<std::vector<std::complex<double> > >& 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<int>(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]; + } + + // Lower matrix main diagonal. + for(int i = 1; i < size; i++) { + matrixL[i][i] = std::complex<double>(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<std::complex<double> > ElectricCalculation::LUEvaluate(std::vector<std::vector<std::complex<double> > > u, + std::vector<std::vector<std::complex<double> > > l, + std::vector<std::complex<double> > b) +{ + int size = static_cast<int>(b.size()); + std::vector<std::complex<double> > x; + std::vector<std::complex<double> > 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; +} |