Multivalued Coefficient Prestorage and Block Parallel Method for Real-Time Simulation of Microgrid on FRTDS

The microgrid containing a large amount of high frequency power switches and nonlinear components has put forward high requirements for power system real-time simulation technology. Multivalued coefficient prestorage can reduce the calculation steps in real-time simulation. In order to reduce the storage pressure of the multivalued coefficients, the whole network is divided into multiple subnetworks that can be simulated in parallel, and only the parameters for computing input variables and internal variables are prestored. The multiport hybrid equivalent is performed to reduce the number of simultaneous network equations. The input variables are tied to state variables of the circuit so that the iterative calculation is limited to the local network. The devised methodology is validated through simulation of a low-voltage microgrid on a field programmable gate array (FPGA)-based real-time digital simulation (FRTDS) platform at a 5 μs time step. Comparison with a power systems computer aided design (PSCAD)/electromagnetic transients including DC (EMTDC) model shows that the proposed method is effective.


Introduction
The energy Internet aims to build a clean, low-carbon, safe and efficient modern energy system, which supports the access of high proportion, distributed clean energy power to improve the efficiency of energy utilization and the utilization of renewable energy [1][2][3].Microgrid is an important part of the new generation energy system as an effective way to play the role of distributed power supply.
A hardware-in-the-loop simulation (HILS) of a microgrid is a kind of semi-physical simulation, and its real-time simulator uses mathematical models to simulate microgrid electrical system, the actual manager or controller to manage or control its models.The HILS technique is used to test the real device and the dynamic performance of a microgrid in different operation modes and control modes, which plays an important role in the construction and development of a microgrid [4][5][6].
At present, the mainstream real-time simulators in the field of electrical engineering are real-time digital simulator (RTDS) and real time lab (RT-LAB).RTDS is a combination of specialized computer hardware and software designed specifically for the solution of power system electromagnetic transients, which simulates power system components to test power control equipment and protection devices in real time with a typically time-step size of about 50 µs [7].RT-LAB runs real-time simulation of a Simulink model on multi-CPU computer with shared memory, which integrates dynamic system mathematical models established by MATLAB/Simulink (R2015b, The MathWorks, Inc., Natick, MA, USA) and is mainly used for a power electronic controller verification test with time-steps

Electromagnetic Transient Simulation Model
In order to carry out the electromagnetic transient simulation of a microgrid, it is necessary to obtain the network model that contains only the power supply and the resistor.

Power Supply and Lumped Resistor
The linear resistor has a determined resistance value by the volt-ampere characteristic, where the voltage across the resistor is proportional to the current through it.In this paper, natural commutation is simulated using the state machine technique illustrated in Figure 1a.The value of the resistive model of the switch (R on /R off ) depends on the voltage across the diode u D .
Energies 2017, 10, 1248 3 of 18 is to improve the speed of real-time simulation by means of multivalued coefficient prestorage, network partitioning, multiport hybrid equivalent and local iterative calculation.The remainder of this paper is organized as follows.The electromagnetic transient simulation model is presented in Section 2. Section 3 describes the subnetwork multiport hybrid equivalent and its internal solution method.The process of writing and solving the network equation is discussed in Section 4. Section 5 presents the local iterative method for nonlinear networks.A low-voltage microgrid is simulated to verify the feasibility of the proposed method in Section 6. Section 7 concludes this paper.

Electromagnetic Transient Simulation Model
In order to carry out the electromagnetic transient simulation of a microgrid, it is necessary to obtain the network model that contains only the power supply and the resistor.

Power Supply and Lumped Resistor
The linear resistor has a determined resistance value by the volt-ampere characteristic, where the voltage across the resistor is proportional to the current through it.In this paper, natural commutation is simulated using the state machine technique illustrated in Figure 1a.The value of the resistive model of the switch (Ron/Roff) depends on the voltage across the diode uD.When the diode is connected in parallel to a forced commutation switch such as an insulated gate bipolar transistor (IGBT), the state machine in Figure 1b is used instead, where the gate signal is accounted for.
When the voltage across the voltage source is independent of the current flowing through it, and the output current of the current source has nothing to do with the terminal voltage.The voltage source or current source is considered to be the ideal power source and its voltage or current output is known at each simulation step.The nonlinear elements such as photovoltaic power supply can be expressed using piecewise linear approximation model in Figure 2a, i.e., + = i i U RI U (1) where i = 1, 2, 3 is the number of segments, Ri is the equivalent resistance of the i-th segment, Ui is the equivalent voltage source for piecewise linearization.Equation (1) can be expressed as an equivalent circuit shown in Figure 2b.When the diode is connected in parallel to a forced commutation switch such as an insulated gate bipolar transistor (IGBT), the state machine in Figure 1b is used instead, where the gate signal is accounted for.
When the voltage across the voltage source is independent of the current flowing through it, and the output current of the current source has nothing to do with the terminal voltage.The voltage source or current source is considered to be the ideal power source and its voltage or current output is known at each simulation step.The nonlinear elements such as photovoltaic power supply can be expressed using piecewise linear approximation model in Figure 2a, i.e., where i = 1, 2, 3 is the number of segments, R i is the equivalent resistance of the i-th segment, U i is the equivalent voltage source for piecewise linearization.Equation (1) can be expressed as an equivalent circuit shown in Figure 2b.
Energies 2017, 10, 1248 3 of 18 is to improve the speed of real-time simulation by means of multivalued coefficient prestorage, network partitioning, multiport hybrid equivalent and local iterative calculation.The remainder of this paper is organized as follows.The electromagnetic transient simulation model is presented in Section 2. Section 3 describes the subnetwork multiport hybrid equivalent and its internal solution method.The process of writing and solving the network equation is discussed in Section 4. Section 5 presents the local iterative method for nonlinear networks.A low-voltage microgrid is simulated to verify the feasibility of the proposed method in Section 6. Section 7 concludes this paper.

Electromagnetic Transient Simulation Model
In order to carry out the electromagnetic transient simulation of a microgrid, it is necessary to obtain the network model that contains only the power supply and the resistor.

Power Supply and Lumped Resistor
The linear resistor has a determined resistance value by the volt-ampere characteristic, where the voltage across the resistor is proportional to the current through it.In this paper, natural commutation is simulated using the state machine technique illustrated in Figure 1a.The value of the resistive model of the switch (Ron/Roff) depends on the voltage across the diode uD.When the diode is connected in parallel to a forced commutation switch such as an insulated gate bipolar transistor (IGBT), the state machine in Figure 1b is used instead, where the gate signal is accounted for.
When the voltage across the voltage source is independent of the current flowing through it, and the output current of the current source has nothing to do with the terminal voltage.The voltage source or current source is considered to be the ideal power source and its voltage or current output is known at each simulation step.The nonlinear elements such as photovoltaic power supply can be expressed using piecewise linear approximation model in Figure 2a, i.e., where i = 1, 2, 3 is the number of segments, Ri is the equivalent resistance of the i-th segment, Ui is the equivalent voltage source for piecewise linearization.Equation (1) can be expressed as an equivalent circuit shown in Figure 2b.

