Saturation initialization implementation start

2 files changed, 68 insertions, 17 deletions

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<double> eq0 = data.terminalVoltage + std::complex<double>(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<double>(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<double>(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;
diff --git a/Project/Project.mk b/Project/Project.mk
index 922d62d..f3bd51d 100644
--- a/Project/Project.mk
+++ b/Project/Project.mk
@@ -13,7 +13,7 @@ CurrentFileName        :=
 CurrentFilePath        :=
 CurrentFileFullPath    :=
 User                   :=NDSE-69
-Date                   :=19/09/2017
+Date                   :=20/09/2017
 CodeLitePath           :="C:/Program Files/CodeLite"
 LinkerName             :=C:/TDM-GCC-64/bin/g++.exe
 SharedObjectLinkerName :=C:/TDM-GCC-64/bin/g++.exe -shared -fPIC