OpenVVVF · saeedamooei · Jul 31, 2026
diff --git a/Assets/NodeTemplates/control.mpcc/inline.cpp b/Assets/NodeTemplates/control.mpcc/inline.cpp
@@ -0,0 +1,129 @@
+/* FCS-MPCC node — Zhang et al. (2017) baseline + improved modes.
+ * Technical reference: Y. Zhang, D. Xu, J. Liu, S. Gao, and W. Xu,
+ * "Performance Improvement of Model-Predictive Current Control of
+ * Permanent Magnet Synchronous Motor Drives," IEEE Transactions on
+ * Industry Applications, vol. 53, no. 4, pp. 3683-3695, July/August 2017.
+ * DOI: 10.1109/TIA.2017.2690998.
+ * Mode: 0=ConventionalOneStep, 1=DelayCompensated, 2=BackEMFCompensated, 3=OptimalDutyCycle
+ * Fixed-size state only; no heap allocation in update. */
+
+static float mpcc_prev_sa = 0.0f;
+static float mpcc_prev_sb = 0.0f;
+static float mpcc_prev_sc = 0.0f;
+static float mpcc_u_alpha_prev = 0.0f;
+static float mpcc_u_beta_prev = 0.0f;
+static float mpcc_id_prev = 0.0f;
+static float mpcc_iq_prev = 0.0f;
+
+const float ts = Ts;
+const float rs = Rs;
+const float ld = Ld;
+const float lq = Lq;
+const float psi_f = PsiF;
+const float i_base = (I_Base > 0.0f) ? I_Base : 10.0f;
+const float i_max = (I_Max > 0.0f) ? I_Max : 30.0f;
+const float vdc = V_Dc.in(au::volts);
+const float id = I_D.in(au::amperes);
+const float iq = I_Q.in(au::amperes);
+const float id_ref = I_D_Ref.in(au::amperes);
+const float iq_ref = I_Q_Ref.in(au::amperes);
+const float theta_e = Theta_E;
+const float omega_e = Omega_E.in(au::radians_per_second);
+const bool enable = Enable > 0.5f;
+
+if (!enable || !(ts > 0.0f) || !(vdc > 0.0f) || ld <= 0.0f || lq <= 0.0f) {
+    S_A = mpcc_prev_sa;
+    S_B = mpcc_prev_sb;
+    S_C = mpcc_prev_sc;
+} else {
+    const float two_thirds = 2.0f / 3.0f;
+    const float inv_sqrt3 = 0.57735026919f;
+    const float active_vectors[6][2] = {
+        {two_thirds, 0.0f},
+        {1.0f / 3.0f, inv_sqrt3},
+        {-1.0f / 3.0f, inv_sqrt3},
+        {-two_thirds, 0.0f},
+        {-1.0f / 3.0f, -inv_sqrt3},
+        {1.0f / 3.0f, -inv_sqrt3}
+    };
+    const int switch_bits[8][3] = {
+        {0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,1,1},{0,0,1},{1,0,1},{1,1,1}
+    };
+
+    float id_base = id;
+    float iq_base = iq;
+    if (Mode >= 1.0f) {
+        const float di_d_sp = (mpcc_u_alpha_prev - rs * id) / lq;
+        const float di_q_sp = (mpcc_u_beta_prev - rs * iq) / lq;
+        const float id_sp = id + ts * di_d_sp;
+        const float iq_sp = iq + ts * di_q_sp;
+        id_base = id_sp + (-rs * (id_sp - id) * ts) / (2.0f * lq);
+        iq_base = iq_sp + (-rs * (iq_sp - iq) * ts) / (2.0f * lq);
+    }
+
+    float best_cost = 1.0e30f;
+    int best_idx = 0;
+    float best_pred_id = id;
+    float best_pred_iq = iq;
+    float best_valpha = 0.0f;
+    float best_vbeta = 0.0f;
+
+    for (int idx = 0; idx < 8; ++idx) {
+        const float sa = static_cast<float>(switch_bits[idx][0]);
+        const float sb = static_cast<float>(switch_bits[idx][1]);
+        const float sc = static_cast<float>(switch_bits[idx][2]);
+        const float valpha = two_thirds * vdc * (sa - 0.5f * (sb + sc));
+        const float vbeta = vdc * inv_sqrt3 * (sb - sc);
+
+        const float cos_t = cosf(theta_e);
+        const float sin_t = sinf(theta_e);
+        const float vd = valpha * cos_t + vbeta * sin_t;
+        const float vq = -valpha * sin_t + vbeta * cos_t;
+
+        const float pred_id = id_base + (ts / ld) * (vd - rs * id_base + omega_e * lq * iq_base);
+        const float pred_iq = iq_base + (ts / lq) * (vq - rs * iq_base - omega_e * ld * id_base - omega_e * psi_f);
+
+        float cost = ((id_ref - pred_id) / i_base) * ((id_ref - pred_id) / i_base)
+                   + ((iq_ref - pred_iq) / i_base) * ((iq_ref - pred_iq) / i_base);
+        const float imag = sqrtf(pred_id * pred_id + pred_iq * pred_iq);
+        if (imag > i_max) cost += 1.0e6f;
+
+        const int trans = ((sa != mpcc_prev_sa) ? 1 : 0) + ((sb != mpcc_prev_sb) ? 1 : 0) + ((sc != mpcc_prev_sc) ? 1 : 0);
+        const int best_trans = ((switch_bits[best_idx][0] != static_cast<int>(mpcc_prev_sa)) ? 1 : 0)
+                             + ((switch_bits[best_idx][1] != static_cast<int>(mpcc_prev_sb)) ? 1 : 0)
+                             + ((switch_bits[best_idx][2] != static_cast<int>(mpcc_prev_sc)) ? 1 : 0);
+
+        if (cost < best_cost - 1.0e-9f
+            || (fabsf(cost - best_cost) <= 1.0e-9f && (trans < best_trans || (trans == best_trans && idx < best_idx)))) {
+            best_cost = cost;
+            best_idx = idx;
+            best_pred_id = pred_id;
+            best_pred_iq = pred_iq;
+            best_valpha = valpha;
+            best_vbeta = vbeta;
+        }
+        (void)active_vectors;
+    }
+
+    S_A = static_cast<float>(switch_bits[best_idx][0]);
+    S_B = static_cast<float>(switch_bits[best_idx][1]);
+    S_C = static_cast<float>(switch_bits[best_idx][2]);
+    State_Index = static_cast<float>(best_idx);
+    Pred_I_D = rte::Amperes(best_pred_id);
+    Pred_I_Q = rte::Amperes(best_pred_iq);
+    Cost = best_cost;
+    V_Alpha = rte::Volts(best_valpha);
+    V_Beta = rte::Volts(best_vbeta);
+    const float cos_t = cosf(theta_e);
+    const float sin_t = sinf(theta_e);
+    V_D = rte::Volts(best_valpha * cos_t + best_vbeta * sin_t);
+    V_Q = rte::Volts(-best_valpha * sin_t + best_vbeta * cos_t);
+
+    mpcc_prev_sa = S_A;
+    mpcc_prev_sb = S_B;
+    mpcc_prev_sc = S_C;
+    mpcc_u_alpha_prev = best_valpha;
+    mpcc_u_beta_prev = best_vbeta;
+    mpcc_id_prev = id;
+    mpcc_iq_prev = iq;
+}
diff --git a/Assets/NodeTemplates/control.mpcc/node.json b/Assets/NodeTemplates/control.mpcc/node.json
@@ -0,0 +1,41 @@
+{
+  "id": "control.mpcc",
+  "displayName": "FCS-MPCC (Three-Phase PMSM)",
+  "defaultName": "Mpcc",
+  "maxInstances": 1,
+  "isEntryPoint": false,
+  "domain": "tim_isr",
+  "inputPorts": [
+    { "name": "I_D", "direction": "input", "type": { "quantity": "current", "frame": "scalar", "dtype": "f32" } },
+    { "name": "I_Q", "direction": "input", "type": { "quantity": "current", "frame": "scalar", "dtype": "f32" } },
+    { "name": "I_D_Ref", "direction": "input", "type": { "quantity": "current", "frame": "scalar", "dtype": "f32" } },
+    { "name": "I_Q_Ref", "direction": "input", "type": { "quantity": "current", "frame": "scalar", "dtype": "f32" } },
+    { "name": "Theta_E", "direction": "input", "type": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" } },
+    { "name": "Omega_E", "direction": "input", "type": { "quantity": "angular_velocity", "frame": "scalar", "dtype": "f32" } },
+    { "name": "V_Dc", "direction": "input", "type": { "quantity": "voltage", "frame": "scalar", "dtype": "f32" } },
+    { "name": "Enable", "direction": "input", "type": { "quantity": "boolean", "frame": "scalar", "dtype": "f32" } }
+  ],
+  "outputPorts": [
+    { "name": "S_A", "direction": "output", "type": { "quantity": "boolean", "frame": "scalar", "dtype": "f32" } },
+    { "name": "S_B", "direction": "output", "type": { "quantity": "boolean", "frame": "scalar", "dtype": "f32" } },
+    { "name": "S_C", "direction": "output", "type": { "quantity": "boolean", "frame": "scalar", "dtype": "f32" } },
+    { "name": "State_Index", "direction": "output", "type": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" } },
+    { "name": "Pred_I_D", "direction": "output", "type": { "quantity": "current", "frame": "scalar", "dtype": "f32" } },
+    { "name": "Pred_I_Q", "direction": "output", "type": { "quantity": "current", "frame": "scalar", "dtype": "f32" } },
+    { "name": "Cost", "direction": "output", "type": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" } },
+    { "name": "V_Alpha", "direction": "output", "type": { "quantity": "voltage", "frame": "scalar", "dtype": "f32" } },
+    { "name": "V_Beta", "direction": "output", "type": { "quantity": "voltage", "frame": "scalar", "dtype": "f32" } },
+    { "name": "V_D", "direction": "output", "type": { "quantity": "voltage", "frame": "scalar", "dtype": "f32" } },
+    { "name": "V_Q", "direction": "output", "type": { "quantity": "voltage", "frame": "scalar", "dtype": "f32" } }
+  ],
+  "parameterTypes": {
+    "Ts": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" },
+    "Rs": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" },
+    "Ld": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" },
+    "Lq": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" },
+    "PsiF": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" },
+    "I_Base": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" },
+    "I_Max": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" },
+    "Mode": { "quantity": "dimensionless", "frame": "scalar", "dtype": "f32" }
+  }
+}
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -15,3 +15,5 @@ add_subdirectory(Lib/InverterProtocol)
 add_subdirectory(Source/RTECodeEmitter)
 add_subdirectory(Source/RTEFirmwareBuilder)
 add_subdirectory(Source/NodeGUI)
+
+add_subdirectory(Lib/Simulation)
diff --git a/Lib/Simulation/CMakeLists.txt b/Lib/Simulation/CMakeLists.txt
@@ -0,0 +1,40 @@
+cmake_minimum_required(VERSION 3.24)
+project(Simulation VERSION 0.1.0 LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+add_library(Simulation INTERFACE)
+target_include_directories(Simulation INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include)
+
+if(TARGET project_warnings)
+    target_link_libraries(Simulation INTERFACE project_warnings)
+endif()
+
+enable_testing()
+include(GoogleTest)
+
+function(add_sim_test name source)
+    add_executable(${name} ${source})
+    target_link_libraries(${name} PRIVATE Simulation GTest::gtest_main)
+    if(TARGET project_warnings)
+        target_link_libraries(${name} PRIVATE project_warnings)
+    endif()
+    gtest_discover_tests(${name})
+endfunction()
+
+add_sim_test(Simulation_inverter_tests tests/test_inverter.cpp)
+add_sim_test(Simulation_pmsm_tests tests/test_pmsm.cpp)
+add_sim_test(Simulation_mpcc_tests tests/test_mpcc.cpp)
+
+add_executable(mpcc_closed_loop simulation/mpcc_closed_loop.cpp)
+target_link_libraries(mpcc_closed_loop PRIVATE Simulation)
+if(TARGET project_warnings)
+    target_link_libraries(mpcc_closed_loop PRIVATE project_warnings)
+endif()
+
+add_executable(compare_mpcc_foc simulation/compare_mpcc_foc.cpp)
+target_link_libraries(compare_mpcc_foc PRIVATE Simulation)
+if(TARGET project_warnings)
+    target_link_libraries(compare_mpcc_foc PRIVATE project_warnings)
+endif()
diff --git a/Lib/Simulation/include/simulation/ClosedLoopSimulator.h b/Lib/Simulation/include/simulation/ClosedLoopSimulator.h
@@ -0,0 +1,176 @@
+#pragma once
+
+#include <chrono>
+#include <cmath>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include "simulation/MpccController.h"
+#include "simulation/PmsmPlant.h"
+#include "simulation/TwoLevelInverter.h"
+#include "simulation/Transforms.h"
+
+namespace simulation {
+
+struct SimulationSample {
+    double time = 0.0;
+    float ia = 0.0f;
+    float ib = 0.0f;
+    float ic = 0.0f;
+    float id = 0.0f;
+    float iq = 0.0f;
+    float id_ref = 0.0f;
+    float iq_ref = 0.0f;
+    float omega_m = 0.0f;
+    float omega_e = 0.0f;
+    float theta_m = 0.0f;
+    float theta_e = 0.0f;
+    float torque_em = 0.0f;
+    float load_torque = 0.0f;
+    float vdc = 0.0f;
+    bool sa = false;
+    bool sb = false;
+    bool sc = false;
+    int switching_state = 0;
+    float valpha = 0.0f;
+    float vbeta = 0.0f;
+    float vd = 0.0f;
+    float vq = 0.0f;
+    float predicted_id = 0.0f;
+    float predicted_iq = 0.0f;
+    float cost = 0.0f;
+    double controller_exec_us = 0.0;
+};
+
+struct Scenario {
+    std::string name;
+    double duration = 0.1;
+    float ts = 100e-6f;
+    float vdc = 540.0f;
+    float load_torque = 0.0f;
+    float id_ref = 0.0f;
+    float iq_ref = 10.0f;
+    bool use_speed_loop = false;
+    float speed_ref = 0.0f;
+    float speed_kp = 0.5f;
+    float speed_ki = 20.0f;
+    MPCCMode mode = MPCCMode::ConventionalOneStep;
+    float param_scale_rs = 1.0f;
+    float param_scale_ld = 1.0f;
+    float param_scale_lq = 1.0f;
+    float param_scale_psi = 1.0f;
+};
+
+class ClosedLoopSimulator {
+public:
+    ClosedLoopSimulator(PmsmParameters plant_params, MpccParameters ctrl_params)
+        : plant_(plant_params), controller_(ctrl_params) {}
+
+    std::vector<SimulationSample> run(const Scenario& scenario) {
+        MpccParameters ctrl = controller_.parameters();
+        ctrl.mode = scenario.mode;
+        ctrl.ts = scenario.ts;
+        controller_.parameters() = ctrl;
+
+        PmsmParameters plant_params = plant_.parameters();
+        plant_params.rs *= scenario.param_scale_rs;
+        plant_params.ld *= scenario.param_scale_ld;
+        plant_params.lq *= scenario.param_scale_lq;
+        plant_params.psi_f *= scenario.param_scale_psi;
+        plant_ = PmsmPlant(plant_params);
+        controller_ = MpccController(ctrl);
+        plant_.reset();
+
+        std::vector<SimulationSample> samples;
+        const int steps = static_cast<int>(scenario.duration / scenario.ts);
+        float speed_integral = 0.0f;
+
+        for (int k = 0; k < steps; ++k) {
+            const double t = k * scenario.ts;
+            const auto& st = plant_.state();
+
+            float id_ref = scenario.id_ref;
+            float iq_ref = scenario.iq_ref;
+
+            if (scenario.use_speed_loop) {
+                const float speed_error = scenario.speed_ref - st.omega_m;
+                speed_integral += speed_error * scenario.ts;
+                iq_ref = scenario.speed_kp * speed_error + scenario.speed_ki * speed_integral;
+            }
+
+            MpccInputs in;
+            in.id = st.id;
+            in.iq = st.iq;
+            in.id_ref = id_ref;
+            in.iq_ref = iq_ref;
+            in.theta_e = st.theta_e;
+            in.omega_e = st.omega_e;
+            in.vdc = scenario.vdc;
+            in.enable = true;
+
+            const auto t0 = std::chrono::steady_clock::now();
+            const MpccOutputs out = controller_.update(in);
+            const auto t1 = std::chrono::steady_clock::now();
+            const double exec_us =
+                std::chrono::duration<double, std::micro>(t1 - t0).count();
+
+            const float valpha = out.valpha;
+            const float vbeta = out.vbeta;
+            plant_.stepAlphaBeta(valpha, vbeta, scenario.load_torque, scenario.ts);
+
+            SimulationSample sample;
+            sample.time = t;
+            sample.ia = plant_.state().ia;
+            sample.ib = plant_.state().ib;
+            sample.ic = plant_.state().ic;
+            sample.id = plant_.state().id;
+            sample.iq = plant_.state().iq;
+            sample.id_ref = id_ref;
+            sample.iq_ref = iq_ref;
+            sample.omega_m = plant_.state().omega_m;
+            sample.omega_e = plant_.state().omega_e;
+            sample.theta_m = plant_.state().theta_m;
+            sample.theta_e = plant_.state().theta_e;
+            sample.torque_em = plant_.state().torque_em;
+            sample.load_torque = scenario.load_torque;
+            sample.vdc = scenario.vdc;
+            sample.sa = out.sa;
+            sample.sb = out.sb;
+            sample.sc = out.sc;
+            sample.switching_state = static_cast<int>(out.switching_state);
+            sample.valpha = out.valpha;
+            sample.vbeta = out.vbeta;
+            sample.vd = out.vd;
+            sample.vq = out.vq;
+            sample.predicted_id = out.predicted_id;
+            sample.predicted_iq = out.predicted_iq;
+            sample.cost = out.min_cost;
+            sample.controller_exec_us = exec_us;
+            samples.push_back(sample);
+        }
+        return samples;
+    }
+
+    static void writeCsv(const std::string& path, const std::vector<SimulationSample>& samples) {
+        std::ofstream out(path);
+        out << "time,ia,ib,ic,id,iq,id_reference,iq_reference,mechanical_speed,electrical_speed,"
+               "mechanical_angle,electrical_angle,electromagnetic_torque,load_torque,dc_link_voltage,"
+               "Sa,Sb,Sc,switching_state,v_alpha,v_beta,v_d,v_q,predicted_id,predicted_iq,cost,"
+               "controller_execution_time\n";
+        for (const auto& s : samples) {
+            out << s.time << ',' << s.ia << ',' << s.ib << ',' << s.ic << ',' << s.id << ',' << s.iq << ','
+                << s.id_ref << ',' << s.iq_ref << ',' << s.omega_m << ',' << s.omega_e << ',' << s.theta_m << ','
+                << s.theta_e << ',' << s.torque_em << ',' << s.load_torque << ',' << s.vdc << ',' << (s.sa ? 1 : 0)
+                << ',' << (s.sb ? 1 : 0) << ',' << (s.sc ? 1 : 0) << ',' << s.switching_state << ',' << s.valpha
+                << ',' << s.vbeta << ',' << s.vd << ',' << s.vq << ',' << s.predicted_id << ',' << s.predicted_iq
+                << ',' << s.cost << ',' << s.controller_exec_us << '\n';
+        }
+    }
+
+private:
+    PmsmPlant plant_;
+    MpccController controller_;
+};
+
+}  // namespace simulation