diff options
Diffstat (limited to 'Project/ElectricCalculation.cpp')
| -rw-r--r-- | Project/ElectricCalculation.cpp | 241 |
1 files changed, 218 insertions, 23 deletions
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<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; @@ -73,7 +77,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > if(load->IsOnline()) { int n = static_cast<Bus*>(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<double>(data.activePower, -data.reactivePower); } } @@ -160,8 +164,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 +185,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 +200,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; } @@ -258,10 +263,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; @@ -287,9 +292,9 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl 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]); @@ -310,8 +315,8 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl 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 + 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); @@ -491,8 +496,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 +537,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 +604,193 @@ 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; +} + +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<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]; + } + + // Upper 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; +} |