Efficient Prototyping of a Field-Programmable Gate Array-Based Real-Time Model of a Modular Multilevel Converter

: Field-programmable gate array (FPGA)-based real-time simulation plays a crucial role in testing power–electronic dominated systems with the formation of controller hardware-in-the-loop (CHIL) or power hardware-in-the-loop (PHIL). This work describes an efficient implementation of computation time and resource usage in the FPGA-based study of a modular multilevel converter (MMC) with detailed electromagnetic transients. The proposed modeling technique can be used in continuous control mode (CCM) and discontinuous control mode (DCM) for high-switching frequency semiconductor technologies. An FPGA-based designed solver structure is also presented to take advantage of the parallel features of FPGAs to achieve an ultra-fast calculation speed. In addition, two different switch modeling techniques are discussed with a five-level MMC case study. Experimental results on the NI PXIe platform show the feasibility of the proposed implementation, and a time step of 100 nanoseconds is achieved.


Introduction
Modular multilevel converters (MMCs) are widely applied in power transmission and conversion in power grids and electrified transportation.However, testing MMCs with high power is a great danger due to the high voltage and current involved.In addition, the prototype is prone to be damaged without precautions.On the other hand, hardwarein-the-loop simulation (HiLs) serves as an effective form of testing and validation of the prototype by simulating the power electronic system using a digital real-time simulator (DRTS) [1].This can minimize the developing period and reduce the danger of damaging and destroying prototypes [2].However, the design of HiLs depends on the mathematical model of the MMC, the calculation time of which will affect the stability of the realtime system.
In a power-electronic-dominated system simulation, the choice of the time step relates to the typical application and the implementation hardware.A relatively low time step is required to maintain numerical stability for studying power systems or electrical machines.A time step ranging from 50 µs to 200 µs is necessary to consider the effect of the third or fifth harmonic in the power system [3].This time step can provide acceptable results for electromagnetic transients up to 1 kHz.This can be implemented on a simulator using CPU/DSP as the core computational engine [4].A power system involving multiple power electronic systems requires a time below 50 µs to perform the transient analysis of power systems.
For time steps below microseconds, FPGAs have to be involved [5].This time step presents a detailed representation of electromagnetic transients and high frequencies inside power electronic systems.A higher-order numerical solver can achieve greater precision.The transient process of the switch devices, similar to the IGBT, can be simulated [6].Compared with processor-based studies, the FPGA-based simulator can affect the system's behavior by sampling the controller's signals with an ultra-high resolution.This gives an acceptable accuracy with an interpolation method compensating for the internal switch events [7].
From this discussion, we can conclude that an FPGA simulator is essential for achieving a time step with hundreds of nanoseconds.In addition, the discrete solver and the modeling of the power switch are two main parts.In discrete solver design, paper [6] presents an FPGA-based DRTS environment for analyzing the electromagnetic transients of power systems, including multiple power electronic converters.Different subsystems are divided by utilizing the Switching-Network Partitioning (SNP) method.Using the parallel calculation of FPGAs enables sub-microsecond simulation time steps.Paper [7] utilizes the LabVIEW Platform 2021 and implements the power converters model using LabVIEW FPGAs.LabVIEW provides a single-cycle time loop structure and achieves a time step of hundreds of nanoseconds.Paper [8] presents a general implementation of the FPGA-based simulation of photovoltaic applications.After representing the system with ODEs using Xilinx system generator tools, a time step below 1 microsecond is obtained.
However, this solver is restricted by sequential calculation, and parallel solving can improve the solving speed [9,10].Paper [11] proposes a similar approach to real-time simulation using a multi-rate structure.By modeling the power electronic system using different time resolutions, the solving is not sequential, and all subsystems can be solved simultaneously.Paper [12] proposes apredictor-corrector parallel solver for FPGA implementation.An independent solving process is obtained after separating the solving relationship between the predictor and corrector calculation.After implementing LabVIEW FPGAs, the time step is significantly reduced to 100 nanoseconds.The other decoupling method is a multi-rate simulation with a different solver [13,14].A parallel process can be achieved by partitioning the whole system using different time steps.However, the latency among different subsystems and the synchrony of the subsystem are two potential problems.
The other key factor that affects the calculation speed is the judgment of the switches.In general, the ideal switch function model [15], the R on /R off model [16], and the associated discrete circuit model [17] are three commonly used models for system-level simulations with a detailed representation of the switch steady-state effect.Although the transient state of the switch can further improve the system accuracy, the model size is largely limited due to the hardware resource of FPGAs [18].In system-level simulations, similar to the DCM condition, an iteration is required to maintain the stability of the real-time model.Paper [19,20] utilizes the iteration process to find the value when the current reaches zero.The maximum iteration times are defined to avoid overrun under the real-time constraint.On the other hand, paper [21] proposes a zero-crossing method seeking the current value when the MMC's bridge is in the blocked mode.The zero-crossing method forces the current value to zero once the transmission from CCM to DCM is captured.This avoids the complex iteration process and can be used for CCM and DCM statutes.
This paper aims to develop a highly efficient MMC model using FPGAs.The rest of this article is organized as follows: Section 2 describes the FPGA solver design, including the ODE-solving and switch judgment processes.Section 3 presents the system model of MMCs using the ideal and R on /R off switch models.Section 4 offers the conclusion.

