Next Article in Journal
Evaluation Method of Naturalistic Driving Behaviour for Shared-Electrical Car
Previous Article in Journal
Simplified Recovery Process for Resistive Solder Bond (RSB) Hotspots Caused by Poor Soldering of Crystalline Silicon Photovoltaic Modules Using Resin
Previous Article in Special Issue
Numbers, Please: Power- and Voltage-Related Indices in Control of a Turbine-Generator Set
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

FPGA-Based Real-Time Simulation of Dual-Port Submodule MMC–HVDC System

1
The Key Laboratory of Smart Grid of Ministry of Education, Tianjin University, Tianjin 300072, China
2
State Grid Handan Electric Power Company, Handan 056035, China
*
Author to whom correspondence should be addressed.
Energies 2022, 15(13), 4624; https://doi.org/10.3390/en15134624
Submission received: 20 May 2022 / Revised: 18 June 2022 / Accepted: 22 June 2022 / Published: 24 June 2022
(This article belongs to the Special Issue Challenge and Research Trends of Power System Simulation)

Abstract

:
Aiming at the problems of the high switch numbers, complex working mechanisms, and complicated real-time simulation of modular multilevel converters (MMCs) composed of dual-port submodules, in this study, we designed a unified equivalent model of the multiple submodule network by analyzing the combination of parallel submodules in the bridge arm. The proposed model decouples the submodules that do not affect each other in the subnetwork calculation process, thereby reducing the number of prestored parameters in the subnetwork simulation. In the Xilinx Virtex-7 FPGA VC709 (Xilinx Corporation, San Jose, CA, USA) development board, we replaced the inline computation combined with the prestorage of parameters with the proposed equivalent model to optimize the execution unit structure and redesigned the FPGA-Based Real-Time Digital Solver (FRTDS). Taking the P-FBSM-based MMC–HVDC system as the simulation object, we performed a real-time simulation with a step size of 10 μs, which verified the effectiveness of the proposed model and the improvement in the hardware. We compared the results with the offline MATLAB/Simulink simulation results to verify the accuracy of the simulation.

1. Introduction

Energy and environmental problems are becoming increasingly severe, and the energy structure is transforming toward clean and low-carbon energy. Flexible DC transmission technology has many advantages in long-distance power transmission and large-scale renewable energy grid integration and thus is playing an important role in energy transformation [1,2,3]. Owing to their scalability, low harmonic content, and strong fault ride-through capability, MMCs have gradually become the mainstream flexible HVDC converter structure. The control and protection system has an important function in the secure and reliable operation of the power grid. The hardware-in-loop (HIL) test based on real-time simulation is used to verify the correctness of the protection devices’ action and the effectiveness of the control strategy in HVDC flexible transmission [4,5,6]. As the transmission voltage level increases, the switching element groups sometimes need to be used in parallel. In this case, a new topology can be constructed with the same components to add new functions to the submodules. As such, researchers have designed an MMC topology that provides capacitor voltage self-balancing and bridge arm current-sharing [7,8,9]. The parallel full-bridge submodule (P-FBSM) topology [10] is a dual-port structure that reduces switching frequency and power loss and performs the sensorless operation of capacitor voltage, providing practicality and advantages in high-voltage MMC applications. The dual-port submodule has four outlets, and its current path is complex and changeable, requiring the Thevenin equivalent model, which is applicable to real-time simulations of the single-port submodule [11]. Previously proposed is the electromagnetic transient simulation model suitable for the topology of the dual-port submodule, but this method requires iterative calculation and involves a large number of matrix operations, so it is not suitable for real-time simulation due to high calculation cost. The many high-frequency interrupted power electronic devices used in MMCs require a small simulation step size for real-time simulation, but a real-time simulation model that can reflect the diverse working modes and complex current paths of dual-port submodules has not been designed. A hardware platform with strong computing power is also required for real-time simulations of the MMC–HVDC system.
At present, the mainstream real-time simulators used in the electrical engineering field include RTDS and RT-LAB [12,13,14]. RT-LAB integrates the dynamic system mathematical model established by MATLAB/Simulink and performs real-time simulations. It is mainly used for the simulation of a single-power electronic device and less so for system-level simulation. The RTDS hardware used for computing tasks includes several racks and multiple processors for parallel computing to increase the simulation scale. However, the underlying computing hardware is still a serial device in essence, which leads to limitations in small-step simulation. In addition, the existing platform has strict hardware configuration requirements, so the equipment cost is too high for most practical situations [15]. The field-programmable logic gate array (FPGA) has a fully configurable parallel hardware structure, with distributed memory and deep pipeline structures, and also is low cost and small. As such, FPGA has gradually become the main hardware used for the real-time simulation of power systems [16]. A specialized calculation module was constructed on an FPGA based on the mathematical model of the power system, which improved the real-time simulation speed [17,18,19]. However, the dedicated modules often remain in an idle state due to limitations of the calculation process, resulting in a waste of FPGA resources. Zhang et al. [20] designed a real-time digital solver (FRTDS) based on FPGA and an instruction stream. FRTDS breaks the design concept of designing specific hardware circuits for each function involved in the simulated objects and provides a new idea for the real-time simulation of electromagnetic transients in power systems. FRTDS was designed for traditional AC power systems; thus, when used for real-time simulation of MMC–HVDC systems, it needs to be modified according to the characteristics of the calculation process. To perform real-time simulations of the dual-port submodule in MMC–HVDC systems, in this study, we optimized both the simulation models and the simulation platform. Section 2 introduces the basic structure of the FPGA-based real-time digital solver (FRTDS), focusing on the operation mechanism of the general computing components; Section 3 describes the method of replacing computation with storage, which is often used for real-time simulations to reduce the computation amount. In Section 4, we propose an equivalent model of the PFBSM subnetwork from the perspective of storage capacity optimization and analyze the computation amount. In Section 5, we describe appropriate improvements to the original FRTDS according to the calculation characteristics of the MMC–HVDC system; in Section 6, we verify the effectiveness and accuracy of the proposed simulation model and design method by describing a study case of the MMC–HVDC system; Section 7 outlines the studies conclusions.

2. FPGA-Based Real-Time Digital Solver (FRTDS)

