diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | Project/ElectricCalculation.cpp | 213 | ||||
| -rw-r--r-- | Project/ElectricCalculation.h | 125 | ||||
| -rw-r--r-- | Project/Fault.cpp | 238 | ||||
| -rw-r--r-- | Project/Fault.h | 37 | ||||
| -rw-r--r-- | Project/PowerFlow.cpp | 55 | ||||
| -rw-r--r-- | Project/PowerFlow.h | 19 |
7 files changed, 620 insertions, 68 deletions
@@ -4,6 +4,7 @@ .codelite/ Project/Release/ Project/Debug/ +doxygen/ # Compiled source # ################### diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp index c78864c..8f8fa35 100644 --- a/Project/ElectricCalculation.cpp +++ b/Project/ElectricCalculation.cpp @@ -37,7 +37,10 @@ void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList) } } -bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, double systemPowerBase) +bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, + double systemPowerBase, + YBusSequence sequence, + bool includeSyncMachines) { if(m_busList.size() == 0) return false; @@ -62,8 +65,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > } // Load - for(auto itlo = m_loadList.begin(); itlo != m_loadList.end(); itlo++) { - Load* load = *itlo; + 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; LoadElectricalData data = load->GetPUElectricalData(systemPowerBase); @@ -73,8 +76,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > } // Capacitor - for(auto itca = m_capacitorList.begin(); itca != m_capacitorList.end(); itca++) { - Capacitor* capacitor = *itca; + 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; CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase); @@ -83,8 +86,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > } // Inductor - for(auto itin = m_inductorList.begin(); itin != m_inductorList.end(); itin++) { - Inductor* inductor = *itin; + 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; InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase); @@ -93,28 +96,43 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > } // Power line - for(auto itl = m_lineList.begin(); itl != m_lineList.end(); itl++) { - Line* line = *itl; + for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) { + Line* line = *it; 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; - yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance); - yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance); + switch(sequence) { + case POSITIVE_SEQ: + case NEGATIVE_SEQ: { + yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance); + yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance); - yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance); - yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance); + yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance); + yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance); - yBus[n1][n1] += std::complex<double>(0.0, data.capSusceptance / 2.0); - yBus[n2][n2] += std::complex<double>(0.0, data.capSusceptance / 2.0); + yBus[n1][n1] += std::complex<double>(0.0, data.capSusceptance / 2.0); + yBus[n2][n2] += std::complex<double>(0.0, data.capSusceptance / 2.0); + } break; + case ZERO_SEQ: { + yBus[n1][n2] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance); + yBus[n2][n1] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance); + + yBus[n1][n1] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance); + yBus[n2][n2] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance); + + yBus[n1][n1] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0); + yBus[n2][n2] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0); + } + } } } // Transformer - for(auto itt = m_transformerList.begin(); itt != m_transformerList.end(); ++itt) { - Transformer* transformer = *itt; + for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) { + Transformer* transformer = *it; if(transformer->IsOnline()) { TransformerElectricalData data = transformer->GetElectricalData(); @@ -124,7 +142,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > // 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) { + if(data.turnsRatio == 1.0 && data.phaseShift == 0.0 && sequence != ZERO_SEQ) { yBus[n1][n2] += -1.0 / std::complex<double>(data.resistance, data.indReactance); yBus[n2][n1] += -1.0 / std::complex<double>(data.resistance, data.indReactance); @@ -136,7 +154,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > //[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 { + else if(sequence != ZERO_SEQ) { // Complex turns ratio double radPhaseShift = wxDegToRad(data.phaseShift); std::complex<double> a = std::complex<double>( @@ -145,10 +163,90 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > // Transformer admitance std::complex<double> y = 1.0 / std::complex<double>(data.resistance, data.indReactance); - 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; + if(sequence == POSITIVE_SEQ) { + 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; + } else if(sequence == NEGATIVE_SEQ) { + yBus[n1][n1] += y / std::pow(std::abs(a), 2.0); + yBus[n1][n2] += -(y / a); + yBus[n2][n1] += -(y / std::conj(a)); + yBus[n2][n2] += y; + } + } else if(sequence == ZERO_SEQ) { + 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)); + std::complex<double> a = std::complex<double>(data.turnsRatio, 0.0); + + yBus[n1][n1] += y / (a * a); + yBus[n1][n2] += -(y / a); + yBus[n2][n1] += -(y / a); + yBus[n2][n2] += y; + } break; + case DELTA_GWYE: { + std::complex<double> y = + 1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.secondaryGrndResistance), + 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)); + yBus[n1][n1] += y; + break; + } + default: + break; + } + } + } + } + + if(includeSyncMachines) { + // 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])->GetEletricalData().number; + SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase); + switch(sequence) { + case POSITIVE_SEQ: { + yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance); + } break; + case NEGATIVE_SEQ: { + yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance); + } break; + case ZERO_SEQ: { + yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroReactance); + } break; + } + } + } + // 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])->GetEletricalData().number; + SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase); + switch(sequence) { + case POSITIVE_SEQ: { + yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance); + } break; + case NEGATIVE_SEQ: { + yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance); + } break; + case ZERO_SEQ: { + yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroReactance); + } break; + } } } } @@ -429,3 +527,72 @@ 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) +{ + int order = static_cast<int>(matrix.size()); + + inverse.clear(); + // Fill the inverse matrix with identity. + for(int i = 0; i < order; ++i) { + std::vector<std::complex<double> > line; + for(int j = 0; j < order; ++j) { + line.push_back(i == j ? std::complex<double>(1.0, 0.0) : std::complex<double>(0.0, 0.0)); + } + inverse.push_back(line); + } + + // Check if a main diagonal value of the matrix is zero, if one is zero, try a linear combination to remove it. + for(int i = 0; i < order; ++i) { + for(int j = 0; j < order; ++j) { + if(i == j && matrix[i][j] == std::complex<double>(0.0, 0.0)) { + int row = 0; + while(row < order) { + if(matrix[row][j] != std::complex<double>(0.0, 0.0)) { + for(int k = 0; k < order; ++k) { + matrix[i][k] += matrix[row][k]; + inverse[i][k] += inverse[row][k]; + } + break; + } + row++; + } + // If all line values are zero, the matrix is singular and the solution is impossible. + if(row == order) return false; + } + } + } + + // Linear combinations are made in both matrices, the goal is the input matrix become the identity. The final result + // have two matrices: the identity and the inverse of the input. + for(int i = 0; i < order; ++i) { + for(int j = 0; j < order; ++j) { + if(i != j) { + if(matrix[i][i] == std::complex<double>(0.0, 0.0)) return false; + + std::complex<double> factor = matrix[j][i] / matrix[i][i]; + for(int k = 0; k < order; ++k) { + matrix[j][k] -= factor * matrix[i][k]; + inverse[j][k] -= factor * inverse[i][k]; + } + } + } + } + // Main diagonal calculation. + for(int i = 0; i < order; ++i) { + for(int j = 0; j < order; ++j) { + if(i == j) { + if(matrix[i][j] == std::complex<double>(0.0, 0.0)) return false; + + std::complex<double> factor = (matrix[i][j] - std::complex<double>(1.0, 0.0)) / matrix[i][j]; + for(int k = 0; k < order; ++k) { + matrix[j][k] -= factor * matrix[i][k]; + inverse[j][k] -= factor * inverse[i][k]; + } + } + } + } + + return true; +} diff --git a/Project/ElectricCalculation.h b/Project/ElectricCalculation.h index e4e4ef7..3eea5dd 100644 --- a/Project/ElectricCalculation.h +++ b/Project/ElectricCalculation.h @@ -18,14 +18,16 @@ enum BusType { BUS_SLACK = 0, BUS_PV, BUS_PQ }; enum ReactiveLimitsType { - RL_UNLIMITED = 0, // The bus can generate any ammount of reactive power. - RL_LIMITED, // The bus reactive power generation is limited. - RL_UNLIMITED_SOURCE, // The bus have at least one source of infinite reative power. - RL_MAX_REACHED, // Max limit reached - RL_MIN_REACHED, // Min limit reached - RL_NONE_REACHED // No limits reached + RL_UNLIMITED = 0, // The bus can generate any ammount of reactive power. + RL_LIMITED, // The bus reactive power generation is limited. + RL_UNLIMITED_SOURCE, // The bus have at least one source of infinite reative power. + RL_MAX_REACHED, // Max limit reached + RL_MIN_REACHED, // Min limit reached + RL_NONE_REACHED // No limits reached }; +enum YBusSequence { POSITIVE_SEQ = 0, NEGATIVE_SEQ, ZERO_SEQ }; + struct ReactiveLimits { double maxLimit = 0.0; double minLimit = 0.0; @@ -34,30 +36,123 @@ struct ReactiveLimits { ReactiveLimitsType limitReached = RL_NONE_REACHED; }; +/** + * @class ElectricCalculation + * @author Thales Lima Oliveira + * @date 09/01/2017 + * @file ElectricCalculation.h + * @brief Base class of electric calculations, with general methods. + */ class ElectricCalculation { - public: +public: + /** + * @brief Constructor. + */ ElectricCalculation(); + + /** + * @brief Destructor. + */ ~ElectricCalculation(); + + /** + * @brief Separate the elements from a generic list. + * @param elementList List of generic elements. + */ virtual void GetElementsFromList(std::vector<Element*> elementList); - virtual bool GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, double systemPowerBase); + + /** + * @brief Get the admittance matrix from the list of elements (use GetElementsFromList first). + * @param yBus Admittance matrix. The previous content will be erased. + * @param systemPowerBase Base power of the system. + * @param sequence Sequence of admittance matrix (positive, negative and zero). + * @param includeSyncMachines Include the synchronous machines on calculation. + * @return Return true if was possible to build the admittance matrix. + */ + virtual bool GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, + double systemPowerBase, + YBusSequence sequence = POSITIVE_SEQ, + bool includeSyncMachines = false); + + /** + * @brief Invert a matrix. + * @param matrix Matrix to invert. + * @param inverse Inverted matrix. The previous content will be erased. + * @return Return true if was possible to invert the matrix. + */ + virtual bool InvertMatrix(std::vector<std::vector<std::complex<double> > > matrix, + std::vector<std::vector<std::complex<double> > >& inverse); + + /** + * @brief Update the elements after the power flow calculation. + * @param voltage Array with the buses voltages. + * @param power Array with the buses injected power. + * @param busType Array with the buses type. + * @param reactiveLimit Array with the reactive limit data. + * @param systemPowerBase Base power of the system. + */ virtual void 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); + /** + * @brief Get the buses of the system (use GetElementsFromList first). + * @return A list of bus elements. + */ const std::vector<Bus*> GetBusList() const { return m_busList; } + + /** + * @brief Get the capacitors of the system (use GetElementsFromList first). + * @return A list of capacitor elements. + */ const std::vector<Capacitor*> GetCapacitorList() const { return m_capacitorList; } + + /** + * @brief Get the induction motors of the system (use GetElementsFromList first). + * @return A list of induction motor elements. + */ const std::vector<IndMotor*> GetIndMotorList() const { return m_indMotorList; } + + /** + * @brief Get the inductors of the system (use GetElementsFromList first). + * @return A list of inductor elements. + */ const std::vector<Inductor*> GetInductorList() const { return m_inductorList; } + + /** + * @brief Get the lines of the system (use GetElementsFromList first). + * @return A list of line elements. + */ const std::vector<Line*> GetLineList() const { return m_lineList; } + + /** + * @brief Get the loads of the system (use GetElementsFromList first). + * @return A list of load elements. + */ const std::vector<Load*> GetLoadList() const { return m_loadList; } + + /** + * @brief Get the synchronous generators of the system (use GetElementsFromList first). + * @return A list of synchronous generator elements. + */ const std::vector<SyncGenerator*> GetSyncGeneratorList() const { return m_syncGeneratorList; } + + /** + * @brief Get the synchronous motors of the system (use GetElementsFromList first). + * @return A list of synchronous motor elements. + */ const std::vector<SyncMotor*> GetSyncMotorList() const { return m_syncMotorList; } - const std::vector<Transformer*> GetTransformerList() const { return m_transformerList; } - protected: + /** + * @brief Get the transformers of the system (use GetElementsFromList first). + * @return A list of transformer elements. + */ + const std::vector<Transformer*> GetTransformerList() const { return m_transformerList; } + +protected: std::vector<Bus*> m_busList; std::vector<Capacitor*> m_capacitorList; std::vector<IndMotor*> m_indMotorList; @@ -69,4 +164,4 @@ class ElectricCalculation std::vector<Transformer*> m_transformerList; }; -#endif // ELECTRICCALCULATION_H +#endif // ELECTRICCALCULATION_H diff --git a/Project/Fault.cpp b/Project/Fault.cpp new file mode 100644 index 0000000..058b4ea --- /dev/null +++ b/Project/Fault.cpp @@ -0,0 +1,238 @@ +#include "Fault.h" + +Fault::Fault() + : ElectricCalculation() +{ +} + +Fault::Fault(std::vector<Element*> elementList) { GetElementsFromList(elementList); } + +Fault::~Fault() {} + +bool Fault::RunFaultCalculation(double systemPowerBase) +{ + int numberOfBuses = static_cast<int>(m_busList.size()); + if(numberOfBuses == 0) { + m_errorMsg = _("There is no buses in the system."); + return false; + } + + // Get adimittance matrices. + std::vector<std::vector<std::complex<double> > > yBusPos; + GetYBus(yBusPos, systemPowerBase, POSITIVE_SEQ, true); + std::vector<std::vector<std::complex<double> > > yBusNeg; + GetYBus(yBusNeg, systemPowerBase, NEGATIVE_SEQ, true); + std::vector<std::vector<std::complex<double> > > yBusZero; + GetYBus(yBusZero, systemPowerBase, ZERO_SEQ, true); + + // Calculate the impedance matrices. + if(!InvertMatrix(yBusPos, m_zBusPos)) { + m_errorMsg = _("Fail to invert the positive sequence admittance matrix."); + return false; + } + if(!InvertMatrix(yBusNeg, m_zBusNeg)) { + m_errorMsg = _("Fail to invert the negative sequence admittance matrix."); + return false; + } + if(!InvertMatrix(yBusZero, m_zBusZero)) { + m_errorMsg = _("Fail to invert the zero sequence admittance matrix."); + return false; + } + + // Pre-fault voltages (power flow solution). + std::vector<std::complex<double> > preFaultVoltages; + + // Get fault parameters. + int fNumber = -1; + FaultData fType = FAULT_THREEPHASE; + FaultData fLocation = FAULT_LINE_A; + std::complex<double> fImpedance = std::complex<double>(0.0, 0.0); + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + BusElectricalData data = bus->GetEletricalData(); + preFaultVoltages.push_back(data.voltage); + if(data.hasFault) { + fNumber = data.number; + fType = data.faultType; + fLocation = data.faultLocation; + fImpedance = std::complex<double>(data.faultResistance, data.faultReactance); + } + } + + if(fNumber == -1) { + m_errorMsg = _("There is no fault in the system."); + return false; + } + + // Fault calculation. + std::complex<double> fCurrentPos = std::complex<double>(0.0, 0.0); + std::complex<double> fCurrentNeg = std::complex<double>(0.0, 0.0); + std::complex<double> fCurrentZero = std::complex<double>(0.0, 0.0); + + std::complex<double> preFaultVoltage = preFaultVoltages[fNumber]; + std::complex<double> a = std::complex<double>(-0.5, 0.866025403784); + std::complex<double> a2 = std::complex<double>(-0.5, -0.866025403784); + + switch(fType) { + case FAULT_THREEPHASE: { + fCurrentPos = preFaultVoltage / (m_zBusPos[fNumber][fNumber] + fImpedance); + } break; + case FAULT_2LINE: { + fCurrentPos = preFaultVoltage / (m_zBusPos[fNumber][fNumber] + m_zBusNeg[fNumber][fNumber] + fImpedance); + + switch(fLocation) { + case FAULT_LINE_A: { + fCurrentNeg = -a2 * fCurrentPos; + } break; + case FAULT_LINE_B: { + fCurrentNeg = -fCurrentPos; + } break; + case FAULT_LINE_C: { + fCurrentNeg = -a * fCurrentPos; + } break; + default: + break; + } + } break; + case FAULT_2LINE_GROUND: { + std::complex<double> z1 = m_zBusPos[fNumber][fNumber]; + std::complex<double> z2 = m_zBusNeg[fNumber][fNumber]; + std::complex<double> z0 = m_zBusZero[fNumber][fNumber]; + std::complex<double> zf_3 = std::complex<double>(3.0, 0.0) * fImpedance; + + fCurrentPos = (preFaultVoltage * (z2 + z0 + zf_3)) / (z1 * z2 + z2 * z0 + z2 * zf_3 + z1 * z0 + z1 * zf_3); + + switch(fLocation) { + case FAULT_LINE_A: { + fCurrentNeg = -a2 * ((preFaultVoltage - z1 * fCurrentPos) / z2); + fCurrentZero = -a * ((preFaultVoltage - z1 * fCurrentPos) / (z0 + zf_3)); + } break; + case FAULT_LINE_B: { + fCurrentNeg = -((preFaultVoltage - z1 * fCurrentPos) / z2); + fCurrentZero = -((preFaultVoltage - z1 * fCurrentPos) / (z0 + zf_3)); + } break; + case FAULT_LINE_C: { + fCurrentNeg = -a * ((preFaultVoltage - z1 * fCurrentPos) / z2); + fCurrentZero = -a2 * ((preFaultVoltage - z1 * fCurrentPos) / (z0 + zf_3)); + } break; + default: + break; + } + } break; + case FAULT_LINE_GROUND: { + fCurrentPos = + preFaultVoltage / (m_zBusPos[fNumber][fNumber] + m_zBusNeg[fNumber][fNumber] + + m_zBusZero[fNumber][fNumber] + std::complex<double>(3.0, 0.0) * fImpedance); + switch(fLocation) { + case FAULT_LINE_A: { + fCurrentNeg = fCurrentPos; + fCurrentZero = fCurrentPos; + } break; + case FAULT_LINE_B: { + fCurrentNeg = a * fCurrentPos; + fCurrentZero = a2 * fCurrentPos; + } break; + case FAULT_LINE_C: { + fCurrentNeg = a2 * fCurrentPos; + fCurrentZero = a * fCurrentPos; + } break; + default: + break; + } + } break; + default: + break; + } + + // Convert sequence currents to ABC. [Iabc] = [A]*[I012] + m_fCurrentA = fCurrentPos + fCurrentNeg + fCurrentZero; + m_fCurrentB = fCurrentPos + a2 * fCurrentNeg + a * fCurrentZero; + m_fCurrentC = fCurrentPos + a * fCurrentNeg + a2 * fCurrentZero; + + // Pos-fault voltages calculation + m_posFaultVoltagePos.clear(); + m_posFaultVoltageNeg.clear(); + m_posFaultVoltageZero.clear(); + m_posFaultVoltageA.clear(); + m_posFaultVoltageB.clear(); + m_posFaultVoltageC.clear(); + + for(int i = 0; i < numberOfBuses; ++i) { + m_posFaultVoltagePos.push_back(preFaultVoltages[i] - m_zBusPos[i][fNumber] * fCurrentPos); + m_posFaultVoltageNeg.push_back(-m_zBusNeg[i][fNumber] * fCurrentNeg); + m_posFaultVoltageZero.push_back(-m_zBusZero[i][fNumber] * fCurrentZero); + + // V012 -> Vabc + m_posFaultVoltageA.push_back(m_posFaultVoltagePos[i] + m_posFaultVoltageNeg[i] + m_posFaultVoltageZero[i]); + m_posFaultVoltageB.push_back( + m_posFaultVoltagePos[i] + a2 * m_posFaultVoltageNeg[i] + a * m_posFaultVoltageZero[i]); + m_posFaultVoltageC.push_back( + m_posFaultVoltagePos[i] + a * m_posFaultVoltageNeg[i] + a2 * m_posFaultVoltageZero[i]); + } + + UpdateElementsFault(systemPowerBase); + return true; +} + +void Fault::UpdateElementsFault(double systemPowerBase) +{ + std::complex<double> a = std::complex<double>(-0.5, 0.866025403784); + std::complex<double> a2 = std::complex<double>(-0.5, -0.866025403784); + + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + auto data = bus->GetEletricalData(); + if(data.hasFault) { + data.faultCurrent[0] = m_fCurrentA; + data.faultCurrent[1] = m_fCurrentB; + data.faultCurrent[2] = m_fCurrentC; + } else { + data.faultCurrent[0] = data.faultCurrent[1] = data.faultCurrent[2] = std::complex<double>(0.0, 0.0); + } + data.faultVoltage[0] = m_posFaultVoltageA[data.number]; + data.faultVoltage[1] = m_posFaultVoltageB[data.number]; + data.faultVoltage[2] = m_posFaultVoltageC[data.number]; + bus->SetElectricalData(data); + } + + for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) { + Line* line = *it; + if(line->IsOnline()) { + int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetEletricalData().number; + int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetEletricalData().number; + auto data = line->GetElectricalData(); + std::complex<double> vPos[2] = { m_posFaultVoltagePos[n1], m_posFaultVoltagePos[n2] }; + std::complex<double> vNeg[2] = { m_posFaultVoltageNeg[n1], m_posFaultVoltageNeg[n2] }; + std::complex<double> vZero[2] = { m_posFaultVoltageZero[n1], m_posFaultVoltageZero[n2] }; + std::complex<double> zPos(data.resistance, data.indReactance); + std::complex<double> bPos(0.0, data.capSusceptance / 2.0); + std::complex<double> zZero(data.zeroResistance, data.zeroIndReactance); + std::complex<double> bZero(0.0, data.zeroCapSusceptance / 2.0); + + std::complex<double> lineCurrentPos[2]; + std::complex<double> lineCurrentNeg[2]; + std::complex<double> lineCurrentZero[2]; + + lineCurrentPos[0] = ((vPos[0] - vPos[1]) / zPos) + (vPos[0] * bPos); + lineCurrentNeg[0] = ((vNeg[0] - vNeg[1]) / zPos) + (vNeg[0] * bPos); + lineCurrentZero[0] = ((vZero[0] - vZero[1]) / zZero) + (vZero[0] * bZero); + lineCurrentPos[1] = ((vPos[1] - vPos[0]) / zPos) + (vPos[1] * bPos); + lineCurrentNeg[1] = ((vNeg[1] - vNeg[0]) / zPos) + (vNeg[1] * bPos); + lineCurrentZero[1] = ((vZero[1] - vZero[0]) / zZero) + (vZero[1] * bZero); + + data.faultCurrent[0][0] = lineCurrentPos[0] + lineCurrentNeg[0] + lineCurrentZero[0]; + data.faultCurrent[1][0] = lineCurrentPos[0] + a2 * lineCurrentNeg[0] + a * lineCurrentZero[0]; + data.faultCurrent[2][0] = lineCurrentPos[0] + a * lineCurrentNeg[0] + a2 * lineCurrentZero[0]; + data.faultCurrent[0][1] = lineCurrentPos[1] + lineCurrentNeg[1] + lineCurrentZero[1]; + data.faultCurrent[1][1] = lineCurrentPos[1] + a2 * lineCurrentNeg[1] + a * lineCurrentZero[1]; + data.faultCurrent[2][1] = lineCurrentPos[1] + a * lineCurrentNeg[1] + a2 * lineCurrentZero[1]; + + line->SetElectricalData(data); + } + } + + for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) { + Transformer* transformer = *it; + + } +} diff --git a/Project/Fault.h b/Project/Fault.h new file mode 100644 index 0000000..dfe0e75 --- /dev/null +++ b/Project/Fault.h @@ -0,0 +1,37 @@ +#ifndef FAULT_H +#define FAULT_H + +#include "ElectricCalculation.h" + +class Fault : public ElectricCalculation +{ +public: + Fault(std::vector<Element*> elementList); + Fault(); + ~Fault(); + + virtual bool RunFaultCalculation(double systemPowerBase); + virtual void UpdateElementsFault(double systemPowerBase); + virtual wxString GetErrorMessage() { return m_errorMsg; } + +protected: + wxString m_errorMsg = ""; + + std::vector<std::vector<std::complex<double> > > m_zBusPos; + std::vector<std::vector<std::complex<double> > > m_zBusNeg; + std::vector<std::vector<std::complex<double> > > m_zBusZero; + + std::vector<std::complex<double> > m_posFaultVoltagePos; + std::vector<std::complex<double> > m_posFaultVoltageNeg; + std::vector<std::complex<double> > m_posFaultVoltageZero; + + std::complex<double> m_fCurrentA; + std::complex<double> m_fCurrentB; + std::complex<double> m_fCurrentC; + + std::vector<std::complex<double> > m_posFaultVoltageA; + std::vector<std::complex<double> > m_posFaultVoltageB; + std::vector<std::complex<double> > m_posFaultVoltageC; +}; + +#endif // FAULT_H diff --git a/Project/PowerFlow.cpp b/Project/PowerFlow.cpp index c94d8d8..bfee2af 100644 --- a/Project/PowerFlow.cpp +++ b/Project/PowerFlow.cpp @@ -1,12 +1,23 @@ #include "PowerFlow.h" -PowerFlow::PowerFlow(std::vector<Element*> elementList) : ElectricCalculation() { GetElementsFromList(elementList); } +PowerFlow::PowerFlow() + : ElectricCalculation() +{ +} + +PowerFlow::PowerFlow(std::vector<Element*> elementList) + : ElectricCalculation() +{ + GetElementsFromList(elementList); +} + PowerFlow::~PowerFlow() {} + bool PowerFlow::RunGaussSeidel(double systemPowerBase, - int maxIteration, - double error, - double initAngle, - double accFactor) + int maxIteration, + double error, + double initAngle, + double accFactor) { // Calculate the Ybus. if(!GetYBus(m_yBus, systemPowerBase)) { @@ -15,13 +26,13 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase, } // Number of buses on the system. - int numberOfBuses = (int)m_busList.size(); + int numberOfBuses = static_cast<int>(m_busList.size()); - std::vector<BusType> busType; // Bus type - std::vector<std::complex<double> > voltage; // Voltage of buses - std::vector<std::complex<double> > power; // Injected power - std::vector<std::complex<double> > loadPower; // Only the load power - std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses + std::vector<BusType> busType; // Bus type + std::vector<std::complex<double> > voltage; // Voltage of buses + std::vector<std::complex<double> > power; // Injected power + std::vector<std::complex<double> > loadPower; // Only the load power + std::vector<ReactiveLimits> reactiveLimit; // Limit of reactive power on PV buses reactiveLimit.resize(numberOfBuses); @@ -61,7 +72,7 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase, } // Fill the power array - power.push_back(std::complex<double>(0.0, 0.0)); // Initial value + power.push_back(std::complex<double>(0.0, 0.0)); // Initial value loadPower.push_back(std::complex<double>(0.0, 0.0)); // Synchronous generator @@ -168,12 +179,12 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase, } // Gauss-Seidel method - std::vector<std::complex<double> > oldVoltage; // Old voltage array. + std::vector<std::complex<double> > oldVoltage; // Old voltage array. oldVoltage.resize(voltage.size()); auto oldBusType = busType; - int iteration = 0; // Current itaration number. + int iteration = 0; // Current itaration number. while(true) { // Reach the max number of iterations. @@ -203,7 +214,7 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase, // Apply the acceleration factor. newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(), - accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag()); + accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag()); voltage[i] = newVolt; } @@ -227,15 +238,15 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase, // Apply the acceleration factor. newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(), - accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag()); + accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag()); // Keep the same voltage magnitude. voltage[i] = std::complex<double>(std::abs(voltage[i]) * std::cos(std::arg(newVolt)), - std::abs(voltage[i]) * std::sin(std::arg(newVolt))); + std::abs(voltage[i]) * std::sin(std::arg(newVolt))); } - double busError = std::max(std::abs(voltage[i].real() - oldVoltage[i].real()), - std::abs(voltage[i].imag() - oldVoltage[i].imag())); + double busError = std::max( + std::abs(voltage[i].real() - oldVoltage[i].real()), std::abs(voltage[i].imag() - oldVoltage[i].imag())); if(busError > iterationError) iterationError = busError; } @@ -246,7 +257,8 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase, if(busType[i] == BUS_PV) { if(reactiveLimit[i].maxLimitType == RL_LIMITED) { if(power[i].imag() - loadPower[i].imag() > reactiveLimit[i].maxLimit) { - power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag()); + power[i] = + std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag()); busType[i] = BUS_PQ; reactiveLimit[i].limitReached = RL_MAX_REACHED; limitReach = true; @@ -254,7 +266,8 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase, } if(reactiveLimit[i].minLimitType == RL_LIMITED) { if(power[i].imag() - loadPower[i].imag() < reactiveLimit[i].minLimit) { - power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag()); + power[i] = + std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag()); busType[i] = BUS_PQ; reactiveLimit[i].limitReached = RL_MIN_REACHED; limitReach = true; diff --git a/Project/PowerFlow.h b/Project/PowerFlow.h index 1e5c621..29e205a 100644 --- a/Project/PowerFlow.h +++ b/Project/PowerFlow.h @@ -4,24 +4,25 @@ #include "ElectricCalculation.h" #include <wx/string.h> -#include <wx/intl.h> //_() -#include <wx/log.h> //temp +#include <wx/intl.h> //_() class PowerFlow : public ElectricCalculation { - public: +public: + PowerFlow(); PowerFlow(std::vector<Element*> elementList); ~PowerFlow(); virtual bool RunGaussSeidel(double systemPowerBase = 100e6, - int maxIteration = 5000, - double error = 1e-6, - double initAngle = 0.0, - double accFactor = 1.0); + int maxIteration = 5000, + double error = 1e-6, + double initAngle = 0.0, + double accFactor = 1.0); virtual wxString GetErrorMessage() { return m_errorMsg; } - protected: + +protected: std::vector<std::vector<std::complex<double> > > m_yBus; wxString m_errorMsg = ""; }; -#endif // POWERFLOW_H +#endif // POWERFLOW_H |