FPGA Solver Design
In this section, the solver of a power electronic system is illustrated.Three factors are involved.One is the discrete integration solver.The second part is the judgment of the switch status.The last is the circuit implementation on FPGAs.

ODE Solver Design
In power electronic simulations, the ordinary differential equation (ODE) form in Equation ( 1) can be used to describe the system behavior.
where u is the input vector, x is the state vector, f (t n , x n ) is the derivative function, and f (t n , x n ) = Ax + Bu, y is the output vector.A, B, C, D is the coefficient matrix containing the system parameters.
In general, the solving of (1) involves the integration method.One commonly used method is the Backward Euler (BE) method.In the n-th time step, (1) is rewritten using BE as Equation ( 2) [22], To solve x(n), two methods have been used in the literature.One is the iteration method.The other is the direct solution [23].The error should be defined in the iteration method to end the calculation loop in each time step.The maximum iteration times are restricted so that real-time constraints can be met.One type of iteration is the improved Euler method, with the formation in Equation ( 3) where x(n) ′ is the value first calculated with the forward Euler method, and it is used to approximate the value of x(n).Then, in the backward Euler formation, x(n) is replaced with x(n) ′ to finish the calculation.This significantly saves the calculation time without much iteration time.The other method is the direct solution.In this method, Equation ( 2) is deduced from Equation (4), Equation ( 3) is a direct solution to Equation (2).Since (1 − h•A) −1 exists in Equation (4), x(n) has value with the condition that (1 − h•A) −1 ̸ = 0. Compared with the iteration method in Equation (2), the direct solution requires more mathematical power and calculation time to solve the inverse matrix of (1 − h•A).So, usually, the iteration process with limited iteration time has the advantage of fast calculation.However, iteration is a sequential process, the calculation time of which can be further accelerated with a parallel process.
As shown in Figure 1, the parallel structure calculates the value xn+2 at time point t n with a time step of 2h.At the same time, x n is calculated with the value xn obtained from the time step of t n−2 .The value of xn is calculated and known at the beginning of time step t n .Thus, the calculation of x n does not need to wait for the calculation of xn at the time step of t n .Since the solving of xn+2 and x n is independent, they can be calculated at the same time.This method can be effectively used when the simulation time step is relatively small.
The parallel solving process can be written as shown in Equation ( 5), The parallel solving process can be written as shown in Equation ( 5), Compared with Equation (4), the solving process is independent.The tion of Equation ( 5) will have a speed advantage.

Switch Status Update and DCM Modeling
One of the essential elements in the power electronic system is the po ductor device.The number of possible topologies in one power electronic s to the number of switches.If N is the switch's total number, the potential pow system's combination is 2 .This time-varying feature of the power electro one of the biggest challenges in modeling.For instance, if  = 10, the siz required to store all the statuses of the circuit elements can be enormous.
For non-controlled switches such as diode elements, its state depends  as input for the current calculation step.In Table 1, the three basic sub of switches are summarized.In the parallel calculation structure, the calcu has to be associated with the inductor connected to it.When the diode is in an insulated gate bipolar transistor (IGBT), its status will also be decided by current as long as the drive signal is equal to zero.The switch state update co value of the predictor or corrector is calculated.

