From 4a72b22d12c0fc71aba4b4df0e5f057a17b8cdf6 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Wed, 20 Sep 2017 21:54:36 -0300 Subject: Saturation initialization implementation start --- Project/Electromechanical.cpp | 83 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 67 insertions(+), 16 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 57a275a..7f93b6d 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -473,21 +473,72 @@ bool Electromechanical::InitializeDynamicElements() } else if(data.syncXq == 0.0) xq = data.syncXd * k; + double sd = 1.0; + double sq = 1.0; + double xp = data.potierReactance; + bool hasSaturation = false; + if(data.satFactor != 0.0) { // Have saturation. + sd = data.satFactor; + if(xp == 0.0) xp = 0.8 * (data.transXd * k); + hasSaturation = true; + } + // Initialize state variables std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; - data.delta = std::arg(eq0); + double delta = std::arg(eq0); double fi0 = std::arg(ia); - double id0, iq0; - // ABCtoDQ0(ia, data.delta - fi0, id0, iq0); - iq0 = std::abs(ia) * std::cos(data.delta - fi0); - id0 = -std::abs(ia) * std::sin(data.delta - fi0); + double id0, iq0, vd0, vq0; + // iq and id + // iq0 = std::abs(ia) * std::cos(delta - fi0); + // id0 = -std::abs(ia) * std::sin(delta - fi0); + // WRONG!!!! Review!! + ABCtoDQ0(ia, delta, id0, iq0); + ABCtoDQ0(data.terminalVoltage, delta, vd0, vq0); + + // Initialize saturation + double xqs = xq; + double xds = xd; + if(hasSaturation) { + double oldDelta = 0; + bool exit = false; + int it = 0; + while(!exit) { + oldDelta = delta; + + ABCtoDQ0(ia, delta, id0, iq0); + ABCtoDQ0(data.terminalVoltage, delta, vd0, vq0); + + // Direct-axis Potier voltage. + double epd = vd0 + ra * id0 - xp * iq0; + + sq = 1.0 + data.satFactor * (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_tolerance) { + exit = true; + } else if(it >= m_maxIterations) { + m_errorMsg = _("Error on initializate the saturation values of \"") + data.name + _("\"."); + return false; + } + it++; + } + // Quadrature-axis Potier voltage. + double epq = vq0 + ra * iq0 - xp * id0; + sd = 1.0 + data.satFactor * std::pow(epq, 6); + xds = (xd - xp) / sd + xp; + } + + double ef0 = vq0 + ra * iq0 + xds * id0; - data.initialFieldVoltage = std::abs(eq0) - (xd - xq) * id0; + // data.initialFieldVoltage = std::abs(eq0) - (xd - xq) * 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)); data.speed = 2.0 * M_PI * m_systemFreq; + data.delta = delta; data.pe = data.pm; data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); @@ -510,7 +561,7 @@ bool Electromechanical::InitializeDynamicElements() case Machines::SM_MODEL_2: { double tranXd = data.transXd * k; - data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; + data.tranEq = data.initialFieldVoltage + (xd - tranXd) * (id0 / sd); data.tranEd = 0.0; data.subEd = 0.0; data.subEq = 0.0; @@ -520,8 +571,8 @@ bool Electromechanical::InitializeDynamicElements() double tranXq = data.transXq * k; if(tranXq == 0.0) tranXq = tranXd; - data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; - data.tranEd = -(xq - tranXq) * iq0; + data.tranEq = data.initialFieldVoltage + (xd - tranXd) * (id0 / sd); + data.tranEd = -(xq - tranXq) * (iq0 / sq); data.subEd = 0.0; data.subEq = 0.0; @@ -533,10 +584,10 @@ bool Electromechanical::InitializeDynamicElements() if(subXd == 0.0) subXd = subXq; if(subXq == 0.0) subXq = subXd; - data.tranEq = data.initialFieldVoltage + (xd - tranXd) * id0; + data.tranEq = data.initialFieldVoltage + (xd - tranXd) * (id0 / sd); data.tranEd = 0.0; - data.subEq = data.tranEq + (tranXd - subXd) * id0; - data.subEd = -(xq - subXq) * iq0; + 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; @@ -546,10 +597,10 @@ bool Electromechanical::InitializeDynamicElements() 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.tranEq = data.initialFieldVoltage + (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; -- cgit From 748b7c8a73c4bc48ff11469e990532e9e673de0e Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 21 Sep 2017 21:48:21 -0300 Subject: Saturation initialization fixed --- Project/Electromechanical.cpp | 38 ++++++++++++++++++-------------------- 1 file changed, 18 insertions(+), 20 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 7f93b6d..aff509b 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -475,10 +475,11 @@ bool Electromechanical::InitializeDynamicElements() double sd = 1.0; double sq = 1.0; - double xp = data.potierReactance; + double satF = 1.0; + double xp = data.potierReactance * k; bool hasSaturation = false; if(data.satFactor != 0.0) { // Have saturation. - sd = data.satFactor; + satF = (data.satFactor - 1.2) / std::pow(1.2, 7); if(xp == 0.0) xp = 0.8 * (data.transXd * k); hasSaturation = true; } @@ -487,12 +488,7 @@ bool Electromechanical::InitializeDynamicElements() std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; double delta = std::arg(eq0); - double fi0 = std::arg(ia); double id0, iq0, vd0, vq0; - // iq and id - // iq0 = std::abs(ia) * std::cos(delta - fi0); - // id0 = -std::abs(ia) * std::sin(delta - fi0); - // WRONG!!!! Review!! ABCtoDQ0(ia, delta, id0, iq0); ABCtoDQ0(data.terminalVoltage, delta, vd0, vq0); @@ -502,7 +498,7 @@ bool Electromechanical::InitializeDynamicElements() if(hasSaturation) { double oldDelta = 0; bool exit = false; - int it = 0; + int numIt = 0; while(!exit) { oldDelta = delta; @@ -510,42 +506,44 @@ bool Electromechanical::InitializeDynamicElements() ABCtoDQ0(data.terminalVoltage, delta, vd0, vq0); // Direct-axis Potier voltage. - double epd = vd0 + ra * id0 - xp * iq0; + double epd = vd0 + ra * id0 + xp * iq0; - sq = 1.0 + data.satFactor * (xq / xd) * std::pow(epd, 6); + 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_tolerance) { exit = true; - } else if(it >= m_maxIterations) { + } else if(numIt >= m_maxIterations) { m_errorMsg = _("Error on initializate the saturation values of \"") + data.name + _("\"."); return false; } - it++; + numIt++; } // Quadrature-axis Potier voltage. double epq = vq0 + ra * iq0 - xp * id0; - sd = 1.0 + data.satFactor * std::pow(epq, 6); + 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 = std::abs(eq0) - (xd - xq) * 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)); data.speed = 2.0 * M_PI * m_systemFreq; - data.delta = delta; data.pe = data.pm; data.electricalPower = std::complex(dataPU.activePower, dataPU.reactivePower); + data.sd = sd; + data.sq = sq; // 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: { @@ -561,7 +559,7 @@ bool Electromechanical::InitializeDynamicElements() case Machines::SM_MODEL_2: { double tranXd = data.transXd * k; - data.tranEq = data.initialFieldVoltage + (xd - tranXd) * (id0 / sd); + data.tranEq = ef0 + (xd - tranXd) * (id0 / sd); data.tranEd = 0.0; data.subEd = 0.0; data.subEq = 0.0; @@ -571,7 +569,7 @@ bool Electromechanical::InitializeDynamicElements() double tranXq = data.transXq * k; if(tranXq == 0.0) tranXq = tranXd; - data.tranEq = data.initialFieldVoltage + (xd - tranXd) * (id0 / sd); + data.tranEq = ef0 + (xd - tranXd) * (id0 / sd); data.tranEd = -(xq - tranXq) * (iq0 / sq); data.subEd = 0.0; @@ -584,7 +582,7 @@ bool Electromechanical::InitializeDynamicElements() if(subXd == 0.0) subXd = subXq; if(subXq == 0.0) subXq = subXd; - data.tranEq = data.initialFieldVoltage + (xd - tranXd) * (id0 / sd); + 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); @@ -597,7 +595,7 @@ bool Electromechanical::InitializeDynamicElements() if(subXd == 0.0) subXd = subXq; if(subXq == 0.0) subXq = subXd; - data.tranEq = data.initialFieldVoltage + (xd - tranXd) * (id0 / sd); + 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); -- cgit From daff46ee3b27b1fc2088de6b49682b4f41aeb080 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Fri, 22 Sep 2017 21:29:37 -0300 Subject: Saturated voltages implemented Need to fix the current and update the saturation factors --- Project/Electromechanical.cpp | 62 ++++++++++++++++++++++++++++++------------- 1 file changed, 43 insertions(+), 19 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index aff509b..b0ab049 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -762,31 +762,41 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat 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 * transTd0 + m_timeStep); - data.icTranEq.c = (1.0f - 2.0 * data.icTranEq.m) * data.tranEq + + // data.icTranEq.c = (1.0f - 2.0 * data.icTranEq.m) * data.tranEq + + // data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id); + data.icTranEq.c = (1.0 - data.icTranEq.m * (1.0 + data.sd)) * data.tranEq + 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 * transTq0 + m_timeStep); - data.icTranEd.c = (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (syncXq - transXq) * iq; + // data.icTranEd.c = (1.0f - 2.0f * data.icTranEd.m) * data.tranEd - data.icTranEd.m * (syncXq - transXq) * iq; + data.icTranEd.c = + (1.0 - data.icTranEd.m * (1.0 + data.sq)) * 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 * subTd0 + m_timeStep); - data.icSubEq.c = - (1.0f - 2.0f * data.icSubEq.m) * data.subEq + data.icSubEq.m * (data.tranEq + (transXd - subXd) * id); + // data.icSubEq.c = + // (1.0f - 2.0f * data.icSubEq.m) * data.subEq + data.icSubEq.m * (data.tranEq + (transXd - subXd) * id); + data.icSubEq.c = (1.0 - data.icSubEq.m * (1.0 + data.sd)) * data.subEq + + data.icSubEq.m * (data.sd * data.tranEq + (transXd - subXd) * id); } // Ed'' if(data.model == Machines::SM_MODEL_4) { 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; + // data.icSubEd.c = (1.0f - 2.0f * data.icSubEd.m) * data.subEd - data.icSubEd.m * (syncXq - subXq) * iq; + data.icSubEd.c = + (1.0f - data.icSubEd.m * (1.0 + data.sq)) * data.subEd - data.icSubEd.m * (syncXq - subXq) * iq; } if(data.model == Machines::SM_MODEL_5) { 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); + // data.icSubEd.c = + // (1.0f - 2.0f * data.icSubEd.m) * data.subEd + data.icSubEd.m * (data.tranEd - (transXq - subXq) * iq); + data.icSubEd.c = (1.0f - data.icSubEd.m * (1.0 + data.sq)) * data.subEd + + data.icSubEd.m * (data.sq * data.tranEd - (transXq - subXq) * iq); } syncGenerator->SetElectricalData(data); @@ -794,7 +804,6 @@ void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerat bool Electromechanical::SolveSynchronousMachines() { - // CalculateMachinesCurrents(); for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { SyncGenerator* syncGenerator = *it; auto data = syncGenerator->GetElectricalData(); @@ -823,6 +832,8 @@ bool Electromechanical::SolveSynchronousMachines() id = 2.0 * id - data.oldId; iq = 2.0 * iq - data.oldIq; pe = 2.0 * pe - data.oldPe; + data.sd = 2.0 * data.sd - data.oldSd; + data.sq = 2.0 * data.sq - data.oldSq; CalculateSyncMachineIntVariables(syncGenerator, id, iq, pe, k); } else { @@ -1028,8 +1039,12 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe // 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); + double syncXd, transXd; + syncXd = data.syncXd * k; + transXd = data.transXd * k; + + double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / + (1.0 + data.icTranEq.m * (data.sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); data.tranEq = tranEq; @@ -1043,10 +1058,12 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe 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); + double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / + (1.0 + data.icTranEq.m * (data.sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); - double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq; + double tranEd = (data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq) / + (1.0 + data.icTranEd.m * (data.sq - 1.0)); error = std::max(error, std::abs(data.tranEd - tranEd)); data.tranEq = tranEq; @@ -1069,13 +1086,16 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe 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); + double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / + (1.0 + data.icTranEq.m * (data.sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); - double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id); + double subEq = (data.icSubEq.c + data.icSubEq.m * (data.sd * tranEq + (transXd - subXd) * id)) / + (1.0 + data.icSubEq.m * (data.sd - 1.0)); error = std::max(error, std::abs(data.subEq - subEq)); - double subEd = data.icSubEd.c - data.icSubEd.m * (syncXq - subXq) * iq; + double subEd = + (data.icSubEd.c - data.icSubEd.m * ((syncXq - subXq) * iq)) / (1.0 + data.icSubEd.m * (data.sq - 1.0)); error = std::max(error, std::abs(data.subEd - subEd)); data.tranEq = tranEq; @@ -1095,16 +1115,20 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe 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); + double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / + (1.0 + data.icTranEq.m * (data.sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); - double tranEd = data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq; + double tranEd = (data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq) / + (1.0 + data.icTranEd.m * (data.sq - 1.0)); error = std::max(error, std::abs(data.tranEd - tranEd)); - double subEq = data.icSubEq.c + data.icSubEq.m * (tranEq + (transXd - subXd) * id); + double subEq = (data.icSubEq.c + data.icSubEq.m * (data.sd * tranEq + (transXd - subXd) * id)) / + (1.0 + data.icSubEq.m * (data.sd - 1.0)); error = std::max(error, std::abs(data.subEq - subEq)); - double subEd = data.icSubEd.c + data.icSubEd.m * (tranEd - (transXq - subXq) * iq); + double subEd = (data.icSubEd.c + data.icSubEd.m * (data.sq * tranEd - (transXq - subXq) * iq)) / + (1.0 + data.icSubEd.m * (data.sq - 1.0)); error = std::max(error, std::abs(data.subEd - subEd)); data.tranEq = tranEq; -- cgit From 54cf77fe93d3862d8718c6c85c83cf8f8317e68d Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Sat, 23 Sep 2017 22:10:08 -0300 Subject: Saturation Implemented Finally! The calculations are converging! Need to be tested. --- Project/Electromechanical.cpp | 354 ++++++++++++++++++++++++------------------ 1 file changed, 204 insertions(+), 150 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index b0ab049..f0614e3 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -7,7 +7,8 @@ Electromechanical::Electromechanical(wxWindow* parent, std::vector ele GetElementsFromList(elementList); SetEventTimeList(); - m_powerSystemBase = GetPowerValue(data.basePower, data.basePowerUnit); + Bus dummyBus; + m_powerSystemBase = dummyBus.GetValueFromUnit(data.basePower, data.basePowerUnit); m_systemFreq = data.stabilityFrequency; m_simTime = data.stabilitySimulationTime; m_timeStep = data.timeStep; @@ -397,42 +398,14 @@ std::complex Electromechanical::GetSyncMachineAdmittance(SyncGenerator* auto data = generator->GetElectricalData(); double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + double oldBase = generator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } - double xd = 0.0; - double xq = 0.0; double ra = data.armResistance * k; - - 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_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; - } + auto smModelData = GetSyncMachineModelData(generator); + double xd = smModelData.xd; + double xq = smModelData.xq; double xdq = 0.5 * (xd + xq); return (std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0)); } @@ -454,7 +427,7 @@ bool Electromechanical::InitializeDynamicElements() if(syncGenerator->IsOnline()) { double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; @@ -501,9 +474,8 @@ bool Electromechanical::InitializeDynamicElements() int numIt = 0; while(!exit) { oldDelta = delta; - - ABCtoDQ0(ia, delta, id0, iq0); - ABCtoDQ0(data.terminalVoltage, delta, vd0, vq0); + // ABCtoDQ0(ia, delta, id0, iq0); + // ABCtoDQ0(data.terminalVoltage, delta, vd0, vq0); // Direct-axis Potier voltage. double epd = vd0 + ra * id0 + xp * iq0; @@ -512,7 +484,7 @@ bool Electromechanical::InitializeDynamicElements() xqs = (xq - xp) / sq + xp; eq0 = data.terminalVoltage + std::complex(ra, xqs) * ia; delta = std::arg(eq0); - if(std::abs(delta - oldDelta) < m_tolerance) { + 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 + _("\"."); @@ -642,7 +614,7 @@ bool Electromechanical::InitializeDynamicElements() return true; } -void Electromechanical::CalculateMachinesCurrents() +bool 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); @@ -653,66 +625,62 @@ void Electromechanical::CalculateMachinesCurrents() if(syncGenerator->IsOnline()) { double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } - double xd = 0.0; - double xq = 0.0; double ra = data.armResistance * k; + + double xp = data.potierReactance * k; + if(xp == 0.0) xp = 0.8 * (data.transXd * 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(data.model) { - 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_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; + auto smModelData = GetSyncMachineModelData(syncGenerator); + DQ0toABC(smModelData.ed, smModelData.eq, data.delta, e); + double xd = smModelData.xd; + double xq = smModelData.xq; + + // Calculate the saturation effect + if(data.satFactor != 0.0) { + if(!CalculateSyncMachineSaturation(syncGenerator, k)) return false; } + + double sd = data.sd; + double sq = data.sq; + + double xdq, xds, xqs, xdqs; xdq = 0.5 * (xd + xq); + xds = (xd - xp) / sd + xp; + xqs = (xq - xp) / sq + xp; + xdqs = 0.5 * (xds + xqs); std::complex y0 = std::complex(ra, -xdq) / std::complex(ra * ra + xd * xq, 0.0); - std::complex iUnadj = y0 * e; + // std::complex iUnadjusted = y0 * e; + std::complex iUnadjusted = y0 * 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)); + // [Ref] Arrillaga, J.; Arnold, C. P.. "Computer Modelling of Electrical Power Systems". Pg. 225-226 + // [Ref] Dommell, H. W.; Sato, N.. "Fast transient stability solutions". IEEE Transactions on Power + // Apparatus and Systems, PAS-91 (4), 1643-1650 + std::complex iSaliency = std::complex(0.0, -((0.5 * (xqs - xds)) / (ra * ra + xds * xqs))) * + (std::conj(e) - std::conj(v)); + iSaliency = iSaliency * std::cos(2.0 * data.delta) + + iSaliency * std::complex(0.0, std::sin(2.0 * data.delta)); - iInj = iUnadj + iAdj; + // [Ref] Arrillaga, J.; Arnold, C. P.; Computer Modelling of Electrical Power Systems. Pg. 258-259 + std::complex y0s = std::complex(ra, -xdqs) / std::complex(ra * ra + xds * xqs, 0.0); + std::complex iSaturation = y0s * (e - v); - m_iBus[n] += iInj; + iInj = iUnadjusted + iSaliency + iSaturation; - std::complex iMachine = iInj - y0 * v; + m_iBus[n] += iInj; + // Remove the current flowing through y0 (i.e. iUnadjusted in this case, y0 is inserted in admittance + // matrix) to calculate the electrical power. + std::complex iMachine = iInj - iUnadjusted; data.electricalPower = v * std::conj(iMachine); } else { data.electricalPower = std::complex(0.0, 0.0); @@ -720,6 +688,7 @@ void Electromechanical::CalculateMachinesCurrents() syncGenerator->SetElectricalData(data); } + return true; } void Electromechanical::CalculateIntegrationConstants(SyncGenerator* syncGenerator, double id, double iq, double k) @@ -810,16 +779,13 @@ bool Electromechanical::SolveSynchronousMachines() if(syncGenerator->IsOnline()) { int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; - double id, iq, pe; - - pe = data.pe; + double id, iq, pe, sd, sq; double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + double oldBase = syncGenerator->GetValueFromUnit(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); @@ -827,15 +793,15 @@ bool Electromechanical::SolveSynchronousMachines() // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq, k); - CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, pe, k); + CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k); // Extrapolate nonintegrable variables. id = 2.0 * id - data.oldId; iq = 2.0 * iq - data.oldIq; pe = 2.0 * pe - data.oldPe; - data.sd = 2.0 * data.sd - data.oldSd; - data.sq = 2.0 * data.sq - data.oldSq; + sd = 2.0 * sd - data.oldSd; + sq = 2.0 * sq - data.oldSq; - CalculateSyncMachineIntVariables(syncGenerator, id, iq, pe, k); + CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); } else { CalculateIntegrationConstants(syncGenerator, 0.0f, 0.0f); } @@ -849,7 +815,7 @@ bool Electromechanical::SolveSynchronousMachines() error = 0.0; // Calculate the injected currents. - CalculateMachinesCurrents(); + if(!CalculateMachinesCurrents()) return false; // Calculate the buses voltages. m_vBus = LUEvaluate(m_yBusU, m_yBusL, m_iBus); @@ -860,16 +826,16 @@ bool Electromechanical::SolveSynchronousMachines() auto data = syncGenerator->GetElectricalData(); - double id, iq, pe; + double id, iq, pe, sd, sq; double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } - CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, pe, k); + CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k); - double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, pe, k); + double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); if(genError > error) error = genError; } @@ -904,46 +870,6 @@ 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; -} - void Electromechanical::SaveData() { for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) { @@ -985,6 +911,8 @@ void Electromechanical::SetSyncMachinesModel() void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncGenerator, double& id, double& iq, + double& sd, + double& sq, double& pe, double k) { @@ -998,6 +926,9 @@ void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG double vd, vq; ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); + sd = data.sd; + sq = data.sq; + if(syncGenerator->IsOnline()) { std::complex iMachine = std::conj(data.electricalPower) / std::conj(m_vBus[n]); ABCtoDQ0(iMachine, data.delta, id, iq); @@ -1007,14 +938,19 @@ void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG pe = id = iq = 0.0f; } data.pe = pe; + data.oldPe = pe; data.oldId = id; data.oldIq = iq; + data.oldSd = sd; + data.oldSq = sq; syncGenerator->SetElectricalData(data); } double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGenerator, double id, double iq, + double sd, + double sq, double pe, double k) { @@ -1044,7 +980,7 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe transXd = data.transXd * k; double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / - (1.0 + data.icTranEq.m * (data.sd - 1.0)); + (1.0 + data.icTranEq.m * (sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); data.tranEq = tranEq; @@ -1059,11 +995,11 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe if(transXq == 0.0) transXq = transXd; double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / - (1.0 + data.icTranEq.m * (data.sd - 1.0)); + (1.0 + data.icTranEq.m * (sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); - double tranEd = (data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq) / - (1.0 + data.icTranEd.m * (data.sq - 1.0)); + double tranEd = + (data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq) / (1.0 + data.icTranEd.m * (sq - 1.0)); error = std::max(error, std::abs(data.tranEd - tranEd)); data.tranEq = tranEq; @@ -1087,11 +1023,11 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe if(subXq == 0.0) subXq = subXd; double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / - (1.0 + data.icTranEq.m * (data.sd - 1.0)); + (1.0 + data.icTranEq.m * (sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); - double subEq = (data.icSubEq.c + data.icSubEq.m * (data.sd * tranEq + (transXd - subXd) * id)) / - (1.0 + data.icSubEq.m * (data.sd - 1.0)); + double subEq = (data.icSubEq.c + data.icSubEq.m * (sd * tranEq + (transXd - subXd) * id)) / + (1.0 + data.icSubEq.m * (sd - 1.0)); error = std::max(error, std::abs(data.subEq - subEq)); double subEd = @@ -1116,19 +1052,19 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe if(subXq == 0.0) subXq = subXd; double tranEq = (data.icTranEq.c + data.icTranEq.m * (data.fieldVoltage + (syncXd - transXd) * id)) / - (1.0 + data.icTranEq.m * (data.sd - 1.0)); + (1.0 + data.icTranEq.m * (sd - 1.0)); error = std::max(error, std::abs(data.tranEq - tranEq)); - double tranEd = (data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq) / - (1.0 + data.icTranEd.m * (data.sq - 1.0)); + double tranEd = + (data.icTranEd.c - data.icTranEd.m * (syncXq - transXq) * iq) / (1.0 + data.icTranEd.m * (sq - 1.0)); error = std::max(error, std::abs(data.tranEd - tranEd)); - double subEq = (data.icSubEq.c + data.icSubEq.m * (data.sd * tranEq + (transXd - subXd) * id)) / - (1.0 + data.icSubEq.m * (data.sd - 1.0)); + double subEq = (data.icSubEq.c + data.icSubEq.m * (sd * tranEq + (transXd - subXd) * id)) / + (1.0 + data.icSubEq.m * (sd - 1.0)); error = std::max(error, std::abs(data.subEq - subEq)); - double subEd = (data.icSubEd.c + data.icSubEd.m * (data.sq * tranEd - (transXq - subXq) * iq)) / - (1.0 + data.icSubEd.m * (data.sq - 1.0)); + double subEd = (data.icSubEd.c + data.icSubEd.m * (sq * tranEd - (transXq - subXq) * iq)) / + (1.0 + data.icSubEd.m * (sq - 1.0)); error = std::max(error, std::abs(data.subEd - subEd)); data.tranEq = tranEq; @@ -1153,7 +1089,7 @@ void Electromechanical::CalculateReferenceSpeed() auto data = syncGenerator->GetElectricalData(); double k = 1.0; // Power base change factor. if(data.useMachineBase) { - double oldBase = GetPowerValue(data.nominalPower, data.nominalPowerUnit); + double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } sumH += data.inertia / k; @@ -1165,3 +1101,121 @@ void Electromechanical::CalculateReferenceSpeed() m_refSpeed = 2.0 * M_PI * m_systemFreq; } } + +bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachine, double k) +{ + // [Ref] Arrillaga, J.; Arnold, C. P.. "Computer Modelling of Electrical Power Systems". Pg. 254-260 + auto data = syncMachine->GetElectricalData(); + auto smDataModel = GetSyncMachineModelData(syncMachine); + + double vd, vq; + ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); + double deltaVd = smDataModel.ed - vd; + double deltaVq = smDataModel.eq - vq; + + double ra = data.armResistance * k; + double xd = smDataModel.xd; + double xq = smDataModel.xq; + + double syncXd = data.syncXd * k; + double syncXq = data.syncXq * k; + + if(data.model == Machines::SM_MODEL_1) { + syncXq = data.transXd * k; + syncXd = syncXq; + } else if(data.syncXq == 0.0) + syncXd = data.syncXd * k; + + double xp = data.potierReactance * k; + if(xp == 0.0) xp = 0.8 * (data.transXd * k); + double sd = data.sd; + double sq = data.sq; + double satFacd = (data.satFactor - 1.2) / std::pow(1.2, 7); + double satFacq = satFacd * (syncXq / syncXd); + + bool exit = false; + int iterations = 0; + while(!exit) { + double oldSd = sd; + double oldSq = sq; + + // Saturated reactances. + double xds = (xd - xp) / sd + xp; + double xqs = (xq - xp) / sq + xp; + // dq currents. + double den = 1.0 / (ra * ra + xds * xqs); + double iq = den * (ra * deltaVq + xds * deltaVd); + double id = den * (-xqs * deltaVq + ra * deltaVd); + // Potier voltages + double epq = vq + ra * iq - xp * id; + double epd = vd + ra * id + xp * iq; + // Saturation factors. + sd = 1.0 + satFacd * std::pow(epq, 6); + sq = 1.0 + satFacq * std::pow(epd, 6); + + double error = std::abs(sd - oldSd) + std::abs(sq - oldSq); + if(error < m_saturationTolerance) exit = true; + + iterations++; + if((iterations >= m_maxIterations) & !exit) { + m_errorMsg = + _("It was not possible to solve the saturation of the synchronous machine \"") + data.name + wxT("\"."); + return false; + } + } + + data.sd = sd; + data.sq = sq; + syncMachine->SetElectricalData(data); + return true; +} + +SyncMachineModelData Electromechanical::GetSyncMachineModelData(SyncGenerator* syncMachine) +{ + SyncMachineModelData smModelData; + + auto data = syncMachine->GetElectricalData(); + double k = 1.0; // Power base change factor. + if(data.useMachineBase) { + double oldBase = syncMachine->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); + k = m_powerSystemBase / oldBase; + } + + switch(data.model) { + case Machines::SM_MODEL_1: { + smModelData.ed = data.tranEd; + smModelData.eq = data.tranEq; + smModelData.xq = data.transXd * k; + smModelData.xd = smModelData.xq; + } break; + case Machines::SM_MODEL_2: { + smModelData.ed = data.tranEd; + smModelData.eq = data.tranEq; + smModelData.xd = data.transXd * k; + smModelData.xq = data.transXq * k; + if(smModelData.xq == 0.0) { + smModelData.xq = data.syncXq * k; + if(smModelData.xq == 0.0) { + smModelData.xq = data.syncXd * k; + } + } + } break; + case Machines::SM_MODEL_3: { + smModelData.ed = data.tranEd; + smModelData.eq = data.tranEq; + smModelData.xd = data.transXd * k; + smModelData.xq = data.transXq * k; + if(smModelData.xq == 0.0) smModelData.xq = smModelData.xd; + } break; + case Machines::SM_MODEL_4: + case Machines::SM_MODEL_5: { + smModelData.ed = data.subEd; + smModelData.eq = data.subEq; + smModelData.xd = data.subXd * k; + smModelData.xq = data.subXq * k; + if(smModelData.xd == 0.0) smModelData.xd = smModelData.xq; + if(smModelData.xq == 0.0) smModelData.xq = smModelData.xd; + } break; + } + return smModelData; +} -- cgit From 5e99db91624cea9308cb296544d81ce1fa4951e5 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Mon, 25 Sep 2017 01:50:23 -0300 Subject: Zero xd bug fixed, some parameters showed The calculation is wrong --- Project/Electromechanical.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index f0614e3..49c1352 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -682,6 +682,8 @@ bool Electromechanical::CalculateMachinesCurrents() // matrix) to calculate the electrical power. std::complex iMachine = iInj - iUnadjusted; data.electricalPower = v * std::conj(iMachine); + data.sd = sd; + data.sq = sq; } else { data.electricalPower = std::complex(0.0, 0.0); } @@ -838,6 +840,8 @@ bool Electromechanical::SolveSynchronousMachines() double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); if(genError > error) error = genError; + m_sdC = data.sd; + m_sqC = data.sq; } ++iterations; @@ -896,6 +900,8 @@ void Electromechanical::SaveData() m_wErrorVector.push_back(m_wError); m_numItVector.push_back(m_numIt); + m_sdCVector.push_back(m_sdC); + m_sqCVector.push_back(m_sqC); } void Electromechanical::SetSyncMachinesModel() @@ -926,6 +932,7 @@ void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG double vd, vq; ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); + CalculateSyncMachineSaturation(syncGenerator, k); sd = data.sd; sq = data.sq; @@ -938,6 +945,8 @@ void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG pe = id = iq = 0.0f; } data.pe = pe; + data.sd = sd; + data.sq = sq; data.oldPe = pe; data.oldId = id; data.oldIq = iq; @@ -1124,7 +1133,7 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin syncXq = data.transXd * k; syncXd = syncXq; } else if(data.syncXq == 0.0) - syncXd = data.syncXd * k; + syncXq = data.syncXd * k; double xp = data.potierReactance * k; if(xp == 0.0) xp = 0.8 * (data.transXd * k); -- cgit From 4f2e475b8de0f19c0cda560c3213809cb9c8cfca Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Mon, 25 Sep 2017 19:56:07 -0300 Subject: Saturation calculated by Newton raphson implmeneted Not working yet --- Project/Electromechanical.cpp | 86 +++++++++++++++++++++++++++++++------------ 1 file changed, 62 insertions(+), 24 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 49c1352..73fd122 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -509,6 +509,8 @@ bool Electromechanical::InitializeDynamicElements() 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; @@ -630,9 +632,8 @@ bool Electromechanical::CalculateMachinesCurrents() } double ra = data.armResistance * k; - double xp = data.potierReactance * k; - if(xp == 0.0) xp = 0.8 * (data.transXd * k); + if(xp == 0.0) xp = 0.8 * data.transXd * k; int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; std::complex e = std::complex(0.0, 0.0); @@ -646,7 +647,7 @@ bool Electromechanical::CalculateMachinesCurrents() // Calculate the saturation effect if(data.satFactor != 0.0) { - if(!CalculateSyncMachineSaturation(syncGenerator, k)) return false; + if(!CalculateSyncMachineSaturation(syncGenerator, false, k)) return false; } double sd = data.sd; @@ -682,6 +683,12 @@ bool Electromechanical::CalculateMachinesCurrents() // matrix) to calculate the electrical power. std::complex iMachine = iInj - iUnadjusted; data.electricalPower = v * std::conj(iMachine); + + double id, iq; + ABCtoDQ0(iMachine, data.delta, id, iq); + + data.id = id; + data.iq = iq; data.sd = sd; data.sq = sq; } else { @@ -780,7 +787,6 @@ bool Electromechanical::SolveSynchronousMachines() auto data = syncGenerator->GetElectricalData(); if(syncGenerator->IsOnline()) { - int n = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().number; double id, iq, pe, sd, sq; double k = 1.0; // Power base change factor. @@ -788,14 +794,13 @@ bool Electromechanical::SolveSynchronousMachines() double oldBase = syncGenerator->GetValueFromUnit(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); + id = data.id; + iq = data.iq; // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq, k); - CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k); + if(!CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k)) return false; // Extrapolate nonintegrable variables. id = 2.0 * id - data.oldId; iq = 2.0 * iq - data.oldIq; @@ -835,13 +840,16 @@ bool Electromechanical::SolveSynchronousMachines() k = m_powerSystemBase / oldBase; } - CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k); + // Calculate id, iq, dq, sd + if(!CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k)) return false; + if(data.satFactor != 0.0) { + m_sdC = sd; + m_sqC = sq; + } double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); if(genError > error) error = genError; - m_sdC = data.sd; - m_sqC = data.sq; } ++iterations; @@ -914,7 +922,7 @@ void Electromechanical::SetSyncMachinesModel() } } -void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncGenerator, +bool Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncGenerator, double& id, double& iq, double& sd, @@ -931,15 +939,15 @@ void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG double vd, vq; ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); - - CalculateSyncMachineSaturation(syncGenerator, k); + if(data.satFactor != 0.0) { + if(!CalculateSyncMachineSaturation(syncGenerator, k)) return false; + } sd = data.sd; sq = data.sq; + id = data.id; + iq = data.iq; 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; @@ -947,12 +955,16 @@ void Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG data.pe = pe; data.sd = sd; data.sq = sq; + data.id = id; + data.iq = iq; data.oldPe = pe; data.oldId = id; data.oldIq = iq; data.oldSd = sd; data.oldSq = sq; syncGenerator->SetElectricalData(data); + + return true; } double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGenerator, @@ -1111,16 +1123,27 @@ void Electromechanical::CalculateReferenceSpeed() } } -bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachine, double k) +bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachine, bool updateCurrents, double k) { // [Ref] Arrillaga, J.; Arnold, C. P.. "Computer Modelling of Electrical Power Systems". Pg. 254-260 auto data = syncMachine->GetElectricalData(); auto smDataModel = GetSyncMachineModelData(syncMachine); + int n = static_cast(syncMachine->GetParentList()[0])->GetElectricalData().number; + if(syncMachine->IsOnline()) { + data.terminalVoltage = m_vBus[n]; + } + double vd, vq; ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); double deltaVd = smDataModel.ed - vd; double deltaVq = smDataModel.eq - vq; + // wxMessageBox(wxString::Format("Vd = %.13f\nEd = %.13f\nVq = %.13f\nEq = %.13f", vd, smDataModel.ed, vq, + // smDataModel.eq)); + + double id, iq; + id = data.id; + iq = data.iq; double ra = data.armResistance * k; double xd = smDataModel.xd; @@ -1128,7 +1151,6 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin double syncXd = data.syncXd * k; double syncXq = data.syncXq * k; - if(data.model == Machines::SM_MODEL_1) { syncXq = data.transXd * k; syncXd = syncXq; @@ -1136,7 +1158,7 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin syncXq = data.syncXd * k; double xp = data.potierReactance * k; - if(xp == 0.0) xp = 0.8 * (data.transXd * k); + if(xp == 0.0) xp = 0.8 * data.transXd * k; double sd = data.sd; double sq = data.sq; double satFacd = (data.satFactor - 1.2) / std::pow(1.2, 7); @@ -1144,6 +1166,7 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin bool exit = false; int iterations = 0; + // wxMessageBox(wxString::Format("%.13f\n%.13f", deltaVd, deltaVq)); while(!exit) { double oldSd = sd; double oldSq = sq; @@ -1153,14 +1176,25 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin double xqs = (xq - xp) / sq + xp; // dq currents. double den = 1.0 / (ra * ra + xds * xqs); - double iq = den * (ra * deltaVq + xds * deltaVd); - double id = den * (-xqs * deltaVq + ra * deltaVd); + iq = den * (ra * deltaVq + xds * deltaVd); + id = den * (-xqs * deltaVq + ra * deltaVd); // Potier voltages double epq = vq + ra * iq - xp * id; double epd = vd + ra * id + xp * iq; // Saturation factors. - sd = 1.0 + satFacd * std::pow(epq, 6); - sq = 1.0 + satFacq * std::pow(epd, 6); + + // Gauss + // sd = 1.0 + satFacd * std::pow(epq, 6); + // sq = 1.0 + satFacq * std::pow(epd, 6); + + // Newton-raphson + double f1 = 1.0 - sd + satFacd * std::pow(epq, 6); + double f2 = 1.0 - sq + satFacq * std::pow(epd, 6); + double dF1dSd = (6.0 * satFacd * std::pow(epq, 5) * xp * (xd - xp) * deltaVq) / ((sd - 1.0) * xp + xd) - 1.0; + double dF2dSq = (6.0 * satFacq * std::pow(epd, 5) * xp * (xq - xp) * deltaVd) / ((sq - 1.0) * xp + xq) - 1.0; + + sd = sd - f1 / dF1dSd; + sq = sq - f2 / dF2dSq; double error = std::abs(sd - oldSd) + std::abs(sq - oldSq); if(error < m_saturationTolerance) exit = true; @@ -1175,6 +1209,10 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin data.sd = sd; data.sq = sq; + if(updateCurrents) { + data.id = id; + data.iq = iq; + } syncMachine->SetElectricalData(data); return true; } -- cgit From 3cd92c1acb9a6fed1067b49cb079e96d98ac35fa Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 28 Sep 2017 20:11:12 -0300 Subject: Some saturation bugs fixed --- Project/Electromechanical.cpp | 124 +++++++++++++++++++++++------------------- 1 file changed, 69 insertions(+), 55 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 73fd122..84e49d5 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -496,6 +496,11 @@ bool Electromechanical::InitializeDynamicElements() double epq = vq0 + ra * iq0 - xp * id0; sd = 1.0 + satF * std::pow(epq, 6); xds = (xd - xp) / sd + xp; + /*CalculateSyncMachineSaturation(syncGenerator, id0, iq0, sq, sd, true, k); + xqs = (xq - xp) / sq + xp; + xds = (xd - xp) / sd + xp; + eq0 = data.terminalVoltage + std::complex(ra, xqs) * ia; + delta = std::arg(eq0);*/ } double ef0 = vq0 + ra * iq0 - xds * id0; @@ -518,6 +523,9 @@ bool Electromechanical::InitializeDynamicElements() data.oldPe = data.pe; data.oldSd = sd; data.oldSq = sq; + + m_sdC = sd; + m_sqC = sq; switch(data.model) { case Machines::SM_MODEL_1: { @@ -645,14 +653,15 @@ bool Electromechanical::CalculateMachinesCurrents() double xd = smModelData.xd; double xq = smModelData.xq; + double sd = data.sd; + double sq = data.sq; + double id, iq; + // Calculate the saturation effect if(data.satFactor != 0.0) { - if(!CalculateSyncMachineSaturation(syncGenerator, false, k)) return false; + if(!CalculateSyncMachineSaturation(syncGenerator, id, iq, sd, sq, false, k)) return false; } - double sd = data.sd; - double sq = data.sq; - double xdq, xds, xqs, xdqs; xdq = 0.5 * (xd + xq); xds = (xd - xp) / sd + xp; @@ -684,7 +693,6 @@ bool Electromechanical::CalculateMachinesCurrents() std::complex iMachine = iInj - iUnadjusted; data.electricalPower = v * std::conj(iMachine); - double id, iq; ABCtoDQ0(iMachine, data.delta, id, iq); data.id = id; @@ -788,14 +796,17 @@ bool Electromechanical::SolveSynchronousMachines() if(syncGenerator->IsOnline()) { double id, iq, pe, sd, sq; + pe = data.pe; + id = data.id; + iq = data.iq; + sd = data.sd; + sq = data.sq; double k = 1.0; // Power base change factor. if(data.useMachineBase) { double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; } - id = data.id; - iq = data.iq; // Calculate integration constants. CalculateIntegrationConstants(syncGenerator, id, iq, k); @@ -807,6 +818,9 @@ bool Electromechanical::SolveSynchronousMachines() pe = 2.0 * pe - data.oldPe; sd = 2.0 * sd - data.oldSd; sq = 2.0 * sq - data.oldSq; + + m_sdC = sd; + m_sqC = sq; CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); } else { @@ -833,8 +847,14 @@ bool Electromechanical::SolveSynchronousMachines() auto data = syncGenerator->GetElectricalData(); - double id, iq, pe, sd, sq; - double k = 1.0; // Power base change factor. + double id = data.id; + double iq = data.iq; + double pe = data.pe; + double sd = data.sd; + double sq = data.sq; + + // Power base change factor. + double k = 1.0; if(data.useMachineBase) { double oldBase = syncGenerator->GetValueFromUnit(data.nominalPower, data.nominalPowerUnit); k = m_powerSystemBase / oldBase; @@ -842,10 +862,6 @@ bool Electromechanical::SolveSynchronousMachines() // Calculate id, iq, dq, sd if(!CalculateSyncMachineNonIntVariables(syncGenerator, id, iq, sd, sq, pe, k)) return false; - if(data.satFactor != 0.0) { - m_sdC = sd; - m_sqC = sq; - } double genError = CalculateSyncMachineIntVariables(syncGenerator, id, iq, sd, sq, pe, k); @@ -939,13 +955,14 @@ bool Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG double vd, vq; ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); + if(data.satFactor != 0.0) { - if(!CalculateSyncMachineSaturation(syncGenerator, k)) return false; + if(!CalculateSyncMachineSaturation(syncGenerator, id, iq, sd, sq, true, k)) return false; + data.sd = sd; + data.sq = sq; + data.oldSd = sd; + data.oldSq = sq; } - sd = data.sd; - sq = data.sq; - id = data.id; - iq = data.iq; if(syncGenerator->IsOnline()) { pe = id * vd + iq * vq + (id * id + iq * iq) * data.armResistance * k; @@ -953,15 +970,11 @@ bool Electromechanical::CalculateSyncMachineNonIntVariables(SyncGenerator* syncG pe = id = iq = 0.0f; } data.pe = pe; - data.sd = sd; - data.sq = sq; data.id = id; data.iq = iq; data.oldPe = pe; data.oldId = id; data.oldIq = iq; - data.oldSd = sd; - data.oldSq = sq; syncGenerator->SetElectricalData(data); return true; @@ -1123,7 +1136,13 @@ void Electromechanical::CalculateReferenceSpeed() } } -bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachine, bool updateCurrents, double k) +bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachine, + double& id, + double& iq, + double& sd, + double& sq, + bool updateCurrents, + double k) { // [Ref] Arrillaga, J.; Arnold, C. P.. "Computer Modelling of Electrical Power Systems". Pg. 254-260 auto data = syncMachine->GetElectricalData(); @@ -1133,17 +1152,15 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin if(syncMachine->IsOnline()) { data.terminalVoltage = m_vBus[n]; } + double idCalc = id; + double iqCalc = iq; + double sdCalc = sd; + double sqCalc = sq; double vd, vq; ABCtoDQ0(data.terminalVoltage, data.delta, vd, vq); double deltaVd = smDataModel.ed - vd; double deltaVq = smDataModel.eq - vq; - // wxMessageBox(wxString::Format("Vd = %.13f\nEd = %.13f\nVq = %.13f\nEq = %.13f", vd, smDataModel.ed, vq, - // smDataModel.eq)); - - double id, iq; - id = data.id; - iq = data.iq; double ra = data.armResistance * k; double xd = smDataModel.xd; @@ -1159,44 +1176,42 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin double xp = data.potierReactance * k; if(xp == 0.0) xp = 0.8 * data.transXd * k; - double sd = data.sd; - double sq = data.sq; double satFacd = (data.satFactor - 1.2) / std::pow(1.2, 7); double satFacq = satFacd * (syncXq / syncXd); bool exit = false; int iterations = 0; - // wxMessageBox(wxString::Format("%.13f\n%.13f", deltaVd, deltaVq)); while(!exit) { - double oldSd = sd; - double oldSq = sq; + double oldSd = sdCalc; + double oldSq = sqCalc; // Saturated reactances. - double xds = (xd - xp) / sd + xp; - double xqs = (xq - xp) / sq + xp; + double xds = (xd - xp) / sdCalc + xp; + double xqs = (xq - xp) / sqCalc + xp; // dq currents. double den = 1.0 / (ra * ra + xds * xqs); - iq = den * (ra * deltaVq + xds * deltaVd); - id = den * (-xqs * deltaVq + ra * deltaVd); + iqCalc = den * (ra * deltaVq + xds * deltaVd); + idCalc = den * (-xqs * deltaVq + ra * deltaVd); // Potier voltages - double epq = vq + ra * iq - xp * id; - double epd = vd + ra * id + xp * iq; + double epq = vq + ra * iqCalc - xp * idCalc; + double epd = vd + ra * idCalc + xp * iqCalc; // Saturation factors. - // Gauss - // sd = 1.0 + satFacd * std::pow(epq, 6); - // sq = 1.0 + satFacq * std::pow(epd, 6); + /*sdCalc = 1.0 + satFacd * std::pow(epq, 6); + sqCalc = 1.0 + satFacq * std::pow(epd, 6);*/ // Newton-raphson - double f1 = 1.0 - sd + satFacd * std::pow(epq, 6); - double f2 = 1.0 - sq + satFacq * std::pow(epd, 6); - double dF1dSd = (6.0 * satFacd * std::pow(epq, 5) * xp * (xd - xp) * deltaVq) / ((sd - 1.0) * xp + xd) - 1.0; - double dF2dSq = (6.0 * satFacq * std::pow(epd, 5) * xp * (xq - xp) * deltaVd) / ((sq - 1.0) * xp + xq) - 1.0; + double f1 = 1.0 - sdCalc + satFacd * std::pow(epq, 6); + double f2 = 1.0 - sqCalc + satFacq * std::pow(epd, 6); + double dF1dSd = + (6.0 * satFacd * std::pow(epq, 5) * xp * (xd - xp) * deltaVq) / ((sdCalc - 1.0) * xp + xd) - 1.0; + double dF2dSq = + (6.0 * satFacq * std::pow(epd, 5) * xp * (xq - xp) * deltaVd) / ((sqCalc - 1.0) * xp + xq) - 1.0; - sd = sd - f1 / dF1dSd; - sq = sq - f2 / dF2dSq; + sdCalc = sdCalc - f1 / dF1dSd; + sqCalc = sqCalc - f2 / dF2dSq; - double error = std::abs(sd - oldSd) + std::abs(sq - oldSq); + double error = std::abs(sdCalc - oldSd) + std::abs(sqCalc - oldSq); if(error < m_saturationTolerance) exit = true; iterations++; @@ -1207,13 +1222,12 @@ bool Electromechanical::CalculateSyncMachineSaturation(SyncGenerator* syncMachin } } - data.sd = sd; - data.sq = sq; + sd = sdCalc; + sq = sqCalc; if(updateCurrents) { - data.id = id; - data.iq = iq; + id = idCalc; + iq = iqCalc; } - syncMachine->SetElectricalData(data); return true; } -- cgit From c7c1963fae6f88a55b3f7e4a3f9ea17e964d23d2 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Wed, 4 Oct 2017 21:06:20 -0300 Subject: Saturation fixed in q axis --- Project/Electromechanical.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 84e49d5..4b40c53 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -433,8 +433,8 @@ bool Electromechanical::InitializeDynamicElements() data.terminalVoltage = static_cast(syncGenerator->GetParentList()[0])->GetElectricalData().voltage; std::complex conjS(dataPU.activePower, -dataPU.reactivePower); - std::complex conjV = std::conj(data.terminalVoltage); - std::complex ia = conjS / conjV; + std::complex vt = data.terminalVoltage; + std::complex ia = conjS / std::conj(vt); double xd = data.syncXd * k; double xq = data.syncXq * k; @@ -458,12 +458,12 @@ bool Electromechanical::InitializeDynamicElements() } // Initialize state variables - std::complex eq0 = data.terminalVoltage + std::complex(ra, xq) * ia; + 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(data.terminalVoltage, delta, vd0, vq0); + ABCtoDQ0(vt, delta, vd0, vq0); // Initialize saturation double xqs = xq; @@ -474,8 +474,8 @@ bool Electromechanical::InitializeDynamicElements() int numIt = 0; while(!exit) { oldDelta = delta; - // ABCtoDQ0(ia, delta, id0, iq0); - // ABCtoDQ0(data.terminalVoltage, delta, vd0, vq0); + ABCtoDQ0(ia, delta, id0, iq0); + ABCtoDQ0(vt, delta, vd0, vq0); // Direct-axis Potier voltage. double epd = vd0 + ra * id0 + xp * iq0; @@ -502,7 +502,7 @@ bool Electromechanical::InitializeDynamicElements() eq0 = data.terminalVoltage + std::complex(ra, xqs) * ia; delta = std::arg(eq0);*/ } - + double ef0 = vq0 + ra * iq0 - xds * id0; data.initialFieldVoltage = ef0 * sd; @@ -1065,7 +1065,7 @@ double Electromechanical::CalculateSyncMachineIntVariables(SyncGenerator* syncGe error = std::max(error, std::abs(data.subEq - subEq)); double subEd = - (data.icSubEd.c - data.icSubEd.m * ((syncXq - subXq) * iq)) / (1.0 + data.icSubEd.m * (data.sq - 1.0)); + (data.icSubEd.c - data.icSubEd.m * ((syncXq - subXq) * iq)) / (1.0 + data.icSubEd.m * (sq - 1.0)); error = std::max(error, std::abs(data.subEd - subEd)); data.tranEq = tranEq; -- cgit From bc5a3e8923cb8efedbbd5b88e212eb0e9009cf87 Mon Sep 17 00:00:00 2001 From: Thales Lima Oliveira Date: Thu, 5 Oct 2017 20:02:59 -0300 Subject: Add copyright on files and documentation update --- Project/Electromechanical.cpp | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) (limited to 'Project/Electromechanical.cpp') diff --git a/Project/Electromechanical.cpp b/Project/Electromechanical.cpp index 4b40c53..4b8a711 100644 --- a/Project/Electromechanical.cpp +++ b/Project/Electromechanical.cpp @@ -1,3 +1,20 @@ +/* + * Copyright (C) 2017 Thales Lima Oliveira + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + #include "Electromechanical.h" #include "ControlElementSolver.h" @@ -502,7 +519,7 @@ bool Electromechanical::InitializeDynamicElements() eq0 = data.terminalVoltage + std::complex(ra, xqs) * ia; delta = std::arg(eq0);*/ } - + double ef0 = vq0 + ra * iq0 - xds * id0; data.initialFieldVoltage = ef0 * sd; @@ -523,7 +540,7 @@ bool Electromechanical::InitializeDynamicElements() data.oldPe = data.pe; data.oldSd = sd; data.oldSq = sq; - + m_sdC = sd; m_sqC = sq; @@ -818,7 +835,7 @@ bool Electromechanical::SolveSynchronousMachines() pe = 2.0 * pe - data.oldPe; sd = 2.0 * sd - data.oldSd; sq = 2.0 * sq - data.oldSq; - + m_sdC = sd; m_sqC = sq; -- cgit