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/Electromechanical.cpp | 59 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 Project/Electromechanical.cpp (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp new file mode 100644 index 0000000..85e10b4 --- /dev/null +++ b/Project/Electromechanical.cpp @@ -0,0 +1,59 @@ +#include "Electromechanical.h" + +Electromechanical::Electromechanical(std::vector elementList) +{ + GetElementsFromList(elementList); + SetEventTimeList(); +} + +Electromechanical::~Electromechanical() {} +void Electromechanical::SetEventTimeList() +{ + // Fault + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + auto data = bus->GetEletricalData(); + if(data.stabHasFault) { + m_eventTimeList.push_back(data.stabFaultTime); + m_eventOccurrenceList.push_back(false); + m_eventTimeList.push_back(data.stabFaultTime + data.stabFaultLength); + m_eventOccurrenceList.push_back(false); + } + } + // Switching + for(auto it = m_powerElementList.begin(), itEnd = m_powerElementList.end(); it != itEnd; ++it) { + PowerElement* element = *it; + SwitchingData swData = element->GetSwitchingData(); + for(unsigned int i = 0; i < swData.swTime.size(); ++i) { + m_eventTimeList.push_back(swData.swTime[i]); + m_eventOccurrenceList.push_back(false); + } + } +} + +bool Electromechanical::HasEvent(double currentTime) +{ + for(unsigned int i = 0; i < m_eventTimeList.size(); ++i) { + if(!m_eventOccurrenceList[i]) { + if((m_eventTimeList[i] - m_timeStep) < currentTime && m_eventTimeList[i] >= currentTime) { + m_eventOccurrenceList[i] = true; + return true; + } + } + } + return false; +} + +void Electromechanical::SetEvent(double currentTime) {} +bool Electromechanical::RunStabilityCalculation() +{ + // test + double simTime = 10.0; + double currentTime = 0.0; + while(currentTime <= simTime) { + if(HasEvent(currentTime)) { + } + currentTime += m_timeStep; + } + return true; +} -- 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/Electromechanical.cpp | 317 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 314 insertions(+), 3 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 85e10b4..191a1da 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -12,7 +12,7 @@ void Electromechanical::SetEventTimeList() // Fault for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { Bus* bus = *it; - auto data = bus->GetEletricalData(); + auto data = bus->GetElectricalData(); if(data.stabHasFault) { m_eventTimeList.push_back(data.stabFaultTime); m_eventOccurrenceList.push_back(false); @@ -35,7 +35,7 @@ bool Electromechanical::HasEvent(double currentTime) { for(unsigned int i = 0; i < m_eventTimeList.size(); ++i) { if(!m_eventOccurrenceList[i]) { - if((m_eventTimeList[i] - m_timeStep) < currentTime && m_eventTimeList[i] >= currentTime) { + if(EventTrigger(m_eventTimeList[i], currentTime)) { m_eventOccurrenceList[i] = true; return true; } @@ -44,16 +44,327 @@ bool Electromechanical::HasEvent(double currentTime) return false; } -void Electromechanical::SetEvent(double currentTime) {} +void Electromechanical::SetEvent(double currentTime) +{ + // Fault + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + auto data = bus->GetElectricalData(); + if(data.stabHasFault) { + int n = data.number; + + // Insert fault + if(EventTrigger(data.stabFaultTime, currentTime)) { + double r, x; + r = data.stabFaultResistance; + x = data.stabFaultReactance; + if(x < 1e-5) x = 1e-5; + m_yBus[n][n] += std::complex(1.0, 0.0) / std::complex(r, x); + } + + // Remove fault + else if(EventTrigger(data.stabFaultTime + data.stabFaultLength, currentTime)) { + double r, x; + r = data.stabFaultResistance; + x = data.stabFaultReactance; + if(x < 1e-5) x = 1e-5; + m_yBus[n][n] -= std::complex(1.0, 0.0) / std::complex(r, x); + } + } + } + + // SyncGenerator switching + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* generator = *it; + auto swData = generator->GetSwitchingData(); + for(unsigned int i = 0; i < swData.swType.size(); ++i) { + if(EventTrigger(swData.swTime[i], currentTime)) { + // Remove machine (only connected machines) + if(swData.swType[i] == SW_REMOVE && generator->IsOnline()) { + generator->SetOnline(false); + auto data = generator->GetPUElectricalData(m_powerSystemBase); + int n = static_cast(generator->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] -= GetSyncMachineAdmittance(generator); + } + + // Insert machine (only disconnected machines) + if(swData.swType[i] == SW_INSERT && !generator->IsOnline() && generator->GetParentList().size() == 1) { + generator->SetOnline(true); + auto data = generator->GetPUElectricalData(m_powerSystemBase); + int n = static_cast(generator->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] += GetSyncMachineAdmittance(generator); + } + } + } + } + + // Load switching + for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) { + Load* load = *it; + auto swData = load->GetSwitchingData(); + for(unsigned int i = 0; i < swData.swType.size(); ++i) { + if(EventTrigger(swData.swTime[i], currentTime)) { + // Remove load (only connected loads) + if(swData.swType[i] == SW_REMOVE && load->IsOnline()) { + load->SetOnline(false); + auto data = load->GetPUElectricalData(m_powerSystemBase); + Bus* parentBus = static_cast(load->GetParentList()[0]); + int n = parentBus->GetElectricalData().number; + std::complex v = parentBus->GetElectricalData().voltage; + m_yBus[n][n] -= std::complex(data.activePower, -data.reactivePower) / (v * v); + } + + // Insert load (only disconnected load) + if(swData.swType[i] == SW_INSERT && !load->IsOnline() && load->GetParentList().size() == 1) { + load->SetOnline(true); + auto data = load->GetPUElectricalData(m_powerSystemBase); + Bus* parentBus = static_cast(load->GetParentList()[0]); + int n = parentBus->GetElectricalData().number; + std::complex v = parentBus->GetElectricalData().voltage; + m_yBus[n][n] += std::complex(data.activePower, -data.reactivePower) / (v * v); + } + } + } + } + + // Line switching + for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) { + Line* line = *it; + auto swData = line->GetSwitchingData(); + for(unsigned int i = 0; i < swData.swType.size(); ++i) { + if(EventTrigger(swData.swTime[i], currentTime)) { + // Remove line (only connected lines) + if(swData.swType[i] == SW_REMOVE && line->IsOnline()) { + line->SetOnline(false); + auto data = line->GetElectricalData(); + + int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(line->GetParentList()[1])->GetElectricalData().number; + + m_yBus[n1][n2] += 1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n1] += 1.0 / std::complex(data.resistance, data.indReactance); + + m_yBus[n1][n1] -= 1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n2] -= 1.0 / std::complex(data.resistance, data.indReactance); + + m_yBus[n1][n1] -= std::complex(0.0, data.capSusceptance / 2.0); + m_yBus[n2][n2] -= std::complex(0.0, data.capSusceptance / 2.0); + } + + // Insert line (only disconnected lines) + if(swData.swType[i] == SW_INSERT && !line->IsOnline() && line->GetParentList().size() == 2) { + if(line->GetParentList()[0] && line->GetParentList()[1]) { + line->SetOnline(true); + auto data = line->GetElectricalData(); + + int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(line->GetParentList()[1])->GetElectricalData().number; + + m_yBus[n1][n2] -= 1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n1] -= 1.0 / std::complex(data.resistance, data.indReactance); + + m_yBus[n1][n1] += 1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n2] += 1.0 / std::complex(data.resistance, data.indReactance); + + m_yBus[n1][n1] += std::complex(0.0, data.capSusceptance / 2.0); + m_yBus[n2][n2] += std::complex(0.0, data.capSusceptance / 2.0); + } + } + } + } + } + + // Transformer switching + for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) { + Transformer* transformer = *it; + auto swData = transformer->GetSwitchingData(); + for(unsigned int i = 0; i < swData.swType.size(); ++i) { + if(EventTrigger(swData.swTime[i], currentTime)) { + // Remove transformer (only connected transformers) + if(swData.swType[i] == SW_REMOVE && transformer->IsOnline()) { + transformer->SetOnline(false); + auto data = transformer->GetElectricalData(); + + int n1 = static_cast(transformer->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(transformer->GetParentList()[1])->GetElectricalData().number; + + if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) { + m_yBus[n1][n2] -= -1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n1] -= -1.0 / std::complex(data.resistance, data.indReactance); + + m_yBus[n1][n1] -= 1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n2] -= 1.0 / std::complex(data.resistance, data.indReactance); + } else { + // Complex turns ratio + double radPhaseShift = wxDegToRad(data.phaseShift); + 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); + m_yBus[n1][n1] -= y / std::pow(std::abs(a), 2.0); + m_yBus[n1][n2] -= -(y / std::conj(a)); + m_yBus[n2][n1] -= -(y / a); + m_yBus[n2][n2] -= y; + } + } + + // Insert transformer (only disconnected transformers) + if(swData.swType[i] == SW_INSERT && !transformer->IsOnline() && + transformer->GetParentList().size() == 2) { + if(transformer->GetParentList()[0] && transformer->GetParentList()[1]) { + transformer->SetOnline(true); + auto data = transformer->GetElectricalData(); + + int n1 = static_cast(transformer->GetParentList()[0])->GetElectricalData().number; + int n2 = static_cast(transformer->GetParentList()[1])->GetElectricalData().number; + + if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) { + m_yBus[n1][n2] += -1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n1] += -1.0 / std::complex(data.resistance, data.indReactance); + + m_yBus[n1][n1] += 1.0 / std::complex(data.resistance, data.indReactance); + m_yBus[n2][n2] += 1.0 / std::complex(data.resistance, data.indReactance); + } else { + // Complex turns ratio + double radPhaseShift = wxDegToRad(data.phaseShift); + 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); + m_yBus[n1][n1] += y / std::pow(std::abs(a), 2.0); + m_yBus[n1][n2] += -(y / std::conj(a)); + m_yBus[n2][n1] += -(y / a); + m_yBus[n2][n2] += y; + } + } + } + } + } + } + + // Capacitor switching + for(auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) { + Capacitor* capacitor = *it; + auto swData = capacitor->GetSwitchingData(); + for(unsigned int i = 0; i < swData.swType.size(); ++i) { + if(EventTrigger(swData.swTime[i], currentTime)) { + // Remove capacitor (only connected capacitors) + if(swData.swType[i] == SW_REMOVE && capacitor->IsOnline()) { + capacitor->SetOnline(false); + auto data = capacitor->GetPUElectricalData(m_powerSystemBase); + int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] += std::complex(0.0, data.reactivePower); + } + + // Insert capacitor (only disconnected capacitors) + if(swData.swType[i] == SW_INSERT && !capacitor->IsOnline() && capacitor->GetParentList().size() == 1) { + if(capacitor->GetParentList()[0]) { + capacitor->SetOnline(true); + auto data = capacitor->GetPUElectricalData(m_powerSystemBase); + int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] -= std::complex(0.0, data.reactivePower); + } + } + } + } + } + + // Capacitor switching + for(auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) { + Inductor* inductor = *it; + auto swData = inductor->GetSwitchingData(); + for(unsigned int i = 0; i < swData.swType.size(); ++i) { + if(EventTrigger(swData.swTime[i], currentTime)) { + // Remove inductor (only connected inductors) + if(swData.swType[i] == SW_REMOVE && inductor->IsOnline()) { + inductor->SetOnline(false); + auto data = inductor->GetPUElectricalData(m_powerSystemBase); + int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] += std::complex(0.0, -data.reactivePower); + } + + // Insert inductor (only disconnected inductors) + if(swData.swType[i] == SW_INSERT && !inductor->IsOnline() && inductor->GetParentList().size() == 1) { + if(inductor->GetParentList()[0]) { + inductor->SetOnline(true); + auto data = inductor->GetPUElectricalData(m_powerSystemBase); + int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] -= std::complex(0.0, -data.reactivePower); + } + } + } + } + } +} + bool Electromechanical::RunStabilityCalculation() { + // Calculate the admittance matrix with the synchronous machines. + if(!GetYBus(m_yBus, m_powerSystemBase, POSITIVE_SEQ, false, true)) { + m_errorMsg = _("It was not possible to build the admittance matrix."); + return false; + } + InsertSyncMachinesOnYBus(); + // test double simTime = 10.0; double currentTime = 0.0; while(currentTime <= simTime) { if(HasEvent(currentTime)) { + SetEvent(currentTime); } currentTime += m_timeStep; } return true; } + +void Electromechanical::InsertSyncMachinesOnYBus() +{ + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* generator = *it; + if(generator->IsOnline()) { + auto data = generator->GetElectricalData(); + int n = static_cast(generator->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] += GetSyncMachineAdmittance(generator); + } + } +} + +bool Electromechanical::EventTrigger(double eventTime, double currentTime) +{ + return (((eventTime - m_timeStep) < currentTime) && (eventTime >= currentTime)); +} + +std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* generator) +{ + auto data = generator->GetPUElectricalData(m_powerSystemBase); + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + + double xd = 0.0; + double xq = 0.0; + double ra = data.armResistance * k; + + switch(GetMachineModel(generator)) { + case SM_MODEL_1: + case SM_MODEL_2: + case SM_MODEL_3: { + xd = data.transXd * k; + xq = data.transXq * k; + break; + } + case SM_MODEL_4: + case SM_MODEL_5: { + xd = data.subXd * k; + xq = data.subXq * k; + break; + } + } + double xdq = 0.5 * (xd + xq); + return std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); +} -- cgit From 2c60c7c13ecc8e45e072346193a820f697fb6fdf Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Tue, 23 May 2017 20:18:04 -0300 Subject: Switching event bug fixed --- Project/Electromechanical.cpp | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 191a1da..895a6f3 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -89,10 +89,11 @@ void Electromechanical::SetEvent(double currentTime) // Insert machine (only disconnected machines) if(swData.swType[i] == SW_INSERT && !generator->IsOnline() && generator->GetParentList().size() == 1) { - generator->SetOnline(true); - auto data = generator->GetPUElectricalData(m_powerSystemBase); - int n = static_cast(generator->GetParentList()[0])->GetElectricalData().number; - m_yBus[n][n] += GetSyncMachineAdmittance(generator); + if(generator->SetOnline(true)) { + auto data = generator->GetPUElectricalData(m_powerSystemBase); + int n = static_cast(generator->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] += GetSyncMachineAdmittance(generator); + } } } } @@ -116,12 +117,13 @@ void Electromechanical::SetEvent(double currentTime) // Insert load (only disconnected load) if(swData.swType[i] == SW_INSERT && !load->IsOnline() && load->GetParentList().size() == 1) { - load->SetOnline(true); - auto data = load->GetPUElectricalData(m_powerSystemBase); - Bus* parentBus = static_cast(load->GetParentList()[0]); - int n = parentBus->GetElectricalData().number; - std::complex v = parentBus->GetElectricalData().voltage; - m_yBus[n][n] += std::complex(data.activePower, -data.reactivePower) / (v * v); + if(load->SetOnline(true)) { + auto data = load->GetPUElectricalData(m_powerSystemBase); + Bus* parentBus = static_cast(load->GetParentList()[0]); + int n = parentBus->GetElectricalData().number; + std::complex v = parentBus->GetElectricalData().voltage; + m_yBus[n][n] += std::complex(data.activePower, -data.reactivePower) / (v * v); + } } } } @@ -153,8 +155,7 @@ void Electromechanical::SetEvent(double currentTime) // Insert line (only disconnected lines) if(swData.swType[i] == SW_INSERT && !line->IsOnline() && line->GetParentList().size() == 2) { - if(line->GetParentList()[0] && line->GetParentList()[1]) { - line->SetOnline(true); + if(line->SetOnline(true)) { auto data = line->GetElectricalData(); int n1 = static_cast(line->GetParentList()[0])->GetElectricalData().number; @@ -212,8 +213,7 @@ void Electromechanical::SetEvent(double currentTime) // Insert transformer (only disconnected transformers) if(swData.swType[i] == SW_INSERT && !transformer->IsOnline() && transformer->GetParentList().size() == 2) { - if(transformer->GetParentList()[0] && transformer->GetParentList()[1]) { - transformer->SetOnline(true); + if(transformer->SetOnline(true)) { auto data = transformer->GetElectricalData(); int n1 = static_cast(transformer->GetParentList()[0])->GetElectricalData().number; @@ -260,8 +260,7 @@ void Electromechanical::SetEvent(double currentTime) // Insert capacitor (only disconnected capacitors) if(swData.swType[i] == SW_INSERT && !capacitor->IsOnline() && capacitor->GetParentList().size() == 1) { - if(capacitor->GetParentList()[0]) { - capacitor->SetOnline(true); + if(capacitor->SetOnline(true)) { auto data = capacitor->GetPUElectricalData(m_powerSystemBase); int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; m_yBus[n][n] -= std::complex(0.0, data.reactivePower); @@ -287,8 +286,7 @@ void Electromechanical::SetEvent(double currentTime) // Insert inductor (only disconnected inductors) if(swData.swType[i] == SW_INSERT && !inductor->IsOnline() && inductor->GetParentList().size() == 1) { - if(inductor->GetParentList()[0]) { - inductor->SetOnline(true); + if(inductor->SetOnline(true)) { auto data = inductor->GetPUElectricalData(m_powerSystemBase); int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; m_yBus[n][n] -= std::complex(0.0, -data.reactivePower); -- cgit From c91e58bb903adeff1e8c0fff1868e80783010e58 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Wed, 24 May 2017 17:50:36 -0300 Subject: Dynamic elements initialization implemented --- Project/Electromechanical.cpp | 164 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 142 insertions(+), 22 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 895a6f3..d2d9ae9 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -1,12 +1,37 @@ #include "Electromechanical.h" +#include "ControlElementSolver.h" -Electromechanical::Electromechanical(std::vector elementList) +Electromechanical::Electromechanical(wxWindow* parent, std::vector elementList) { + m_parent = parent; GetElementsFromList(elementList); SetEventTimeList(); } Electromechanical::~Electromechanical() {} +bool Electromechanical::RunStabilityCalculation() +{ + // Calculate the admittance matrix with the synchronous machines. + if(!GetYBus(m_yBus, m_powerSystemBase, POSITIVE_SEQ, false, true)) { + m_errorMsg = _("It was not possible to build the admittance matrix."); + return false; + } + InsertSyncMachinesOnYBus(); + + if(!InitializeDynamicElements()) return false; + + // test + double simTime = 10.0; + double currentTime = 0.0; + while(currentTime <= simTime) { + if(HasEvent(currentTime)) { + SetEvent(currentTime); + } + currentTime += m_timeStep; + } + return true; +} + void Electromechanical::SetEventTimeList() { // Fault @@ -297,27 +322,6 @@ void Electromechanical::SetEvent(double currentTime) } } -bool Electromechanical::RunStabilityCalculation() -{ - // Calculate the admittance matrix with the synchronous machines. - if(!GetYBus(m_yBus, m_powerSystemBase, POSITIVE_SEQ, false, true)) { - m_errorMsg = _("It was not possible to build the admittance matrix."); - return false; - } - InsertSyncMachinesOnYBus(); - - // test - double simTime = 10.0; - double currentTime = 0.0; - while(currentTime <= simTime) { - if(HasEvent(currentTime)) { - SetEvent(currentTime); - } - currentTime += m_timeStep; - } - return true; -} - void Electromechanical::InsertSyncMachinesOnYBus() { for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { @@ -366,3 +370,119 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* double xdq = 0.5 * (xd + xq); return std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); } + +bool Electromechanical::InitializeDynamicElements() +{ + // Synchronous generators + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + if(syncGenerator->IsOnline()) { + auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; + + std::complex conjS(data.activePower, -data.reactivePower); + std::complex conjV = std::conj(data.terminalVoltage); + std::complex ia = conjS / conjV; + + double xd = data.syncXd * k; + double xq = data.syncXq * k; + double ra = data.armResistance * k; + + // Initialize state variables + std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; + data.delta = std::arg(eq0); + + double teta0 = std::arg(data.terminalVoltage); + double vd0, vq0; + ABCtoDQ0(data.terminalVoltage, data.delta - teta0, vd0, vq0); + + double fi0 = std::arg(ia); + double id0, iq0; + ABCtoDQ0(ia, data.delta - fi0, id0, iq0); + + data.initialFieldVoltage = std::abs(eq0) - (xd - xq) * id0; + data.pm = std::real((data.terminalVoltage * std::conj(ia)) + (std::abs(ia) * std::abs(ia) * ra)); + data.speed = 2.0 * M_PI * 60.0; + + switch(GetMachineModel(syncGenerator)) { + case SM_MODEL_1: { + data.tranEq = 0.0; + data.tranEd = 0.0; + data.subEq = 0.0; + data.subEd = 0.0; + } break; + case SM_MODEL_2: { + double tranXd = data.transXd * k; + + data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; + data.tranEd = 0.0; + data.subEd = 0.0; + data.subEq = 0.0; + } break; + case SM_MODEL_3: { + double tranXd = data.transXd * k; + double tranXq = data.transXq * k; + + data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; + data.tranEd = -(xq - tranXq) * iq0; + data.subEd = 0.0; + data.subEq = 0.0; + } break; + case SM_MODEL_4: { + double tranXd = data.transXd * k; + double subXd = data.subXd * k; + double subXq = data.subXq * k; + + data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; + data.tranEd = 0.0; + data.subEq = data.tranEq + (tranXd - subXd) * id0; + data.subEd = -(xq - subXq) * iq0; + } break; + case SM_MODEL_5: { + double tranXd = data.transXd * k; + double tranXq = data.transXq * k; + double subXd = data.subXd * k; + double subXq = data.subXq * k; + + data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; + data.tranEd = -(xq - tranXq) * iq0; + data.subEq = data.tranEq + (tranXd - subXd) * id0; + data.subEd = data.tranEd + (tranXq - subXq) * iq0; + } break; + default: { + break; + } + } + + // Initialize controllers + if(data.useAVR) { + if(data.avrSolver) delete data.avrSolver; + data.avrSolver = new ControlElementSolver(data.avr, m_timeStep, 1e-3, false, + std::abs(data.terminalVoltage), m_parent); + if(!data.avrSolver->IsOK()) { + m_errorMsg = _("Error on initializate the AVR of \"") + data.name + _("\"."); + syncGenerator->SetElectricalData(data); + return false; + } + } + if(data.useSpeedGovernor) { + if(data.speedGovSolver) delete data.speedGovSolver; + data.speedGovSolver = + new ControlElementSolver(data.speedGov, m_timeStep, 1e-3, false, data.speed, m_parent); + if(!data.speedGovSolver->IsOK()) { + m_errorMsg = _("Error on initializate the speed governor of \"") + data.name + _("\"."); + syncGenerator->SetElectricalData(data); + return false; + } + } + + syncGenerator->SetElectricalData(data); + } + } + return true; +} -- cgit From 993288a099a4ba08c40cfb5ff79620257193131c Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 25 May 2017 18:50:28 -0300 Subject: Sync machines solver almost finished --- Project/Electromechanical.cpp | 252 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 251 insertions(+), 1 deletion(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index d2d9ae9..eceb197 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -17,6 +17,22 @@ bool Electromechanical::RunStabilityCalculation() return false; } InsertSyncMachinesOnYBus(); + GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); + + // Get buses voltages. + m_vBus.clear(); + m_vBus.resize(m_busList.size()); + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + auto data = bus->GetElectricalData(); + m_vBus[data.number] = data.voltage; + } + + // Calculate injected currents + m_iBus = ComplexMatrixTimesVector(m_yBus, m_vBus); + for(unsigned int i = 0; i < m_iBus.size(); ++i) { + if(std::abs(m_iBus[i]) < 1e-5) m_iBus[i] = std::complex(0.0, 0.0); + } if(!InitializeDynamicElements()) return false; @@ -26,7 +42,15 @@ bool Electromechanical::RunStabilityCalculation() while(currentTime <= simTime) { if(HasEvent(currentTime)) { SetEvent(currentTime); + GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); + } + + // Solve synchronous machines. + if(!SolveSynchronousMachines()){ + wxMessageBox(wxString::Format("%f", currentTime)); + return false; } + currentTime += m_timeStep; } return true; @@ -406,8 +430,9 @@ bool Electromechanical::InitializeDynamicElements() ABCtoDQ0(ia, data.delta - fi0, id0, iq0); data.initialFieldVoltage = std::abs(eq0) - (xd - xq) * id0; + data.fieldVoltage = data.initialFieldVoltage; data.pm = std::real((data.terminalVoltage * std::conj(ia)) + (std::abs(ia) * std::abs(ia) * ra)); - data.speed = 2.0 * M_PI * 60.0; + data.speed = 2.0 * M_PI * m_systemFreq; switch(GetMachineModel(syncGenerator)) { case SM_MODEL_1: { @@ -486,3 +511,228 @@ bool Electromechanical::InitializeDynamicElements() } return true; } + +void Electromechanical::CalculateMachinesCurrents() +{ + // Reset injected currents vector + for(unsigned int i = 0; i < m_iBus.size(); ++i) m_iBus[i] = std::complex(0.0, 0.0); + + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); + if(syncGenerator->IsOnline()) { + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + + double xd = 0.0; + double xq = 0.0; + double ra = data.armResistance * k; + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + std::complex e = std::complex(0.0, 0.0); + std::complex v = m_vBus[n]; + std::complex iInj = m_iBus[n]; + double xdq = 0.0; + + switch(GetMachineModel(syncGenerator)) { + case SM_MODEL_1: + case SM_MODEL_2: + case SM_MODEL_3: { + DQ0toABC(data.tranEd, data.tranEq, data.delta, e); + xd = data.transXd * k; + xq = data.transXq * k; + } break; + case SM_MODEL_4: + case SM_MODEL_5: { + DQ0toABC(data.subEd, data.subEq, data.delta, e); + xd = data.subXd * k; + xq = data.subXq * k; + } break; + } + xdq = 0.5 * (xd + xq); + + std::complex y0 = std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); + std::complex iUnadj = y0 * e; + + double dVR = std::real(e) - std::real(v); + double dVI = std::imag(e) - std::imag(v); + + double iAdjR = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * + (-std::sin(2.0 * data.delta) * dVR + std::cos(2.0 * data.delta) * dVI); + double iAdjI = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * + (std::cos(2.0 * data.delta) * dVR + std::sin(2.0 * data.delta) * dVI); + + iInj = iUnadj + std::complex(iAdjR, iAdjI); + m_iBus[n] += iInj; + + std::complex iMachine = iInj - y0 * v; + + data.electricalPower = v * std::conj(iMachine); + } else { + data.electricalPower = std::complex(0.0, 0.0); + } + + syncGenerator->SetElectricalData(data); + } +} + +void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq) +{ + auto data = syncGenerator->GetElectricalData(); + + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + double w0 = 2.0f * M_PI * m_systemFreq; + + // Speed + data.icSpeed.m = m_timeStep / ((4.0f * data.inertia / w0) / k + m_timeStep * data.damping * k); + data.icSpeed.c = (1.0f - 2.0f * data.icSpeed.m * data.damping * k) * data.speed + + data.icSpeed.m * (data.pm - data.pe + 2.0f * w0 * data.damping * k); + + // Delta + data.icDelta.m = 0.5f * m_timeStep; + data.icDelta.c = data.delta + data.icDelta.m * (data.speed - 2.0f * w0); + + // Eq' + data.icTranEq.m = m_timeStep / (2.0f * data.transTd0 + m_timeStep); + data.icTranEq.c = (1.0f - 2.0 * data.icTranEq.m) * data.tranEq + + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); + + // Ed' + data.icTranEd.m = m_timeStep / (2.0f * data.transTq0 + m_timeStep); + data.icTranEd.c = + (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; + + // Eq'' + data.icSubEq.m = m_timeStep / (2.0f * data.subTd0 + m_timeStep); + data.icSubEq.c = (1.0f - 2.0f * data.icSubEq.m) * data.subEq + + data.icSubEq.m * (data.tranEq + (data.transXd * k - data.subXd * k) * id); + + // Ed'' + data.icSubEd.m = m_timeStep / (2.0f * data.subTq0 + m_timeStep); + data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd + + data.icSubEd.m * (data.tranEd - (data.transXq * k - data.subXq * k) * iq); + + syncGenerator->SetElectricalData(data); +} + +bool Electromechanical::SolveSynchronousMachines() +{ + double w0 = 2.0 * M_PI * m_systemFreq; + + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + if(syncGenerator->IsOnline()) { + auto data = syncGenerator->GetElectricalData(); + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + double id, iq; + + std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); + + ABCtoDQ0(iMachine, data.delta, id, iq); + + // Calculate integration constants. + CalculateIntegrationConstants(syncGenerator, id, iq); + } + } + + double error = 1.0; + int iterations = 0; + while(error > m_tolerance) { + error = 0.0; + + // Calculate the injected currents. + CalculateMachinesCurrents(); + + // Calculate the buses voltages. + m_vBus = LUEvaluate(m_yBusU, m_yBusL, m_iBus); + + // Solve machine equations. + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + if(syncGenerator->IsOnline()) { + auto data = syncGenerator->GetElectricalData(); + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + + data.terminalVoltage = m_vBus[n]; + + // Mechanical differential equations. + double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe); + error = std::max(error, std::abs(data.speed - w) / w0); + + double delta = data.icDelta.c + data.icDelta.m * w; + error = std::max(error, std::abs(data.delta - delta)); + + data.speed = w; + data.delta = delta; + + // Electric power. + double id, iq, vd, vq; + std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); + + ABCtoDQ0(iMachine, data.delta, id, iq); + ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); + + double pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; + // data.pe = (2 * pe - data.pe); // Extrapolating Pe. + data.pe = pe; // Don't extrapolating Pe. + + // Solve controllers. + if(data.useAVR && data.avrSolver) { + data.avrSolver->SolveNextStep(std::abs(data.terminalVoltage)); + data.fieldVoltage = data.initialFieldVoltage + data.avrSolver->GetLastSolution(); + } + + if(data.useSpeedGovernor && data.speedGovSolver) { + data.speedGovSolver->SolveNextStep(data.speed); + data.pm = data.speedGovSolver->GetLastSolution(); + } + + // Electrical differential equations + switch(GetMachineModel(syncGenerator)) { + case SM_MODEL_1: { + // There is no differential equations. + } break; + case SM_MODEL_2: { + } break; + case SM_MODEL_3: { + double tranEq = + data.icTranEq.c + + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); + error = std::max(error, std::abs(data.tranEq - tranEq)); + + double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; + error = std::max(error, std::abs(data.tranEd - tranEd)); + + data.tranEq = tranEq; + data.tranEd = tranEd; + } break; + case SM_MODEL_4: { + } break; + case SM_MODEL_5: { + } break; + } + syncGenerator->SetElectricalData(data); + } + } + + ++iterations; + if(iterations > m_maxIterations) { + m_errorMsg = _("Impossible to solve the synchronous generators.\nCheck the system parameters and/or " + "decrease the time step."); + return false; + } + } + + return true; +} -- cgit From 734476bea0a2b85b131f5a4d2c9e2b219af7be41 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Fri, 26 May 2017 02:21:15 -0300 Subject: Sync generator plot implemented Electromechanical calc not working --- Project/Electromechanical.cpp | 105 +++++++++++++++++++++++++++++++++++------- 1 file changed, 88 insertions(+), 17 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index eceb197..09d7791 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -38,20 +38,26 @@ bool Electromechanical::RunStabilityCalculation() // test double simTime = 10.0; + double printTime = 0.01; double currentTime = 0.0; + double currentPrintTime = 0.0; while(currentTime <= simTime) { if(HasEvent(currentTime)) { SetEvent(currentTime); GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); } - // Solve synchronous machines. - if(!SolveSynchronousMachines()){ - wxMessageBox(wxString::Format("%f", currentTime)); - return false; + bool saveData = false; + if(currentPrintTime >= printTime) { + m_timeVector.push_back(currentTime); + currentPrintTime = 0.0; + saveData = true; } + if(!SolveSynchronousMachines(saveData)) return false; + currentTime += m_timeStep; + currentPrintTime += m_timeStep; } return true; } @@ -368,7 +374,7 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* auto data = generator->GetPUElectricalData(m_powerSystemBase); double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } @@ -400,11 +406,11 @@ bool Electromechanical::InitializeDynamicElements() // Synchronous generators for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); if(syncGenerator->IsOnline()) { - auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; @@ -434,6 +440,9 @@ bool Electromechanical::InitializeDynamicElements() data.pm = std::real((data.terminalVoltage * std::conj(ia)) + (std::abs(ia) * std::abs(ia) * ra)); data.speed = 2.0 * M_PI * m_systemFreq; + //data.pe = id0 * vd0 + iq0 * vq0 + (id0 * id0 + iq0 * iq0) * data.armResistance * k; + //data.electricalPower = std::complex(data.activePower, data.reactivePower); + switch(GetMachineModel(syncGenerator)) { case SM_MODEL_1: { data.tranEq = 0.0; @@ -505,9 +514,18 @@ bool Electromechanical::InitializeDynamicElements() return false; } } - - syncGenerator->SetElectricalData(data); + } else { + // Initialize open circuit machine. } + // Reset plot data + data.terminalVoltageVector.clear(); + data.electricalPowerVector.clear(); + data.mechanicalPowerVector.clear(); + data.freqVector.clear(); + data.fieldVoltageVector.clear(); + data.deltaVector.clear(); + + syncGenerator->SetElectricalData(data); } return true; } @@ -523,7 +541,7 @@ void Electromechanical::CalculateMachinesCurrents() if(syncGenerator->IsOnline()) { double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } @@ -584,7 +602,7 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } double w0 = 2.0f * M_PI * m_systemFreq; @@ -617,18 +635,18 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat data.icSubEd.m = m_timeStep / (2.0f * data.subTq0 + m_timeStep); data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd + data.icSubEd.m * (data.tranEd - (data.transXq * k - data.subXq * k) * iq); - + syncGenerator->SetElectricalData(data); } -bool Electromechanical::SolveSynchronousMachines() +bool Electromechanical::SolveSynchronousMachines(bool saveValues) { double w0 = 2.0 * M_PI * m_systemFreq; for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; if(syncGenerator->IsOnline()) { - auto data = syncGenerator->GetElectricalData(); + auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; double id, iq; @@ -638,6 +656,16 @@ bool Electromechanical::SolveSynchronousMachines() // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq); + + if(saveValues && data.plotSyncMachine) { + data.terminalVoltageVector.push_back(data.terminalVoltage); + data.electricalPowerVector.push_back(data.pe); + data.mechanicalPowerVector.push_back(data.pm); + data.freqVector.push_back(data.speed / (2.0f * M_PI)); + data.fieldVoltageVector.push_back(data.fieldVoltage); + data.deltaVector.push_back(data.delta); + } + syncGenerator->SetElectricalData(data); } } @@ -655,11 +683,12 @@ bool Electromechanical::SolveSynchronousMachines() // Solve machine equations. for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetElectricalData(); + if(syncGenerator->IsOnline()) { - auto data = syncGenerator->GetElectricalData(); double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = data.nominalPower * std::pow(1000.0f, data.nominalPowerUnit); + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; @@ -669,6 +698,7 @@ bool Electromechanical::SolveSynchronousMachines() // Mechanical differential equations. double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe); error = std::max(error, std::abs(data.speed - w) / w0); + //wxMessageBox(wxString::Format("%f", std::abs(w))); double delta = data.icDelta.c + data.icDelta.m * w; error = std::max(error, std::abs(data.delta - delta)); @@ -722,7 +752,8 @@ bool Electromechanical::SolveSynchronousMachines() case SM_MODEL_5: { } break; } - syncGenerator->SetElectricalData(data); + } else { + // Set values to open circuit machine. } } @@ -736,3 +767,43 @@ bool Electromechanical::SolveSynchronousMachines() return true; } + +double Electromechanical::GetPowerValue(double value, ElectricalUnit unit) +{ + switch(unit) { + case UNIT_PU: { + return value; + } break; + case UNIT_VA: { + return value; + } break; + case UNIT_kVA: { + return value * 1e3; + } break; + case UNIT_MVA: { + return value * 1e6; + } break; + case UNIT_W: { + return value; + } break; + case UNIT_kW: { + return value * 1e3; + } break; + case UNIT_MW: { + return value * 1e6; + } break; + case UNIT_VAr: { + return value; + } break; + case UNIT_kVAr: { + return value * 1e3; + } break; + case UNIT_MVAr: { + return value * 1e6; + } break; + default: { + return 0.0; + } break; + } + return 0.0; +} -- 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/Electromechanical.cpp | 86 +++++++++++++++++++++++++++++-------------- 1 file changed, 58 insertions(+), 28 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 09d7791..4988382 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -20,19 +20,25 @@ bool Electromechanical::RunStabilityCalculation() GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); // Get buses voltages. + wxString str = ""; m_vBus.clear(); m_vBus.resize(m_busList.size()); for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { Bus* bus = *it; auto data = bus->GetElectricalData(); m_vBus[data.number] = data.voltage; + str += wxString::Format("%f /_ %f\n", std::abs(m_vBus[data.number]), std::arg(m_vBus[data.number]) * 180.0f / M_PI); } + wxMessageBox(str); // Calculate injected currents + str = ""; m_iBus = ComplexMatrixTimesVector(m_yBus, m_vBus); for(unsigned int i = 0; i < m_iBus.size(); ++i) { if(std::abs(m_iBus[i]) < 1e-5) m_iBus[i] = std::complex(0.0, 0.0); + str += wxString::Format("%f /_ %f\n", std::abs(m_iBus[i]), std::arg(m_iBus[i]) * 180.0f / M_PI); } + wxMessageBox(str); if(!InitializeDynamicElements()) return false; @@ -47,14 +53,16 @@ bool Electromechanical::RunStabilityCalculation() GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); } - bool saveData = false; if(currentPrintTime >= printTime) { m_timeVector.push_back(currentTime); + SaveData(); currentPrintTime = 0.0; - saveData = true; } - if(!SolveSynchronousMachines(saveData)) return false; + if(!SolveSynchronousMachines()) { + wxMessageBox(wxString::Format("%f", currentTime)); + return false; + } currentTime += m_timeStep; currentPrintTime += m_timeStep; @@ -325,7 +333,7 @@ void Electromechanical::SetEvent(double currentTime) } } - // Capacitor switching + // Inductor switching for(auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) { Inductor* inductor = *it; auto swData = inductor->GetSwitchingData(); @@ -398,7 +406,7 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* } } double xdq = 0.5 * (xd + xq); - return std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); + return (std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0)); } bool Electromechanical::InitializeDynamicElements() @@ -427,21 +435,28 @@ bool Electromechanical::InitializeDynamicElements() std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; data.delta = std::arg(eq0); - double teta0 = std::arg(data.terminalVoltage); - double vd0, vq0; - ABCtoDQ0(data.terminalVoltage, data.delta - teta0, vd0, vq0); + //double teta0 = std::arg(data.terminalVoltage); + //double vd0, vq0; + // ABCtoDQ0(data.terminalVoltage, data.delta - teta0, vd0, vq0); + //vq0 = std::abs(data.terminalVoltage) * cos(data.delta - teta0); + //vd0 = -std::abs(data.terminalVoltage) * sin(data.delta - teta0); double fi0 = std::arg(ia); double id0, iq0; - ABCtoDQ0(ia, data.delta - fi0, id0, iq0); + //ABCtoDQ0(ia, data.delta - fi0, id0, iq0); + iq0 = std::abs(ia) * cos(data.delta - fi0); + id0 = -std::abs(ia) * sin(data.delta - fi0); data.initialFieldVoltage = std::abs(eq0) - (xd - xq) * id0; data.fieldVoltage = data.initialFieldVoltage; data.pm = std::real((data.terminalVoltage * std::conj(ia)) + (std::abs(ia) * std::abs(ia) * ra)); data.speed = 2.0 * M_PI * m_systemFreq; - //data.pe = id0 * vd0 + iq0 * vq0 + (id0 * id0 + iq0 * iq0) * data.armResistance * k; - //data.electricalPower = std::complex(data.activePower, data.reactivePower); + data.pe = data.pm; + data.electricalPower = std::complex(data.activePower, data.reactivePower); + + // data.pe = id0 * vd0 + iq0 * vq0 + (id0 * id0 + iq0 * iq0) * data.armResistance * k; + // data.electricalPower = std::complex(data.activePower, data.reactivePower); switch(GetMachineModel(syncGenerator)) { case SM_MODEL_1: { @@ -577,9 +592,9 @@ void Electromechanical::CalculateMachinesCurrents() double dVR = std::real(e) - std::real(v); double dVI = std::imag(e) - std::imag(v); - double iAdjR = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * + double iAdjR = ((xd - xq) / (ra * ra + xd * xq)) * (-std::sin(2.0 * data.delta) * dVR + std::cos(2.0 * data.delta) * dVI); - double iAdjI = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * + double iAdjI = ((xd - xq) / (ra * ra + xd * xq)) * (std::cos(2.0 * data.delta) * dVR + std::sin(2.0 * data.delta) * dVI); iInj = iUnadj + std::complex(iAdjR, iAdjI); @@ -598,7 +613,7 @@ void Electromechanical::CalculateMachinesCurrents() void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq) { - auto data = syncGenerator->GetElectricalData(); + auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); double k = 1.0; // Power base change factor. if(data.useMachineBase) { @@ -639,14 +654,15 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat syncGenerator->SetElectricalData(data); } -bool Electromechanical::SolveSynchronousMachines(bool saveValues) +bool Electromechanical::SolveSynchronousMachines() { double w0 = 2.0 * M_PI * m_systemFreq; for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetElectricalData(); + if(syncGenerator->IsOnline()) { - auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; double id, iq; @@ -656,16 +672,6 @@ bool Electromechanical::SolveSynchronousMachines(bool saveValues) // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq); - - if(saveValues && data.plotSyncMachine) { - data.terminalVoltageVector.push_back(data.terminalVoltage); - data.electricalPowerVector.push_back(data.pe); - data.mechanicalPowerVector.push_back(data.pm); - data.freqVector.push_back(data.speed / (2.0f * M_PI)); - data.fieldVoltageVector.push_back(data.fieldVoltage); - data.deltaVector.push_back(data.delta); - } - syncGenerator->SetElectricalData(data); } } @@ -675,10 +681,13 @@ bool Electromechanical::SolveSynchronousMachines(bool saveValues) error = 0.0; // Calculate the injected currents. + //wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); CalculateMachinesCurrents(); - + //wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); + // Calculate the buses voltages. m_vBus = LUEvaluate(m_yBusU, m_yBusL, m_iBus); + // Solve machine equations. for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { @@ -698,7 +707,6 @@ bool Electromechanical::SolveSynchronousMachines(bool saveValues) // Mechanical differential equations. double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe); error = std::max(error, std::abs(data.speed - w) / w0); - //wxMessageBox(wxString::Format("%f", std::abs(w))); double delta = data.icDelta.c + data.icDelta.m * w; error = std::max(error, std::abs(data.delta - delta)); @@ -742,10 +750,13 @@ bool Electromechanical::SolveSynchronousMachines(bool saveValues) error = std::max(error, std::abs(data.tranEq - tranEq)); double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; + // double tranEd = data.icTranEq.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * + // id; error = std::max(error, std::abs(data.tranEd - tranEd)); data.tranEq = tranEq; data.tranEd = tranEd; + } break; case SM_MODEL_4: { } break; @@ -755,6 +766,8 @@ bool Electromechanical::SolveSynchronousMachines(bool saveValues) } else { // Set values to open circuit machine. } + + syncGenerator->SetElectricalData(data); } ++iterations; @@ -807,3 +820,20 @@ double Electromechanical::GetPowerValue(double value, ElectricalUnit unit) } return 0.0; } + +void Electromechanical::SaveData() +{ + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetElectricalData(); + if(data.plotSyncMachine) { + data.terminalVoltageVector.push_back(data.terminalVoltage); + data.electricalPowerVector.push_back(data.pe); + data.mechanicalPowerVector.push_back(data.pm); + data.freqVector.push_back(data.speed / (2.0f * M_PI)); + data.fieldVoltageVector.push_back(data.fieldVoltage); + data.deltaVector.push_back(data.delta); + syncGenerator->SetElectricalData(data); + } + } +} -- cgit From 5df0b53a01acfa3e41b0e106426fa874f6af81f7 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Sat, 27 May 2017 01:32:26 -0300 Subject: Some minor changes --- Project/Electromechanical.cpp | 54 +++++++++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 23 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 4988382..4c8d3c1 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -27,7 +27,8 @@ bool Electromechanical::RunStabilityCalculation() Bus* bus = *it; auto data = bus->GetElectricalData(); m_vBus[data.number] = data.voltage; - str += wxString::Format("%f /_ %f\n", std::abs(m_vBus[data.number]), std::arg(m_vBus[data.number]) * 180.0f / M_PI); + str += wxString::Format("%f /_ %f\n", std::abs(m_vBus[data.number]), + std::arg(m_vBus[data.number]) * 180.0f / M_PI); } wxMessageBox(str); @@ -435,15 +436,15 @@ bool Electromechanical::InitializeDynamicElements() std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; data.delta = std::arg(eq0); - //double teta0 = std::arg(data.terminalVoltage); - //double vd0, vq0; + // double teta0 = std::arg(data.terminalVoltage); + // double vd0, vq0; // ABCtoDQ0(data.terminalVoltage, data.delta - teta0, vd0, vq0); - //vq0 = std::abs(data.terminalVoltage) * cos(data.delta - teta0); - //vd0 = -std::abs(data.terminalVoltage) * sin(data.delta - teta0); + // vq0 = std::abs(data.terminalVoltage) * cos(data.delta - teta0); + // vd0 = -std::abs(data.terminalVoltage) * sin(data.delta - teta0); double fi0 = std::arg(ia); double id0, iq0; - //ABCtoDQ0(ia, data.delta - fi0, id0, iq0); + // ABCtoDQ0(ia, data.delta - fi0, id0, iq0); iq0 = std::abs(ia) * cos(data.delta - fi0); id0 = -std::abs(ia) * sin(data.delta - fi0); @@ -588,6 +589,7 @@ void Electromechanical::CalculateMachinesCurrents() std::complex y0 = std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); std::complex iUnadj = y0 * e; + wxMessageBox(wxString::Format("%f / %f\n", std::real(iUnadj), std::imag(iUnadj))); double dVR = std::real(e) - std::real(v); double dVI = std::imag(e) - std::imag(v); @@ -602,7 +604,9 @@ void Electromechanical::CalculateMachinesCurrents() std::complex iMachine = iInj - y0 * v; + //wxMessageBox(wxString::Format("%f /_ %f\n", std::real(data.electricalPower), std::imag(data.electricalPower))); data.electricalPower = v * std::conj(iMachine); + //wxMessageBox(wxString::Format("%f /_ %f\n", std::real(data.electricalPower), std::imag(data.electricalPower))); } else { data.electricalPower = std::complex(0.0, 0.0); } @@ -681,13 +685,12 @@ bool Electromechanical::SolveSynchronousMachines() error = 0.0; // Calculate the injected currents. - //wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); + // wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); CalculateMachinesCurrents(); - //wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); - + // wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); + // Calculate the buses voltages. m_vBus = LUEvaluate(m_yBusU, m_yBusL, m_iBus); - // Solve machine equations. for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { @@ -722,19 +725,8 @@ bool Electromechanical::SolveSynchronousMachines() ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); double pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; - // data.pe = (2 * pe - data.pe); // Extrapolating Pe. - data.pe = pe; // Don't extrapolating Pe. - - // Solve controllers. - if(data.useAVR && data.avrSolver) { - data.avrSolver->SolveNextStep(std::abs(data.terminalVoltage)); - data.fieldVoltage = data.initialFieldVoltage + data.avrSolver->GetLastSolution(); - } - - if(data.useSpeedGovernor && data.speedGovSolver) { - data.speedGovSolver->SolveNextStep(data.speed); - data.pm = data.speedGovSolver->GetLastSolution(); - } + //data.pe = (2 * pe - data.pe); // Extrapolating Pe. + data.pe = pe; // Don't extrapolating Pe // Electrical differential equations switch(GetMachineModel(syncGenerator)) { @@ -778,6 +770,22 @@ bool Electromechanical::SolveSynchronousMachines() } } + // Solve controllers. + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetElectricalData(); + if(data.useAVR && data.avrSolver) { + data.avrSolver->SolveNextStep(std::abs(data.terminalVoltage)); + data.fieldVoltage = data.initialFieldVoltage + data.avrSolver->GetLastSolution(); + } + + if(data.useSpeedGovernor && data.speedGovSolver) { + data.speedGovSolver->SolveNextStep(data.speed); + data.pm = data.speedGovSolver->GetLastSolution(); + } + syncGenerator->SetElectricalData(data); + } + return true; } -- cgit From bee32c0b8432c0a849baa57a1897eab936894af9 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Sat, 27 May 2017 17:25:43 -0300 Subject: More bugfixes If xd' is different of xq' the calculations fail --- Project/Electromechanical.cpp | 43 +++++++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 20 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 4c8d3c1..654f5d5 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -42,7 +42,7 @@ bool Electromechanical::RunStabilityCalculation() wxMessageBox(str); if(!InitializeDynamicElements()) return false; - + // test double simTime = 10.0; double printTime = 0.01; @@ -397,14 +397,12 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* case SM_MODEL_3: { xd = data.transXd * k; xq = data.transXq * k; - break; - } + } break; case SM_MODEL_4: case SM_MODEL_5: { xd = data.subXd * k; xq = data.subXq * k; - break; - } + } break; } double xdq = 0.5 * (xd + xq); return (std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0)); @@ -436,17 +434,17 @@ bool Electromechanical::InitializeDynamicElements() std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; data.delta = std::arg(eq0); - // double teta0 = std::arg(data.terminalVoltage); - // double vd0, vq0; - // ABCtoDQ0(data.terminalVoltage, data.delta - teta0, vd0, vq0); - // vq0 = std::abs(data.terminalVoltage) * cos(data.delta - teta0); - // vd0 = -std::abs(data.terminalVoltage) * sin(data.delta - teta0); + double teta0 = std::arg(data.terminalVoltage); + double vd0, vq0; + ABCtoDQ0(data.terminalVoltage, data.delta - teta0, vd0, vq0); + vq0 = std::abs(data.terminalVoltage) * std::cos(data.delta - teta0); + vd0 = -std::abs(data.terminalVoltage) * std::sin(data.delta - teta0); double fi0 = std::arg(ia); double id0, iq0; // ABCtoDQ0(ia, data.delta - fi0, id0, iq0); - iq0 = std::abs(ia) * cos(data.delta - fi0); - id0 = -std::abs(ia) * sin(data.delta - fi0); + iq0 = std::abs(ia) * std::cos(data.delta - fi0); + id0 = -std::abs(ia) * std::sin(data.delta - fi0); data.initialFieldVoltage = std::abs(eq0) - (xd - xq) * id0; data.fieldVoltage = data.initialFieldVoltage; @@ -480,6 +478,13 @@ bool Electromechanical::InitializeDynamicElements() data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; data.tranEd = -(xq - tranXq) * iq0; + + // data.tranEq = vq0 + ra * iq0 - tranXd * id0; + // data.tranEd = vd0 + ra * id0 - tranXq * iq0; + + // data.tranEq = data.initialFieldVoltage; + // data.tranEd = + data.subEd = 0.0; data.subEq = 0.0; } break; @@ -504,9 +509,8 @@ bool Electromechanical::InitializeDynamicElements() data.subEq = data.tranEq + (tranXd - subXd) * id0; data.subEd = data.tranEd + (tranXq - subXq) * iq0; } break; - default: { + default: break; - } } // Initialize controllers @@ -589,7 +593,7 @@ void Electromechanical::CalculateMachinesCurrents() std::complex y0 = std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); std::complex iUnadj = y0 * e; - wxMessageBox(wxString::Format("%f / %f\n", std::real(iUnadj), std::imag(iUnadj))); + // wxMessageBox(wxString::Format("%f / %f\n", data.tranEd, data.tranEq)); double dVR = std::real(e) - std::real(v); double dVI = std::imag(e) - std::imag(v); @@ -603,10 +607,9 @@ void Electromechanical::CalculateMachinesCurrents() m_iBus[n] += iInj; std::complex iMachine = iInj - y0 * v; + //std::complex iMachine = iUnadj - y0 * v; - //wxMessageBox(wxString::Format("%f /_ %f\n", std::real(data.electricalPower), std::imag(data.electricalPower))); data.electricalPower = v * std::conj(iMachine); - //wxMessageBox(wxString::Format("%f /_ %f\n", std::real(data.electricalPower), std::imag(data.electricalPower))); } else { data.electricalPower = std::complex(0.0, 0.0); } @@ -725,7 +728,7 @@ bool Electromechanical::SolveSynchronousMachines() ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); double pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; - //data.pe = (2 * pe - data.pe); // Extrapolating Pe. + // data.pe = (2 * pe - data.pe); // Extrapolating Pe. data.pe = pe; // Don't extrapolating Pe // Electrical differential equations @@ -836,11 +839,11 @@ void Electromechanical::SaveData() auto data = syncGenerator->GetElectricalData(); if(data.plotSyncMachine) { data.terminalVoltageVector.push_back(data.terminalVoltage); - data.electricalPowerVector.push_back(data.pe); + data.electricalPowerVector.push_back(data.electricalPower); data.mechanicalPowerVector.push_back(data.pm); data.freqVector.push_back(data.speed / (2.0f * M_PI)); data.fieldVoltageVector.push_back(data.fieldVoltage); - data.deltaVector.push_back(data.delta); + data.deltaVector.push_back(wxRadToDeg(data.delta)); syncGenerator->SetElectricalData(data); } } -- cgit From e1a11643e0245676b04d6c9fce5eb35d68163121 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Sun, 28 May 2017 02:37:17 -0300 Subject: Type 3 machine working properly now --- Project/Electromechanical.cpp | 56 ++++++++++++++++--------------------------- 1 file changed, 21 insertions(+), 35 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 654f5d5..6cf5d9e 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -11,6 +11,9 @@ Electromechanical::Electromechanical(wxWindow* parent, std::vector ele Electromechanical::~Electromechanical() {} bool Electromechanical::RunStabilityCalculation() { + wxProgressDialog pbd(_("Running simulation"), _("Initializing..."), 100, m_parent, + wxPD_APP_MODAL | wxPD_AUTO_HIDE | wxPD_CAN_ABORT | wxPD_SMOOTH); + // Calculate the admittance matrix with the synchronous machines. if(!GetYBus(m_yBus, m_powerSystemBase, POSITIVE_SEQ, false, true)) { m_errorMsg = _("It was not possible to build the admittance matrix."); @@ -20,34 +23,29 @@ bool Electromechanical::RunStabilityCalculation() GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); // Get buses voltages. - wxString str = ""; m_vBus.clear(); m_vBus.resize(m_busList.size()); for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { Bus* bus = *it; auto data = bus->GetElectricalData(); m_vBus[data.number] = data.voltage; - str += wxString::Format("%f /_ %f\n", std::abs(m_vBus[data.number]), - std::arg(m_vBus[data.number]) * 180.0f / M_PI); } - wxMessageBox(str); // Calculate injected currents - str = ""; m_iBus = ComplexMatrixTimesVector(m_yBus, m_vBus); for(unsigned int i = 0; i < m_iBus.size(); ++i) { if(std::abs(m_iBus[i]) < 1e-5) m_iBus[i] = std::complex(0.0, 0.0); - str += wxString::Format("%f /_ %f\n", std::abs(m_iBus[i]), std::arg(m_iBus[i]) * 180.0f / M_PI); } - wxMessageBox(str); if(!InitializeDynamicElements()) return false; - + // test - double simTime = 10.0; + double simTime = 20.0; double printTime = 0.01; + double pbdTime = 0.01; double currentTime = 0.0; double currentPrintTime = 0.0; + double currentPbdTime = 0.0; while(currentTime <= simTime) { if(HasEvent(currentTime)) { SetEvent(currentTime); @@ -60,13 +58,22 @@ bool Electromechanical::RunStabilityCalculation() currentPrintTime = 0.0; } - if(!SolveSynchronousMachines()) { - wxMessageBox(wxString::Format("%f", currentTime)); - return false; + if(currentPbdTime > pbdTime) { + if(!pbd.Update((currentTime / simTime) * 100, wxString::Format("Time = %.2fs", currentTime))) { + pbd.Update(100); + currentTime = simTime; + + m_errorMsg = _("Simulation cancelled."); + return false; + } + currentPbdTime = 0.0; } + if(!SolveSynchronousMachines()) return false; + currentTime += m_timeStep; currentPrintTime += m_timeStep; + currentPbdTime += m_timeStep; } return true; } @@ -434,12 +441,6 @@ bool Electromechanical::InitializeDynamicElements() std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; data.delta = std::arg(eq0); - double teta0 = std::arg(data.terminalVoltage); - double vd0, vq0; - ABCtoDQ0(data.terminalVoltage, data.delta - teta0, vd0, vq0); - vq0 = std::abs(data.terminalVoltage) * std::cos(data.delta - teta0); - vd0 = -std::abs(data.terminalVoltage) * std::sin(data.delta - teta0); - double fi0 = std::arg(ia); double id0, iq0; // ABCtoDQ0(ia, data.delta - fi0, id0, iq0); @@ -454,9 +455,6 @@ bool Electromechanical::InitializeDynamicElements() data.pe = data.pm; data.electricalPower = std::complex(data.activePower, data.reactivePower); - // data.pe = id0 * vd0 + iq0 * vq0 + (id0 * id0 + iq0 * iq0) * data.armResistance * k; - // data.electricalPower = std::complex(data.activePower, data.reactivePower); - switch(GetMachineModel(syncGenerator)) { case SM_MODEL_1: { data.tranEq = 0.0; @@ -479,12 +477,6 @@ bool Electromechanical::InitializeDynamicElements() data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; data.tranEd = -(xq - tranXq) * iq0; - // data.tranEq = vq0 + ra * iq0 - tranXd * id0; - // data.tranEd = vd0 + ra * id0 - tranXq * iq0; - - // data.tranEq = data.initialFieldVoltage; - // data.tranEd = - data.subEd = 0.0; data.subEq = 0.0; } break; @@ -593,21 +585,19 @@ void Electromechanical::CalculateMachinesCurrents() std::complex y0 = std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); std::complex iUnadj = y0 * e; - // wxMessageBox(wxString::Format("%f / %f\n", data.tranEd, data.tranEq)); double dVR = std::real(e) - std::real(v); double dVI = std::imag(e) - std::imag(v); - double iAdjR = ((xd - xq) / (ra * ra + xd * xq)) * + double iAdjR = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * (-std::sin(2.0 * data.delta) * dVR + std::cos(2.0 * data.delta) * dVI); - double iAdjI = ((xd - xq) / (ra * ra + xd * xq)) * + double iAdjI = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * (std::cos(2.0 * data.delta) * dVR + std::sin(2.0 * data.delta) * dVI); iInj = iUnadj + std::complex(iAdjR, iAdjI); m_iBus[n] += iInj; std::complex iMachine = iInj - y0 * v; - //std::complex iMachine = iUnadj - y0 * v; data.electricalPower = v * std::conj(iMachine); } else { @@ -688,9 +678,7 @@ bool Electromechanical::SolveSynchronousMachines() error = 0.0; // Calculate the injected currents. - // wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); CalculateMachinesCurrents(); - // wxMessageBox(wxString::Format("%f %f", m_iBus[0].real(), m_iBus[0].imag())); // Calculate the buses voltages. m_vBus = LUEvaluate(m_yBusU, m_yBusL, m_iBus); @@ -745,8 +733,6 @@ bool Electromechanical::SolveSynchronousMachines() error = std::max(error, std::abs(data.tranEq - tranEq)); double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; - // double tranEd = data.icTranEq.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * - // id; error = std::max(error, std::abs(data.tranEd - tranEd)); data.tranEq = tranEq; -- cgit From 41c6ab0cac47046db7b7a3faf360c60944fd39b5 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Mon, 29 May 2017 00:40:46 -0300 Subject: Removing sync generator is now working, bus plot implemented --- Project/Electromechanical.cpp | 147 +++++++++++++++++++++++++++--------------- 1 file changed, 95 insertions(+), 52 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 6cf5d9e..a81912f 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -60,10 +60,8 @@ bool Electromechanical::RunStabilityCalculation() if(currentPbdTime > pbdTime) { if(!pbd.Update((currentTime / simTime) * 100, wxString::Format("Time = %.2fs", currentTime))) { + m_errorMsg = wxString::Format(_("Simulation cancelled at %.2fs."), currentTime); pbd.Update(100); - currentTime = simTime; - - m_errorMsg = _("Simulation cancelled."); return false; } currentPbdTime = 0.0; @@ -417,6 +415,13 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* bool Electromechanical::InitializeDynamicElements() { + // Buses + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + auto data = bus->GetElectricalData(); + data.stabVoltageVector.clear(); + bus->SetElectricalData(data); + } // Synchronous generators for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; @@ -670,8 +675,16 @@ bool Electromechanical::SolveSynchronousMachines() // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq); } + else { + CalculateIntegrationConstants(syncGenerator, 0.0f, 0.0f); + } } + m_wError = 0; + m_deltaError = 0; + m_transEdError = 0; + m_transEqError = 0; + double error = 1.0; int iterations = 0; while(error > m_tolerance) { @@ -688,64 +701,79 @@ bool Electromechanical::SolveSynchronousMachines() SyncGenerator* syncGenerator = *it; auto data = syncGenerator->GetElectricalData(); - if(syncGenerator->IsOnline()) { - double k = 1.0; // Power base change factor. - if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); - k = m_powerSystemBase / oldBase; - } - int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + if(syncGenerator->IsOnline()) { data.terminalVoltage = m_vBus[n]; + } - // Mechanical differential equations. - double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe); - error = std::max(error, std::abs(data.speed - w) / w0); + // Mechanical differential equations. + double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe); + error = std::max(error, std::abs(data.speed - w) / w0); - double delta = data.icDelta.c + data.icDelta.m * w; - error = std::max(error, std::abs(data.delta - delta)); + m_wError += std::abs(data.speed - w) / w0; - data.speed = w; - data.delta = delta; + double delta = data.icDelta.c + data.icDelta.m * w; + error = std::max(error, std::abs(data.delta - delta)); - // Electric power. - double id, iq, vd, vq; - std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); + data.speed = w; + data.delta = delta; + + m_deltaError += std::abs(data.delta - delta); + // Electric power. + double id, iq, vd, vq, pe; + ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); + + if(syncGenerator->IsOnline()) { + std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); ABCtoDQ0(iMachine, data.delta, id, iq); - ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); - - double pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; - // data.pe = (2 * pe - data.pe); // Extrapolating Pe. - data.pe = pe; // Don't extrapolating Pe - - // Electrical differential equations - switch(GetMachineModel(syncGenerator)) { - case SM_MODEL_1: { - // There is no differential equations. - } break; - case SM_MODEL_2: { - } break; - case SM_MODEL_3: { - double tranEq = - data.icTranEq.c + - data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); - error = std::max(error, std::abs(data.tranEq - tranEq)); - - double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; - error = std::max(error, std::abs(data.tranEd - tranEd)); - - data.tranEq = tranEq; - data.tranEd = tranEd; - - } break; - case SM_MODEL_4: { - } break; - case SM_MODEL_5: { - } break; - } + + pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; + // pe = (2 * pe - data.pe); // Extrapolating Pe. } else { - // Set values to open circuit machine. + pe = id = iq = 0.0f; + } + + data.pe = pe; + + // Electrical differential equations + switch(GetMachineModel(syncGenerator)) { + case SM_MODEL_1: { + // There is no differential equations. + } break; + case SM_MODEL_2: { + } break; + case SM_MODEL_3: { + double tranEq = data.icTranEq.c + + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); + error = std::max(error, std::abs(data.tranEq - tranEq)); + + m_transEqError += std::abs(data.tranEq - tranEq); + + double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; + error = std::max(error, std::abs(data.tranEd - tranEd)); + + m_transEdError += std::abs(data.tranEd - tranEd); + + data.tranEq = tranEq; + data.tranEd = tranEd; + + if(!syncGenerator->IsOnline()) { + std::complex e; + DQ0toABC(data.tranEd, data.tranEq, data.delta, e); + data.terminalVoltage = e; + } + } break; + case SM_MODEL_4: { + } break; + case SM_MODEL_5: { + } break; } syncGenerator->SetElectricalData(data); @@ -758,6 +786,7 @@ bool Electromechanical::SolveSynchronousMachines() return false; } } + m_numIt = iterations; // Solve controllers. for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { @@ -833,4 +862,18 @@ void Electromechanical::SaveData() syncGenerator->SetElectricalData(data); } } + for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { + Bus* bus = *it; + auto data = bus->GetElectricalData(); + if(data.plotBus) { + data.stabVoltageVector.push_back(m_vBus[data.number]); + bus->SetElectricalData(data); + } + } + + m_wErrorVector.push_back(m_wError); + m_deltaErrorVector.push_back(m_deltaError); + m_transEdErrorVector.push_back(m_transEdError); + m_transEqErrorVector.push_back(m_transEqError); + m_numItVector.push_back(m_numIt); } -- cgit From 73eacaf08f695bb3261f072d82ba9fe88366f1c1 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Mon, 29 May 2017 20:51:26 -0300 Subject: Unit data and control import bugs fixed --- Project/Electromechanical.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index a81912f..cdfac9c 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -151,7 +151,6 @@ void Electromechanical::SetEvent(double currentTime) // Remove machine (only connected machines) if(swData.swType[i] == SW_REMOVE && generator->IsOnline()) { generator->SetOnline(false); - auto data = generator->GetPUElectricalData(m_powerSystemBase); int n = static_cast(generator->GetParentList()[0])->GetElectricalData().number; m_yBus[n][n] -= GetSyncMachineAdmittance(generator); } @@ -159,7 +158,6 @@ void Electromechanical::SetEvent(double currentTime) // Insert machine (only disconnected machines) if(swData.swType[i] == SW_INSERT && !generator->IsOnline() && generator->GetParentList().size() == 1) { if(generator->SetOnline(true)) { - auto data = generator->GetPUElectricalData(m_powerSystemBase); int n = static_cast(generator->GetParentList()[0])->GetElectricalData().number; m_yBus[n][n] += GetSyncMachineAdmittance(generator); } @@ -385,7 +383,7 @@ bool Electromechanical::EventTrigger(double eventTime, double currentTime) std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* generator) { - auto data = generator->GetPUElectricalData(m_powerSystemBase); + auto data = generator->GetElectricalData(); double k = 1.0; // Power base change factor. if(data.useMachineBase) { double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); @@ -425,7 +423,8 @@ bool Electromechanical::InitializeDynamicElements() // Synchronous generators for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; - auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); + auto dataPU = syncGenerator->GetPUElectricalData(m_powerSystemBase); + auto data = syncGenerator->GetElectricalData(); if(syncGenerator->IsOnline()) { double k = 1.0; // Power base change factor. if(data.useMachineBase) { @@ -434,7 +433,7 @@ bool Electromechanical::InitializeDynamicElements() } data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; - std::complex conjS(data.activePower, -data.reactivePower); + std::complex conjS(dataPU.activePower, -dataPU.reactivePower); std::complex conjV = std::conj(data.terminalVoltage); std::complex ia = conjS / conjV; @@ -458,7 +457,7 @@ bool Electromechanical::InitializeDynamicElements() data.speed = 2.0 * M_PI * m_systemFreq; data.pe = data.pm; - data.electricalPower = std::complex(data.activePower, data.reactivePower); + data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); switch(GetMachineModel(syncGenerator)) { case SM_MODEL_1: { @@ -554,7 +553,7 @@ void Electromechanical::CalculateMachinesCurrents() for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; - auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); + auto data = syncGenerator->GetElectricalData(); if(syncGenerator->IsOnline()) { double k = 1.0; // Power base change factor. if(data.useMachineBase) { @@ -615,7 +614,7 @@ void Electromechanical::CalculateMachinesCurrents() void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq) { - auto data = syncGenerator->GetPUElectricalData(m_powerSystemBase); + auto data = syncGenerator->GetElectricalData(); double k = 1.0; // Power base change factor. if(data.useMachineBase) { -- cgit From 0e0a956edd98d71b22d5be2c85b0fc2049a77c4a Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Tue, 30 May 2017 17:56:04 -0300 Subject: Control step different from machines implemented --- Project/Electromechanical.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index cdfac9c..c08e4fd 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -512,8 +512,8 @@ bool Electromechanical::InitializeDynamicElements() // Initialize controllers if(data.useAVR) { if(data.avrSolver) delete data.avrSolver; - data.avrSolver = new ControlElementSolver(data.avr, m_timeStep, 1e-3, false, - std::abs(data.terminalVoltage), m_parent); + data.avrSolver = new ControlElementSolver(data.avr, m_timeStep * m_ctrlTimeStepMultiplier, m_tolerance, + false, std::abs(data.terminalVoltage), m_parent); if(!data.avrSolver->IsOK()) { m_errorMsg = _("Error on initializate the AVR of \"") + data.name + _("\"."); syncGenerator->SetElectricalData(data); @@ -522,8 +522,8 @@ bool Electromechanical::InitializeDynamicElements() } if(data.useSpeedGovernor) { if(data.speedGovSolver) delete data.speedGovSolver; - data.speedGovSolver = - new ControlElementSolver(data.speedGov, m_timeStep, 1e-3, false, data.speed, m_parent); + data.speedGovSolver = new ControlElementSolver(data.speedGov, m_timeStep * m_ctrlTimeStepMultiplier, + m_tolerance, false, data.speed, m_parent); if(!data.speedGovSolver->IsOK()) { m_errorMsg = _("Error on initializate the speed governor of \"") + data.name + _("\"."); syncGenerator->SetElectricalData(data); @@ -673,8 +673,7 @@ bool Electromechanical::SolveSynchronousMachines() // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq); - } - else { + } else { CalculateIntegrationConstants(syncGenerator, 0.0f, 0.0f); } } @@ -762,7 +761,7 @@ bool Electromechanical::SolveSynchronousMachines() data.tranEq = tranEq; data.tranEd = tranEd; - + if(!syncGenerator->IsOnline()) { std::complex e; DQ0toABC(data.tranEd, data.tranEq, data.delta, e); @@ -788,16 +787,17 @@ bool Electromechanical::SolveSynchronousMachines() m_numIt = iterations; // Solve controllers. + int ctrlRatio = static_cast(1 / m_ctrlTimeStepMultiplier); for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; auto data = syncGenerator->GetElectricalData(); if(data.useAVR && data.avrSolver) { - data.avrSolver->SolveNextStep(std::abs(data.terminalVoltage)); + for(int i = 0; i < ctrlRatio; ++i) data.avrSolver->SolveNextStep(std::abs(data.terminalVoltage)); data.fieldVoltage = data.initialFieldVoltage + data.avrSolver->GetLastSolution(); } if(data.useSpeedGovernor && data.speedGovSolver) { - data.speedGovSolver->SolveNextStep(data.speed); + for(int i = 0; i < ctrlRatio; ++i) data.speedGovSolver->SolveNextStep(data.speed); data.pm = data.speedGovSolver->GetLastSolution(); } syncGenerator->SetElectricalData(data); -- cgit From 54291220edb747c2a059ac5d316f6f7b14445404 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Fri, 2 Jun 2017 20:08:33 -0300 Subject: Insert capacitor/inductor fixed --- Project/Electromechanical.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index c08e4fd..8d683c4 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -322,7 +322,7 @@ void Electromechanical::SetEvent(double currentTime) capacitor->SetOnline(false); auto data = capacitor->GetPUElectricalData(m_powerSystemBase); int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; - m_yBus[n][n] += std::complex(0.0, data.reactivePower); + m_yBus[n][n] -= std::complex(0.0, data.reactivePower); } // Insert capacitor (only disconnected capacitors) @@ -330,7 +330,7 @@ void Electromechanical::SetEvent(double currentTime) if(capacitor->SetOnline(true)) { auto data = capacitor->GetPUElectricalData(m_powerSystemBase); int n = static_cast(capacitor->GetParentList()[0])->GetElectricalData().number; - m_yBus[n][n] -= std::complex(0.0, data.reactivePower); + m_yBus[n][n] += std::complex(0.0, data.reactivePower); } } } @@ -348,7 +348,7 @@ void Electromechanical::SetEvent(double currentTime) inductor->SetOnline(false); auto data = inductor->GetPUElectricalData(m_powerSystemBase); int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; - m_yBus[n][n] += std::complex(0.0, -data.reactivePower); + m_yBus[n][n] -= std::complex(0.0, -data.reactivePower); } // Insert inductor (only disconnected inductors) @@ -356,7 +356,7 @@ void Electromechanical::SetEvent(double currentTime) if(inductor->SetOnline(true)) { auto data = inductor->GetPUElectricalData(m_powerSystemBase); int n = static_cast(inductor->GetParentList()[0])->GetElectricalData().number; - m_yBus[n][n] -= std::complex(0.0, -data.reactivePower); + m_yBus[n][n] += std::complex(0.0, -data.reactivePower); } } } -- 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/Electromechanical.cpp | 105 ++++++++++++++++++++++++++---------------- 1 file changed, 65 insertions(+), 40 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 8d683c4..9830de9 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -394,15 +394,17 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* double xq = 0.0; double ra = data.armResistance * k; - switch(GetMachineModel(generator)) { - case SM_MODEL_1: - case SM_MODEL_2: - case SM_MODEL_3: { + switch(data.model) { + case Machines::SM_MODEL_1: { + xq = data.transXd * k; + } break; + case Machines::SM_MODEL_2: + case Machines::SM_MODEL_3: { xd = data.transXd * k; xq = data.transXq * k; } break; - case SM_MODEL_4: - case SM_MODEL_5: { + case Machines::SM_MODEL_4: + case Machines::SM_MODEL_5: { xd = data.subXd * k; xq = data.subXq * k; } break; @@ -432,6 +434,8 @@ bool Electromechanical::InitializeDynamicElements() k = m_powerSystemBase / oldBase; } data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; + + data.model = GetMachineModel(syncGenerator); std::complex conjS(dataPU.activePower, -dataPU.reactivePower); std::complex conjV = std::conj(data.terminalVoltage); @@ -440,6 +444,8 @@ bool Electromechanical::InitializeDynamicElements() double xd = data.syncXd * k; double xq = data.syncXq * k; double ra = data.armResistance * k; + + if(data.model == Machines::SM_MODEL_1) xq = data.transXd * k; // Initialize state variables std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; @@ -459,14 +465,17 @@ bool Electromechanical::InitializeDynamicElements() data.pe = data.pm; data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); - switch(GetMachineModel(syncGenerator)) { - case SM_MODEL_1: { - data.tranEq = 0.0; + switch(data.model) { + case Machines::SM_MODEL_1: { + double tranXd = data.transXd * k; + + data.tranEq = data.initialFieldVoltage - tranXd * id0; + wxMessageBox(wxString::Format("%.5f", data.tranEq)); data.tranEd = 0.0; data.subEq = 0.0; data.subEd = 0.0; } break; - case SM_MODEL_2: { + case Machines::SM_MODEL_2: { double tranXd = data.transXd * k; data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; @@ -474,7 +483,7 @@ bool Electromechanical::InitializeDynamicElements() data.subEd = 0.0; data.subEq = 0.0; } break; - case SM_MODEL_3: { + case Machines::SM_MODEL_3: { double tranXd = data.transXd * k; double tranXq = data.transXq * k; @@ -484,7 +493,7 @@ bool Electromechanical::InitializeDynamicElements() data.subEd = 0.0; data.subEq = 0.0; } break; - case SM_MODEL_4: { + case Machines::SM_MODEL_4: { double tranXd = data.transXd * k; double subXd = data.subXd * k; double subXq = data.subXq * k; @@ -494,7 +503,7 @@ bool Electromechanical::InitializeDynamicElements() data.subEq = data.tranEq + (tranXd - subXd) * id0; data.subEd = -(xq - subXq) * iq0; } break; - case SM_MODEL_5: { + case Machines::SM_MODEL_5: { double tranXd = data.transXd * k; double tranXq = data.transXq * k; double subXd = data.subXd * k; @@ -570,16 +579,19 @@ void Electromechanical::CalculateMachinesCurrents() std::complex iInj = m_iBus[n]; double xdq = 0.0; - switch(GetMachineModel(syncGenerator)) { - case SM_MODEL_1: - case SM_MODEL_2: - case SM_MODEL_3: { + switch(data.model) { + case Machines::SM_MODEL_1: { + DQ0toABC(data.tranEd, data.tranEq, data.delta, e); + xq = data.transXd * k; + } break; + case Machines::SM_MODEL_2: + case Machines::SM_MODEL_3: { DQ0toABC(data.tranEd, data.tranEq, data.delta, e); xd = data.transXd * k; xq = data.transXq * k; } break; - case SM_MODEL_4: - case SM_MODEL_5: { + case Machines::SM_MODEL_4: + case Machines::SM_MODEL_5: { DQ0toABC(data.subEd, data.subEq, data.delta, e); xd = data.subXd * k; xq = data.subXq * k; @@ -633,24 +645,37 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat data.icDelta.c = data.delta + data.icDelta.m * (data.speed - 2.0f * w0); // Eq' - data.icTranEq.m = m_timeStep / (2.0f * data.transTd0 + m_timeStep); - data.icTranEq.c = (1.0f - 2.0 * data.icTranEq.m) * data.tranEq + - data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); + if(data.model == Machines::SM_MODEL_2 || data.model == Machines::SM_MODEL_3 || data.model == Machines::SM_MODEL_4 || + data.model == Machines::SM_MODEL_5) { + data.icTranEq.m = m_timeStep / (2.0f * data.transTd0 + m_timeStep); + data.icTranEq.c = (1.0f - 2.0 * data.icTranEq.m) * data.tranEq + + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); + } // Ed' - data.icTranEd.m = m_timeStep / (2.0f * data.transTq0 + m_timeStep); - data.icTranEd.c = - (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; + if(data.model == Machines::SM_MODEL_3 || data.model == Machines::SM_MODEL_4 || data.model == Machines::SM_MODEL_5) { + data.icTranEd.m = m_timeStep / (2.0f * data.transTq0 + m_timeStep); + data.icTranEd.c = + (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; + } // Eq'' - data.icSubEq.m = m_timeStep / (2.0f * data.subTd0 + m_timeStep); - data.icSubEq.c = (1.0f - 2.0f * data.icSubEq.m) * data.subEq + - data.icSubEq.m * (data.tranEq + (data.transXd * k - data.subXd * k) * id); - + if(data.model == Machines::SM_MODEL_4 || data.model == Machines::SM_MODEL_5) { + data.icSubEq.m = m_timeStep / (2.0f * data.subTd0 + m_timeStep); + data.icSubEq.c = (1.0f - 2.0f * data.icSubEq.m) * data.subEq + + data.icSubEq.m * (data.tranEq + (data.transXd * k - data.subXd * k) * id); + } // Ed'' - data.icSubEd.m = m_timeStep / (2.0f * data.subTq0 + m_timeStep); - data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd + - data.icSubEd.m * (data.tranEd - (data.transXq * k - data.subXq * k) * iq); + if(data.model == Machines::SM_MODEL_4) { + data.icSubEd.m = m_timeStep / (2.0f * data.subTq0 + m_timeStep); + data.icSubEd.c = + (1.0f - 2.0f * data.icSubEd.m) * data.subEd - data.icSubEd.m * (data.syncXq * k - data.subXq * k) * iq; + } + if(data.model == Machines::SM_MODEL_5) { + data.icSubEd.m = m_timeStep / (2.0f * data.subTq0 + m_timeStep); + data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd + + data.icSubEd.m * (data.tranEd - (data.transXq * k - data.subXq * k) * iq); + } syncGenerator->SetElectricalData(data); } @@ -719,11 +744,11 @@ bool Electromechanical::SolveSynchronousMachines() double delta = data.icDelta.c + data.icDelta.m * w; error = std::max(error, std::abs(data.delta - delta)); + m_deltaError += std::abs(data.delta - delta); + data.speed = w; data.delta = delta; - m_deltaError += std::abs(data.delta - delta); - // Electric power. double id, iq, vd, vq, pe; ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); @@ -741,13 +766,13 @@ bool Electromechanical::SolveSynchronousMachines() data.pe = pe; // Electrical differential equations - switch(GetMachineModel(syncGenerator)) { - case SM_MODEL_1: { + switch(data.model) { + case Machines::SM_MODEL_1: { // There is no differential equations. } break; - case SM_MODEL_2: { + case Machines::SM_MODEL_2: { } break; - case SM_MODEL_3: { + case Machines::SM_MODEL_3: { double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); error = std::max(error, std::abs(data.tranEq - tranEq)); @@ -768,9 +793,9 @@ bool Electromechanical::SolveSynchronousMachines() data.terminalVoltage = e; } } break; - case SM_MODEL_4: { + case Machines::SM_MODEL_4: { } break; - case SM_MODEL_5: { + case Machines::SM_MODEL_5: { } break; } -- 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/Electromechanical.cpp | 400 ++++++++++++++++++++++++++++++------------ 1 file changed, 289 insertions(+), 111 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 9830de9..0a959c2 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -14,6 +14,8 @@ bool Electromechanical::RunStabilityCalculation() wxProgressDialog pbd(_("Running simulation"), _("Initializing..."), 100, m_parent, wxPD_APP_MODAL | wxPD_AUTO_HIDE | wxPD_CAN_ABORT | wxPD_SMOOTH); + SetSyncMachinesModel(); + // Calculate the admittance matrix with the synchronous machines. if(!GetYBus(m_yBus, m_powerSystemBase, POSITIVE_SEQ, false, true)) { m_errorMsg = _("It was not possible to build the admittance matrix."); @@ -397,16 +399,29 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* switch(data.model) { case Machines::SM_MODEL_1: { xq = data.transXd * k; + xd = xq; + } break; + case Machines::SM_MODEL_2: { + xd = data.transXd * k; + xq = data.transXq * k; + if(xq == 0.0) { + xq = data.syncXq * k; + if(xq == 0.0) { + xq = data.syncXd * k; + } + } } break; - case Machines::SM_MODEL_2: case Machines::SM_MODEL_3: { xd = data.transXd * k; xq = data.transXq * k; + if(xq == 0.0) xq = xd; } break; case Machines::SM_MODEL_4: case Machines::SM_MODEL_5: { xd = data.subXd * k; xq = data.subXq * k; + if(xd == 0.0) xd = xq; + if(xq == 0.0) xq = xd; } break; } double xdq = 0.5 * (xd + xq); @@ -434,8 +449,6 @@ bool Electromechanical::InitializeDynamicElements() k = m_powerSystemBase / oldBase; } data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; - - data.model = GetMachineModel(syncGenerator); std::complex conjS(dataPU.activePower, -dataPU.reactivePower); std::complex conjV = std::conj(data.terminalVoltage); @@ -444,8 +457,12 @@ bool Electromechanical::InitializeDynamicElements() double xd = data.syncXd * k; double xq = data.syncXq * k; double ra = data.armResistance * k; - - if(data.model == Machines::SM_MODEL_1) xq = data.transXd * k; + + if(data.model == Machines::SM_MODEL_1) { + xq = data.transXd * k; + xd = xq; + } + if(data.syncXq == 0.0) xq = data.syncXd * k; // Initialize state variables std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; @@ -465,12 +482,18 @@ bool Electromechanical::InitializeDynamicElements() data.pe = data.pm; data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); + // Variables to extrapolate. + data.oldIq = iq0; + data.oldId = id0; + data.oldPe = data.pe; + switch(data.model) { case Machines::SM_MODEL_1: { - double tranXd = data.transXd * k; - - data.tranEq = data.initialFieldVoltage - tranXd * id0; - wxMessageBox(wxString::Format("%.5f", data.tranEq)); + // double tranXd = data.transXd * k; + + // data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; + data.tranEq = std::abs(eq0); + data.tranEd = 0.0; data.subEq = 0.0; data.subEd = 0.0; @@ -486,6 +509,7 @@ bool Electromechanical::InitializeDynamicElements() case Machines::SM_MODEL_3: { double tranXd = data.transXd * k; double tranXq = data.transXq * k; + if(tranXq == 0.0) tranXq = tranXd; data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; data.tranEd = -(xq - tranXq) * iq0; @@ -497,6 +521,8 @@ bool Electromechanical::InitializeDynamicElements() double tranXd = data.transXd * k; double subXd = data.subXd * k; double subXq = data.subXq * k; + if(subXd == 0.0) subXd = subXq; + if(subXq == 0.0) subXq = subXd; data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; data.tranEd = 0.0; @@ -508,11 +534,13 @@ bool Electromechanical::InitializeDynamicElements() double tranXq = data.transXq * k; double subXd = data.subXd * k; double subXq = data.subXq * k; + if(subXd == 0.0) subXd = subXq; + if(subXq == 0.0) subXq = subXd; data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; data.tranEd = -(xq - tranXq) * iq0; data.subEq = data.tranEq + (tranXd - subXd) * id0; - data.subEd = data.tranEd + (tranXq - subXq) * iq0; + data.subEd = data.tranEd - (tranXq - subXq) * iq0; } break; default: break; @@ -583,18 +611,32 @@ void Electromechanical::CalculateMachinesCurrents() case Machines::SM_MODEL_1: { DQ0toABC(data.tranEd, data.tranEq, data.delta, e); xq = data.transXd * k; + xd = xq; + } break; + case Machines::SM_MODEL_2: { + DQ0toABC(data.tranEd, data.tranEq, data.delta, e); + xd = data.transXd * k; + xq = data.transXq * k; + if(xq == 0.0) { + xq = data.syncXq * k; + if(xq == 0.0) { + xq = data.syncXd * k; + } + } } break; - case Machines::SM_MODEL_2: case Machines::SM_MODEL_3: { DQ0toABC(data.tranEd, data.tranEq, data.delta, e); xd = data.transXd * k; xq = data.transXq * k; + if(xq == 0.0) xq = xd; } break; case Machines::SM_MODEL_4: case Machines::SM_MODEL_5: { DQ0toABC(data.subEd, data.subEq, data.delta, e); xd = data.subXd * k; xq = data.subXq * k; + if(xd == 0.0) xd = xq; + if(xq == 0.0) xq = xd; } break; } xdq = 0.5 * (xd + xq); @@ -602,15 +644,12 @@ void Electromechanical::CalculateMachinesCurrents() std::complex y0 = std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); std::complex iUnadj = y0 * e; - double dVR = std::real(e) - std::real(v); - double dVI = std::imag(e) - std::imag(v); + std::complex iAdj = + std::complex(0.0, -((0.5 * (xq - xd)) / (ra * ra + xd * xq))) * (std::conj(e) - std::conj(v)); + iAdj = iAdj * std::cos(2.0 * data.delta) + iAdj * std::complex(0.0, std::sin(2.0 * data.delta)); - double iAdjR = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * - (-std::sin(2.0 * data.delta) * dVR + std::cos(2.0 * data.delta) * dVI); - double iAdjI = ((0.5 * (xd - xq)) / (ra * ra + xd * xq)) * - (std::cos(2.0 * data.delta) * dVR + std::sin(2.0 * data.delta) * dVI); + iInj = iUnadj + iAdj; - iInj = iUnadj + std::complex(iAdjR, iAdjI); m_iBus[n] += iInj; std::complex iMachine = iInj - y0 * v; @@ -624,17 +663,33 @@ void Electromechanical::CalculateMachinesCurrents() } } -void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq) +void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq, double k) { auto data = syncGenerator->GetElectricalData(); - - double k = 1.0; // Power base change factor. - if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); - k = m_powerSystemBase / oldBase; - } double w0 = 2.0f * M_PI * m_systemFreq; + double syncXd, syncXq, transXd, transXq, subXd, subXq; + syncXd = data.syncXd * k; + syncXq = data.syncXq * k; + transXd = data.transXd * k; + transXq = data.transXq * k; + subXd = data.subXd * k; + subXq = data.subXq * k; + + if(syncXq == 0.0) syncXq = syncXd; + if(transXq == 0.0) transXq = transXd; + if(subXd == 0.0) subXd = subXq; + if(subXq == 0.0) subXq = subXd; + + double transTd0, transTq0, subTd0, subTq0; + transTd0 = data.transTd0; + transTq0 = data.transTq0; + subTd0 = data.subTd0; + subTq0 = data.subTq0; + + if(subTd0 == 0.0) subTd0 = subTq0; + if(subTq0 == 0.0) subTq0 = subTd0; + // Speed data.icSpeed.m = m_timeStep / ((4.0f * data.inertia / w0) / k + m_timeStep * data.damping * k); data.icSpeed.c = (1.0f - 2.0f * data.icSpeed.m * data.damping * k) * data.speed + @@ -647,34 +702,32 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat // Eq' if(data.model == Machines::SM_MODEL_2 || data.model == Machines::SM_MODEL_3 || data.model == Machines::SM_MODEL_4 || data.model == Machines::SM_MODEL_5) { - data.icTranEq.m = m_timeStep / (2.0f * data.transTd0 + m_timeStep); + data.icTranEq.m = m_timeStep / (2.0f * transTd0 + m_timeStep); data.icTranEq.c = (1.0f - 2.0 * data.icTranEq.m) * data.tranEq + - data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); } // Ed' if(data.model == Machines::SM_MODEL_3 || data.model == Machines::SM_MODEL_4 || data.model == Machines::SM_MODEL_5) { - data.icTranEd.m = m_timeStep / (2.0f * data.transTq0 + m_timeStep); - data.icTranEd.c = - (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; + data.icTranEd.m = m_timeStep / (2.0f * transTq0 + m_timeStep); + data.icTranEd.c = (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (syncXq - transXq) * iq; } // Eq'' if(data.model == Machines::SM_MODEL_4 || data.model == Machines::SM_MODEL_5) { - data.icSubEq.m = m_timeStep / (2.0f * data.subTd0 + m_timeStep); - data.icSubEq.c = (1.0f - 2.0f * data.icSubEq.m) * data.subEq + - data.icSubEq.m * (data.tranEq + (data.transXd * k - data.subXd * k) * id); + data.icSubEq.m = m_timeStep / (2.0f * subTd0 + m_timeStep); + data.icSubEq.c = + (1.0f - 2.0f * data.icSubEq.m) * data.subEq + data.icSubEq.m * (data.tranEq + (transXd - subXd) * id); } // Ed'' if(data.model == Machines::SM_MODEL_4) { - data.icSubEd.m = m_timeStep / (2.0f * data.subTq0 + m_timeStep); - data.icSubEd.c = - (1.0f - 2.0f * data.icSubEd.m) * data.subEd - data.icSubEd.m * (data.syncXq * k - data.subXq * k) * iq; + data.icSubEd.m = m_timeStep / (2.0f * subTq0 + m_timeStep); + data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd - data.icSubEd.m * (syncXq - subXq) * iq; } if(data.model == Machines::SM_MODEL_5) { - data.icSubEd.m = m_timeStep / (2.0f * data.subTq0 + m_timeStep); - data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd + - data.icSubEd.m * (data.tranEd - (data.transXq * k - data.subXq * k) * iq); + data.icSubEd.m = m_timeStep / (2.0f * subTq0 + m_timeStep); + data.icSubEd.c = + (1.0f - 2.0f * data.icSubEd.m) * data.subEd + data.icSubEd.m * (data.tranEd - (transXq - subXq) * iq); } syncGenerator->SetElectricalData(data); @@ -682,22 +735,37 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat bool Electromechanical::SolveSynchronousMachines() { - double w0 = 2.0 * M_PI * m_systemFreq; - + // CalculateMachinesCurrents(); for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; auto data = syncGenerator->GetElectricalData(); if(syncGenerator->IsOnline()) { int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; - double id, iq; + double id, iq, pe; + + pe = data.pe; + + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); ABCtoDQ0(iMachine, data.delta, id, iq); // Calculate integration constants. - CalculateIntegrationConstants(syncGenerator, id, iq); + CalculateIntegrationConstants(syncGenerator, id, iq, k); + + CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, pe, k); + // Extrapolate nonintegrable variables. + id = 2.0 * id - data.oldId; + iq = 2.0 * iq - data.oldIq; + pe = 2.0 * pe - data.oldPe; + + CalculateSyncMachineIntVariables(syncGenerator, id, iq, pe, k); } else { CalculateIntegrationConstants(syncGenerator, 0.0f, 0.0f); } @@ -722,87 +790,25 @@ bool Electromechanical::SolveSynchronousMachines() // Solve machine equations. for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetElectricalData(); + double id, iq, pe; double k = 1.0; // Power base change factor. if(data.useMachineBase) { double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } - int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; - - if(syncGenerator->IsOnline()) { - data.terminalVoltage = m_vBus[n]; - } - - // Mechanical differential equations. - double w = data.icSpeed.c + data.icSpeed.m * (data.pm - data.pe); - error = std::max(error, std::abs(data.speed - w) / w0); - - m_wError += std::abs(data.speed - w) / w0; - - double delta = data.icDelta.c + data.icDelta.m * w; - error = std::max(error, std::abs(data.delta - delta)); - - m_deltaError += std::abs(data.delta - delta); - - data.speed = w; - data.delta = delta; - - // Electric power. - double id, iq, vd, vq, pe; - ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); - - if(syncGenerator->IsOnline()) { - std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); - ABCtoDQ0(iMachine, data.delta, id, iq); - - pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; - // pe = (2 * pe - data.pe); // Extrapolating Pe. - } else { - pe = id = iq = 0.0f; - } - - data.pe = pe; - - // Electrical differential equations - switch(data.model) { - case Machines::SM_MODEL_1: { - // There is no differential equations. - } break; - case Machines::SM_MODEL_2: { - } break; - case Machines::SM_MODEL_3: { - double tranEq = data.icTranEq.c + - data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); - error = std::max(error, std::abs(data.tranEq - tranEq)); - m_transEqError += std::abs(data.tranEq - tranEq); + CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, pe, k); - double tranEd = data.icTranEd.c - data.icTranEd.m * (data.syncXq * k - data.transXq * k) * iq; - error = std::max(error, std::abs(data.tranEd - tranEd)); + double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, pe, k); - m_transEdError += std::abs(data.tranEd - tranEd); - - data.tranEq = tranEq; - data.tranEd = tranEd; - - if(!syncGenerator->IsOnline()) { - std::complex e; - DQ0toABC(data.tranEd, data.tranEq, data.delta, e); - data.terminalVoltage = e; - } - } break; - case Machines::SM_MODEL_4: { - } break; - case Machines::SM_MODEL_5: { - } break; - } - - syncGenerator->SetElectricalData(data); + if(genError > error) error = genError; } ++iterations; + if(iterations > m_maxIterations) { m_errorMsg = _("Impossible to solve the synchronous generators.\nCheck the system parameters and/or " "decrease the time step."); @@ -901,3 +907,175 @@ void Electromechanical::SaveData() m_transEqErrorVector.push_back(m_transEqError); m_numItVector.push_back(m_numIt); } + +void Electromechanical::SetSyncMachinesModel() +{ + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + auto data = syncGenerator->GetElectricalData(); + data.model = GetMachineModel(syncGenerator); + syncGenerator->SetElectricalData(data); + } +} + +void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncGenerator, + double& id, + double& iq, + double& pe, + double k) +{ + auto data = syncGenerator->GetElectricalData(); + int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; + + if(syncGenerator->IsOnline()) { + data.terminalVoltage = m_vBus[n]; + } + + double vd, vq; + ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); + + if(syncGenerator->IsOnline()) { + std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); + ABCtoDQ0(iMachine, data.delta, id, iq); + + pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; + } else { + pe = id = iq = 0.0f; + } + data.pe = pe; + data.oldId = id; + data.oldIq = iq; + syncGenerator->SetElectricalData(data); +} + +double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGenerator, + double id, + double iq, + double pe, + double k) +{ + double w0 = 2.0 * M_PI * m_systemFreq; + double error = 0.0; + auto data = syncGenerator->GetElectricalData(); + + // Mechanical differential equations. + double w = data.icSpeed.c + data.icSpeed.m * (data.pm - pe); + error = std::max(error, std::abs(data.speed - w) / w0); + + m_wError += std::abs(data.speed - w) / w0; + + double delta = data.icDelta.c + data.icDelta.m * w; + error = std::max(error, std::abs(data.delta - delta)); + + m_deltaError += std::abs(data.delta - delta); + + data.speed = w; + data.delta = delta; + + // Electrical differential equations + switch(data.model) { + case Machines::SM_MODEL_1: { + // There is no differential equations. + } break; + case Machines::SM_MODEL_2: { + double tranEq = + data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); + error = std::max(error, std::abs(data.tranEq - tranEq)); + + m_transEqError += std::abs(data.tranEq - tranEq); + + data.tranEq = tranEq; + } break; + case Machines::SM_MODEL_3: { + double syncXd, syncXq, transXd, transXq; + syncXd = data.syncXd * k; + syncXq = data.syncXq * k; + transXd = data.transXd * k; + transXq = data.transXq * k; + if(syncXq == 0.0) syncXq = syncXd; + if(transXq == 0.0) transXq = transXd; + + double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); + error = std::max(error, std::abs(data.tranEq - tranEq)); + + m_transEqError += std::abs(data.tranEq - tranEq); + + double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq; + error = std::max(error, std::abs(data.tranEd - tranEd)); + + m_transEdError += std::abs(data.tranEd - tranEd); + + data.tranEq = tranEq; + data.tranEd = tranEd; + + if(!syncGenerator->IsOnline()) { + std::complex e; + DQ0toABC(data.tranEd, data.tranEq, data.delta, e); + data.terminalVoltage = e; + } + } break; + case Machines::SM_MODEL_4: { + double syncXd, syncXq, transXd, subXd, subXq; + syncXd = data.syncXd * k; + syncXq = data.syncXq * k; + transXd = data.transXd * k; + subXd = data.subXd * k; + subXq = data.subXq * k; + if(syncXq == 0.0) syncXq = syncXd; + if(subXd == 0.0) subXd = subXq; + if(subXq == 0.0) subXq = subXd; + + double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); + error = std::max(error, std::abs(data.tranEq - tranEq)); + + m_transEqError += std::abs(data.tranEq - tranEq); + + double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id); + error = std::max(error, std::abs(data.subEq - subEq)); + + double subEd = data.icSubEd.c - data.icSubEd.m * (syncXq - subXq) * iq; + error = std::max(error, std::abs(data.subEd - subEd)); + + data.tranEq = tranEq; + data.subEq = subEq; + data.subEd = subEd; + } break; + case Machines::SM_MODEL_5: { + double syncXd, syncXq, transXd, transXq, subXd, subXq; + syncXd = data.syncXd * k; + syncXq = data.syncXq * k; + transXd = data.transXd * k; + transXq = data.transXq * k; + subXd = data.subXd * k; + subXq = data.subXq * k; + if(syncXq == 0.0) syncXq = syncXd; + if(transXq == 0.0) transXq = transXd; + if(subXd == 0.0) subXd = subXq; + if(subXq == 0.0) subXq = subXd; + + double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); + error = std::max(error, std::abs(data.tranEq - tranEq)); + + m_transEqError += std::abs(data.tranEq - tranEq); + + double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq; + error = std::max(error, std::abs(data.tranEd - tranEd)); + + m_transEdError += std::abs(data.tranEd - tranEd); + + double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id); + error = std::max(error, std::abs(data.subEq - subEq)); + + double subEd = data.icSubEd.c + data.icSubEd.m * (tranEd - (transXq - subXq) * iq); + error = std::max(error, std::abs(data.subEd - subEd)); + + data.tranEq = tranEq; + data.tranEd = tranEd; + data.subEq = subEq; + data.subEd = subEd; + } break; + } + + syncGenerator->SetElectricalData(data); + return error; +} -- cgit From 8a1ffbc01135a1466ad7400518e4c56a4bdc3af5 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Wed, 9 Aug 2017 21:01:41 -0300 Subject: Stability properties implemented --- Project/Electromechanical.cpp | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 0a959c2..e26a21d 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -1,11 +1,22 @@ #include "Electromechanical.h" #include "ControlElementSolver.h" -Electromechanical::Electromechanical(wxWindow* parent, std::vector elementList) +Electromechanical::Electromechanical(wxWindow* parent, std::vector elementList, SimulationData data) { m_parent = parent; GetElementsFromList(elementList); SetEventTimeList(); + + m_powerSystemBase = GetPowerValue(data.basePower, data.basePowerUnit); + m_systemFreq = data.stabilityFrequency; + m_simTime = data.stabilitySimulationTime; + m_timeStep = data.timeStep; + m_tolerance = data.stabilityTolerance; + m_maxIterations = data.stabilityMaxIterations; + + m_ctrlTimeStepMultiplier = 1.0 / static_cast(data.controlTimeStepRatio); + + m_plotTime = data.plotTime; } Electromechanical::~Electromechanical() {} @@ -41,27 +52,24 @@ bool Electromechanical::RunStabilityCalculation() if(!InitializeDynamicElements()) return false; - // test - double simTime = 20.0; - double printTime = 0.01; - double pbdTime = 0.01; + double pbdTime = m_plotTime; double currentTime = 0.0; - double currentPrintTime = 0.0; + double currentPlotTime = 0.0; double currentPbdTime = 0.0; - while(currentTime <= simTime) { + while(currentTime < m_simTime) { if(HasEvent(currentTime)) { SetEvent(currentTime); GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); } - if(currentPrintTime >= printTime) { + if(currentPlotTime >= m_plotTime || currentTime == 0.0) { m_timeVector.push_back(currentTime); SaveData(); - currentPrintTime = 0.0; + currentPlotTime = 0.0; } if(currentPbdTime > pbdTime) { - if(!pbd.Update((currentTime / simTime) * 100, wxString::Format("Time = %.2fs", currentTime))) { + if(!pbd.Update((currentTime / m_simTime) * 100, wxString::Format("Time = %.2fs", currentTime))) { m_errorMsg = wxString::Format(_("Simulation cancelled at %.2fs."), currentTime); pbd.Update(100); return false; @@ -72,7 +80,7 @@ bool Electromechanical::RunStabilityCalculation() if(!SolveSynchronousMachines()) return false; currentTime += m_timeStep; - currentPrintTime += m_timeStep; + currentPlotTime += m_timeStep; currentPbdTime += m_timeStep; } return true; -- cgit From e5f19af9a5829a307951785e850e53fc4fcfac90 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 10 Aug 2017 18:07:02 -0300 Subject: Model 1 machine fixed, init file under implementation --- Project/Electromechanical.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index e26a21d..446ade2 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -470,7 +470,7 @@ bool Electromechanical::InitializeDynamicElements() xq = data.transXd * k; xd = xq; } - if(data.syncXq == 0.0) xq = data.syncXd * k; + else if(data.syncXq == 0.0) xq = data.syncXd * k; // Initialize state variables std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; -- cgit From 516cdb72d3ff99a1ee786d3ea24c9b579272fe76 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Wed, 30 Aug 2017 20:42:27 -0300 Subject: COI (center of inertia) implemented See Milano's book p. 342 --- Project/Electromechanical.cpp | 43 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 9 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 446ade2..f72171e 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -17,6 +17,7 @@ Electromechanical::Electromechanical(wxWindow* parent, std::vector ele m_ctrlTimeStepMultiplier = 1.0 / static_cast(data.controlTimeStepRatio); m_plotTime = data.plotTime; + m_useCOI = data.useCOI; } Electromechanical::~Electromechanical() {} @@ -469,8 +470,8 @@ bool Electromechanical::InitializeDynamicElements() if(data.model == Machines::SM_MODEL_1) { xq = data.transXd * k; xd = xq; - } - else if(data.syncXq == 0.0) xq = data.syncXd * k; + } else if(data.syncXq == 0.0) + xq = data.syncXd * k; // Initialize state variables std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; @@ -588,6 +589,7 @@ bool Electromechanical::InitializeDynamicElements() syncGenerator->SetElectricalData(data); } + CalculateReferenceSpeed(); return true; } @@ -673,8 +675,8 @@ void Electromechanical::CalculateMachinesCurrents() void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq, double k) { + CalculateReferenceSpeed(); auto data = syncGenerator->GetElectricalData(); - double w0 = 2.0f * M_PI * m_systemFreq; double syncXd, syncXq, transXd, transXq, subXd, subXq; syncXd = data.syncXd * k; @@ -699,13 +701,13 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat if(subTq0 == 0.0) subTq0 = subTd0; // Speed - data.icSpeed.m = m_timeStep / ((4.0f * data.inertia / w0) / k + m_timeStep * data.damping * k); + data.icSpeed.m = m_timeStep / ((4.0f * data.inertia / m_refSpeed) / k + m_timeStep * data.damping * k); data.icSpeed.c = (1.0f - 2.0f * data.icSpeed.m * data.damping * k) * data.speed + - data.icSpeed.m * (data.pm - data.pe + 2.0f * w0 * data.damping * k); + data.icSpeed.m * (data.pm - data.pe + 2.0f * m_refSpeed * data.damping * k); // Delta data.icDelta.m = 0.5f * m_timeStep; - data.icDelta.c = data.delta + data.icDelta.m * (data.speed - 2.0f * w0); + data.icDelta.c = data.delta + data.icDelta.m * (data.speed - 2.0f * m_refSpeed); // Eq' if(data.model == Machines::SM_MODEL_2 || data.model == Machines::SM_MODEL_3 || data.model == Machines::SM_MODEL_4 || @@ -962,15 +964,14 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe double pe, double k) { - double w0 = 2.0 * M_PI * m_systemFreq; double error = 0.0; auto data = syncGenerator->GetElectricalData(); // Mechanical differential equations. double w = data.icSpeed.c + data.icSpeed.m * (data.pm - pe); - error = std::max(error, std::abs(data.speed - w) / w0); + error = std::max(error, std::abs(data.speed - w) / m_refSpeed); - m_wError += std::abs(data.speed - w) / w0; + m_wError += std::abs(data.speed - w) / m_refSpeed; double delta = data.icDelta.c + data.icDelta.m * w; error = std::max(error, std::abs(data.delta - delta)); @@ -1087,3 +1088,27 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe syncGenerator->SetElectricalData(data); return error; } + +void Electromechanical::CalculateReferenceSpeed() +{ + if(m_useCOI) { + double sumHW = 0.0; + double sumH = 0.0; + for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { + SyncGenerator* syncGenerator = *it; + if(syncGenerator->IsOnline()) { + auto data = syncGenerator->GetElectricalData(); + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + sumH += data.inertia / k; + sumHW += data.inertia * data.speed / k; + } + } + m_refSpeed = sumHW / sumH; + } else { + m_refSpeed = 2.0 * M_PI * m_systemFreq; + } +} -- cgit From 6f3421c4150e49af026432a2a2be0171d741ad03 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Mon, 4 Sep 2017 20:04:42 -0300 Subject: Some bugfixes --- Project/Electromechanical.cpp | 20 -------------------- 1 file changed, 20 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index f72171e..57a275a 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -782,9 +782,6 @@ bool Electromechanical::SolveSynchronousMachines() } m_wError = 0; - m_deltaError = 0; - m_transEdError = 0; - m_transEqError = 0; double error = 1.0; int iterations = 0; @@ -912,9 +909,6 @@ void Electromechanical::SaveData() } m_wErrorVector.push_back(m_wError); - m_deltaErrorVector.push_back(m_deltaError); - m_transEdErrorVector.push_back(m_transEdError); - m_transEqErrorVector.push_back(m_transEqError); m_numItVector.push_back(m_numIt); } @@ -976,8 +970,6 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe double delta = data.icDelta.c + data.icDelta.m * w; error = std::max(error, std::abs(data.delta - delta)); - m_deltaError += std::abs(data.delta - delta); - data.speed = w; data.delta = delta; @@ -991,8 +983,6 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (data.syncXd * k - data.transXd * k) * id); error = std::max(error, std::abs(data.tranEq - tranEq)); - m_transEqError += std::abs(data.tranEq - tranEq); - data.tranEq = tranEq; } break; case Machines::SM_MODEL_3: { @@ -1007,13 +997,9 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); error = std::max(error, std::abs(data.tranEq - tranEq)); - m_transEqError += std::abs(data.tranEq - tranEq); - double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq; error = std::max(error, std::abs(data.tranEd - tranEd)); - m_transEdError += std::abs(data.tranEd - tranEd); - data.tranEq = tranEq; data.tranEd = tranEd; @@ -1037,8 +1023,6 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); error = std::max(error, std::abs(data.tranEq - tranEq)); - m_transEqError += std::abs(data.tranEq - tranEq); - double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id); error = std::max(error, std::abs(data.subEq - subEq)); @@ -1065,13 +1049,9 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe double tranEq = data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); error = std::max(error, std::abs(data.tranEq - tranEq)); - m_transEqError += std::abs(data.tranEq - tranEq); - double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq; error = std::max(error, std::abs(data.tranEd - tranEd)); - m_transEdError += std::abs(data.tranEd - tranEd); - double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id); error = std::max(error, std::abs(data.subEq - subEq)); -- cgit