Switch Status Update and DCM Modeling
One of the essential elements in the power electronic system is the power semiconductor device.The number of possible topologies in one power electronic system relates to the number of switches.If N is the switch's total number, the potential power electronic system's combination is 2 N .This time-varying feature of the power electronic system is one of the biggest challenges in modeling.For instance, if N = 10, the size of the data required to store all the statuses of the circuit elements can be enormous.
For non-controlled switches such as diode elements, its state depends on i switch or V switch as input for the current calculation step.In Table 1, the three basic subsystem types of switches are summarized.In the parallel calculation structure, the calculation of I switch has to be associated with the inductor connected to it.When the diode is in parallel with an insulated gate bipolar transistor (IGBT), its status will also be decided by the injection current as long as the drive signal is equal to zero.The switch state update comes after the value of the predictor or corrector is calculated.The parallel solving process can be written as shown in Equation ( 5), Compared with Equation (4), the solving process is independent.The implementation of Equation ( 5) will have a speed advantage.

Switch Status Update and DCM Modeling
One of the essential elements in the power electronic system is the power semiconductor device.The number of possible topologies in one power electronic system relates to the number of switches.If N is the switch's total number, the potential power electronic system's combination is 2 .This time-varying feature of the power electronic system is one of the biggest challenges in modeling.For instance, if  = 10, the size of the data required to store all the statuses of the circuit elements can be enormous.
For non-controlled switches such as diode elements, its state depends on  or  as input for the current calculation step.In Table 1, the three basic subsystem types of switches are summarized.In the parallel calculation structure, the calculation of Iswitch has to be associated with the inductor connected to it.When the diode is in parallel with an insulated gate bipolar transistor (IGBT), its status will also be decided by the injection current as long as the drive signal is equal to zero.The switch state update comes after the value of the predictor or corrector is calculated.

Subsystem
Type The system modeling separates the switch's status from the inductor or the capacitor.The circuit partitioning with the switch is shown in Figure 2, with the formation of a towport subsystem.In one time step, since the current flowing through the inductance or the voltage between the capacitance is constant, the connected capacitor is treated as the voltage source, and the inductance is treated as the current source.After obtaining the switch x n+3 x n x n+2 x n+1 x n+3 The parallel solving process can be written as shown in Equation ( 5), Compared with Equation ( 4), the solving process is independent.The implementation of Equation ( 5) will have a speed advantage.

Switch Status Update and DCM Modeling
One of the essential elements in the power electronic system is the power semiconductor device.The number of possible topologies in one power electronic system relates to the number of switches.If N is the switch's total number, the potential power electronic system's combination is 2 .This time-varying feature of the power electronic system is one of the biggest challenges in modeling.For instance, if  = 10, the size of the data required to store all the statuses of the circuit elements can be enormous.
For non-controlled switches such as diode elements, its state depends on  or  as input for the current calculation step.In Table 1, the three basic subsystem types of switches are summarized.In the parallel calculation structure, the calculation of Iswitch has to be associated with the inductor connected to it.When the diode is in parallel with an insulated gate bipolar transistor (IGBT), its status will also be decided by the injection current as long as the drive signal is equal to zero.The switch state update comes after the value of the predictor or corrector is calculated.

Subsystem
Type The system modeling separates the switch's status from the inductor or the capacitor.The circuit partitioning with the switch is shown in Figure 2, with the formation of a towport subsystem.In one time step, since the current flowing through the inductance or the voltage between the capacitance is constant, the connected capacitor is treated as the voltage source, and the inductance is treated as the current source.After obtaining the switch The parallel solving process can be written as shown in Equation ( 5), Compared with Equation ( 4), the solving process is independent.The implementation of Equation ( 5) will have a speed advantage.

Switch Status Update and DCM Modeling
One of the essential elements in the power electronic system is the power semiconductor device.The number of possible topologies in one power electronic system relates to the number of switches.If N is the switch's total number, the potential power electronic system's combination is 2 .This time-varying feature of the power electronic system is one of the biggest challenges in modeling.For instance, if  = 10, the size of the data required to store all the statuses of the circuit elements can be enormous.
For non-controlled switches such as diode elements, its state depends on  or  as input for the current calculation step.In Table 1, the three basic subsystem types of switches are summarized.In the parallel calculation structure, the calculation of Iswitch has to be associated with the inductor connected to it.When the diode is in parallel with an insulated gate bipolar transistor (IGBT), its status will also be decided by the injection current as long as the drive signal is equal to zero.The switch state update comes after the value of the predictor or corrector is calculated.

Subsystem
Type The system modeling separates the switch's status from the inductor or the capacitor.The circuit partitioning with the switch is shown in Figure 2, with the formation of a towport subsystem.In one time step, since the current flowing through the inductance or the voltage between the capacitance is constant, the connected capacitor is treated as the voltage source, and the inductance is treated as the current source.After obtaining the switch The system modeling separates the switch's status from the inductor or the capacitor.The circuit partitioning with the switch is shown in Figure 2, with the formation of a tow-port subsystem.In one time step, since the current flowing through the inductance or the voltage between the capacitance is constant, the connected capacitor is treated as the voltage source, and the inductance is treated as the current source.After obtaining the switch statutes from Table 1, the calculation of the unknown value [V out , I in ] is deduced from Equation (6).8) and ( 5), the calculation of the corrector  and Equations ( 9) and (10), respectively.
Once  ≤  is met, the next step is to obtain the n  ( ,  ).However,  ( ,  ) and  ( ,  ) are different.Different will be involved, which will require more hardware resources.
For a bridge with two switches,  ( ,  ) indicates that the switc is ,  or [,  ],  ( ,  ) is the DCM indicating that the ,  .With the defined boundary  , here, we define the poin  ( ,  ) to  ( ,  ).From Table 1, the value of I switch decides the switch statutes and acts as the guard function describing the time-varying feature in one bridge.The value of I switch is calculated using Equation (5).Three different derivative functions can be used to describe the timevarying topology in both DCM and CCM.A basic model is shown in Equation (7).
where g(t n−1 , x n−1 ) is the guard function, f 1 (t n , x n ) is the derivative function of CCM, and f 2 (t n , x n ) are the derivative functions of DCM.When the power electronic system translates from CCM statues to DCM statues, the zero-crossing function has to find the zero point t r meeting g(t r , x n−1 ) = 0. To reduce the calculation burden in simulation, g(t r , x n−1 ) = 0 can be approximated with a defined neighborhood |g(t n−1 , x n−1 )| ≤ α.Thus, (3) can be rewritten as Combining ( 8) and ( 5), the calculation of the corrector x n and x p t+h are shown in Equations ( 9) and (10), respectively.
Once |I switch p n | ≤ α is met, the next step is to obtain the numerical value of f 2 (t n , x n ).However, f 2 (t n , x n ) and f 1 (t n , x n ) are different.Different iteration equations will be in- volved, which will require more hardware resources.
For a bridge with two switches, With the defined boundary α, here, we define the point translating from When switch [S1, S2] is [OFF, OFF], I switch ≈ 0. Here, the value of I switch during the DCM conditions is limited to a small value β that can be neglected.
Energies 2024, 17, 591 6 of 15 At the time point t n , ( 12) and ( 13) are the calculations dealing with the predictor and the corrector.The function f lag(t n ) is used to monitor whether the system is in a chattering situation.Once chattering is detected, the value x n and x p n+1 in the chattering zone will be assigned to β or −β using ( 16) and (17).Thus, in the next simulation step t n+1 , since the sign of x n and x p n+1 are different, and the DCM will be recognized as the point crossing zero.Therefore, the following point will continue to be recognized as the DCM points in both the predictor and the corrector, but the value of x n+1 and x p n+2 are set to −β and β, which are negligible.

FPGA Implementation Structure
Figure 3 shows the FPGA implementation structure using the proposed zero-crossing and predictor-corrector methods.
At the time point   , ( 12) and ( 13) are the calculations dealing with the predictor and the corrector.The function (  ) is used to monitor whether the system is in a chattering situation.Once chattering is detected, the value   and  +1  in the chattering zone will be assigned to β or −β using ( 16) and (17).Thus, in the next simulation step  +1 , since the sign of   and  +1  are different, and the DCM will be recognized as the point crossing zero.Therefore, the following point will continue to be recognized as the DCM points in both the predictor and the corrector, but the value of  +1 and  +2  are set to −β and β, which are negligible.

