From a1c362b8a5e989cf9e75be14383fe8f7d3625422 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Tue, 22 Jan 2019 14:49:35 -0200 Subject: Merge error removal --- Project/Electromechanical.cpp | 8 -------- 1 file changed, 8 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index de30184..66d5bb2 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -235,11 +235,7 @@ void Electromechanical::SetEvent(double currentTime) Bus* parentBus = static_cast(load->GetParentList()[0]); int n = parentBus->GetElectricalData().number; std::complex v = parentBus->GetElectricalData().voltage; -<<<<<<< HEAD - m_yBus[n][n] -= std::complex(data.activePower, -data.reactivePower) / (abs(v) * abs(v)); -======= m_yBus[n][n] -= std::complex(data.activePower, -data.reactivePower) / (std::abs(v) * std::abs(v)); ->>>>>>> aad89bf4d16d45c0790bd2fc010d9ec06cc35430 } // Insert load (only disconnected load) @@ -249,11 +245,7 @@ void Electromechanical::SetEvent(double currentTime) Bus* parentBus = static_cast(load->GetParentList()[0]); int n = parentBus->GetElectricalData().number; std::complex v = parentBus->GetElectricalData().voltage; -<<<<<<< HEAD - m_yBus[n][n] += std::complex(data.activePower, -data.reactivePower) / (abs(v) * abs(v)); -======= m_yBus[n][n] += std::complex(data.activePower, -data.reactivePower) / (std::abs(v) * std::abs(v)); ->>>>>>> aad89bf4d16d45c0790bd2fc010d9ec06cc35430 } } } -- cgit From 59c28b89bb2a4eb05a87bb28154442266c35b17d Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Tue, 22 Jan 2019 15:06:09 -0200 Subject: Some code stability optimization --- Project/Electromechanical.cpp | 42 +++++++++++++++++++++++------------------- 1 file changed, 23 insertions(+), 19 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 50ce4e5..f74f546 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -70,7 +70,6 @@ bool Electromechanical::RunStabilityCalculation() wxProgressDialog pbd(_("Running simulation"), _("Initializing..."), 100, m_parent, wxPD_APP_MODAL | wxPD_AUTO_HIDE | wxPD_CAN_ABORT | wxPD_SMOOTH); - PreallocateVectors(); // Reserve the vectors' memory with a estimated size, this can optimize the simulation. SetSyncMachinesModel(); // Calculate the admittance matrix with the synchronous machines. @@ -83,6 +82,7 @@ bool Electromechanical::RunStabilityCalculation() // Get buses voltages. m_vBus.clear(); + m_vBus.shrink_to_fit(); m_vBus.resize(m_busList.size()); for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { Bus* bus = *it; @@ -97,6 +97,7 @@ bool Electromechanical::RunStabilityCalculation() } if(!InitializeDynamicElements()) return false; + PreallocateVectors(); // Reserve the vectors' memory with a estimated size, this can optimize the simulation. double pbdTime = m_plotTime; m_currentTime = 0.0; @@ -109,7 +110,7 @@ bool Electromechanical::RunStabilityCalculation() } if(currentPlotTime >= m_plotTime || m_currentTime == 0.0) { - m_timeVector.push_back(m_currentTime); + m_timeVector.emplace_back(m_currentTime); SaveData(); currentPlotTime = 0.0; } @@ -139,10 +140,10 @@ void Electromechanical::SetEventTimeList() Bus* bus = *it; auto data = bus->GetElectricalData(); 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); + m_eventTimeList.emplace_back(data.stabFaultTime); + m_eventOccurrenceList.emplace_back(false); + m_eventTimeList.emplace_back(data.stabFaultTime + data.stabFaultLength); + m_eventOccurrenceList.emplace_back(false); } } // Switching @@ -150,8 +151,8 @@ void Electromechanical::SetEventTimeList() 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); + m_eventTimeList.emplace_back(swData.swTime[i]); + m_eventOccurrenceList.emplace_back(false); } } } @@ -462,6 +463,7 @@ bool Electromechanical::InitializeDynamicElements() Bus* bus = *it; auto data = bus->GetElectricalData(); data.stabVoltageVector.clear(); + data.stabVoltageVector.shrink_to_fit(); bus->SetElectricalData(data); } // Loads @@ -492,7 +494,9 @@ bool Electromechanical::InitializeDynamicElements() data.qp0 = (data.constPowerReactive / 100.0) * reactivePower; data.voltageVector.clear(); + data.voltageVector.shrink_to_fit(); data.electricalPowerVector.clear(); + data.electricalPowerVector.shrink_to_fit(); if(load->IsOnline()) data.electricalPower = std::complex(activePower, reactivePower); @@ -704,11 +708,17 @@ bool Electromechanical::InitializeDynamicElements() } // Reset plot data data.terminalVoltageVector.clear(); + data.terminalVoltageVector.shrink_to_fit(); data.electricalPowerVector.clear(); + data.electricalPowerVector.shrink_to_fit(); data.mechanicalPowerVector.clear(); + data.mechanicalPowerVector.shrink_to_fit(); data.freqVector.clear(); + data.freqVector.shrink_to_fit(); data.fieldVoltageVector.clear(); + data.fieldVoltageVector.shrink_to_fit(); data.deltaVector.clear(); + data.deltaVector.shrink_to_fit(); syncGenerator->SetElectricalData(data); } @@ -1044,20 +1054,14 @@ void Electromechanical::SaveData() SyncGenerator* syncGenerator = *it; auto data = syncGenerator->GetElectricalData(); if(data.plotSyncMachine) { - data.terminalVoltageVector.push_back(data.terminalVoltage); - 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(wxRadToDeg(data.delta)); - syncGenerator->SetElectricalData(data); + syncGenerator->SavePlotData(); } } 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]); + data.stabVoltageVector.emplace_back(m_vBus[data.number]); bus->SetElectricalData(data); } } @@ -1065,12 +1069,12 @@ void Electromechanical::SaveData() Load* load = *it; auto data = load->GetElectricalData(); if(data.plotLoad) { - data.voltageVector.push_back(data.voltage); - data.electricalPowerVector.push_back(data.electricalPower); + data.voltageVector.emplace_back(data.voltage); + data.electricalPowerVector.emplace_back(data.electricalPower); load->SetElectricalData(data); } } - m_iterationsNumVector.push_back(m_iterationsNum); + m_iterationsNumVector.emplace_back(m_iterationsNum); } void Electromechanical::SetSyncMachinesModel() -- cgit From 2771fff79ac9c3c09b70f4668e7142b2e944d1f2 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 25 Apr 2019 01:25:41 -0300 Subject: Matpower Importer and power quality calculation Power quality in implementation --- 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 72b2e84..92aa88f 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -511,7 +511,7 @@ bool Electromechanical::InitializeDynamicElements() SyncGenerator* syncGenerator = *it; auto dataPU = syncGenerator->GetPUElectricalData(m_powerSystemBase); auto data = syncGenerator->GetElectricalData(); - if(syncGenerator->IsOnline()) { + //if(syncGenerator->IsOnline()) { double k = 1.0; // Power base change factor. if(data.useMachineBase) { double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); @@ -590,7 +590,7 @@ bool Electromechanical::InitializeDynamicElements() data.initialFieldVoltage = ef0 * sd; 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; + syncGenerator->IsOnline() ? data.speed = 2.0 * M_PI * m_systemFreq : data.speed = 2.0 * M_PI * data.ocFrequency; data.delta = delta; data.pe = data.pm; data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); @@ -702,9 +702,9 @@ bool Electromechanical::InitializeDynamicElements() return false; } } - } else { + //} else { // Initialize open circuit machine. - } + //} // Reset plot data data.terminalVoltageVector.clear(); data.terminalVoltageVector.shrink_to_fit(); -- cgit From f54297e08079fe1954920ca2742b0bed19f86181 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Wed, 10 Jul 2019 14:14:30 -0300 Subject: Induction motor implementation start Machine initialization implemented. It seems that it's working. Check a OMIB with the motor. --- Project/Electromechanical.cpp | 511 ++++++++++++++++++++++++++---------------- 1 file changed, 322 insertions(+), 189 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 92aa88f..0f32ee8 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -78,6 +78,7 @@ bool Electromechanical::RunStabilityCalculation() return false; } InsertSyncMachinesOnYBus(); + if(!InsertIndMachinesOnYBus()) return false; GetLUDecomposition(m_yBus, m_yBusL, m_yBusU); // Get buses voltages. @@ -236,7 +237,8 @@ void Electromechanical::SetEvent(double currentTime) 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) / (std::abs(v) * std::abs(v)); + m_yBus[n][n] -= + std::complex(data.activePower, -data.reactivePower) / (std::abs(v) * std::abs(v)); } // Insert load (only disconnected load) @@ -246,7 +248,8 @@ void Electromechanical::SetEvent(double currentTime) 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) / (std::abs(v) * std::abs(v)); + m_yBus[n][n] += + std::complex(data.activePower, -data.reactivePower) / (std::abs(v) * std::abs(v)); } } } @@ -511,199 +514,199 @@ bool Electromechanical::InitializeDynamicElements() SyncGenerator* syncGenerator = *it; auto dataPU = syncGenerator->GetPUElectricalData(m_powerSystemBase); auto data = syncGenerator->GetElectricalData(); - //if(syncGenerator->IsOnline()) { - double k = 1.0; // Power base change factor. - if(data.useMachineBase) { - double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); - k = m_powerSystemBase / oldBase; - } - data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; - - std::complex conjS(dataPU.activePower, -dataPU.reactivePower); - std::complex vt = data.terminalVoltage; - std::complex ia = conjS / std::conj(vt); - - 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; - xd = xq; - } else if(data.syncXq == 0.0) - xq = data.syncXd * k; - - double sd = 1.0; - double sq = 1.0; - double satF = 1.0; - double xp = data.potierReactance * k; - bool hasSaturation = false; - if(data.satFactor != 0.0) { // Have saturation. - satF = (data.satFactor - 1.2) / std::pow(1.2, 7); - if(xp == 0.0) xp = 0.8 * (data.transXd * k); - hasSaturation = true; - } + // if(syncGenerator->IsOnline()) { + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; + + std::complex conjS(dataPU.activePower, -dataPU.reactivePower); + std::complex vt = data.terminalVoltage; + std::complex ia = conjS / std::conj(vt); + + 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; + xd = xq; + } else if(data.syncXq == 0.0) + xq = data.syncXd * k; + + double sd = 1.0; + double sq = 1.0; + double satF = 1.0; + double xp = data.potierReactance * k; + bool hasSaturation = false; + if(data.satFactor != 0.0) { // Have saturation. + satF = (data.satFactor - 1.2) / std::pow(1.2, 7); + if(xp == 0.0) xp = 0.8 * (data.transXd * k); + hasSaturation = true; + } - // Initialize state variables - std::complex eq0 = vt + std::complex(ra, xq) * ia; - double delta = std::arg(eq0); - - double id0, iq0, vd0, vq0; - ABCtoDQ0(ia, delta, id0, iq0); - ABCtoDQ0(vt, delta, vd0, vq0); - - // Initialize saturation - double xqs = xq; - double xds = xd; - if(hasSaturation) { - double oldDelta = 0; - bool exit = false; - int numIt = 0; - while(!exit) { - oldDelta = delta; - ABCtoDQ0(ia, delta, id0, iq0); - ABCtoDQ0(vt, delta, vd0, vq0); - - // Direct-axis Potier voltage. - double epd = vd0 + ra * id0 + xp * iq0; - - sq = 1.0 + satF * (xq / xd) * std::pow(epd, 6); - xqs = (xq - xp) / sq + xp; - eq0 = data.terminalVoltage + std::complex(ra, xqs) * ia; - delta = std::arg(eq0); - if(std::abs(delta - oldDelta) < m_saturationTolerance) { - exit = true; - } else if(numIt >= m_maxIterations) { - m_errorMsg = _("Error on initializate the saturation values of \"") + data.name + _("\"."); - return false; - } - numIt++; + // Initialize state variables + std::complex eq0 = vt + std::complex(ra, xq) * ia; + double delta = std::arg(eq0); + + double id0, iq0, vd0, vq0; + ABCtoDQ0(ia, delta, id0, iq0); + ABCtoDQ0(vt, delta, vd0, vq0); + + // Initialize saturation + double xqs = xq; + double xds = xd; + if(hasSaturation) { + double oldDelta = 0; + bool exit = false; + int numIt = 0; + while(!exit) { + oldDelta = delta; + ABCtoDQ0(ia, delta, id0, iq0); + ABCtoDQ0(vt, delta, vd0, vq0); + + // Direct-axis Potier voltage. + double epd = vd0 + ra * id0 + xp * iq0; + + sq = 1.0 + satF * (xq / xd) * std::pow(epd, 6); + xqs = (xq - xp) / sq + xp; + eq0 = data.terminalVoltage + std::complex(ra, xqs) * ia; + delta = std::arg(eq0); + if(std::abs(delta - oldDelta) < m_saturationTolerance) { + exit = true; + } else if(numIt >= m_maxIterations) { + m_errorMsg = _("Error on initializate the saturation values of \"") + data.name + _("\"."); + return false; } - // Quadrature-axis Potier voltage. - double epq = vq0 + ra * iq0 - xp * id0; - sd = 1.0 + satF * std::pow(epq, 6); - xds = (xd - xp) / sd + xp; + numIt++; } + // Quadrature-axis Potier voltage. + double epq = vq0 + ra * iq0 - xp * id0; + sd = 1.0 + satF * std::pow(epq, 6); + xds = (xd - xp) / sd + xp; + } - double ef0 = vq0 + ra * iq0 - xds * id0; + double ef0 = vq0 + ra * iq0 - xds * id0; - data.initialFieldVoltage = ef0 * sd; - data.fieldVoltage = data.initialFieldVoltage; - data.pm = std::real((data.terminalVoltage * std::conj(ia)) + (std::abs(ia) * std::abs(ia) * ra)); - syncGenerator->IsOnline() ? data.speed = 2.0 * M_PI * m_systemFreq : data.speed = 2.0 * M_PI * data.ocFrequency; - data.delta = delta; - data.pe = data.pm; - data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); - data.sd = sd; - data.sq = sq; - data.id = id0; - data.iq = iq0; - - // Variables to extrapolate. - data.oldIq = iq0; - data.oldId = id0; - data.oldPe = data.pe; - data.oldSd = sd; - data.oldSq = sq; - - switch(data.model) { - case Machines::SM_MODEL_1: { - data.tranEq = std::abs(eq0); - - data.tranEd = 0.0; - data.subEq = 0.0; - data.subEd = 0.0; - } break; - case Machines::SM_MODEL_2: { - double tranXd = data.transXd * k; - - data.tranEq = ef0 + (xd - tranXd) * (id0 / sd); - data.tranEd = 0.0; - data.subEd = 0.0; - data.subEq = 0.0; - } break; - case Machines::SM_MODEL_3: { - double tranXd = data.transXd * k; - double tranXq = data.transXq * k; - if(tranXq == 0.0) tranXq = tranXd; - - data.tranEq = ef0 + (xd - tranXd) * (id0 / sd); - data.tranEd = -(xq - tranXq) * (iq0 / sq); - - data.subEd = 0.0; - data.subEq = 0.0; - } break; - case Machines::SM_MODEL_4: { - 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 = ef0 + (xd - tranXd) * (id0 / sd); - data.tranEd = 0.0; - data.subEq = data.tranEq + (tranXd - subXd) * (id0 / sd); - data.subEd = -(xq - subXq) * (iq0 / sq); - } break; - case Machines::SM_MODEL_5: { - double tranXd = data.transXd * k; - 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 = ef0 + (xd - tranXd) * (id0 / sd); - data.tranEd = -(xq - tranXq) * (iq0 / sq); - data.subEq = data.tranEq + (tranXd - subXd) * (id0 / sd); - data.subEd = data.tranEd - (tranXq - subXq) * (iq0 / sq); - } break; - default: - break; - } + data.initialFieldVoltage = ef0 * sd; + data.fieldVoltage = data.initialFieldVoltage; + data.pm = std::real((data.terminalVoltage * std::conj(ia)) + (std::abs(ia) * std::abs(ia) * ra)); + syncGenerator->IsOnline() ? data.speed = 2.0 * M_PI* m_systemFreq : data.speed = 2.0 * M_PI * data.ocFrequency; + data.delta = delta; + data.pe = data.pm; + data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); + data.sd = sd; + data.sq = sq; + data.id = id0; + data.iq = iq0; - // Initialize controllers - if(data.useAVR) { - if(data.avrSolver) delete data.avrSolver; - data.avrSolver = - new ControlElementSolver(data.avr, m_timeStep * m_ctrlTimeStepMultiplier, m_tolerance, m_parent); - data.avrSolver->SetSwitchStatus(syncGenerator->IsOnline()); - data.avrSolver->SetCurrentTime(m_currentTime); - data.avrSolver->SetTerminalVoltage(std::abs(data.terminalVoltage)); - data.avrSolver->SetInitialTerminalVoltage(std::abs(data.terminalVoltage)); - data.avrSolver->SetActivePower(dataPU.activePower); - data.avrSolver->SetReactivePower(dataPU.reactivePower); - data.avrSolver->SetVelocity(data.speed); - data.avrSolver->SetInitialVelocity(data.speed); - data.avrSolver->InitializeValues(false); - if(!data.avrSolver->IsOK()) { - m_errorMsg = _("Error on initializate the AVR of \"") + data.name + wxT("\".\n") + - data.avrSolver->GetErrorMessage(); - syncGenerator->SetElectricalData(data); - return false; - } + // Variables to extrapolate. + data.oldIq = iq0; + data.oldId = id0; + data.oldPe = data.pe; + data.oldSd = sd; + data.oldSq = sq; + + switch(data.model) { + case Machines::SM_MODEL_1: { + data.tranEq = std::abs(eq0); + + data.tranEd = 0.0; + data.subEq = 0.0; + data.subEd = 0.0; + } break; + case Machines::SM_MODEL_2: { + double tranXd = data.transXd * k; + + data.tranEq = ef0 + (xd - tranXd) * (id0 / sd); + data.tranEd = 0.0; + data.subEd = 0.0; + data.subEq = 0.0; + } break; + case Machines::SM_MODEL_3: { + double tranXd = data.transXd * k; + double tranXq = data.transXq * k; + if(tranXq == 0.0) tranXq = tranXd; + + data.tranEq = ef0 + (xd - tranXd) * (id0 / sd); + data.tranEd = -(xq - tranXq) * (iq0 / sq); + + data.subEd = 0.0; + data.subEq = 0.0; + } break; + case Machines::SM_MODEL_4: { + 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 = ef0 + (xd - tranXd) * (id0 / sd); + data.tranEd = 0.0; + data.subEq = data.tranEq + (tranXd - subXd) * (id0 / sd); + data.subEd = -(xq - subXq) * (iq0 / sq); + } break; + case Machines::SM_MODEL_5: { + double tranXd = data.transXd * k; + 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 = ef0 + (xd - tranXd) * (id0 / sd); + data.tranEd = -(xq - tranXq) * (iq0 / sq); + data.subEq = data.tranEq + (tranXd - subXd) * (id0 / sd); + data.subEd = data.tranEd - (tranXq - subXq) * (iq0 / sq); + } break; + default: + break; + } + + // Initialize controllers + if(data.useAVR) { + if(data.avrSolver) delete data.avrSolver; + data.avrSolver = + new ControlElementSolver(data.avr, m_timeStep * m_ctrlTimeStepMultiplier, m_tolerance, m_parent); + data.avrSolver->SetSwitchStatus(syncGenerator->IsOnline()); + data.avrSolver->SetCurrentTime(m_currentTime); + data.avrSolver->SetTerminalVoltage(std::abs(data.terminalVoltage)); + data.avrSolver->SetInitialTerminalVoltage(std::abs(data.terminalVoltage)); + data.avrSolver->SetActivePower(dataPU.activePower); + data.avrSolver->SetReactivePower(dataPU.reactivePower); + data.avrSolver->SetVelocity(data.speed); + data.avrSolver->SetInitialVelocity(data.speed); + data.avrSolver->InitializeValues(false); + if(!data.avrSolver->IsOK()) { + m_errorMsg = _("Error on initializate the AVR of \"") + data.name + wxT("\".\n") + + data.avrSolver->GetErrorMessage(); + syncGenerator->SetElectricalData(data); + return false; } - if(data.useSpeedGovernor) { - if(data.speedGovSolver) delete data.speedGovSolver; - data.speedGovSolver = new ControlElementSolver(data.speedGov, m_timeStep * m_ctrlTimeStepMultiplier, - m_tolerance, m_parent); - data.speedGovSolver->SetSwitchStatus(syncGenerator->IsOnline()); - data.speedGovSolver->SetCurrentTime(m_currentTime); - data.speedGovSolver->SetActivePower(dataPU.activePower); - data.speedGovSolver->SetReactivePower(dataPU.reactivePower); - data.speedGovSolver->SetVelocity(data.speed); - data.speedGovSolver->SetInitialVelocity(data.speed); - data.speedGovSolver->SetInitialMecPower(data.pm); - data.speedGovSolver->InitializeValues(false); - if(!data.speedGovSolver->IsOK()) { - m_errorMsg = _("Error on initializate the speed governor of \"") + data.name + wxT("\".\n") + - data.speedGovSolver->GetErrorMessage(); - syncGenerator->SetElectricalData(data); - return false; - } + } + if(data.useSpeedGovernor) { + if(data.speedGovSolver) delete data.speedGovSolver; + data.speedGovSolver = + new ControlElementSolver(data.speedGov, m_timeStep * m_ctrlTimeStepMultiplier, m_tolerance, m_parent); + data.speedGovSolver->SetSwitchStatus(syncGenerator->IsOnline()); + data.speedGovSolver->SetCurrentTime(m_currentTime); + data.speedGovSolver->SetActivePower(dataPU.activePower); + data.speedGovSolver->SetReactivePower(dataPU.reactivePower); + data.speedGovSolver->SetVelocity(data.speed); + data.speedGovSolver->SetInitialVelocity(data.speed); + data.speedGovSolver->SetInitialMecPower(data.pm); + data.speedGovSolver->InitializeValues(false); + if(!data.speedGovSolver->IsOK()) { + m_errorMsg = _("Error on initializate the speed governor of \"") + data.name + wxT("\".\n") + + data.speedGovSolver->GetErrorMessage(); + syncGenerator->SetElectricalData(data); + return false; } + } //} else { - // Initialize open circuit machine. + // Initialize open circuit machine. //} // Reset plot data data.terminalVoltageVector.clear(); @@ -721,6 +724,40 @@ bool Electromechanical::InitializeDynamicElements() syncGenerator->SetElectricalData(data); } + // Induction motors + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* indMotor = *it; + auto dataPU = indMotor->GetPUElectricalData(m_powerSystemBase); + auto data = indMotor->GetElectricalData(); + + double w0 = 2.0 * M_PI * m_systemFreq; + + std::complex i1 = std::complex(dataPU.activePower, -data.q0) / std::conj(data.terminalVoltage); + data.ir = std::real(i1); + data.im = std::imag(i1); + std::complex e = data.terminalVoltage - std::complex(data.r1t, data.x1t) * i1; + data.te = std::real(e * std::conj(i1)) / w0; + + double wi = w0 * (1 - data.s0); // Initial velocity + data.as = data.te * (data.aw + (data.bw * w0) / wi + (data.cw * w0 * w0) / (wi * wi)); + data.bs = data.te * ((data.bw * w0) / wi + (2.0 * data.cw * w0 * w0) / (wi * wi)); + data.cs = (data.te * data.cw * w0 * w0) / (wi * wi); + + std::complex tranE = + (std::complex(0, data.x0 - data.xt) * i1) / std::complex(1.0, w0 * data.s0 * data.t0); + + data.tranEr = std::real(tranE); + data.tranEm = std::imag(tranE); + + /* + double pEr = w0 * data.s0 * data.tranEm - (data.tranEr + (data.x0 - data.xt) * data.im) / (data.t0); + double pEm = -w0 * data.s0 * data.tranEr - (data.tranEm - (data.x0 - data.xt) * data.ir) / (data.t0); + wxMessageBox(wxString::Format("%.2f\%\n%f\n%f\n%f\n\n%f\n\n%f\n%f\n\n%f\n%f", data.s0 * 100, std::abs(i1), + std::abs(e), data.te, data.as - data.bs * data.s0 + data.cs * data.s0 * data.s0, + data.tranEr, data.tranEm, pEr, pEm)); + */ + indMotor->SetElectricalData(data); + } CalculateReferenceSpeed(); return true; } @@ -806,6 +843,17 @@ bool Electromechanical::CalculateInjectedCurrents() syncGenerator->SetElectricalData(data); } + // Induction motors + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* motor = *it; + if(motor->IsOnline()) { + auto data = motor->GetElectricalData(); + int n = static_cast(motor->GetParentList()[0])->GetElectricalData().number; + std::complex y0 = std::complex(1.0, 0.0) / std::complex(data.r1t, data.xt); + std::complex e = std::complex(data.tranEr, data.tranEm); + m_iBus[n] += y0 * e; + } + } // Loads for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) { @@ -1052,9 +1100,7 @@ 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) { - syncGenerator->SavePlotData(); - } + if(data.plotSyncMachine) { syncGenerator->SavePlotData(); } } for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) { Bus* bus = *it; @@ -1457,3 +1503,90 @@ void Electromechanical::PreallocateVectors() } } } + +std::complex Electromechanical::GetIndMachineAdmittance(IndMotor* motor) +{ + auto data = motor->GetElectricalData(); + auto dataPU = motor->GetPUElectricalData(m_powerSystemBase); + std::complex y0 = (std::complex(1, 0) / std::complex(data.r1t, data.xt)); + // The difference between calculated and defined reactive power + std::complex yq = std::complex(0.0, data.q0 - dataPU.reactivePower); + return y0 + yq; +} + +bool Electromechanical::InsertIndMachinesOnYBus() +{ + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* motor = *it; + if(motor->IsOnline()) { + if(!CalculateIndMachinesTransientValues(motor)) return false; + int n = static_cast(motor->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] += GetIndMachineAdmittance(motor); + } + } + return true; +} + +bool Electromechanical::CalculateIndMachinesTransientValues(IndMotor* motor) +{ + auto data = motor->GetElectricalData(); + auto dataPU = motor->GetPUElectricalData(m_powerSystemBase); + double k = 1.0; // Power base change factor. + if(data.useMachinePowerAsBase) { + double oldBase = motor->GetValueFromUnit(data.ratedPower, data.ratedPowerUnit); + k = m_powerSystemBase / oldBase; + } + + // Calculate the induction machine transient constants at the machine base + double r1t = data.r1 * k; + double r2t = data.r2 * k; + double x1t = data.x1 * k; + double x2t = data.x2 * k; + double xmt = data.xm * k; + data.r1t = r1t; + data.r2t = r2t; + data.x1t = x1t; + data.x2t = x2t; + data.xmt = xmt; + + double xt = x1t + (x2t * xmt) / (x2t + xmt); + double t0 = (x2t + xmt) / (2.0 * M_PI * m_systemFreq * r2t); + double x0 = x1t + xmt; + data.xt = xt; + data.t0 = t0; + data.x0 = x0; + + data.terminalVoltage = static_cast(motor->GetParentList()[0])->GetElectricalData().voltage; + + //[Ref.] Induction Motor Static Models for Power Flow and Voltage Stability Studies + double p = dataPU.activePower; + double v = std::abs(data.terminalVoltage); + double r1 = data.r1t; + double r2 = data.r2t; + double x1 = data.x1t; + double x2 = data.x2t; + double xm = data.xmt; + double k1 = x2 + xm; + double k2 = -x1 * k1 - x2 * xm; + double k3 = xm + x1; + double k4 = r1 * k1; + double a = p * (r1 * r1 + k3 * k3) - v * v * r1; + double b = 2.0 * p * (r1 * k2 + k3 * k4) - v * v * (k2 + k1 * k3); + double c = p * (k2 * k2 + k4 * k4) - v * v * k1 * k4; + double d = (b * b - 4.0 * a * c); + if(d < 0.0) { + m_errorMsg = _("Error on initializate the slip of \"") + data.name + _("\"."); + return false; + } + double r2_s = (-b + std::sqrt(d)) / (2.0 * a); + data.s0 = r2 / r2_s; + + double qa = k1 * (r2_s * r1 - x1 * k1 - x2 * xm); + double qb = r2_s * (r2_s * (xm + x1) + r1 * k1); + double qc = r2_s * r1 - x1 * k1 - x2 * xm; + double qd = r2_s * (xm + x1) + r1 * k1; + data.q0 = (-v * v * (qa - qb)) / (qc * qc + qd * qd); + + motor->SetElectricalData(data); + return true; +} -- cgit From ec20e2ac95f8a5013f2293a7a2a96016fb805581 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 11 Jul 2019 18:06:08 -0300 Subject: Induction motor DAE implemented Need some testing Check inertia --- Project/Electromechanical.cpp | 250 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 219 insertions(+), 31 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 0f32ee8..1a8a0d0 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -125,7 +125,7 @@ bool Electromechanical::RunStabilityCalculation() currentPbdTime = 0.0; } - if(!SolveSynchronousMachines()) return false; + if(!SolveMachines()) return false; m_currentTime += m_timeStep; currentPlotTime += m_timeStep; @@ -749,13 +749,16 @@ bool Electromechanical::InitializeDynamicElements() data.tranEr = std::real(tranE); data.tranEm = std::imag(tranE); - /* - double pEr = w0 * data.s0 * data.tranEm - (data.tranEr + (data.x0 - data.xt) * data.im) / (data.t0); - double pEm = -w0 * data.s0 * data.tranEr - (data.tranEm - (data.x0 - data.xt) * data.ir) / (data.t0); - wxMessageBox(wxString::Format("%.2f\%\n%f\n%f\n%f\n\n%f\n\n%f\n%f\n\n%f\n%f", data.s0 * 100, std::abs(i1), - std::abs(e), data.te, data.as - data.bs * data.s0 + data.cs * data.s0 * data.s0, - data.tranEr, data.tranEm, pEr, pEm)); - */ + data.slip = data.s0; + + // Reset plot data + data.slipVector.clear(); + data.slipVector.shrink_to_fit(); + data.electricalTorqueVector.clear(); + data.electricalTorqueVector.shrink_to_fit(); + data.mechanicalTorqueVector.clear(); + data.mechanicalTorqueVector.shrink_to_fit(); + indMotor->SetElectricalData(data); } CalculateReferenceSpeed(); @@ -846,13 +849,31 @@ bool Electromechanical::CalculateInjectedCurrents() // Induction motors for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { IndMotor* motor = *it; + auto data = motor->GetElectricalData(); if(motor->IsOnline()) { - auto data = motor->GetElectricalData(); int n = static_cast(motor->GetParentList()[0])->GetElectricalData().number; std::complex y0 = std::complex(1.0, 0.0) / std::complex(data.r1t, data.xt); + std::complex v = m_vBus[n]; std::complex e = std::complex(data.tranEr, data.tranEm); - m_iBus[n] += y0 * e; + std::complex iInj = y0 * e; + m_iBus[n] += iInj; + + std::complex iMachine = y0 * (v - e); + + // wxMessageBox(wxString::Format("%f\n%f", data.ir, data.im)); + + data.ir = std::real(iMachine); + data.im = std::imag(iMachine); + // double y = (1.0 / (data.r1t * data.r1t + data.xt * data.xt)); + // data.ir = y * (data.r1t * data.tranEr + data.xt * data.tranEm); + // data.im = y * (-data.xt * data.tranEr + data.r1t * data.tranEm); + + // wxMessageBox(wxString::Format("%f\n%f", data.ir, data.im)); + } else { + data.ir = 0.0; + data.im = 0.0; } + motor->SetElectricalData(data); } // Loads @@ -972,8 +993,37 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat syncGenerator->SetElectricalData(data); } -bool Electromechanical::SolveSynchronousMachines() +void Electromechanical::CalculateIntegrationConstants(IndMotor* indMotor, double ir, double im, double k) { + double w0 = 2.0 * M_PI * m_systemFreq; + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* motor = *it; + auto data = motor->GetElectricalData(); + + // Mechanical equations + // s + data.icSlip.m = m_timeStep / (4.0f * data.inertia / k + data.bs * m_timeStep); + data.icSlip.c = data.slip * (1.0 - 2.0 * data.bs * data.icSlip.m) + + data.icSlip.m * (2.0 * data.as + data.cs * data.slip * data.slip - data.te); + + // Electrical equations + // Er + data.icTranEr.m = m_timeStep / (2.0 * data.t0 + m_timeStep); + data.icTranEr.c = data.tranEr * (1.0 - 2.0 * data.icTranEr.m) + + data.icTranEr.m * (w0 * data.t0 * data.slip * data.tranEm - (data.x0 - data.xt) * im); + // Em + data.icTranEm.m = m_timeStep / (2.0 * data.t0 + m_timeStep); + data.icTranEm.c = data.tranEm * (1.0 - 2.0 * data.icTranEm.m) - + data.icTranEm.m * (w0 * data.t0 * data.slip * data.tranEr - (data.x0 - data.xt) * ir); + + motor->SetElectricalData(data); + } +} + +bool Electromechanical::SolveMachines() +{ + // First interation values and extrapolations + // Synchronous machines for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; auto data = syncGenerator->GetElectricalData(); @@ -995,7 +1045,7 @@ bool Electromechanical::SolveSynchronousMachines() // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq, k); - if(!CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k)) return false; + if(!CalculateNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k)) return false; // Extrapolate nonintegrable variables. id = 2.0 * id - data.oldId; iq = 2.0 * iq - data.oldIq; @@ -1003,11 +1053,42 @@ bool Electromechanical::SolveSynchronousMachines() sd = 2.0 * sd - data.oldSd; sq = 2.0 * sq - data.oldSq; - CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); + CalculateIntVariables(syncGenerator, id, iq, sd, sq, pe, k); } else { CalculateIntegrationConstants(syncGenerator, 0.0f, 0.0f); } } + // Induction machines + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* indMotor = *it; + auto data = indMotor->GetElectricalData(); + + if(indMotor->IsOnline()) { + double ir, im, te; + te = data.te; + ir = data.ir; + im = data.im; + + double k = 1.0; // Power base change factor. + if(data.useMachinePowerAsBase) { + double oldBase = indMotor->GetValueFromUnit(data.ratedPower, data.ratedPowerUnit); + k = m_powerSystemBase / oldBase; + } + + // Calculate integration constants. + CalculateIntegrationConstants(indMotor, ir, im, k); + + if(!CalculateNonIntVariables(indMotor, ir, im, te, k)) return false; + // Extrapolate nonintegrable variables. + // ir = 2.0 * ir - data.oldIr; + // im = 2.0 * im - data.oldIm; + // te = 2.0 * te - data.oldTe; + + CalculateIntVariables(indMotor, ir, im, te, k); + } else { + CalculateIntegrationConstants(indMotor, 0.0f, 0.0f); + } + } double error = 1.0; int iterations = 0; @@ -1020,7 +1101,7 @@ bool Electromechanical::SolveSynchronousMachines() // Calculate the buses voltages. m_vBus = LUEvaluate(m_yBusU, m_yBusL, m_iBus); - // Solve machine equations. + // Solve synchronous machine equations. for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; @@ -1040,17 +1121,42 @@ bool Electromechanical::SolveSynchronousMachines() } // Calculate id, iq, dq, sd - if(!CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k)) return false; + if(!CalculateNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k)) return false; - double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); + double genError = CalculateIntVariables(syncGenerator, id, iq, sd, sq, pe, k); if(genError > error) error = genError; } + // Solve induction machine equation + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* indMotor = *it; + + auto data = indMotor->GetElectricalData(); + + double ir = data.ir; + double im = data.im; + double te = data.te; + + // Power base change factor. + double k = 1.0; + if(data.useMachinePowerAsBase) { + double oldBase = indMotor->GetValueFromUnit(data.ratedPower, data.ratedPowerUnit); + k = m_powerSystemBase / oldBase; + } + + // Calculate te + if(!CalculateNonIntVariables(indMotor, ir, im, te, k)) return false; + + double motorError = CalculateIntVariables(indMotor, ir, im, te, k); + + if(motorError > error) error = motorError; + } + ++iterations; if(iterations > m_maxIterations) { - m_errorMsg = _("Impossible to solve the synchronous generators.\nCheck the system parameters and/or " + m_errorMsg = _("Impossible to solve the system machines.\nCheck the system parameters and/or " "decrease the time step."); return false; } @@ -1119,6 +1225,19 @@ void Electromechanical::SaveData() load->SetElectricalData(data); } } + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* motor = *it; + auto data = motor->GetElectricalData(); + if(data.plotIndMachine) { + data.slipVector.emplace_back(data.slip * 100.0); + double te = std::real(std::complex(data.tranEr, data.tranEm) * std::complex(data.ir, -data.im)); + data.electricalTorqueVector.emplace_back(te); + //data.electricalTorqueVector.emplace_back(data.te * 2.0 * M_PI * m_systemFreq); + double tm = (data.as - data.bs * data.slip + data.cs * data.slip * data.slip) * 2.0 * M_PI * m_systemFreq; + data.mechanicalTorqueVector.emplace_back(tm); + motor->SetElectricalData(data); + } + } m_iterationsNumVector.emplace_back(m_iterationsNum); } @@ -1132,13 +1251,13 @@ void Electromechanical::SetSyncMachinesModel() } } -bool Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncGenerator, - double& id, - double& iq, - double& sd, - double& sq, - double& pe, - double k) +bool Electromechanical::CalculateNonIntVariables(SyncGenerator* syncGenerator, + double& id, + double& iq, + double& sd, + double& sq, + double& pe, + double k) { auto data = syncGenerator->GetElectricalData(); int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; @@ -1172,13 +1291,33 @@ bool Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG return true; } -double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGenerator, - double id, - double iq, - double sd, - double sq, - double pe, - double k) +bool Electromechanical::CalculateNonIntVariables(IndMotor* indMotor, double& ir, double& im, double& te, double k) +{ + auto data = indMotor->GetElectricalData(); + if(indMotor->IsOnline()) { + double w0 = 2.0 * M_PI * m_systemFreq; + te = (data.tranEr * ir + data.tranEm * im) / w0; + } else { + te = ir = im = 0.0; + } + data.te = te; + data.ir = ir; + data.im = im; + data.oldTe = te; + data.oldIr = ir; + data.oldIm = im; + indMotor->SetElectricalData(data); + + return true; +} + +double Electromechanical::CalculateIntVariables(SyncGenerator* syncGenerator, + double id, + double iq, + double sd, + double sq, + double pe, + double k) { double error = 0.0; auto data = syncGenerator->GetElectricalData(); @@ -1302,6 +1441,46 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe return error; } +double Electromechanical::CalculateIntVariables(IndMotor* indMotor, double ir, double im, double te, double k) +{ + double error = 0.0; + auto data = indMotor->GetElectricalData(); + double w0 = 2.0 * M_PI * m_systemFreq; + + // Mechanical differential equations + // Using Newton method to solve the non-linear slip equation: s = Cs + Ms * (C * s^2 - Te): + double slip = data.slip; // Initial value + double ds = (data.icSlip.c + data.icSlip.m * (data.cs * slip * slip - te) - slip) / + (1.0 - 2.0 * data.icSlip.m * data.cs * slip * slip); + int iteration = 0; + while(std::abs(ds) > 1e-8) { + slip += ds; + ds = (data.icSlip.c + data.icSlip.m * (data.cs * slip * slip - te) - slip) / + (1.0 - 2.0 * data.icSlip.m * data.cs * slip * slip); + iteration++; + if(iteration > m_maxIterations) break; + } + + error = std::max(error, std::abs(data.slip - slip)); + data.slip = slip; + + // Electrical differential equations + + double tranEr = data.icTranEr.c + data.icTranEr.m * (w0 * data.t0 * slip * data.tranEm - (data.x0 - data.xt) * im); + error = std::max(error, std::abs(data.tranEr - tranEr)); + + double tranEm = data.icTranEm.c - data.icTranEm.m * (w0 * data.t0 * slip * data.tranEr - (data.x0 - data.xt) * ir); + error = std::max(error, std::abs(data.tranEm - tranEm)); + + // wxMessageBox(wxString::Format("%f\n%f\n\n%f\n%f", data.tranEr, data.tranEm, tranEr, tranEm)); + + data.tranEr = tranEr; + data.tranEm = tranEm; + + indMotor->SetElectricalData(data); + return error; +} + void Electromechanical::CalculateReferenceSpeed() { if(m_useCOI) { @@ -1502,6 +1681,15 @@ void Electromechanical::PreallocateVectors() load->SetElectricalData(data); } } + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* motor = *it; + auto data = motor->GetElectricalData(); + if(data.plotIndMachine) { + data.slipVector.reserve(numPoints); + data.electricalTorqueVector.reserve(numPoints); + motor->SetElectricalData(data); + } + } } std::complex Electromechanical::GetIndMachineAdmittance(IndMotor* motor) -- cgit From fed838330b39f69a682b64b1b885fec17ff901d7 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Fri, 12 Jul 2019 02:03:28 -0300 Subject: Motor start implemented and bugfixes --- Project/Electromechanical.cpp | 108 +++++++++++++++++++++++++++++------------- 1 file changed, 74 insertions(+), 34 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 1a8a0d0..64fdb9a 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -224,6 +224,32 @@ void Electromechanical::SetEvent(double currentTime) } } + // Induction motor switching + for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { + IndMotor* motor = *it; + auto swData = motor->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 && motor->IsOnline() && motor->GetParentList().size() == 1) { + auto data = motor->GetElectricalData(); + motor->SetOnline(false); + int n = static_cast(motor->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] -= (std::complex(1, 0) / std::complex(data.r1t, data.xt)); + } + + // Insert machine (only disconnected machines) + if(swData.swType[i] == SW_INSERT && !motor->IsOnline() && motor->GetParentList().size() == 1) { + auto data = motor->GetElectricalData(); + if(motor->SetOnline(true)) { + int n = static_cast(motor->GetParentList()[0])->GetElectricalData().number; + m_yBus[n][n] += (std::complex(1, 0) / std::complex(data.r1t, data.xt)); + } + } + } + } + } + // Load switching for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) { Load* load = *it; @@ -231,7 +257,7 @@ void Electromechanical::SetEvent(double currentTime) 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()) { + if(swData.swType[i] == SW_REMOVE && load->IsOnline() && load->GetParentList().size() == 1) { load->SetOnline(false); auto data = load->GetPUElectricalData(m_powerSystemBase); Bus* parentBus = static_cast(load->GetParentList()[0]); @@ -731,25 +757,36 @@ bool Electromechanical::InitializeDynamicElements() auto data = indMotor->GetElectricalData(); double w0 = 2.0 * M_PI * m_systemFreq; - std::complex i1 = std::complex(dataPU.activePower, -data.q0) / std::conj(data.terminalVoltage); - data.ir = std::real(i1); - data.im = std::imag(i1); + double ir = std::real(i1); + double im = std::imag(i1); std::complex e = data.terminalVoltage - std::complex(data.r1t, data.x1t) * i1; - data.te = std::real(e * std::conj(i1)) / w0; + double te = std::real(e * std::conj(i1)) / w0; double wi = w0 * (1 - data.s0); // Initial velocity - data.as = data.te * (data.aw + (data.bw * w0) / wi + (data.cw * w0 * w0) / (wi * wi)); - data.bs = data.te * ((data.bw * w0) / wi + (2.0 * data.cw * w0 * w0) / (wi * wi)); - data.cs = (data.te * data.cw * w0 * w0) / (wi * wi); + data.as = te * (data.aw + (data.bw * w0) / wi + (data.cw * w0 * w0) / (wi * wi)); + data.bs = te * ((data.bw * w0) / wi + (2.0 * data.cw * w0 * w0) / (wi * wi)); + data.cs = (te * data.cw * w0 * w0) / (wi * wi); - std::complex tranE = - (std::complex(0, data.x0 - data.xt) * i1) / std::complex(1.0, w0 * data.s0 * data.t0); - - data.tranEr = std::real(tranE); - data.tranEm = std::imag(tranE); - - data.slip = data.s0; + if(indMotor->IsOnline()) { + std::complex tranE = + (std::complex(0, data.x0 - data.xt) * i1) / std::complex(1.0, w0 * data.s0 * data.t0); + + data.tranEr = std::real(tranE); + data.tranEm = std::imag(tranE); + + data.slip = data.s0; + data.ir = ir; + data.im = im; + data.te = te; + } else { // Offline + data.tranEr = 0.0; + data.tranEm = 0.0; + data.slip = 1.0 - 1e-7; + data.ir = 0.0; + data.im = 0.0; + data.te = 0.0; + } // Reset plot data data.slipVector.clear(); @@ -757,7 +794,11 @@ bool Electromechanical::InitializeDynamicElements() data.electricalTorqueVector.clear(); data.electricalTorqueVector.shrink_to_fit(); data.mechanicalTorqueVector.clear(); - data.mechanicalTorqueVector.shrink_to_fit(); + data.mechanicalTorqueVector.shrink_to_fit(); + data.velocityVector.clear(); + data.velocityVector.shrink_to_fit(); + data.currentVector.clear(); + data.currentVector.shrink_to_fit(); indMotor->SetElectricalData(data); } @@ -860,15 +901,8 @@ bool Electromechanical::CalculateInjectedCurrents() std::complex iMachine = y0 * (v - e); - // wxMessageBox(wxString::Format("%f\n%f", data.ir, data.im)); - data.ir = std::real(iMachine); data.im = std::imag(iMachine); - // double y = (1.0 / (data.r1t * data.r1t + data.xt * data.xt)); - // data.ir = y * (data.r1t * data.tranEr + data.xt * data.tranEm); - // data.im = y * (-data.xt * data.tranEr + data.r1t * data.tranEm); - - // wxMessageBox(wxString::Format("%f\n%f", data.ir, data.im)); } else { data.ir = 0.0; data.im = 0.0; @@ -1002,7 +1036,7 @@ void Electromechanical::CalculateIntegrationConstants(IndMotor* indMotor, double // Mechanical equations // s - data.icSlip.m = m_timeStep / (4.0f * data.inertia / k + data.bs * m_timeStep); + data.icSlip.m = m_timeStep / ((4.0f * data.inertia / w0) / k + data.bs * m_timeStep); data.icSlip.c = data.slip * (1.0 - 2.0 * data.bs * data.icSlip.m) + data.icSlip.m * (2.0 * data.as + data.cs * data.slip * data.slip - data.te); @@ -1080,9 +1114,9 @@ bool Electromechanical::SolveMachines() if(!CalculateNonIntVariables(indMotor, ir, im, te, k)) return false; // Extrapolate nonintegrable variables. - // ir = 2.0 * ir - data.oldIr; - // im = 2.0 * im - data.oldIm; - // te = 2.0 * te - data.oldTe; + ir = 2.0 * ir - data.oldIr; + im = 2.0 * im - data.oldIm; + te = 2.0 * te - data.oldTe; CalculateIntVariables(indMotor, ir, im, te, k); } else { @@ -1230,11 +1264,13 @@ void Electromechanical::SaveData() auto data = motor->GetElectricalData(); if(data.plotIndMachine) { data.slipVector.emplace_back(data.slip * 100.0); - double te = std::real(std::complex(data.tranEr, data.tranEm) * std::complex(data.ir, -data.im)); - data.electricalTorqueVector.emplace_back(te); - //data.electricalTorqueVector.emplace_back(data.te * 2.0 * M_PI * m_systemFreq); + data.electricalTorqueVector.emplace_back(data.te * 2.0 * M_PI * m_systemFreq); double tm = (data.as - data.bs * data.slip + data.cs * data.slip * data.slip) * 2.0 * M_PI * m_systemFreq; data.mechanicalTorqueVector.emplace_back(tm); + double w = (1.0 - data.slip) * 2.0 * M_PI * m_systemFreq; + data.velocityVector.emplace_back(w); + double i1 = std::sqrt(data.ir * data.ir + data.im * data.im); + data.currentVector.emplace_back(i1); motor->SetElectricalData(data); } } @@ -1449,7 +1485,7 @@ double Electromechanical::CalculateIntVariables(IndMotor* indMotor, double ir, d // Mechanical differential equations // Using Newton method to solve the non-linear slip equation: s = Cs + Ms * (C * s^2 - Te): - double slip = data.slip; // Initial value + double slip = data.slip; // Initial value. CAN BE THE PROBLEM ON MOTOR START! double ds = (data.icSlip.c + data.icSlip.m * (data.cs * slip * slip - te) - slip) / (1.0 - 2.0 * data.icSlip.m * data.cs * slip * slip); int iteration = 0; @@ -1460,7 +1496,8 @@ double Electromechanical::CalculateIntVariables(IndMotor* indMotor, double ir, d iteration++; if(iteration > m_maxIterations) break; } - + + if(!indMotor->IsOnline()) slip = 1.0 - 1e-7; error = std::max(error, std::abs(data.slip - slip)); data.slip = slip; @@ -1696,9 +1733,10 @@ std::complex Electromechanical::GetIndMachineAdmittance(IndMotor* motor) { auto data = motor->GetElectricalData(); auto dataPU = motor->GetPUElectricalData(m_powerSystemBase); + std::complex v = static_cast(motor->GetParentList()[0])->GetElectricalData().voltage; std::complex y0 = (std::complex(1, 0) / std::complex(data.r1t, data.xt)); // The difference between calculated and defined reactive power - std::complex yq = std::complex(0.0, data.q0 - dataPU.reactivePower); + std::complex yq = std::complex(0.0, data.q0 - dataPU.reactivePower) / (std::abs(v) * std::abs(v)); return y0 + yq; } @@ -1706,8 +1744,8 @@ bool Electromechanical::InsertIndMachinesOnYBus() { for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { IndMotor* motor = *it; + if(!CalculateIndMachinesTransientValues(motor)) return false; if(motor->IsOnline()) { - if(!CalculateIndMachinesTransientValues(motor)) return false; int n = static_cast(motor->GetParentList()[0])->GetElectricalData().number; m_yBus[n][n] += GetIndMachineAdmittance(motor); } @@ -1749,6 +1787,8 @@ bool Electromechanical::CalculateIndMachinesTransientValues(IndMotor* motor) //[Ref.] Induction Motor Static Models for Power Flow and Voltage Stability Studies double p = dataPU.activePower; double v = std::abs(data.terminalVoltage); + // If the motor is offline, calculate the nominal slip to user-defined power input and 1.0 p.u. voltage + if(!motor->IsOnline()) v = 1.0; double r1 = data.r1t; double r2 = data.r2t; double x1 = data.x1t; -- cgit From 4d4f39cb101bb0ccb1d336587a9e143b392dbefc Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Sat, 13 Jul 2019 01:14:14 -0300 Subject: Multiple motor initialization implemented --- Project/Electromechanical.cpp | 82 ++++++++++++++++++++++++++----------------- 1 file changed, 50 insertions(+), 32 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 64fdb9a..bd95fbe 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -788,6 +788,11 @@ bool Electromechanical::InitializeDynamicElements() data.te = 0.0; } + // Variables to extrapolate. + data.oldTe = data.te; + data.oldIr = data.ir; + data.oldIm = data.im; + // Reset plot data data.slipVector.clear(); data.slipVector.shrink_to_fit(); @@ -799,6 +804,12 @@ bool Electromechanical::InitializeDynamicElements() data.velocityVector.shrink_to_fit(); data.currentVector.clear(); data.currentVector.shrink_to_fit(); + data.terminalVoltageVector.clear(); + data.terminalVoltageVector.shrink_to_fit(); + data.activePowerVector.clear(); + data.activePowerVector.shrink_to_fit(); + data.reactivePowerVector.clear(); + data.reactivePowerVector.shrink_to_fit(); indMotor->SetElectricalData(data); } @@ -1030,28 +1041,26 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat void Electromechanical::CalculateIntegrationConstants(IndMotor* indMotor, double ir, double im, double k) { double w0 = 2.0 * M_PI * m_systemFreq; - for(auto it = m_indMotorList.begin(), itEnd = m_indMotorList.end(); it != itEnd; ++it) { - IndMotor* motor = *it; - auto data = motor->GetElectricalData(); - // Mechanical equations - // s - data.icSlip.m = m_timeStep / ((4.0f * data.inertia / w0) / k + data.bs * m_timeStep); - data.icSlip.c = data.slip * (1.0 - 2.0 * data.bs * data.icSlip.m) + - data.icSlip.m * (2.0 * data.as + data.cs * data.slip * data.slip - data.te); - - // Electrical equations - // Er - data.icTranEr.m = m_timeStep / (2.0 * data.t0 + m_timeStep); - data.icTranEr.c = data.tranEr * (1.0 - 2.0 * data.icTranEr.m) + - data.icTranEr.m * (w0 * data.t0 * data.slip * data.tranEm - (data.x0 - data.xt) * im); - // Em - data.icTranEm.m = m_timeStep / (2.0 * data.t0 + m_timeStep); - data.icTranEm.c = data.tranEm * (1.0 - 2.0 * data.icTranEm.m) - - data.icTranEm.m * (w0 * data.t0 * data.slip * data.tranEr - (data.x0 - data.xt) * ir); + auto data = indMotor->GetElectricalData(); - motor->SetElectricalData(data); - } + // Mechanical equations + // s + data.icSlip.m = m_timeStep / ((4.0f * data.inertia / w0) / k + data.bs * m_timeStep); + data.icSlip.c = data.slip * (1.0 - 2.0 * data.bs * data.icSlip.m) + + data.icSlip.m * (2.0 * data.as + data.cs * data.slip * data.slip - data.te); + + // Electrical equations + // Er + data.icTranEr.m = m_timeStep / (2.0 * data.t0 + m_timeStep); + data.icTranEr.c = data.tranEr * (1.0 - 2.0 * data.icTranEr.m) + + data.icTranEr.m * (w0 * data.t0 * data.slip * data.tranEm - (data.x0 - data.xt) * im); + // Em + data.icTranEm.m = m_timeStep / (2.0 * data.t0 + m_timeStep); + data.icTranEm.c = data.tranEm * (1.0 - 2.0 * data.icTranEm.m) - + data.icTranEm.m * (w0 * data.t0 * data.slip * data.tranEr - (data.x0 - data.xt) * ir); + + indMotor->SetElectricalData(data); } bool Electromechanical::SolveMachines() @@ -1269,8 +1278,14 @@ void Electromechanical::SaveData() data.mechanicalTorqueVector.emplace_back(tm); double w = (1.0 - data.slip) * 2.0 * M_PI * m_systemFreq; data.velocityVector.emplace_back(w); - double i1 = std::sqrt(data.ir * data.ir + data.im * data.im); - data.currentVector.emplace_back(i1); + std::complex i1 = std::complex(data.ir, data.im); + data.currentVector.emplace_back(std::abs(i1)); + int n = static_cast(motor->GetParentList()[0])->GetElectricalData().number; + std::complex v = m_vBus[n]; + data.terminalVoltageVector.emplace_back(std::abs(v)); + std::complex s = v * std::conj(i1); + data.activePowerVector.emplace_back(std::real(s)); + data.reactivePowerVector.emplace_back(std::imag(s)); motor->SetElectricalData(data); } } @@ -1485,7 +1500,7 @@ double Electromechanical::CalculateIntVariables(IndMotor* indMotor, double ir, d // Mechanical differential equations // Using Newton method to solve the non-linear slip equation: s = Cs + Ms * (C * s^2 - Te): - double slip = data.slip; // Initial value. CAN BE THE PROBLEM ON MOTOR START! + double slip = 0.0; // Initial value. CAN BE THE PROBLEM ON MOTOR START! double ds = (data.icSlip.c + data.icSlip.m * (data.cs * slip * slip - te) - slip) / (1.0 - 2.0 * data.icSlip.m * data.cs * slip * slip); int iteration = 0; @@ -1496,21 +1511,22 @@ double Electromechanical::CalculateIntVariables(IndMotor* indMotor, double ir, d iteration++; if(iteration > m_maxIterations) break; } - + if(!indMotor->IsOnline()) slip = 1.0 - 1e-7; error = std::max(error, std::abs(data.slip - slip)); data.slip = slip; - // Electrical differential equations + // Change T'0 with the cage factor + if(data.useKf) + data.t0 = (data.x2t + data.xmt) / (2.0 * M_PI * m_systemFreq * data.r2t * (1.0 + data.kf * data.r2t)); + // Electrical differential equations double tranEr = data.icTranEr.c + data.icTranEr.m * (w0 * data.t0 * slip * data.tranEm - (data.x0 - data.xt) * im); error = std::max(error, std::abs(data.tranEr - tranEr)); double tranEm = data.icTranEm.c - data.icTranEm.m * (w0 * data.t0 * slip * data.tranEr - (data.x0 - data.xt) * ir); error = std::max(error, std::abs(data.tranEm - tranEm)); - // wxMessageBox(wxString::Format("%f\n%f\n\n%f\n%f", data.tranEr, data.tranEm, tranEr, tranEm)); - data.tranEr = tranEr; data.tranEm = tranEm; @@ -1763,6 +1779,8 @@ bool Electromechanical::CalculateIndMachinesTransientValues(IndMotor* motor) k = m_powerSystemBase / oldBase; } + data.terminalVoltage = static_cast(motor->GetParentList()[0])->GetElectricalData().voltage; + // Calculate the induction machine transient constants at the machine base double r1t = data.r1 * k; double r2t = data.r2 * k; @@ -1776,21 +1794,19 @@ bool Electromechanical::CalculateIndMachinesTransientValues(IndMotor* motor) data.xmt = xmt; double xt = x1t + (x2t * xmt) / (x2t + xmt); - double t0 = (x2t + xmt) / (2.0 * M_PI * m_systemFreq * r2t); double x0 = x1t + xmt; data.xt = xt; - data.t0 = t0; data.x0 = x0; - data.terminalVoltage = static_cast(motor->GetParentList()[0])->GetElectricalData().voltage; - - //[Ref.] Induction Motor Static Models for Power Flow and Voltage Stability Studies double p = dataPU.activePower; double v = std::abs(data.terminalVoltage); + + //[Ref.] Induction Motor Static Models for Power Flow and Voltage Stability Studies // If the motor is offline, calculate the nominal slip to user-defined power input and 1.0 p.u. voltage if(!motor->IsOnline()) v = 1.0; double r1 = data.r1t; double r2 = data.r2t; + if(data.useKf) r2 *= (1.0 + data.kf * r2t); double x1 = data.x1t; double x2 = data.x2t; double xm = data.xmt; @@ -1815,6 +1831,8 @@ bool Electromechanical::CalculateIndMachinesTransientValues(IndMotor* motor) double qd = r2_s * (xm + x1) + r1 * k1; data.q0 = (-v * v * (qa - qb)) / (qc * qc + qd * qd); + data.t0 = (x2t + xmt) / (2.0 * M_PI * m_systemFreq * r2); + motor->SetElectricalData(data); return true; } -- cgit