The FPGA-Based Real-Time Digital Solver (FRTDS), which uses instruction stream-driven architecture, is to design a high reusability calculation unit based on the simulation method, the data address and the operation type of the calculation unit are given by instruction to realize the calculation. FRTDS was designed to simulate traditional AC power systems, adopting the node voltage analysis method. Its calculation process and commonly used arithmetic operations are shown in Figure 1. If each calculation process is designed as a specific function module, it not only wastes FPGA resources but also increases the simulation calculation time. According to the typical arithmetic operations of the calculation process, when the operating components of an FRTDS are designed, the following seven arithmetic operations are defined as basic operations: Y = A × B, Y = A + B, Y = A/B, Y = A × B + C × D, Y = A/B, Y = A × B/C, and Y = A × B/C + D; all the other arithmetic operations can be split into these basic operations. For example, Y = A × B/C + D/E + F can be split into two basic operations: Y = A × B/C + Y1 and Y1 = 1 × D/E + F. To save the DSP resources in an FPGA, all basic operations are implemented with two adders, two multipliers, and one divider in the general operation component.
FRTDS describes the data memory address and the selected word of the basic operation involved in completing a certain operation in the specified instruction format and saves them in the program’s memory. The arithmetic component takes the instructions from the program memory at the specified operating frequency, and transmits them to the read, write, and selection controllers. The read controller obtains the operation data from the data memory through the multiport read-and-write circuit, the arithmetic unit performs calculation according to the operation formula determined by the select controller, and the write controller stores the calculation result in the data memory through the multiport read-and-write circuit. In each step, the arithmetic components continuously and cyclically execute the calculation process according to the instructions; the specific structure of the universal computing component is shown in Figure 2. To enable the arithmetic components to work in a pipeline, an instruction buffer queue is added to the read, write, and selection controllers to implement a delayed output of data memory addresses and selection words. Several arithmetic components can be formed into a microprocessor core with stronger computing power. The arithmetic components in the core communicate through shared memory, and the microprocessor cores communicate through the straight-through wires of the arithmetic components via message passing. When allowed by available hardware resources, as many microprocessor cores as possible should be built so that the FRTDS will have powerful parallel computing capabilities. The data exchange between the microprocessor cores adopts the way of “hand in hand and data pipeline”, which means that two adjacent microprocessor cores share data storage units while two non-adjacent microprocessor cores exchange data through the data pipeline.
When used for HIL testing of real equipment, the time when the control signals of the switching elements such as IGBTs sent by the controller reach the FRTDS is uncertain; therefore, these signals are first stored in the backup data memory. At the start of each simulation step, the backup data memory becomes the running data memory of the arithmetic components, and the external information is provided to the arithmetic components. The output data of the arithmetic components cannot be directly provided to the peripheral equipment. First, the output data are saved in the backup data memory, and at the beginning of each simulation step, the backup data memory is changed to the running data memory of the peripheral equipment to ensure that the voltage and current peripherals are used at the same time. The ping-pong operation of the data memory improves the working efficiency of the real-time simulation platform.
During the simulation calculations, the switching element is modeled as a binary resistance, and the nonlinear elements are calculated in a piecewise linear manner. The existence of a switching element and a nonlinear element leads to many simulation parameters having various values. To facilitate the selection of the current value of a multivalued parameter, all the values of a multivalued parameter are stored as an array, in which each element represents the value of a multivalued parameter for a certain case. The information describing the multivalued parameters is stored in the boot word, which includes the starting address of the multivalued parameter array, the address of the affect word, and the decoding method. After obtaining the state of the affect word, FRTDS calculates the address offset according to a certain decoding method and adds this offset to the starting address to obtain the real address of the current value of the multivalued parameter. The specific process is shown in Figure 3. The affect word is divided into internal and external affect words. The state of the internal affect word is determined by the calculation result of the arithmetic component, and the state of the external influence word is determined by the input of the external device.
The graphical simulation software provided by the FRTDS includes the graphical modeling of simulation objects and instruction-level assignment of computing tasks. The graphical modeling of the simulation object automatically provides the specific simulation calculation process according to the characteristics of the electrical equipment and their connection relationships. The instruction-level assignment of computing tasks automatically generates a directed acyclic graph describing the dependencies of computing tasks according to the positions of variables in the computing expression. For each single clock cycle, the task with the least amount of time slack among the ready computing tasks is found, and the cost to arrange it in every arithmetic component is calculated. Then, the task is arranged in the arithmetic component with the least cost until all computing tasks are scheduled, or the arithmetic component cannot take on more computing tasks.

3. Method of Replacing Inline Computation with Parameters Prestorage