FPGA Implementation Structure
Figure 3 shows the FPGA implementation structure using the proposed zero-crossing and predictor-corrector methods.

Gate Signals
Zero-crossing

Register
x n(t+h)

Register
x n(t+2h)  As shown in Figure 3, the system can be divided into independent subsystems by using a current source, voltage source, and two ports as a division.In each predictor calculation or corrector calculation process, the switch calculation and circuit solving have a serial connection relationship.At the time step   , the solver first uses Equation (5) to calculate  ̂(+2) and  () after sampling the information of gate signals. ̂(+2ℎ) is a value after two steps, and  (+ℎ) is the calculated system status.Because each variable in  ̂(+2) and  () is independent,  ̂(+2) and  () can be calculated at the same time.Compared with the forward Euler Method, the order in the parallel structure is improved to 2 [12,24].

ˆ+ n X
These two values are further used to decide the system's status in CCM or DCM using Equations ( 11)- (13).Furthermore, in Equations ( 11)-( 13), the FPGA resources utilized in the zero-crossing process only involve the comparison unit, and no math calculation is used.As a result, the calculation speed is almost unaffected by the zero-crossing unit insertion.As shown in Figure 3, the system can be divided into independent subsystems by using a current source, voltage source, and two ports as a division.In each predictor calculation or corrector calculation process, the switch calculation and circuit solving have a serial connection relationship.At the time step t n , the solver first uses Equation (5) to calculate x(n+2) and x (n) after sampling the information of gate signals.x(n+2h) is a value after two steps, and x (n+h) is the calculated system status.Because each variable in x(n+2) and x (n) is independent, x(n+2) and x (n) can be calculated at the same time.Compared with the forward Euler Method, the order in the parallel structure is improved to 2 [12,24].
These two values are further used to decide the system's status in CCM or DCM using Equations ( 11)- (13).Furthermore, in Equations ( 11)-( 13), the FPGA resources utilized in the zero-crossing process only involve the comparison unit, and no math calculation is used.As a result, the calculation speed is almost unaffected by the zero-crossing unit insertion.

System Model and Experimental Results
In this section, we show the effectiveness and the general results of the proposed method.The other case study of a five-level MMC is shown in Figure 4. Figure 4 divides the system into two parts: N1 and N2.N1 is the subsystem with inductance and capacitance.The N2 subsystem consists of the switches, which are made of a large number of power electronic submodules (SMs).The related simulation parameters are shown in Table 2.
gies 2024, 17, x FOR PEER REVIEW capacitance.The N2 subsystem consists of the switches, which are m of power electronic submodules (SMs).The related simulation par Table 2.

Real-Time Simulation with ISF Model
The switch device of the SM in the upper bridge is shown in F S2 operate in the dead zone, the switch function can be obtained w

Real-Time Simulation with ISF Model
The switch device of the SM in the upper bridge is shown in Figure 5a.When S1 and S2 operate in the dead zone, the switch function can be obtained with (14).
where p = A1, B1, and C1, which represent phases A, B, and C, respectively, and j is the number of SM in the upper bridge.The output nodal voltage   and the output current   can be calculated (15).
Thus, the current and voltage relationship in the upper bridge can be calculated ( 16) where  _0 =  1 .With the same method, the current and voltage relationship in the lower bridge B2, C2) can be calculated with ( 17)- (19).
where  = 2, 2, and 2, which represent phases A, B, and C, respectively, j is the n ber of SM in the upper bridge, and  _0 =  2 .After obtaining  1_ and  2_ , the voltage relationship inside the N1 subsys can be expressed with KVL as shown in (20), Using (20),   ,  1 , and  2 can be obtained.Then, the status of the inducta The output nodal voltage V out and the output current I out can be calculated with (15).
Thus, the current and voltage relationship in the upper bridge can be calculated with ( 16) where With the same method, the current and voltage relationship in the lower bridge (A2, B2, C2) can be calculated with ( 17)- (19).
where q = A2, B2, and C2, which represent phases A, B, and C, respectively, j is the number of SM in the upper bridge, and I out_0 = i Lp2 .
After obtaining V p1_N and V p2_N , the voltage relationship inside the N1 subsystem can be expressed with KVL as shown in (20), Using (20), V NO , V Lp1 , and V Lp2 can be obtained.Then, the status of the inductances can be calculated with (21): Its discretization formulation with the predictor-corrector method can be expressed as (22), For the capacitance in each SM, the mathematical model can be described as (23): Its discretization formulation with the predictor-corrector method can be expressed as (24), With an FPGA Kintex-7 XC7K410T embedded in the National Instrument (NI, Austin, TX, USA) PXIe-7975R FlexRIO PX, the model is implemented in the Express FPGA module with the structure shown in Figure 6.Open Loop Dead-time Eq. ( 18) Eq. ( 28)-( Eq. ( 21)-( 23) Ideal Switch Function Model Eq. ( 25)-( 26) Ron/Roff Model The Implementation of MMC model on 7975R After building on the FPGAs, the hardware resource utilization which compares the P-C method with the zero-crossing method, th the zero-crossing method, and the forward Euler method.Altho (with/without the zero-crossing unit) has the same calculation spee The model is implemented using the NI platform, as shown in Figure 6.The FlexRIO PXI platform contains a Kintex-7 XC7K410T FPGA.The single-cycle time loop (SCTL) allows all functions inside the loop to execute within a single tick [25,26].The implementation of FPGAs has four sequential steps.Firstly, based on ( 14) and ( 17), the driving signals and the predictor values i p li1(n+1) and i p li2(n+1) determine the switch status in the correction process.Meanwhile, the driving signals and the corrector values i li1(n+1) and i li2(n+1) determine the switch status in the prediction calculation process.Secondly, we substitute the corresponding switch function into (15)-( 16) and ( 18)-( 19), so the nodal voltages and the branch currents are calculated.Then, the port voltage and current are updated using (20).Finally, ( 22) and ( 24) are computed in parallel.
After building on the FPGAs, the hardware resource utilization is shown in Figure 7, which compares the P-C method with the zero-crossing method, the P-C method without the zero-crossing method, and the forward Euler method.Although the P-C method (with/without the zero-crossing unit) has the same calculation speed (−100 ns) as the forward Euler method, the resource of the DSP48s is doubled in the P-C method.Furthermore, it can be noted that the zero-crossing process only increases by 0.5% in the Slice L.U.T.s and 0.3% in the total slices.In the meantime, the zero-crossing unit does not cause an increase in the time step or the DSP48s utilization.Open Loop Dead-time Eq. ( 18) Eq. ( 29)-( Eq. ( 21)-( 23) Ideal Switch Function Model Eq. ( 29)-(32) Ron/Roff Model After building on the FPGAs, the hardware resource utilization is shown in Figure 7, which compares the P-C method with the zero-crossing method, the P-C method without the zero-crossing method, and the forward Euler method.Although the P-C method (with/without the zero-crossing unit) has the same calculation speed (−100 ns) as the forward Euler method, the resource of the DSP48s is doubled in the P-C method.Furthermore, it can be noted that the zero-crossing process only increases by 0.5% in the Slice L.U.T.s and 0.3% in the total slices.In the meantime, the zero-crossing unit does not cause an increase in the time step or the DSP48s utilization.With a simulation step of 100 ns, the simulation results of  1 are shown in Figure 8. Due to a long dead time of 30 µs in the PWM pulses, the current  1 is not continuous.In the zero-crossing setup, the value of β is 1 × 10 −7 , and the value of α is 2 × 10 −7 .Figure 8a compares the Simulink results with the proposed zero-crossing method.In the Simulink results, when the current is close to zero and the MMC is in the blocked mode, the current is calculated to a value that relates to the open voltage of the IGBT device.However, the calculation of this current is a complex process.The Simulink model utilizes iteration and zero-crossing methods to seek this point.This is an offline process requiring much calculation time to finish.For a real-time simulation, the calculation burden is too heavy to seek the point crossing zero.With a simulation step of 100 ns, the simulation results of i La1 are shown in Figure 8. Due to a long dead time of 30 µs in the PWM pulses, the current i La1 is not continuous.In the zero-crossing setup, the value of β is 1 × 10 −7 , and the value of α is 2 × 10 −7 .Figure 8a compares the Simulink results with the proposed zero-crossing method.In the Simulink results, when the current is close to zero and the MMC is in the blocked mode, the current is calculated to a value that relates to the open voltage of the IGBT device.However, the calculation of this current is a complex process.The Simulink model utilizes iteration and zero-crossing methods to seek this point.This is an offline process requiring much calculation time to finish.For a real-time simulation, the calculation burden is too heavy to seek the point crossing zero.
Figure 8b shows the results of the method without using the zero-crossing method.Although the system can work normally, high oscillation exists when i La1 is close to zero, which causes a potential unstable problem in the real-time simulation.Furthermore, the high current oscillation also causes the blocked model's virtual losses.However, in the proposed MMC model, when the bridge is in the blocked mode and the current is around zero, the proposed method can regulate the current oscillation with −1 × 10 −7 and 1 × 10 −7 .It can also achieve a real-time simulation step of 100 nanoseconds.The proposed zero-crossing method can reduce the error and improve the accuracy when the current is around zero. Figure 8b shows the results of the method without usin Although the system can work normally, high oscillation exis which causes a potential unstable problem in the real-time s high current oscillation also causes the blocked model's virt proposed MMC model, when the bridge is in the blocked mo zero, the proposed method can regulate the current oscillatio It can also achieve a real-time simulation step of 100 nanose crossing method can reduce the error and improve the ac around zero.

Real-Time Simulation with the Ron/Roff Model
In the Ron/Roff model, the current flowing through the sw with the value of the current   .

Real-Time Simulation with the R on /R off Model
In the R on /R off model, the current flowing through the switch is calculated using (25) with the value of the current I in .
where the status of S1 p_j is updated by the direction of the current As shown in Figure 9, based on Thevetin's theorem, the SM is substituted with a voltage source V eq in a series connection with a resistance R eq .The value of V eq and R eq can be calculated with (29).
where p = A, B, and C, which represent phases A, B, and C, respectively, and j is the number of SM in the upper bridge.
rgies 2024, 17, x FOR PEER REVIEW After obtaining  1_ and  2_ , ( 25)-( 27) can be further used to age/current status of the system.The implementation structure on the LabVIEW is also based on Figure 10, where the SM update unit is calcu (29).The final FPGA resource utilization is shown in Figure 10.With R eq_p_j and V eq_p_j , the equivalent voltage V upper eq and V lower eq in the upper phase and the equivalent resistance R upper eq and R lower eq can be calculated with (27).
The current and voltage relationship in the upper and lower bridges can be calculated with ( 27) and (28). V After obtaining V p1_N and V p2_N , ( 25)-( 27) can be further used to calculate the voltage/current status of the system.The implementation structure on the FPGA board using LabVIEW is also based on Figure 10, where the SM update unit is calculated with (28) and (29).The final FPGA resource utilization is shown in Figure 10.With  __ and  __ , the equivalent voltage    and    in the upper phase and the equivalent resistance    and    can be calculated with (27).
The current and voltage relationship in the upper and lower bridges can be calculated with ( 27) and (28).After obtaining  1_ and  2_ , ( 25)-( 27) can be further used to calculate the voltage/current status of the system.The implementation structure on the FPGA board using LabVIEW is also based on Figure 10, where the SM update unit is calculated with (28) and (29).The final FPGA resource utilization is shown in Figure 10.Compared with Figures 7 and 10, the ISF and R on /R off models can achieve a calculation speed of 100 nanoseconds.And, because the calculation of (26) in the R on /R off model requires more math operations, the R on /R off model occupies more FPGA resources than the ISF model.
Figure 11 shows the simulation results of the current i La1 in the R on /R off model.In Figure 11a, the Simulink results reference the proposed zero-crossing method.The proposed method can regulate the current oscillation with −1 × 10 −7 and 1 × 10 −7 when the current is around zero.Compared with the P-C method without the zero-crossing method (Figure 11b), the proposed zero-crossing method increases the stability and the accuracy of the whole system simulation.

024, 17, x FOR PEER REVIEW
Figure 11 shows the simulation results of the current i L Figure 11a, the Simulink results reference the proposed zeroposed method can regulate the current oscillation with −1 × current is around zero.Compared with the P-C method withou (Figure 11b), the proposed zero-crossing method increases the of the whole system simulation.

Conclusions
This paper proposes an efficient real-time FPGA model o sentation of transient with a time step of 100 ns.Performanc switch models (ISF and Ron/Roff models) using a five-level M pared with the traditional FPGA solver, the proposed solving

Conclusions
This paper proposes an efficient real-time FPGA model of MMC to give a full representation of transient with a time step of 100 ns.Performance evaluations on different switch models (ISF and R on /R off models) using a five-level MMC are carried out.Compared with the traditional FPGA solver, the proposed solving structure keeps the same speed Energies 2024, 17, 591 14 of 15 as the FE method, and the total calculation time is only 100 ns under the STCL of the LabVIEW FPGA.When the PWM has a dead-time zone, and the switch operates in the blocked mode, the proposed method can regulate the current to −1 × 10 −7 and 1 × 10 −7 .Furthermore, the proposed method has a strong generality, which can be used in the ISF and R on /R off models.
The proposed solver for the power electronic system contains a parallel ODE-solving structure and the zero-crossing method.The solver can enhance the robustness by allowing the simulation of DCM and CCM with the same model.Implementing the FPGA will further reduce the simulation step under 100 ns with the detailed representation of the MMC electromagnetically transient.The ultra-fast speed advantage also allows for predicting the system behavior as fast as possible.Thus, it can be further used in the digital twin setup [27] and the model predictive control.

Figure 3 .
Figure 3. Power electronic simulation with parallel calculation.

Figure 3 .
Figure 3. Power electronic simulation with parallel calculation.

Figure 4 .
Figure 4. Topology of a five-level MMC.

Figure 4 .
Figure 4. Topology of a five-level MMC.

Figure 5 .
Figure 5.The system partitioning.(a) SM is in the upper bridge; (b) SM is in the lower bridge (c) the N1 and N2 subsystems.

Figure 5 .
Figure 5.The system partitioning.(a) SM is in the upper bridge; (b) SM is in the lower bridge; and (c) the N1 and N2 subsystems.

2 Figure 8 .
Figure 8.Current  1 in the ISF model: (a) comparison with the Sim ison with the non-zero-crossing method.

Figure 8 .
Figure 8.Current i La1 in the ISF model: (a) comparison with the Simulink results and (b) comparison with the non-zero-crossing method.

Figure 9 . 4 𝑗=1
Figure 9.The system partitioning: (a) SM is in the upper bridge; (b) SM is in (c) the N1 subsystem and N2 subsystem.

Figure 9 .
Figure 9.The system partitioning: (a) SM is in the upper bridge; (b) SM is in the lower bridge; and (c) the N1 subsystem and N2 subsystem.

Figure 9 .
Figure 9.The system partitioning: (a) SM is in the upper bridge; (b) SM is in the lower bridge; and (c) the N1 subsystem and N2 subsystem.

2 Figure 11 .
Figure 11.Current  1 in the Ron/Roff model: (a) comparison with comparison with the non-zero-crossing method.

Figure 11 .
Figure 11.Current i La1 in the R on /R off model: (a) comparison with the Simulink results and (b) comparison with the non-zero-crossing method.

Table 1 ,
(5) value of Iswitch decides the switch statutes and acts tion describing the time-varying feature in one bridge.The value of Isw ing Equation(5).Three different derivative functions can be used to varying topology in both DCM and CCM.A basic model is shown in E where ( ,  ) is the guard function,  ( ,  ) is the derivative and  ( ,  ) are the derivative functions of DCM.When the powe translates from CCM statues to DCM statues, the zero-crossing funct zero point  meeting ( ,  ) = 0 .To reduce the calculation bur ( ,  ) = 0 can be approximated with a defined neighborhood Thus, (3) can be rewritten as   =  ( ,  )  |( ,  )| > α  ( ,  )  |( ,  )| ≤ α Combining (