Lumped Elements (L, C)
The lumped elements (L, C) are represented by an equivalent conductance in parallel with a current source representing the historical values.The values for the equivalent conductance and the historical current source are determined by the element type together with the adopted numerical integration method.The excitation inductor of power transformer and the inductor of reactor are nonlinear because of saturation.The characteristic equation of the nonlinear inductor can be described as: u = dψ/dt, by using the implicit trapezoidal rule, the following equation is derived: where the latter two are the equivalent flux linkage.The process of determining the current by the flux linkage can be done using the piecewise linearization method.Figure 3a presents the piecewise linear representation of the nonlinear inductor.

Lumped Elements (L, C)
The lumped elements (L, C) are represented by an equivalent conductance in parallel with a current source representing the historical values.The values for the equivalent conductance and the historical current source are determined by the element type together with the adopted numerical integration method.The excitation inductor of power transformer and the inductor of reactor are nonlinear because of saturation.The characteristic equation of the nonlinear inductor can be described as: u = dΨ/dt, by using the implicit trapezoidal rule, the following equation is derived: where the latter two are the equivalent flux linkage.The process of determining the current by the flux linkage can be done using the piecewise linearization method.Figure 3a presents the piecewise linear representation of the nonlinear inductor.The piecewise linear representation of the nonlinear inductor is reduced to where k = 1, 2, 3 is the number of segments, Lk is the equivalent inductance of the k-th segment, and Ψk is the equivalent flux linkage generated by piecewise linearization.Substituting Label (3) in Label (2), the relation between voltage and current of the nonlinear inductor is derived: Equation ( 4) can be reduced to (5), which can be expressed as an equivalent circuit shown in Figure 3b: where k R is the equivalent resistance of the k-th segment, kLine u is the equivalent voltage source generated by piecewise linearization, and kdif u is the equivalent historical voltage source because of the discretization.The equivalent method for the nonlinear capacitor can be derived similarly.First, discretize the characteristic equation of the nonlinear capacitor i = dQ/dt and then use the piecewise linearization method to obtain the relation between voltage and current of the nonlinear capacitor.

Fault Model and Circuit Breaker Model
To simulate the various types of faults that occur in microgrid better, the fault model is shown in Figure 4a.When a single phase(a)-to-ground short-circuit fault occurs, let Ra be a small resistance and Rb, Rc be large resistances, and Rg, Lg be small values.When an interphase (b, c) short-circuit fault The piecewise linear representation of the nonlinear inductor is reduced to where k = 1, 2, 3 is the number of segments, L k is the equivalent inductance of the k-th segment, and ψ k is the equivalent flux linkage generated by piecewise linearization.Substituting Label (3) in Label (2), the relation between voltage and current of the nonlinear inductor is derived: Equation ( 4) can be reduced to (5), which can be expressed as an equivalent circuit shown in Figure 3b: where R k is the equivalent resistance of the k-th segment, u kLine is the equivalent voltage source generated by piecewise linearization, and u kdi f is the equivalent historical voltage source because of the discretization.The equivalent method for the nonlinear capacitor can be derived similarly.First, discretize the characteristic equation of the nonlinear capacitor i = dQ/dt and then use the piecewise linearization method to obtain the relation between voltage and current of the nonlinear capacitor.

Fault Model and Circuit Breaker Model
To simulate the various types of faults that occur in microgrid better, the fault model is shown in Figure 4a.When a single phase(a)-to-ground short-circuit fault occurs, let R a be a small resistance and R b , R c be large resistances, and R g , L g be small values.When an interphase (b, c) short-circuit fault occurs, let R a , R g be large resistances and R b , R c be small resistances.In order to obtain the equivalent circuit model for transient simulation, discretize the fault model by applying the implicit trapezoidal rule as shown in Figure 4b.
Energies 2017, 10, 1248 5 of 18 occurs, let Ra, Rg be large resistances and Rb, Rc be small resistances.In order to obtain the equivalent circuit model for transient simulation, discretize the fault model by applying the implicit trapezoidal rule as shown in Figure 4b.Circuit breaker model is illustrated in Figure 5. Ra, Rb, Rc are, respectively, the resistance of the three phases.Ra, Rb, Rc are equal to a small resistance in the normal case (circuit breakers closed) to indicate the circuit breakers are switched on.When the circuit breakers are switched off, the resistances of the circuit breakers take a large value.When the resistance of fault model and circuit breaker model changes from a small value to a large value, the zero-cross detection to the current should be applied.Because of the complexity of the arc extinguishing process, it is hereby specified that the model can only trip when the current passes through zero.

Subnetwork Multiport Hybrid Equivalent and Internal Solution Method
After electromagnetic transient simulation model, there are only two kinds of linear components, power supply and resistance in the network.Any one-port network with linear source can be reduced to an equivalent Norton/Thevenin circuit.
For the current port in the equivalent Thevenin circuit illustrated in Figure 6a, iA is the input variable, uA is the output variable, i.e., where Ueq is the open-circuit voltage of the original network port, Req is equal to the voltage value on the port when the input current is the unit value after all the independent sources are set to zero in the original network.
For the voltage port in the equivalent Norton circuit illustrated in Figure 6b, uB is the input variable, iB is the output variable, i.e., B B

= − +
eq eq i G u I (7) Circuit breaker model is illustrated in Figure 5. R a , R b , R c are, respectively, the resistance of the three phases.R a , R b , R c are equal to a small resistance in the normal case (circuit breakers closed) to indicate the circuit breakers are switched on.When the circuit breakers are switched off, the resistances of the circuit breakers take a large value.
Energies 2017, 10, 1248 5 of 18 occurs, let Ra, Rg be large resistances and Rb, Rc be small resistances.In order to obtain the equivalent circuit model for transient simulation, discretize the fault model by applying the implicit trapezoidal rule as shown in Figure 4b.Circuit breaker model is illustrated in Figure 5. Ra, Rb, Rc are, respectively, the resistance of the three phases.Ra, Rb, Rc are equal to a small resistance in the normal case (circuit breakers closed) to indicate the circuit breakers are switched on.When the circuit breakers are switched off, the resistances of the circuit breakers take a large value.When the resistance of fault model and circuit breaker model changes from a small value to a large value, the zero-cross detection to the current should be applied.Because of the complexity of the arc extinguishing process, it is hereby specified that the model can only trip when the current passes through zero.

Subnetwork Multiport Hybrid Equivalent and Internal Solution Method
After electromagnetic transient simulation model, there are only two kinds of linear components, power supply and resistance in the network.Any one-port network with linear source can be reduced to an equivalent Norton/Thevenin circuit.
For the current port in the equivalent Thevenin circuit illustrated in Figure 6a, iA is the input variable, uA is the output variable, i.e., where Ueq is the open-circuit voltage of the original network port, Req is equal to the voltage value on the port when the input current is the unit value after all the independent sources are set to zero in the original network.
For the voltage port in the equivalent Norton circuit illustrated in Figure 6b, uB is the input variable, iB is the output variable, i.e., B B

= − +
eq eq i G u I (7) When the resistance of fault model and circuit breaker model changes from a small value to a large value, the zero-cross detection to the current should be applied.Because of the complexity of the arc extinguishing process, it is hereby specified that the model can only trip when the current passes through zero.

Subnetwork Multiport Hybrid Equivalent and Internal Solution Method
After electromagnetic transient simulation model, there are only two kinds of linear components, power supply and resistance in the network.Any one-port network with linear source can be reduced to an equivalent Norton/Thevenin circuit.
For the current port in the equivalent Thevenin circuit illustrated in Figure 6a, i A is the input variable, u A is the output variable, i.e., u A = −R eq i A + U eq (6) where U eq is the open-circuit voltage of the original network port, R eq is equal to the voltage value on the port when the input current is the unit value after all the independent sources are set to zero in the original network.
Energies 2017, 10, 1248 6 of 18 For the voltage port in the equivalent Norton circuit illustrated in Figure 6b, u B is the input variable, i B is the output variable, i.e., i B = −G eq u B + I eq (7) where I eq is the short-circuit current of the original network port, and G eq is equal to the current value on the port when the input voltage is the unit value after all the independent sources are set to zero in the original network.The sign of the ports is also shown: black circles stand for current ports and white circles for voltage ports.
Energies 2017, 10, 1248 6 of 18 where Ieq is the short-circuit current of the original network port, and Geq is equal to the current value on the port when the input voltage is the unit value after all the independent sources are set to zero in the original network.The sign of the ports is also shown: black circles stand for current ports and white circles for voltage ports.The multiport hybrid equivalent is a general formulation of the Norton/Thevenin equivalents for the complex network.The hybrid nature of the results from the fact that the equivalent circuit encompasses voltage ports and current ports.For a circuit with k current ports (input variable iA = [i1……ik] T , output variable uA = [u1……uk] T ) and m voltage ports (input variable uB = [u1……um] T , output variable iB = [i1……im] T ), and the multiport hybrid equivalent of the circuit is given by the exterior characteristic equation: where [Ueq Ieq] T is the vector of the equivalent voltages/currents in the subnetwork, Ui(1 in Ieq are, respectively, the voltage and current values of the i-th port in the case all current ports are open and all voltage ports are shorted.
in HBA are, respectively, the voltage and current values of the i-th port in the case the voltage ports are all short and the current ports are open except for the j-th port with a unit current source imposed.
in GBB are, respectively, the voltage and current values of the i-th port in the case the current ports are all open and the voltage ports are shorted except for the j-th port with a unit voltage source imposed.The elements of matrices mentioned above need to be prestored as the multivalued coefficients, which are constants after the state of the circuit element is determined.Assume that the numbers of historical voltage source Uh, historical current source Ih, independent voltage source Us, independent current source Is in the subnetwork are a, b, c, d.The [Ueq Ieq] T in Equation ( 8) are associated with the historical source vector [Uh Ih] T and independent source vector [Us Is] T , i.e., where in the condition that Us and Is are set to zero, Kij(1 in Gh are, respectively, the voltage and current values of the i-th equivalent source in the case Ih are set to zero and Uh are set to zero except for the j-th historical voltage source with a unit voltage source imposed.
in Hh are, respectively, the voltage and current values of the i-th equivalent source in the case Uh are set to zero and Ih are set to zero except for the j-th historical current source with a unit current source imposed.In the condition that Uh and Ih are set to zero, the elements of Ks, Rs, Gs, Hs are derived similarly.The elements of matrices mentioned above need to be prestored as the multivalued coefficients, which are constants after the state of the circuit element is determined.The multiport hybrid equivalent is a general formulation of the Norton/Thevenin equivalents for the complex network.The hybrid nature of the results from the fact that the equivalent circuit encompasses voltage ports and current ports.For a circuit with k current ports (input variable , and the multiport hybrid equivalent of the circuit is given by the exterior characteristic equation: where [U eq I eq ] T is the vector of the equivalent voltages/currents in the subnetwork, U i (1 ≤ i ≤ k) in U eq and I i (1 ≤ i ≤ m) in I eq are, respectively, the voltage and current values of the i-th port in the case all current ports are open and all voltage ports are shorted.
, respectively, the voltage and current values of the i-th port in the case the voltage ports are all short and the current ports are open except for the j-th port with a unit current source imposed.
respectively, the voltage and current values of the i-th port in the case the current ports are all open and the voltage ports are shorted except for the j-th port with a unit voltage source imposed.The elements of matrices mentioned above need to be prestored as the multivalued coefficients, which are constants after the state of the circuit element is determined.Assume that the numbers of historical voltage source U h , historical current source I h , independent voltage source U s , independent current source I s in the subnetwork are a, b, c, d.The [U eq I eq ] T in Equation ( 8) are associated with the historical source vector [U h I h ] T and independent source vector [U s I s ] T , i.e., where in the condition that U s and I s are set to zero, in G h are, respectively, the voltage and current values of the i-th equivalent source in the case I h are set to zero and U h are set to zero except for the j-th historical voltage source with a unit voltage source imposed.
in H h are, respectively, the voltage and current values of the i-th equivalent source in the case U h are set to zero and I h are set to zero except for the j-th historical current source with a unit current source imposed.In the condition that U h and I h are set to zero, the elements of K s , R s , G s , H s are derived similarly.The elements of matrices mentioned above need to be prestored as the multivalued coefficients, which are constants after the state of the circuit element is determined.The unknown internal variables w within each subnetwork can be solved after the input variables of each subnetwork are determined, i.e., where

Write the Network Equation
Two two-port hybrid equivalents connected to solve interface input currents/voltages shown in Figure 7 are used as an example to introduce the process of writing the network equation.There are two distinct subcircuits identified by a circled number, and subcircuit 1 is reduced to whereas subcircuit 2 is reduced to Consider the Kirchhoff Voltage Law (KVL) and Kirchhoff Current Law (KCL) constraints of the input and output variables in Figure 7, and i 1 A1 , u 1 B2 are selected as the independent input variables.
Energies 2017, 10, 1248 7 of 18 The unknown internal variables w within each subnetwork can be solved after the input variables of each subnetwork are determined, i.e., where [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.The elements of [A1 A2], [B1 B2], and [C1 C2] need to be prestored as the multivalued coefficients, which are constants after the state of the circuit element is determined.

Write the Network Equation
Two two-port hybrid equivalents connected to solve interface input currents/voltages shown in Figure 7 are used as an example to introduce the process of writing the network equation.There are two distinct subcircuits identified by a circled number, and subcircuit 1 is reduced to whereas subcircuit 2 is reduced to Consider the Kirchhoff Voltage Law (KVL) and Kirchhoff Current Law (KCL) constraints of the input and output variables in Figure 7, and 1 A1 i , 1 B2 u are selected as the independent input variables.For the independent input current 1

A1
i , the loop equation can be reduced to For the independent input voltage 1 B2 u , the node equation can be reduced to Equations ( 13) and ( 14) constitute the network equations of Figure 7, and the input variables of the whole network can be obtained by solving the network equations simultaneously.
As shown in Figure 8a, there is a KCL constraint:  For the independent input current i 1 A1 , the loop equation can be reduced to For the independent input voltage u 1 B2 , the node equation can be reduced to Energies 2017, 10, 1248 8 of 18 Equations ( 13) and ( 14) constitute the network equations of Figure 7, and the input variables of the whole network can be obtained by solving the network equations simultaneously.
As shown in Figure 8a, there is a KCL constraint: A , and the ports are connected in series with current as input variables.There is KVL constraint: , and the ports are connected in parallel with voltage as input variables as shown in Figure 8b.All the port input variables can be obtained by solving only one independent current or voltage variable.
As justified and mentioned above, equivalent Thevenin circuit with the input current variables should be selected when the ports are connected in series, and an equivalent Norton circuit with the input voltage variables should be chosen when the ports are connected in parallel.
Energies 2017, 10, 1248 8 of 18 As justified and mentioned above, equivalent Thevenin circuit with the input current variables should be selected when the ports are connected in series, and an equivalent Norton circuit with the input voltage variables should be chosen when the ports are connected in parallel.

Solve the Network Equation
The input variables can be solved by the elimination method after forming the network equation.The basic idea of the elimination method is to convert the coefficient matrix of the network equation into a form that an equation only contains one input variable through the equivalent transformation, including two steps: the matrix is converted to echelon form and a back substitution process sequentially solves for each unknown.
In order to arrange the elimination order of the independent input variables and reduce the computation time, each desired input variable is regarded as a node, and the association of the nodes is represented by the undirected graph G = (V, E).Equation ( 13) is used as an example to illustrate the forming process of the undirected graph, and the independent input variable 1 A1 i is relevant to 1 B2 u .The undirected graph of the network equation can be obtained by the relation between each independent input variable searched in turn according to the network equations.The number of connections between the node and associated nodes is called the degree of the node.In order to reduce the computational burden, the node with the least degree should be eliminated preferentially.The degree of the remaining nodes needs to be recalculated after each process of elimination since the degree of the remaining nodes may change.The so-called minimum degree method means that the node with the least degree in the remaining nodes is given the priority of being eliminated.
The degree of the nodes adjacent to the eliminated node may change after each process of elimination, which has nothing to do with the rest of the node.To improve the parallelism of the elimination, it is necessary to find all nodes that can be simultaneously eliminated, i.e., the maximum node independent set.
In order to realize the two points mentioned above, the minimum degree maximum independent set method is adopted to arrange the elimination order of the network nodes.The node with the least degree is selected as the starting node to be eliminated, the nodes adjacent to which are marked as visited, and the node with the least degree is selected as the next node to be eliminated that has not been visited until the nodes left are already visited.These nodes form a maximum node independent set.The same method is applied to find the next maximum independent set until all nodes are arranged.

Solve the Network Equation
The input variables can be solved by the elimination method after forming the network equation.The basic idea of the elimination method is to convert the coefficient matrix of the network equation into a form that an equation only contains one input variable through the equivalent transformation, including two steps: the matrix is converted to echelon form and a back substitution process sequentially solves for each unknown.
In order to arrange the elimination order of the independent input variables and reduce the computation time, each desired input variable is regarded as a node, and the association of the nodes is represented by the undirected graph G = (V, E).Equation ( 13) is used as an example to illustrate the forming process of the undirected graph, and the independent input variable i 1 A1 is relevant to u 1 B2 .The undirected graph of the network equation can be obtained by the relation between each independent input variable searched in turn according to the network equations.The number of connections between the node and associated nodes is called the degree of the node.In order to reduce the computational burden, the node with the least degree should be eliminated preferentially.The degree of the remaining nodes needs to be recalculated after each process of elimination since the degree of the remaining nodes may change.The so-called minimum degree method means that the node with the least degree in the remaining nodes is given the priority of being eliminated.
The degree of the nodes adjacent to the eliminated node may change after each process of elimination, which has nothing to do with the rest of the node.To improve the parallelism of the elimination, it is necessary to find all nodes that can be simultaneously eliminated, i.e., the maximum node independent set.
In order to realize the two points mentioned above, the minimum degree maximum independent set method is adopted to arrange the elimination order of the network nodes.The node with the least Energies 2017, 10, 1248 9 of 18 degree is selected as the starting node to be eliminated, the nodes adjacent to which are marked as visited, and the node with the least degree is selected as the next node to be eliminated that has not been visited until the nodes left are already visited.These nodes form a maximum node independent set.The same method is applied to find the next maximum independent set until all nodes are arranged.

Local Iterative Method for Nonlinear Networks
The multivalued coefficients in the coefficient matrixes of the exterior characteristic Equation ( 8), the equivalent source Equation ( 9) and the unknown internal variables Equation (10) need to be stored in advance.The analysis above suggests that there are three steps of the block parallel real-time simulation algorithm.(1) Determine the parameters of the coefficient matrix in Labels ( 8), ( 9) and (10) according to the state of the components in each subnetwork, and calculate the equivalent voltage source vector U eq and the equivalent current source vector I eq ; (2) Compute the parameters of the coefficient matrix and the equivalent power vector in the network equations, and calculate the input variables of each subnetwork; (3) Solve for the unknown internal variables w within each subnetwork.
The voltage across each nonlinear element in the subnetwork is related to independent source vector [U s I s ] T , historical source vector [U h I h ] T , as well as input variables i A and u B .On the other hand, the nonlinear element status affects the selected multivalued coefficients, which, in turn, affects the input variables.When updating the status of nonlinear elements, it is necessary to modify the network equation repeatedly and compute the input variables in each subnetwork to determine the state of the nonlinear element.This constraint can be time consuming.In certain situations, this requirement can be relaxed.This is the case when input variables are tied to state variables (inductor current or capacitor voltage), which can not change all at once.It can be considered that the change of the nonlinear element status has very little influence on the external characteristic equation, so that the whole network equation needn't be solved again.Under this assumption, the local iteration method shown in Figure 9 can be used to describe the simulation algorithm.
Energies 2017, 10, 1248 9 of 18 The multivalued coefficients in the coefficient matrixes of the exterior characteristic Equation ( 8), the equivalent source Equation ( 9) and the unknown internal variables Equation (10) need to be stored in advance.The analysis above suggests that there are three steps of the block parallel realtime simulation algorithm.(1) Determine the parameters of the coefficient matrix in Labels ( 8), ( 9) and (10) according to the state of the components in each subnetwork, and calculate the equivalent voltage source vector Ueq and the equivalent current source vector Ieq; (2) Compute the parameters of the coefficient matrix and the equivalent power vector in the network equations, and calculate the input variables of each subnetwork; (3) Solve for the unknown internal variables w within each subnetwork.
The voltage across each nonlinear element in the subnetwork is related to independent source vector [Us Is] T , historical source vector [Uh Ih] T , as well as input variables iA and uB.On the other hand, the nonlinear element status affects the selected multivalued coefficients, which, in turn, affects the input variables.When updating the status of nonlinear elements, it is necessary to modify the network equation repeatedly and compute the input variables in each subnetwork to determine the state of the nonlinear element.This constraint can be time consuming.In certain situations, this requirement can be relaxed.This is the case when input variables are tied to state variables (inductor current or capacitor voltage), which can not change all at once.It can be considered that the change of the nonlinear element status has very little influence on the external characteristic equation, so that the whole network equation needn't be solved again.Under this assumption, the local iteration method shown in Figure 9 can be used to describe the simulation algorithm.
Note that when dealing with the nonlinear inductor and capacitor, it is necessary to apply an integration formula first and then use piecewise linear approximation model to determine the relationship between voltage and current according to the phenomenon that the flux linkage and the quantity of electric charge do not change all at once.It is worth mentioning that, when dividing the network into equivalent block, the following principles should be observed: (1) assemble the switching elements and nonlinear elements into different subnetworks as far as possible and ensure the input variables of the subnetwork containing Note that when dealing with the nonlinear inductor and capacitor, it is necessary to apply an integration formula first and then use piecewise linear approximation model to determine the relationship between voltage and current according to the phenomenon that the flux linkage and the quantity of electric charge do not change all at once.It is worth mentioning that, when dividing the network into equivalent block, the following principles should be observed: (1) assemble the switching elements and nonlinear elements into different subnetworks as far as possible and ensure the input variables of the subnetwork containing the switching element are tied to state variables, which limits the iterative calculation of determining the switching status to the local network; (2) reduce the storage pressure of multivalued coefficients by network partitioning as much as possible.

Simulation Platform
In the FRTDS real-time simulation platform, as shown in Figure 10, a Xilinx Virtex-7 VC709 FPGA development board is adopted for calculation and simulation as a simulator with the XC7 VX690 T-2 FFG1761 as the FPGA chip.It has the following features: 693,120 logic cells, 108,300 configurable logic elements, 3600 digital signal processing (DSP) blocks, 36 KB of 1470 dual-port BRAMs.
Energies 2017, 10, 1248 10 of 18 the switching element are tied to state variables, which limits the iterative calculation of determining the switching status to the local network; (2) reduce the storage pressure of multivalued coefficients by network partitioning as much as possible.

Simulation Platform
In the FRTDS real-time simulation platform, as shown in Figure 10, a Xilinx Virtex-7 VC709 FPGA development board is adopted for calculation and simulation as a simulator with the XC7 VX690 T-2 FFG1761 as the FPGA chip.It has the following features: 693,120 logic cells, 108,300 configurable logic elements, 3600 digital signal processing (DSP) blocks, 36 KB of 1470 dual-port BRAMs.The data interchange station is composed of four sets of dual-port memory, the number of which is the same as the microprocessor core.The data interchange station and the microprocessor core constitutes a hand-in-hand data transmission channel where data interchange stations are embedded in microprocessor cores.The entrance of the shared pipeline is connected to the microprocessor core via first-in-first-out (FIFO), and the exit end is connected to the data exchange station.
The step size for the electromagnetic transient simulation of the power system including commutators is generally in the microsecond range, which is usually achieved by reducing the calculation steps, and the multivalued coefficient prestorage is an effective method.When the impact values, these multivalued coefficients are treated as a whole and stored in DDR3, and the current values are quickly sent to BRAM when they are used.The value guide of the multivalued coefficients includes impact factor (15 bits), the offset word (5 bits) and the base address (12 bits).The impact factor consists of external contributors (settings, inputs) and internal contributors (the calculated information), wherein the data memory of the external contributors is a ping-pong mechanism.The offset word shows the method of calculating the offset by the impact factor.The hardware circuit is reconfigurable and the user can customize the offset word.The solver has 32 K × 64 B total BRAM bits and 128 M × 64 B total DDR3 bits for storing multivalued coefficients and provides a detection circuit to guarantee that external contributors are changed under specified conditions.The operation formula that the solver can handle is regarded as a task, and the parallelism and dependence between these tasks are described by the undirected graph.The influence function is defined to express the relation between the multivalued coefficients and the internal impact factors, thus reflecting the role of piecewise linearization process in electromagnetic transient simulation of the power system.The tasks that are ready for the independent variables are assigned to the microprocessor core according to the data access convenience criteria.The microprocessor core prioritizes the task with the earliest of the latest scheduled time (later than this time will extend the program execution time) that can be performed.The scheduled tasks that are not arranged are queued from ahead to behind in the order of the latest scheduled time, and then assigned to the most convenient microprocessor core that has the ability to execute data.The task can be executed with three conditions: the instruction length does not exceed the specified value, the arithmetic component has resources, and the read and write memory does not conflict.When data memory addresses of the two operation component ports are the same, the last data memory address is ignored by the autolatch function.When certain port data storage addresses satisfy the auto-increment feature, the data storage addresses are replaced by the status of corresponding function select line.When the read and write memory data are in conflict, the memory block stored the arguments can be changed or the arguments data can be sent to the arithmetic component port in advance, or the independent variables can be copied to other blocks in advance (included the shared area and give priority to the use of internal resources in the microprocessor core).The average allocation method and limited use strategy are used to arrange the storage of temporary data, and the overdue temporary data is used to avoid unnecessary waste of arithmetic components and instruction resources.In order to shorten The step size for the electromagnetic transient simulation of the power system including commutators is generally in the microsecond range, which is usually achieved by reducing the calculation steps, and the multivalued coefficient prestorage is an effective method.When the impact factors of the multivalued coefficient are the analog, several typical multivalued coefficients can be prestored or calculated in accordance with the set value of the analog.The multivalued coefficients can be stored in an industrial control machine or BRAM/DDR3 of the solver.When the impact factors of the multivalued coefficient are only the hypothetical information (such as line fault location and fault types), multivalued coefficients possible values are stored in the industrial control machine and only the current values are sent to BRAM of the solver when used.When the impact factors are the calculated information (such as the voltage on freewheeling diode) and the input from the real device (e.g., the IGBT switching pulse), the relevant multivalued coefficient should be stored in the BRAM.When the impact factors of some multivalued coefficients are the same and have many possible values, these multivalued coefficients are treated as a whole and stored in DDR3, and the current values are quickly sent to BRAM when they are used.The value guide of the multivalued coefficients includes impact factor (15 bits), the offset word (5 bits) and the base address (12 bits).The impact factor consists of external contributors (settings, inputs) and internal contributors (the calculated information), wherein the data memory of the external contributors is a ping-pong mechanism.The offset word shows the method of calculating the offset by the impact factor.The hardware circuit is reconfigurable and the user can customize the offset word.The solver has 32 K × 64 B total BRAM bits and 128 M × 64 B total DDR3 bits for storing multivalued coefficients and provides a detection circuit to guarantee that external contributors are changed under specified conditions.
The operation formula that the solver can handle is regarded as a task, and the parallelism and dependence between these tasks are described by the undirected graph.The influence function is defined to express the relation between the multivalued coefficients and the internal impact factors, thus reflecting the role of piecewise linearization process in electromagnetic transient simulation of the power system.The tasks that are ready for the independent variables are assigned to the microprocessor core according to the data access convenience criteria.The microprocessor core prioritizes the task with the earliest of the latest scheduled time (later than this time will extend the program execution time) that can be performed.The scheduled tasks that are not arranged are queued from ahead to behind in the order of the latest scheduled time, and then assigned to the most convenient microprocessor core that has the ability to execute data.The task can be executed with three conditions: the instruction length does not exceed the specified value, the arithmetic component has resources, and the read and write memory does not conflict.When data memory addresses of the two operation component ports are the same, the last data memory address is ignored by the auto-latch function.When certain port data storage addresses satisfy the auto-increment feature, the data storage addresses are replaced by the status of corresponding function select line.When the read and write memory data are in conflict, the memory block stored the arguments can be changed or the arguments data can be sent to the arithmetic component port in advance, or the independent variables can be copied to other blocks in advance (included the shared area and give priority to the use of internal resources in the microprocessor core).The average allocation method and limited use strategy are used to arrange the storage of temporary data, and the overdue temporary data is used to avoid unnecessary waste of arithmetic components and instruction resources.In order to shorten the simulation step as much as possible, a genetic algorithm is used to optimize the storage scheme of multivalued coefficients and cyclic data, and the execution time of the program is reduced.

Simulation and Analysis
The considered case study refers to a benchmark low voltage network developed by the European Union (EU) project "Microgrids".The system can be configured with a variety of line and load models, as well as various forms of distributed power generation systems, which fully embodies the microgrid structure and operation complexity.As shown in Figure 12, a photovoltaic system, a photovoltaic-battery system, a wind turbine system, a micro-gas turbine system, and a fuel cell system are arranged in the example microgrid.The photovoltaic array is simulated by a single-diode equivalent circuit model consisting of a photogenerated current source and a nonlinear diode in parallel to consider the internal loss of the battery [23].Universal model equivalent circuit is selected for the battery, composed of internal resistance and controlled voltage source connected in series [24].The wind turbine is modeled as double-fed wind power generation system, with the stator winding directly connected to the network and the rotor winding connected to the converter [25].The micro-gas turbine adopts single shaft structure with high-frequency alternating current (AC) at the motor outlet [26].The fuel cell uses the mid-term dynamic model, with a fixed impedance to indicate the internal overvoltage of the fuel cell [27].The detailed parameters of other components of the low-voltage microgrid are shown in [28].
The microgrid HIL test system is built by a multi-functional converter controller developed by Key Laboratory of Intelligent Power Grid of Ministry of Education, Tianjin University and relay protection device process control system (PCS)-931 GM-D.The multi-function converter controller uses Altera's FPGA EP3 C25, which integrates the maximum power point tracking control, the proportional-integral (PI) control for the chopper circuit, the constant power control of the inverter, inverter constant DC voltage and constant reactive power control, and other programs.It is capable of handling control tasks of distributed power generation systems in a low-voltage microgrid.In order to test the performance of the converter controller and the protection device, various short-circuit faults and ground faults that may occur at the nodes 14, 16, 18, 20, 22, 24, 26, and 28 are described by the fault model.
The low-voltage microgrid is divided at the nodes 3, 8, 30, 31, 32, 33, 35, 36, 38, 39, 41, 42, 44, 45, and 46, the DC-AC inverters and AC-AC converters are, respectively, divided into three and six subnetworks, so that the whole low-voltage microgrid system is divided into 36 subnetworks.After applying multiport hybrid equivalent method to the 36 subnetworks, the input variables can be obtained by solving a 38-dimensional network equations that are simplified by the constraints of the input and output variables, thus the low-voltage microgrid was decomposed into 36 subnetworks of independent parallel computing.
The photovoltaic-battery hybrid generation system is used as an example to introduce the proposed block parallel method as shown in Figure 13.The network is divided into nine subnetworks according to the principles of network partitioning as illustrated in Figure 14.
The wind turbine is modeled as double-fed wind power generation system, with the stator winding directly connected to the network and the rotor winding connected to the converter [25].The microgas turbine adopts single shaft structure with high-frequency alternating current (AC) at the motor outlet [26].The fuel cell uses the mid-term dynamic model, with a fixed impedance to indicate the internal overvoltage of the fuel cell [27].The detailed parameters of other components of the lowvoltage microgrid are shown in [28].The microgrid HIL test system is built by a multi-functional converter controller developed by Key Laboratory of Intelligent Power Grid of Ministry of Education, Tianjin University and relay protection device process control system (PCS)-931 GM-D.The multi-function converter controller uses Altera's FPGA EP3 C25, which integrates the maximum power point tracking control, the proportional-integral (PI) control for the chopper circuit, the constant power control of the inverter, inverter constant DC voltage and constant reactive power control, and other programs.It is capable of handling control tasks of distributed power generation systems in a low-voltage microgrid.In order to test the performance of the converter controller and the protection device, various short-   and 46, the DC-AC inverters and AC-AC converters are, respectively, divided into three and six subnetworks, so that the whole low-voltage microgrid system is divided into 36 subnetworks.After applying multiport hybrid equivalent method to the 36 subnetworks, the input variables can be obtained by solving a 38-dimensional network equations that are simplified by the constraints of the input and output variables, thus the low-voltage microgrid was decomposed into 36 subnetworks of independent parallel computing.
The photovoltaic-battery hybrid generation system is used as an example to introduce the proposed block parallel method as shown in Figure 13.The network is divided into nine subnetworks according to the principles of network partitioning as illustrated in Figure 14. Figure 15 presents the port connection relationship after applying multiport hybrid equivalent method to the nine subnetworks.Consider the constraints of the input and output variables: 3  u and 5 i , 6  i , 8  i , 8  i are selected as the independent input variables.Six node   and 46, the DC-AC inverters and AC-AC converters are, respectively, divided into three and six subnetworks, so that the whole low-voltage microgrid system is divided into 36 subnetworks.After applying multiport hybrid equivalent method to the 36 subnetworks, the input variables can be obtained by solving a 38-dimensional network equations that are simplified by the constraints of the input and output variables, thus the low-voltage microgrid was decomposed into 36 subnetworks of independent parallel computing.
The photovoltaic-battery hybrid generation system is used as an example to introduce the proposed block parallel method as shown in Figure 13.The network is divided into nine subnetworks according to the principles of network partitioning as illustrated in Figure 14. Figure 15 presents the port connection relationship after applying multiport hybrid equivalent method to the nine subnetworks.Consider the constraints of the input and output variables: Figure 15 presents the port connection relationship after applying multiport hybrid equivalent method to the nine subnetworks.Consider the constraints of the input and output variables: B3 , u 3 B4 and i 5 A3 , i 6 A3 , i 8 A3 , i 8 A4 are selected as the independent input variables.Six node equations and four loop equations are written for the ten independent input variables, which can be obtained by solving the 10-dimensional network equations.
Energies 2017, 10, 1248 14 of 18 Taking subnetwork 3 as an example, the number of ports is 4, the number of L and C is 3, the number of independent power supplies is 0, the number of unknown internal variables w is 1, the number of nonlinear components is 0, and the number of switches is 3.The numbers of the multivalued coefficients to be stored are, respectively, 16, 12 and 7 according to the exterior characteristic Equation ( 8), the equivalent source Equation ( 9), unknown internal variables Equation (10), and the total number x is 35.Each multivalued coefficient variable has eight possible values, and it can be seen that the total amount y of multivalued coefficients that need to be prestored is 280.Table 1 shows the numbers of the relevant data with the multivalued coefficients of the nine subnetworks.In the real-time simulation calculation of the whole low-voltage microgrid system, the numbers of ports, inductors and capacitors, independent power supplies, unknown internal variables w, nonlinear components and switches need to be known, thus the number x of the multivalued coefficients is 1196 and the total amount y of the multivalued coefficients is 17,135.At a step size of 5 μs, each microprocessor core can execute up to 1000 instructions.The equations of each subnetwork include the network equation to solve the input variables and the unknown internal variables equation, and the control instruction contains the operation instruction of the direct conversion from the equation and the communication command among the microprocessor cores.The real-time simulation operation can be carried out after initialization in the FRTDS platform.The numbers of the operations and the valid instructions for each microprocessor core are shown in Table 2.It can be seen that each microprocessor core can complete the real-time simulation operation of a 5 μs simulation step.Taking subnetwork 3 as an example, the number of ports is 4, the number of L and C is 3, the number of independent power supplies is 0, the number of unknown internal variables w is 1, the number of nonlinear components is 0, and the number of switches is 3.The numbers of the multivalued coefficients to be stored are, respectively, 16, 12 and 7 according to the exterior characteristic Equation ( 8), the equivalent source Equation ( 9), unknown internal variables Equation (10), and the total number x is 35.Each multivalued coefficient variable has eight possible values, and it can be seen that the total amount y of multivalued coefficients that need to be prestored is 280.Table 1 shows the numbers of the relevant data with the multivalued coefficients of the nine subnetworks.In the real-time simulation calculation of the whole low-voltage microgrid system, the numbers of ports, inductors and capacitors, independent power supplies, unknown internal variables w, nonlinear components and switches need to be known, thus the number x of the multivalued coefficients is 1196 and the total amount y of the multivalued coefficients is 17,135.At a step size of 5 µs, each microprocessor core can execute up to 1000 instructions.The equations of each subnetwork include the network equation to solve the input variables and the unknown internal variables equation, and the control instruction contains the operation instruction of the direct conversion from the equation and the communication command among the microprocessor cores.The real-time simulation operation can be carried out after initialization in the FRTDS platform.The numbers of the operations and the valid instructions for each microprocessor core are shown in Table 2.It can be seen that each microprocessor core can complete the real-time simulation operation of a 5 µs simulation step.In order to verify the accuracy of the proposed method, the same low-voltage microgrid simulation model with the same parameters is set up in a power systems computer aided design (PSCAD)/electromagnetic transients including DC (EMTDC) offline simulation tool.The simulation step is also set to 5 µs. Figure 16 shows the comparison of FRTDS and PSCAD simulation results of the photovoltaic-battery hybrid generation system after the change of light intensity and inverter active power command.
Energies 2017, 10, 1248 15 of 18 (PSCAD)/electromagnetic transients including DC (EMTDC) offline simulation tool.The simulation step is also set to 5 μs. Figure 16 shows the comparison of FRTDS and PSCAD simulation results of the photovoltaic-battery hybrid generation system after the change of light intensity and inverter active power command.
The initial light intensity of the system is 1000 W/m 2 , the active power of the inverter is 9.8 kW.After the system reaches steady state, the light intensity decreases to 800 W/m 2 at t = 1.5 s, and the active power is reduced to 6 kW at t = 3 s.As shown in Figure 17, the inverter output power is about 9.8 kW, and photovoltaic (PV) array output power is also 9.8 kW, so the battery does not output active power after the system reaches steady state, the best operation voltage of the PV array is about 356 V through the maximum power point tracking (MPPT) control.When the light intensity is reduced to 800 W/m 2 , the optimum operating point of the PV array is reduced to 343 V from 356 V because of the MPPT control, and the PV array output power is reduced to 7.6 kW.The output power of the inverter is still 9.8 kW.Thus, the battery is in the state of power generation and issued 2.2 kW active power.At t = 3 s, the inverter's active power is reduced to 6 kW, the current of the inverter is decreased accordingly, the battery is charged by the PV array, and the charging power is 1.6 kW.The DC bus voltage is continuously controlled to about 400 V in this process.
The waveforms of the filter outlet a phase current and inverter DC side voltage are presented in Figure 17 when a three-phase short-circuit fault at node 24 occurs in a photovoltaic-battery hybrid generation system.The initial light intensity of the system is 1000 W/m 2 , the active power of the inverter is 9.8 kW.After the system reaches steady state, the light intensity decreases to 800 W/m 2 at t = 1.5 s, and the active power is reduced to 6 kW at t = 3 s.As shown in Figure 17, the inverter output power is about 9.8 kW, and photovoltaic (PV) array output power is also 9.8 kW, so the battery does not output active power after the system reaches steady state, the best operation voltage of the PV array is about 356 V through the maximum power point tracking (MPPT) control.When the light intensity is reduced to 800 W/m 2 , the optimum operating point of the PV array is reduced to 343 V from 356 V because of the MPPT control, and the PV array output power is reduced to 7.6 kW.The output power of the inverter is still 9.8 kW.Thus, the battery is in the state of power generation and issued 2.2 kW active power.At t = 3 s, the inverter's active power is reduced to 6 kW, the current of the inverter is decreased accordingly, the battery is charged by the PV array, and the charging power is 1.6 kW.The DC bus voltage is continuously controlled to about 400 V in this process.As can be seen from Figure 17, there is a good agreement between the results of FRTDS and PSCAD, and the feasibility of the multivalued coefficient prestorage and block parallel method is verified.In the partial magnification of the waveforms, we can see that some of the transient features are not completely coincident, which is mainly due to the fact that: (1) the selected input variables are state variables, but they might change slightly sometimes; and (2) there are some errors because the state of the nonlinear element is determined by the local iterative algorithm.Considering that the whole waveform is consistent and the error is within 1%, it is shown that the multivalued coefficient prestorage and block parallel method proposed in this paper has high simulation accuracy.

Conclusions
In the paper, a novel multivalued coefficient prestorage and block parallel method for real-time simulation of a microgrid on FRTDS is proposed.It alleviates memory requirements, and reduces the computational burden.Neither of these results is obtained by compromising the accuracy of the simulation.Equivalent Thevenin circuit with the input current variables should be selected when the ports are connected in series, and an equivalent Norton circuit with the input voltage variables should be chosen when the ports are connected in parallel.It can reduce the number of simultaneous network equations.The input variables are tied to state variables of the circuit, which can limit the iterative calculation of determining the nonlinear elements status to the local network.It allows smaller time steps to be reached for the real-time simulation of a microgrid.The multivalued coefficient prestorage and block parallel method solves the contradiction in real-time simulation between data storage and computing speed, which meets the requirements for large-scale real-time simulation of a microgrid.The waveforms of the filter outlet a phase current and inverter DC side voltage are presented in Figure 17 when a three-phase short-circuit fault at node 24 occurs in a photovoltaic-battery hybrid generation system.
As can be seen from Figure 17, there is a good agreement between the results of FRTDS and PSCAD, and the feasibility of the multivalued coefficient prestorage and block parallel method is verified.In the partial magnification of the waveforms, we can see that some of the transient features are not completely coincident, which is mainly due to the fact that: (1) the selected input variables are state variables, but they might change slightly sometimes; and (2) there are some errors because the state of the nonlinear element is determined by the local iterative algorithm.Considering that the whole waveform is consistent and the error is within 1%, it is shown that the multivalued coefficient prestorage and block parallel method proposed in this paper has high simulation accuracy.

Conclusions
In the paper, a novel multivalued coefficient prestorage and block parallel method for real-time simulation of a microgrid on FRTDS is proposed.It alleviates memory requirements, and reduces the computational burden.Neither of these results is obtained by compromising the accuracy of the simulation.Equivalent Thevenin circuit with the input current variables should be selected when the ports are connected in series, and an equivalent Norton circuit with the input voltage variables should be chosen when the ports are connected in parallel.It can reduce the number of simultaneous network equations.The input variables are tied to state variables of the circuit, which can limit the iterative calculation of determining the nonlinear elements status to the local network.It allows smaller time steps to be reached for the real-time simulation of a microgrid.The multivalued coefficient prestorage and block parallel method solves the contradiction in real-time simulation between data storage and computing speed, which meets the requirements for large-scale real-time simulation of a microgrid.

Figure 2 .
Figure 2. Piecewise linearization and equivalent circuit of the nonlinear source: (a) current-voltage (I-U) curve and piecewise linearization; (b) equivalent circuit.

Figure 2 .
Figure 2. Piecewise linearization and equivalent circuit of the nonlinear source: (a) current-voltage (I-U) curve and piecewise linearization; (b) equivalent circuit.

Figure 2 .
Figure 2. Piecewise linearization and equivalent circuit of the nonlinear source: (a) current-voltage (I-U) curve and piecewise linearization; (b) equivalent circuit.

Figure 3 .
Figure 3. Piecewise linear representation of the nonlinear inductor and its equivalent circuit: (a) piecewise linearization; (b) equivalent circuit.

Figure 3 .
Figure 3. Piecewise linear representation of the nonlinear inductor and its equivalent circuit: (a) piecewise linearization; (b) equivalent circuit.

Figure 4 .
Figure 4. Fault model and its equivalent circuit model for transient simulation: (a) fault model; (b) equivalent circuit model for transient simulation.

Figure 4 .
Figure 4. Fault model and its equivalent circuit model for transient simulation: (a) fault model; (b) equivalent circuit model for transient simulation.

Figure 4 .
Figure 4. Fault model and its equivalent circuit model for transient simulation: (a) fault model; (b) equivalent circuit model for transient simulation.

Figure 7 .
Figure 7. Two two-port hybrid equivalents connected to solve interface input currents/voltages.
and the ports are connected in series with current as input variables.There is KVL constraint: u u u u , and the ports are connected in parallel with voltage as input variables as shown in Figure8b.All the port input variables can be obtained by solving only one independent current or voltage variable.

Figure 7 .
Figure 7. Two two-port hybrid equivalents connected to solve interface input currents/voltages.

Figure 8 .
Figure 8. Series and parallel port equivalents connected to solve interface input currents/voltages: (a) series port equivalents connected to solve interface input currents; (b) parallel port equivalents connected to solve interface input voltages.

Figure 8 .
Figure 8. Series and parallel port equivalents connected to solve interface input currents/voltages: (a) series port equivalents connected to solve interface input currents; (b) parallel port equivalents connected to solve interface input voltages.

Figure 13 .
Figure 13.Configuration of the photovoltaic-battery hybrid generation system.

Figure 14 .
Figure 14.The photovoltaic-battery hybrid generation system with partitioning into nine subnetworks.

Figure 15 .
Figure 15.The nine subnetworks connected to solve interface input currents/voltages.

Figure 15 .
Figure 15.The nine subnetworks connected to solve interface input currents/voltages.

Figure 16 .Figure 16 .
Figure 16.Simulation results of FPGA-based real-time digital simulation (FRTDS) and power systems computer aided design (PSCAD): (a) inverter output power; (b) photovoltaic array output power; (c) Figure 16.Simulation results of FPGA-based real-time digital simulation (FRTDS) and power systems computer aided design (PSCAD): (a) inverter output power; (b) photovoltaic array output power; (c) battery output power; (d) photovoltaic output voltage; (e) inverter Direct Current (DC) side voltage; (f) The filter output a phase current.

Figure 17 .
Figure 17.Simulation results of FRTDS and PSCAD: (a) the filter output a phase current; (b) inverter DC side voltage.

Figure 17 .
Figure 17.Simulation results of FRTDS and PSCAD: (a) the filter output a phase current; (b) inverter DC side voltage.
[A 1 A 2 ] is the coefficient matrix associated with the input variables, [B 1 B 2 ] is the coefficient matrix associated with the independent sources, and [C 1 C 2 ] is the coefficient matrix associated with the historical sources.The elements of [A 1 A 2 ], [B 1 B 2 ], and [C 1 C 2 ] need to be prestored as the multivalued coefficients, which are constants after the state of the circuit element is determined.

Table 1 .
The numbers of the relevant data with multivalued coefficients.

Table 2 .
The numbers of the operations and the valid instructions.

Table 1 .
The numbers of the relevant data with multivalued coefficients.

Table 2 .
The numbers of the operations and the valid instructions.