Real-time simulations require the calculation time to be strictly synchronized with real-time; therefore, the computing power per unit time of the simulation platform is determined, so an algorithm with a small calculation amount must be selected to ensure that the hardware has sufficient ability to complete the calculation task within the expected time. Some calculation tasks can be calculated offline and then stored in the simulation platform. As such, this method of replacing computation with storage reduces the amount of calculation. Large-scale simulation objects are usually divided into multiple subnetworks, and an upper-layer network can be created after the external equivalent circuit of each subnetwork is obtained. After solving the upper-layer network, which has fewer nodes, the input variables of each subnetwork port are obtained, and the internal subnetworks are solved. Some calculation tasks during the equivalence and solution process of each subnetwork are appropriate for being replaced with prestored parameters. A discrete linear multiport network is equivalent to k Thevenin equivalent ports and m Norton equivalent ports. The input variable of the Thevenin ports is the current vector iA = [i1ik]T, and the input variable of the Norton equivalent ports is the current vector uB = [i1im]T. The equivalent voltage source vector Ueq and equivalent current source vector Ieq in the external equivalent circuit of the subnetwork can be expressed as a linear combination of the historical power vector [Uh Ih] and the independent power vector [Us Is] inside the subnetwork. The solution formula is:
[ U eq I e q ] = [ K h R h G h H h ] [ U h I h ] + [ K s R s G s H s ] [ U s I s ]
According to the network connection relationship, after simultaneously solving the port input variables of each subnetwork, the x to be solved in each subnetwork can be further solved as follows:
x = [ A 1 A 2 ] [ i A u B ] + [ B 1 B 2 ] [ U s I s ] + [ C 1 C 2 ] [ U h I h ]
In formula (2) [A1 A2] is the coefficient matrix associated with the input variables, [B1 B2] is the coefficient matrix associated with the independent sources, and [C1 C2] is the coefficient matrix associated with the historical sources. Kh, Rh, Gh, Hh, Ks, Rs, Gs, and Hs in (1), and A1, A2, B1, B2, C1, and C2 in (2) are sets of known constants after the element states within the subnetwork are determined, which can be precalculated and stored in the hardware of the simulation platform. When calculating in each step, the prestored coefficients are removed and used in the vector multiplication and addition operations of Equations (1) and (2) to complete the equivalence and solution of the subnetwork. The state variable method can also be used in the subnetwork equivalence and solution process. The voltage at the port of the subnetwork is taken as the external input vector of the subnetwork, the current at the port of the subnetwork is taken as the output vector of the subnetwork, and the independent voltage source or current source inside the subnetwork is taken as the internal input vector of the subnetwork. Then, the state and output equations of the subnetwork can be written as:
{ x a = A a x a + B a w a + E a u i a = C a x a + D a w a + F a u
where i is the current vector at the port of the subnetwork, u is the voltage vector at the port of the subnetwork, s is the independent power vector inside the subnetwork, and x is the state variable vector of the subnetwork. Use the implicit Euler method to differentiate Equation (3); then, the following formula is derived:
{ x a ( t ) = A a * x a ( t Δ t ) + B a * w a ( t ) + E a * u ( t ) i a ( t ) = C a * x a ( t Δ t ) + D a * w a ( t ) + F a * u ( t )
Equation (4) shows that the state variable and output vectors of the subnetwork at time t can also be expressed as the linear combination of the state vector of the subnetwork at time t − Δt, the independent power vector in the subnetwork at time t, and the voltage vector at the port of the subnetwork at time t. Combining the first two terms of the right-hand side of the output equation in Equation (4), the Norton equivalent circuit expression of the subnetwork is obtained:
i a ( t ) = P a ( t ) + F a * u ( t )
where P a ( t ) = C a * x a ( t Δ t ) + D a * w a ( t ) is the same as the multiport equivalent method based on a discrete network. The coefficient matrices A*, B*, C*, D*, E*, and F* are precalculated and stored in the simulation hardware, and in each step, the vector multiplication and addition operations are performed according to Equations (4) and (5) to complete the equivalence of the subnetwork and update the state variables.

4. P-FBSM Subnetwork Model

The MMC topology based on P-FBSM is shown in Figure 4. Each bridge arm of the converter is composed of a bridge arm reactor Larm and N parallel full-bridge submodules (PFBSM) connected in series. The upper and lower bridge arms of each phase form a phase unit, and the MMC bridge arm is still a single-port structure by shorting the two terminals of sub-modules at the head and end of the bridge arm. The P-FBSM is obtained by symmetrically inverting the full-bridge submodule, which is a typical dual-port submodule. The binary resistor is modeled, the capacitor is then differentiated by the Euler method, and the accompanying PFBSM circuit shown in Figure 4 is obtained. Using a binary resistor to model the entire switch group of the IGBT and diode in parallel with the PFBSM and differentiating the capacitor by the implicit Euler method, the companion circuit shown in Figure 4 can be obtained. For a submodule, the three conduction modes for the switch groups on the left and right sides of the capacitor are: positive (G1, G2 on), negative (G3, G4 on), and parallel (G1, G4 or G2, G3 on).
When simulating the MMC converter composed of a P-FBSM, if the node voltage analysis is directly used to solve the problem, owing to the large number of neutron modules in the bridge arm, the dimensions of the node voltage equation are large, which makes it difficult to quickly solve the problem. Therefore, the bridge arm needs to be divided into multiple subnetworks, and the method of replacing computation with storage needs to be used to perform real-time simulation. The currently used external equivalent circuit of the PFBSM subnetwork still has a dual-port structure, as shown in Figure 5. Each subnetwork contains n submodules; uex1, uex2, uex1, and uex2 are the input voltages of the subnetwork ports; G12, G13, G14, G23, G24, and G34 are the six equivalent conductances in the equivalent model; and Jin1, Jin2, Jin3, and Jin4 are equivalent current sources. The equivalent current source vector J = [Jin1 Jin2 Jin3 Jin4]T can be obtained by the linear combination of the historical value of each capacitor voltage in the subnetwork according to coefficient matrix C:
J ( t ) = C 4 × n · u c ( t Δ t )
By connecting the equivalent circuits of each subnetwork and eliminating the internal nodes, the equivalent circuit of the entire bridge arm is obtained, then the node voltage equation of the whole simulation object can be formed and solved. After finding the input voltage at the interface of each subnetwork, the voltage of each capacitor in the subnetwork can be calculated and updated. Each capacitor voltage can be expressed as the linear combination of the voltage vector uex = [uex1 uex2 uex1 uex2]T and the historical value of each capacitor voltage in the subnetwork:
u c ( t ) = A n × n · u c ( t Δ t ) + B 4 × n · u ex ( t )
The coefficient matrices An×n, B4×n, C4×n and the six equivalent conductances of the subnetwork in Equations (6) and (7) can be prestored. The coefficients in these matrices are related to the conduction mode of the switch groups on the left and right sides of the capacitor in each submodule. However, an analysis of the operation mode shows that if the on-off of the IGBT is controlled by the control pulse during normal operation or if only the on-off of the diode is controlled during the blockade of the whole station after startup or failure, only five connection modes will exist between two adjacent submodules, as shown in Figure 6. The number of parameters that each subnet needs to prestore is:
M e m = 3 2 × 5 n 1 × ( n 2 + 8 n + 6 )
The number of prestored parameters exponentially increases as the number of submodules in the subnetwork increases. Practically applying these parameters in an FRTDS without optimizing them is difficult. Given the parallel state of the newly added submodules in PFBSM–MMC, the concept of a submodule segment was proposed [21], in which the parallel submodule combination is regarded as a submodule segment. The segmented structure of a bridge arm is shown in Figure 7a. The submodule segment in the bridge arm has four different working states, as shown in Figure 7b: positive voltage input, negative voltage input, and two bypass states. A single submodule is also regarded as a submodule segment, and the number of submodules in the segment dynamically changes with the control signal.
By considering the submodule segment as a subnetwork, it is equivalent to a single-port Thevenin circuit. Then, the entire bridge arm is composed of the Thevenin equivalent circuits of n submodule segments in series. The resistance Rceq of a Thevenin equivalent circuit is related to the working mode and length of the submodule segment. The voltage source Uceq in a Thevenin equivalent circuit can be expressed as the linear combination of the historical value of the capacitor voltage of each submodule in the segment:
U ceq ( t ) = C 1 × l · u c ( t Δ t )
The value of element ci in coefficient matrix Cl of the historical capacitance voltage is related to the working mode of the submodule segment, the length of the submodule segment, and the position of the ith submodule in the segment. The Thevenin equivalent circuit of the whole bridge arm can be obtained by summing the equivalent voltage source and the equivalent resistance of each submodule segment in the bridge arm. Then, the upper-layer network node equation of the whole system can be formed and solved. After calculating the bridge arm current iarm, it is regarded as the input variable of each subnetwork. The submodule capacitor voltage in each segment can be updated by the linear combination of the current submodule capacitor voltage in each segment and iarm:
u c ( t ) = A l × l · u c ( t Δ t ) + B l × 1 · i a r m ( t )
The value of element aij in coefficient matrix Al×l of the historical capacitance voltage is determined by four factors: the working mode of the submodule segment, the length of the submodule segment, the position of the ith submodule in the segment, and the position of the jth submodule in the segment. The value of element bi in coefficient vector Bl of the bridge arm current is determined by three factors: the working mode of the submodule segment, the length of the submodule segment, and the position of the ith submodule in the segment. When the length of the module group is l, Al×l, Bl, Cl, and Rceq occupy 4l2, 4l, 4l, and 4 storage spaces, respectively, where l ranges from 1 to the total number of submodules in bridge arm m. When the submodule segment is used as the subnetwork, the number of parameters that need to be prestored is:
M e m = l = 1 m ( 4 l 2 + 4 l + 4 l + 4 ) = 4 3 m 3 + 6 m 2 + 26 3 m
At this point, the number of prestored subnetwork parameters quadratically increases with the number of submodules in the subnetwork, which are substantially fewer compared to the number of subnetwork parameters calculated according to Equation (3). However, if the submodule segment is directly used as the subnetwork, the calculation process changes after the segmentation situation changes, so the control instructions of the FRTDS involve conditional branching. The submodule segmentation status is determined by the control signal from the controller, and the user cannot obtain the segmentation status in real-time, which prevents switching the control instructions in real-time. To ensure that the instructions in FRTDS remain unchanged, only a fixed number of adjacent submodules in the topology can be delineated as a subnetwork, as shown in Figure 8a. Depending on the switch state, a subnetwork may contain several submodule segments or all submodules in the subnetwork may be connected in parallel and belong to one submodule segment. The parallel combination of submodules in the subnetwork, in addition to the four submodule segment forms in Figure 7b, has five added forms, as shown in Figure 8b. The equivalent circuit for the parallel combination of three-terminal submodules appears at the head or end of the subnetwork. When all the submodules in the subnetwork are connected in parallel, the entire subnetwork is a dual-port structure. The equivalent circuit for the parallel combination of three-terminal submodules appearing at the head or end is shown in Figure 9a. Combining the single-ended Thevenin equivalent circuit of the submodule segment and the dual-port equivalent circuit of all submodules in parallel, the P-FBSM subnet equivalent model, shown in Figure 9b, is finally obtained. When solving the model, first, the values of each parameter in the figure are obtained by looking up the table or calculating according to the current operation of the subnetwork, then internal nodes 5 and 6 can be eliminated; finally, the model is simultaneously solved with other subnetworks. When the subnetwork does not contain a certain type of submodule parallel structure at the current time, the unified model of six nodes remains unchanged, and the part of the model that reflects the parallel combination of this type takes the maximum value of the conductance and the minimum value of the resistance.
The data storage capacity of the equivalent model in Figure 9b is the sum of the prestored parameters in the case of a single-port Thevenin circuit, two three-terminal equivalent circuits, and a dual-port equivalent circuit with all submodules connected in parallel. The storage capacity of a single-port Thevenin circuit can be calculated by Equation (11). For a subnetwork with n submodules, when all the submodules are connected in parallel, the equivalent current sources J1, J2, J5, and J6 are calculated by Equation (6), and the voltage of each capacitor in the subnetwork is updated by Equation (7). When the number of submodules in the subnetwork is determined, each element in the coefficient matrix An × n, B4 × n, C4 × n, and the six equivalent conductances G1(G12), G13, G14, G23, G24, and G6(G34) in the subnetwork are all definite values; the storage capacity is n2 + 8n + 6.
For a three-terminal structure, taking the head as an example, the update of the capacitor voltage of the submodule in the parallel combination can be calculated by:
u c ( t ) = A l × l · u c ( t Δ t ) + B l × 2 · [ u ex 1 ( t ) u ex 2 ( t ) ] T + C l × 1 i a r m
which needs to store the historical capacitor voltage value coefficient matrix A l × l , the port input voltage coefficient matrix B l × 2 , and the bridge-arm current coefficient matrix C l × 1 . The equivalent current source vector J = [J1 J2 J3] is obtained by the linear combination of the capacitor voltages of the submodule in parallel combination, for which the voltage coefficient matrix D 3 × l and the three conductances G1, G2, and G3 need to be stored. The value of l may range from 1 to n, and two parallel combinations of submodules corresponding to this model are possible. In summary, the total storage capacity of the three terminals at the head is 2 3 n 3 + 7 n 2 + 37 3 n , and the total storage capacity of the three terminals at the end is the same as that of the head end. The number of prestored parameters required for the equivalent model of the six-node P-FBSM subnetwork is:
M e m = ( 4 3 n 3 + 6 n 2 + 26 3 n ) + 2 × ( 2 3 n 3 + 7 n 2 + 37 3 n ) + ( n 2 + 8 n + 6 ) = 8 3 n 3 + 21 n 2 + 124 3 n + 6
When the total number of submodules in the bridge arm is N, suppose the subnetwork contains n submodules, and the entire bridge arm is divided into N/n subnetworks. The amount of calculation for elimination and back-substitution between two adjacent subnetworks is fixed as E1, which needs to be performed (N/n − 1) times. Inside the P-FBSM model, each conductance value can be obtained by looking up the table. However, the six equivalent injection current sources J1J6, the equivalent voltage source, and the resistance of the single-port Thevenin equivalent circuit all need to be calculated according to the prestored parameters. To ensure that the calculation formula remains unchanged under the model’s various working states, the single-port Thevenin equivalent circuit part of the model is obtained by adding the single-port Thevenin equivalent circuits to the n submodule segments. Only when all the submodules in the subnetwork are connected in series will n submodule segments be formed. In other cases, the number of submodules in each segment is determined by the actual segmentation situation; then, the resistance of the Thevenin equivalent circuit is obtained by looking up the table and the power supply of the Thevenin equivalent circuit is calculated according to the prestored parameters. The resistance of the Thevenin equivalent circuit in the nonexistent segment is set to the minimum value of 10−7 Ω. The calculation formula of this part of the model is:
R ceq = i = 1 n R ceq i U ceq = i = 1 n C i " u i ( t Δ t )
According to the two operational formulas of A × B + C × D and A + B + C + D, the calculation amount is determined. The calculation amount of Equation (13) is 0.75n. The calculation formulas for the six equivalent injection current sources are the same:
J h   = i = 1 n D i " x i ( t Δ t ) , h = 1 , 2 , 6
The amount of calculation for Formula (14) is 3n, and the updated formula of the internal capacitor voltage of the subnetwork is:
x i ( t )   = j = 1 n A j " × x j ( t Δ t ) + B ( u ex 1 ) × u ex 1 ( t ) + B ( u ex 2 ) × u ex 2 ( t ) + B ( u ex 3 ) × u ex 3 ( t )   + B LT ( u ex 4 ) × u ex 4 ( t ) + B ( i arm ) × i arm ( t ) , i = 1 , 2 , , n
The amount of calculation of Formula (15) is approximately 0.5n2 + 3n. In addition to the parameter calculation in the model, the subnetwork six-node model needs to eliminate the two internal nodes, 5 and 6, with a fixed amount of calculation E2. In summary, the total calculation amount for a bridge arm is:
C a l = ( 0.5 n 2 + 6.75 n + E 2 ) × N n + ( N n 1 ) × E 1 = N × ( 0.5 n + E 1 + E 2 n ) + 6.75 N E 1
Equation (17) shows that when the number of submodules n in each subnetwork is determined, the amount of calculation linearly increases with the total number of submodules in the bridge arm. When the total number N of submodules in the bridge arm is determined by taking the derivative of n in Formula (17), its value is smallest when n = 2 ( E 1 + E 2 ) . In actual simulations, a close integer number of submodules can be used to divide the subnetwork of the bridge arm. After division, the number of prestored parameters of the subnetwork needs to be calculated. If they cannot be stored, then the number of submodules in the subnetwork should be reduced to the amount that can be stored.

5. Adaptive Transformation of FRTDS Platform

The original FRTDS is mainly oriented to traditional AC systems, and the amount of calculation is mainly due to the formation and solution of the node voltage equation. The data dependencies between the calculation tasks are complex, and the data flow is heavily coupled. For example, when forming the injection current source for the node equation, the injection current source of a node needs to sum the historical current sources on all its related branches. When a node is eliminated, the self-admittance and mutual-admittance of multiple adjacent nodes are affected. This results in the diverse arithmetic operations and the data reading and writing processes being irregular when FRTDS performs simulation, necessitating read-and-write address instructions and operation type selection instructions at every moment to control the data flow and arithmetic operation selection. When using FRTDS for PFBSM–MMC simulation, the unified subnetwork model proposed in Section 2 is used, in which the external equivalent and internal solutions of the bridge arm adopt the matrix–vector multiplication operation. The execution process of this calculation is shown in Figure 10. The data flow direction in the inner and outer loops is determined, and only the fixed operation formula Y = Σ(A × B + C × D) is executed. By storing the data involved in the operation in the storage area in order, sequential reading and writing can be achieved.
The multiply accumulating operation used by matrix–vector multiplication accounts for a large proportion of the whole calculation process, and a dedicated execution unit was designed to perform this part of the calculation. The floating-point adder needs five clock cycles to obtain the calculation result. To ensure the data input involved in the accumulation is continuous, the floating-point number is converted into a fixed-point number, and the fixed-point number adder is used to complete the accumulation calculation in one clock cycle. The completion of each accumulation is determined by device S1 with a 1-bit control signal. The dedicated execution unit does not need a select controller to select the operation type of each clock. The instruction format used to drive the dedicated execution unit to run is addr(A, B, C, D, Y) + ctrl. The input and output data are arranged in an orderly manner according to the inner and outer loops, and the coefficient address continuously increases from the first address during the entire calculation process. In the inner loop, the data address is continuously increased from the first address, and the value is refetched from the first address after each inner loop. The output address is increased from the first address after each inner loop; the control command outputs 1 in the inner loop and 0 after each inner loop completes. Therefore, the instruction storage area dedicated for executing unit instruction can store only one start instruction, start(AS,BS,CS,DS,YS) +m,n; that is, the start address of each port and the number of inner loops n and the number of outer loops m. When the coefficients involved in matrix–vector multiplication are multivalued parameters, each instruction only needs to calculate the offset once. An instruction loop unit, shown in Figure 11, was designed to generate the instructions required in the calculation process of the dedicated execution unit. The +1 and −1 operations in the instruction loop unit consume one clock; thus, a delay of one clock needs to be set for the data selected in synchronization. Channels A and C are coefficient channels, and channels B and D are data channels. The instruction loop unit reduces the storage resource consumption of the instruction storage area of the dedicated execution unit, leaving FRTDS more storage resources for prestored parameters. When forming the microprocessor core, several universal arithmetic components are replaced with dedicated execution units, thereby reducing the resource consumption of each microprocessor core and enabling more microprocessor cores to be placed in the FPGA. The computing tasks performed by the general-purpose computing components and the dedicated execution units in the same microprocessor core still communicate in a shared memory manner. The group of multiply accumulating operations is not split and is handled by a dedicated execution unit, which does not involve inter-core communication. After the calculation of each step starts, according to all the switch states of the current bridge arm, combined with the current direction of the bridge arm, the parallel combination mode of every submodule segment in each subnetwork and the number of submodules in the segment are determined.

6. Example Verification

6.1. Simulation Platform

FRTDS uses the Virtex-7 FPGA VC709 as the development board for the solver, and the FPGA chip on the board is XC7VX690T-2FFG1761. The chip contains 693,120 logic units, 108,300 configurable logic modules, 3600 DSP slices, and 1470 36KB dual-port BRAMs. The original FRTDS contains four universal arithmetic components in each microprocessor core. The newly designed FRTDS replaces two universal arithmetic components in each computing core with dedicated execution units. Table 1 shows the resources consumed by placing different numbers of microprocessor cores in the FRTDS in VC709. Limited by the configurable logic block (CLB), the original scheme only holds five microprocessor cores at most. The new design scheme has six microprocessor cores. Because the dedicated execution unit consumes less instruction storage space, in the new design scheme, the storage resources used for instruction storage in the original design can be saved for storing the parameters by replacing online computation with parameter prestorage. Finally, we designed the FRTDS according to the new design scheme with six microprocessor cores, working at a clock frequency of 150 MHz, and each microprocessor core is able to store 8000 prestored parameters in the form of double-precision floating points.
The graphical simulation software and task scheduling software are developed by the QT C++ platform. The graphical modeling of the simulation object automatically provides the specific simulation calculation process according to the characteristics of the electrical equipment and their connection relationships. The instruction-level assignment of computing tasks automatically generates a directed acyclic graph (DAG) describing the dependencies of computing tasks according to the positions of variables in the computing expression. The task scheduling software schedules the DAG graph to get an instruction stream, which can be stored or downloaded to FRTDS by the upper computer. FRTDS performs the corresponding operations according to each instruction to complete the simulation calculation. We connected the designed FRTDS to the digital relay protection device through a common ethernet switch to form a hardware-in-loop real-time simulation platform, as shown in Figure 12. We adopted the UDP communication protocol between the FRTDS and the upper computer responsible for command download and system operation monitoring, and we used the IEC61850 protocol to connect the real digital protection device. In particular, the IP address of FRTDS and the upper computer should be set on the same LAN. The message analyzer communicated with the FRTDS using SV messages and is able to draw and display real-time simulation waveforms. For another development board, VC707, we simulated the control system of the flexible DC transmission system to control the Aurora protocol of Xilinx Company between the controller and the FRTDS instead of the real controller.

6.2. Simulation Objects and Results

As the simulation object, we adopted the double-terminal MMC–HVDC system shown in Figure 13, where the AC grid voltage was 230 kV, and the DC bus voltage was 400 kV. The topology of the converter submodule was P-FBSM, in which each bridge arm contained 80 submodules, the submodule capacitance was 2 mF, the bridge arm reactor was 50 mH, and the converter capacity was 350 MVA. For the inverter-side MMC1, we adopted the control mode of fixed active power and fixed reactive power; for the rectifier-side MMC2, we adopted the control mode of fixed DC voltage and fixed reactive power and adopted the previously described low-switching-frequency modulation strategy [22]. The simulation step was set to 10 μs. If using the previously proposed calculation method of eliminating the internal nodes of all dual-port submodules in the MMC bridge arm and the interconnection nodes between modules in cyclic iteration [11], the calculation of one step would cost 64.8 μs, which is insufficient for real-time simulation. In the newly designed FRTDS, the calculation amount E1 of elimination and back-substitution between subnetworks is 38 unit quantity, and the calculation amount E2 of eliminating nodes 5 and 6 in the six-node subnetwork model is 16 unit quantity. Therefore, we divided the bridge arm into subnetworks, with each subnetwork containing 10 submodules. According to Equation (13), the number of parameters that the subnetwork needed to prestore was 5186, and the FRTDS was capable of storing these prestored parameters. After the subnetwork was divided, we used the unified model of the 6-node PFBSM subnetwork and performed the real-time simulation in FRTDS with a step of 9.62 μs using the method of replacing inline computation with the parameters’ prestorage. We built the same simulation object using MATLAB/Simulink, and we compared the offline running results with the real-time calculation results produced by FRTDS:
(1)
Trouble-free operation
Figure 14a–d shows the A-phase AC current Isa on the inverter side, the DC-side voltage Udc, the capacitor current Ic of the No. 1 submodule on the upper bridge arm of the rectifier side in A phase, and the capacitor voltage of each submodule of the upper bridge arm in A phase when the MMC–HVDC system operates in a steady-state, respectively. Figure 14a–c shows that when the MMC simulation system ran without fault, the simulation waveform in FRTDS using the proposed modeling method had a high degree of coincidence with the simulation waveform in the detailed model in MATLAB. The error was mainly reflected in the high-frequency irregular burrs, and the maximum relative error was about 0.89% in the DC-side voltage. The result verified that the simulation accuracy of the proposed modeling method is high. From the waveform of the bridge arm capacitor voltage of the FRTDS shown in Figure 14d, we found that the capacitor voltage was balanced over time, and some capacitors had equal voltage at some moments, which demonstrated the parallel voltage equalization characteristics of P-FBSM.
(2)
DC bipolar fault
We set a DC bipolar short-circuit fault at t = 1.295 s, the whole station was blocked 5 ms after the fault occurred, which unblocked at t = 1.365 s. Figure 15a–d shows the AC current Isa of the inverter side in phase A under the fault, the capacitor voltage Uc of the first submodule of the upper bridge arm in phase A on the rectifier side, the DC-side voltage Udc, and the DC side current Idc. Figure 15 shows that the simulation accuracy of FRTDS was still high during the bipolar fault occurrence in the MMC simulation system before and after the blocking, and the maximum relative error was about 1.09% on the DC-side voltage.

7. Conclusions

(1) We found that the proposed unified equivalent model of the P-FBSM subnetwork, after analyzing the working mechanism of the submodules, requires less computation by replacing inline computation with parameter prestorage and reduces the storage capacity required for prestored parameters, making it suitable for application in FRTDS.
(2) After dividing the subnetwork by a fixed number of submodules, we found a relationship between the amount of simulation calculation and the number of submodules in the subnetwork. Accordingly, the subnetwork division method with the least amount of calculation could be determined.
(3) After improving the structure of the computing components, we found that the FRTDS real-time simulation platform has improved resource use, requires less storage space for prestored parameters, and is more suitable for real-time simulations of dual-port submodule MMC–HVDC systems by replacing inline computation with the prestorage of parameters.
When subnetwork multivalued parameters are prestored, for the convenience of hardware query design, the current method does not consider the problem where different prestored parameters may have the same value. In future research, the amounts of computation and storage should be further optimized, a matching multivalued parameters query circuit should be designed in FRTDS, and the improved method and platform should be applied to the real-time simulation of MMC–HVDC using other types of dual-port submodule topologies.

Author Contributions

Conceptualization, Z.J.; methodology, B.Z. and Z.J.; software, Y.W.; validation, S.W. and Z.J.; formal analysis, B.Z.; investigation, B.Z. and Z.J.; resources, B.Z.; data curation, Y.W. and S.W.; writing—original draft preparation, S.W.; writing—review and editing, B.Z. and S.W.; visualization, Y.W.; supervision, B.Z.; project administration, B.Z.; funding acquisition, B.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (No. 51477114).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We are grateful for the experimental site provided by the Key Laboratory of the Ministry of Education for Smart Grid of Tianjin University.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bresesti, P.; Kling, W.L.; Hendriks, R.L.; Vailati, R. HVDC Connection of Offshore Wind Farms to the Transmission System. IEEE Trans. Energy Convers. 2007, 22, 37–43. [Google Scholar] [CrossRef]
  2. Gomis-Bellmunt, O.; Sau-Bassols, J.; Prieto-Araujo, E.; Cheah-Mane, M. Flexible Converters for Meshed HVDC Grids: From Flexible AC Transmission Systems (FACTS) to Flexible DC Grids. IEEE Trans. Power Deliv. 2020, 35, 2–15. [Google Scholar] [CrossRef]
  3. Tang, L.; Dong, X.; Shi, S.; Qiu, Y. A high-speed protection scheme for the DC transmission line of a MMC–HVDC grid. Electr. Power Syst. Res. 2019, 168, 81–91. [Google Scholar] [CrossRef]
  4. Resch, S.; Friedrich, J.; Wagner, T.; Mehlmann, G.; Luther, M. Stability Analysis of Power Hardware-in-the-Loop Simulations for Grid Applications. Electronics 2022, 11, 7. [Google Scholar] [CrossRef]
  5. Jeon, J.H.; Kim, J.Y.; Kim, H.M.; Kim, S.K.; Cho, C.; Kim, J.M.; Ahn, J.B.; Nam, K.Y. Development of Hardware In-the-Loop Simulation System for Testing Operation and Control Functions of Microgrid. IEEE Trans. Power Electron. 2010, 25, 2919–2929. [Google Scholar] [CrossRef]
  6. Bruno, S.; Giannoccaro, G.; Iurlaro, C.; La Scala, M.; Rodio, C. Power Hardware-in-the-Loop Test of a Low-Cost Synthetic Inertia Controller for Battery Energy Storage System. Energies 2022, 15, 3016. [Google Scholar] [CrossRef]
  7. Ilves, K.; Taffner, F.; Norrga, S.; Antonopoulos, A.; Harnefors, L.; Nee, H.P. A submodule implementation for parallel connection of capacitors in modular multilevel converters. IEEE Trans. Power Electron. 2015, 30, 3518–3527. [Google Scholar] [CrossRef]
  8. Li, Z.; Lizana, R.; Sha, S.; Yu, Z.; Peterchev, A.V.; Goetz, S.M. Module implementation and modulation strategy for sensor-less balancing in modular multilevel converters. IEEE Trans. Power Electron. 2019, 34, 8405–8416. [Google Scholar] [CrossRef] [Green Version]
  9. Xu, J.; Zhang, J.; Li, J.; Shi, L.; Jia, X.; Zhao, C. Series-parallel HBSM and two-port FBSM based hybrid MMC with local capacitor voltage self-balancing capability. Int. J. Electr. Power Energy Syst. 2018, 103, 203–211. [Google Scholar] [CrossRef]
  10. Goetz, S.M.; Peterchev, A.V.; Weyh, T. Modular multilevel converter with series and parallel module connectivity: Topology and control. IEEE Trans. Power Electron. 2015, 30, 203–215. [Google Scholar] [CrossRef]
  11. Xu, J.; Fan, S.; Zhao, C.; Gole, A.M. High-Speed EMT Modeling of MMCs With Arbitrary Multiport Submodule Structures Using Generalized Norton Equivalents. IEEE Trans. Power Deliv. 2018, 33, 1299–1307. [Google Scholar] [CrossRef]
  12. Sharma, A.; Srivastava, S.C.; Chakrabarti, S. Testing and validation of power system dynamic state estimators using real time digital simulator (RTDS). IEEE Trans. Power Syst. 2016, 31, 2338–2347. [Google Scholar] [CrossRef]
  13. Yu, D.; Zhou, N.; Luo, Y.; Dong, L.; Jia, Z. Research on the Cross-Platform Co-Simulation Strategy of Power Systems Based on the Model-Segmentation Algorithm. Electronics 2021, 10, 3185. [Google Scholar] [CrossRef]
  14. Yuan, J.; Samson, S.Y.; Zhang, G.; Lim, C.P.; Trinh, H.; Zhang, Y. Design and HIL Realization of an Online Adaptive Dynamic Programming Approach for Real-Time Economic Operations of Household Energy Systems. IEEE Trans. Smart Grid 2022, 13, 330–341. [Google Scholar] [CrossRef]
  15. Zi, P.; Li, Y.; Tan, B.; Chen, X.; Jia, L.; Wang, M. Current situation of largescale power grid simulation tools and their popularization and application in North China Power Grid. Electr. Power Autom. Equip. 2019, 39, 199–205. [Google Scholar]
  16. Chen, Y.; Dinavahi, V. Digital hardware emulation of universal machine and universal line models for real-time electro-magnetic transient simulation. IEEE Trans. Ind. Electron. 2012, 59, 1300–1309. [Google Scholar] [CrossRef]
  17. Liu, J.; Dinavahi, V. A Real-Time Nonlinear Hysteretic Power Transformer Transient Model on FPGA. IEEE Trans. Ind. Electron. 2013, 61, 3587–3597. [Google Scholar] [CrossRef]
  18. Bai, H.; Luo, H.; Liu, C.; Paire, D.; Gao, F. A Device-Level Transient Modeling Approach for the FPGA-Based Real-Time Simulation of Power Converters. IEEE Trans. Power Electron. 2020, 35, 1282–1292. [Google Scholar] [CrossRef]
  19. Matar, M.; Iravani, R. Massively Parallel Implementation of AC Machine Models for FPGA-Based Real-Time Simulation of Electromagnetic Transients. IEEE Trans. Power Deliv. 2010, 26, 830–840. [Google Scholar] [CrossRef]
  20. Zhang, B.; Jin, X.; Tu, S.; Jin, Z.; Zhang, J. A New FPGA-Based Real-Time Digital Solver for Power System Simulation. Energies 2019, 12, 4666. [Google Scholar] [CrossRef] [Green Version]
  21. Shi, L.; Zhao, C.; Xu, J.; Xu, Y. Research on voltage self-balancing operating characteristics of paralleled full-bridge submodules based MMC. Proc. CSEE 2018, 38, 2419–2428+2551. [Google Scholar]
  22. Xu, J.; Li, J.; Zhang, J.; Shi, L.; Jia, X.; Zhao, C. Open-loop voltage balancing algorithm for two-port full-bridge MMC–HVDC system. Int. J. Electr. Power Energy Syst. 2019, 109, 259–268. [Google Scholar] [CrossRef]
Figure 1. Computational features of node analysis.
Figure 1. Computational features of node analysis.
Energies 15 04624 g001
Figure 2. Structure of universal arithmetic components.
Figure 2. Structure of universal arithmetic components.
Energies 15 04624 g002
Figure 3. Indirect addressing workflow.
Figure 3. Indirect addressing workflow.
Energies 15 04624 g003
Figure 4. Topology of MMC based on P-FBSM.
Figure 4. Topology of MMC based on P-FBSM.
Energies 15 04624 g004
Figure 5. P-FBSM dual-port equivalent circuit.
Figure 5. P-FBSM dual-port equivalent circuit.
Energies 15 04624 g005
Figure 6. Connection (mode) of adjacent submodules.
Figure 6. Connection (mode) of adjacent submodules.
Energies 15 04624 g006
Figure 7. Submodule segment in PFBSM–MMC bridge arm: (a) segmented structure and equivalent model of bridge arm; (b) four working states of the submodule segment.
Figure 7. Submodule segment in PFBSM–MMC bridge arm: (a) segmented structure and equivalent model of bridge arm; (b) four working states of the submodule segment.
Energies 15 04624 g007
Figure 8. Subnetwork divided by a fixed number of submodules: (a) connection of submodules in the subnetwork; (b) added submodule parallel combination forms.
Figure 8. Subnetwork divided by a fixed number of submodules: (a) connection of submodules in the subnetwork; (b) added submodule parallel combination forms.
Energies 15 04624 g008
Figure 9. Equivalent model considering the parallel combination of submodules: (a) equivalent circuit of three-terminal submodule parallel combination; (b) unified equivalent model of PFBSM subnetworks.
Figure 9. Equivalent model considering the parallel combination of submodules: (a) equivalent circuit of three-terminal submodule parallel combination; (b) unified equivalent model of PFBSM subnetworks.
Energies 15 04624 g009
Figure 10. Computational features of matrix–vector multiplication.
Figure 10. Computational features of matrix–vector multiplication.
Energies 15 04624 g010
Figure 11. Structure of dedicated execution units and instruction loop unit.
Figure 11. Structure of dedicated execution units and instruction loop unit.
Energies 15 04624 g011
Figure 12. Hardware-in-loop system based on FRTDS.
Figure 12. Hardware-in-loop system based on FRTDS.
Energies 15 04624 g012
Figure 13. Double-ended MMC–HVDC system.
Figure 13. Double-ended MMC–HVDC system.
Energies 15 04624 g013
Figure 14. Steady-state operation waveform of the MMC simulation system: (a) A-phase AC current waveform of inverter side; (b) DC-side voltage waveform; (c) waveform of the capacitor current of the upper bridge arm in phase A on the rectifier side; (d) capacitor voltage of the upper bridge arm in phase A on the rectifier side.
Figure 14. Steady-state operation waveform of the MMC simulation system: (a) A-phase AC current waveform of inverter side; (b) DC-side voltage waveform; (c) waveform of the capacitor current of the upper bridge arm in phase A on the rectifier side; (d) capacitor voltage of the upper bridge arm in phase A on the rectifier side.
Energies 15 04624 g014
Figure 15. Bipolar fault transient waveform diagram of the MMC simulation system: (a) A-phase AC current waveform of the inverter side; (b) waveform diagram of the capacitor voltage of the upper bridge arm of phase A on the rectifier side; (c) DC-side voltage waveform; (d) DC-side current waveform.
Figure 15. Bipolar fault transient waveform diagram of the MMC simulation system: (a) A-phase AC current waveform of the inverter side; (b) waveform diagram of the capacitor voltage of the upper bridge arm of phase A on the rectifier side; (c) DC-side voltage waveform; (d) DC-side current waveform.
Energies 15 04624 g015aEnergies 15 04624 g015b
Table 1. Resource consumption of the improved FRTDS design.
Table 1. Resource consumption of the improved FRTDS design.
Number of
Microprocessing Cores
Resource Consumption
CLBDSPRAM
New FRTDSOriginal FRTDSNew FRTDSOriginal FRTDS
342.63%51.25%18.7%25.86%50.03%
453.83%70.94%24.56%34.11%64.18%
567.14%92.62%30.42%42.36%78.35%
684.6%/36.28%/92.53%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jin, Z.; Wu, Y.; Wang, S.; Zhang, B. FPGA-Based Real-Time Simulation of Dual-Port Submodule MMC–HVDC System. Energies 2022, 15, 4624. https://doi.org/10.3390/en15134624

AMA Style

Jin Z, Wu Y, Wang S, Zhang B. FPGA-Based Real-Time Simulation of Dual-Port Submodule MMC–HVDC System. Energies. 2022; 15(13):4624. https://doi.org/10.3390/en15134624

Chicago/Turabian Style

Jin, Zhao, Yanjie Wu, Shuyuan Wang, and Bingda Zhang. 2022. "FPGA-Based Real-Time Simulation of Dual-Port Submodule MMC–HVDC System" Energies 15, no. 13: 4624. https://doi.org/10.3390/en15134624

APA Style

Jin, Z., Wu, Y., Wang, S., & Zhang, B. (2022). FPGA-Based Real-Time Simulation of Dual-Port Submodule MMC–HVDC System. Energies, 15(13), 4624. https://doi.org/10.3390/en15134624

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop