2.2.1. PV and Battery Power Generation Systems

Due to how the topology and control of PV and battery systems are similar, this section will describe and establish the impedance model of the power generation system simultaneously. As shown in Figure 3, in the Boost circuit, *u*PV is the output voltage of the PV side; *i*PV is the output current of the PV side; *i*LPV is the inductor current of the PV side; *u*Bus and *i*Bus are the DC bus voltage and current; *L*PV is the energy storage inductor of the PV side; *C*PV is the support capacitor of the PV side; *C*Bus1 is the DC-side capacitor of the RPC; *d*PV is the duty cycle of the Boost converter's switching tube, *d*\*PV = 1 − *d*PV; and *u*\*Bus and *i*\*LPV are the reference signals for their corresponding variables. In the DC/DC circuit, *u*Bat is the output voltage of the lithium battery side, *i*Bat is the output current of the lithium battery, *L*Bat is the energy storage inductor, *C*Bus2 is the DC-side capacitor of the RPC, *d*Bat is the duty cycle of the *T*<sup>1</sup> switching tube, and the *T*<sup>2</sup> switching tube is complementary to the *T*<sup>1</sup> switching tube. In addition, *u*\*Bus and *i*\*Bat are the reference signals for their corresponding variables.

**Figure 3.** Topology and control of DC converters under mode three. (**a**) PV boost converter; (**b**) battery DC/DC converter.

The state equations for the PV and battery systems are given by Equations (1) and (2), respectively.

$$\begin{cases} & L\_{\rm PV} \frac{\text{di}\_{\rm PV}}{\text{d}t} + R\_{\rm PV} \cdot i\_{\rm LPV} = \mu\_{\rm PV} - d\_{\rm PV}^{\*} \cdot u\_{\rm Bus} \\ & C\_{\rm Bus1} \frac{\text{du}\_{\rm Bus}}{\text{d}t} = d\_{\rm PV}^{\*} \cdot i\_{\rm LPV} - i\_{\rm Bus} \\ & i\_{\rm PV} = i\_{\rm LPV} + C\_{\rm PV} \frac{\text{d}u\_{\rm PV}}{\text{d}t} \end{cases} \tag{1}$$
 
$$\begin{cases} & L\_{\rm Bat} \frac{\text{di}\_{\rm Bat}}{\text{d}t} + R\_{\rm Bat} \cdot i\_{\rm Bat} = u\_{\rm Bat} - d\_{\rm Bat} u\_{\rm Bus} \\ & C\_{\rm Bus2} \frac{\text{d}u\_{\rm Bat}}{\text{d}t} = d\_{\rm Bat} i\_{\rm Bat} - i\_{\rm Bus} \end{cases} \tag{2}$$

With small signal linearizing Equations (1) and (2), we obtain Equations (3) and (4):

$$\begin{cases} \left(s\text{L}\_{\text{PV}} + R\_{\text{PV}}\right) \cdot \Delta i\_{\text{LPV}} = -D\_{\text{PV}}^{\*} \cdot \Delta u\_{\text{Bus}} + \mathcal{U}\_{\text{Bus}} \cdot \Delta d\_{\text{PV}} + \Delta u\_{\text{PV}}\\\ s\text{C}\_{\text{Bus}1} \cdot \Delta u\_{\text{Bus}} = D\_{\text{PV}}^{\*} \cdot \Delta i\_{\text{LPV}} - I\_{\text{LPV}} \cdot \Delta d\_{\text{PV}} - \Delta i\_{\text{Bus}}\\\ \Delta i\_{\text{PV}} = \Delta i\_{\text{LPV}} + s\text{C}\_{\text{PV}} \cdot \Delta u\_{\text{PV}} \end{cases} \tag{3}$$

$$\begin{cases} \left( sL\_{\rm Bat} + R\_{\rm Bat} \right) \Delta i\_{\rm Bat} = -D\_{\rm Bat} \cdot \Delta \mu\_{\rm Bus} - \mathcal{U}\_{\rm Bus} \cdot \Delta d\_{\rm Bat} + \Delta \mu\_{\rm Bat} \\\ s \mathcal{C}\_{\rm Bus2} \Delta \mu\_{\rm Bus} = D\_{\rm Bat} \cdot \Delta i\_{\rm Bat} + I\_{\rm Bat} \cdot \Delta d\_{\rm Bat} - \Delta i\_{\rm Bus} \end{cases} \tag{4}$$

According to the superposition theorem, the transfer functions between each variable of the PV and battery systems can be obtained from Equations (3) and (4), and the specific expressions are shown in Appendix A, where Equation (A1) includes the transfer function expansion of the PV system, and Equation (A2) includes the transfer function expansion of the battery system. Then, based on the control strategies of the PV and battery systems, draw the closed-loop block diagrams of the small signal transfer functions as shown in Figure 4.

**Figure 4.** Small-signal closed-loop control block diagram. (**a**) PV transfer function; (**b**) battery transfer function.

Using the Mason's formula, the output impedances of the PV and battery systems are, respectively:

$$Z\_{\rm PV}(\mathbf{s}) = \frac{\Delta u\_{\rm Bus}}{\Delta i\_{\rm Bus}} = \frac{G\_{\rm PVBB} + G\_{\rm PVBB} G\_{\rm PV\bar{i}} G\_{\rm PV\perp} - G\_{\rm PV\perp B} G\_{\rm PV\bar{i}} G\_{\rm PV\bar{B}\bar{P}}}{1 + G\_{\rm PV\bar{i}} G\_{\rm PV\perp P} + G\_{\rm PV\bar{u}} G\_{\rm PV\bar{i}} G\_{\rm PV\bar{B}\bar{P}}} \tag{5}$$

$$Z\_{\rm Bat}(\mathbf{s}) = \frac{\Delta\mu\_{\rm Bus}}{\Delta i\_{\rm Bus}} = \frac{G\_{\rm BiB} + G\_{\rm BiB}G\_{\rm Bi}G\_{\rm BiB} - G\_{\rm Bii}G\_{\rm BdB}G\_{\rm Bi}}{1 + G\_{\rm Bi}G\_{\rm iBd} + G\_{\rm Bu}G\_{\rm Bi}G\_{\rm BdB}}\tag{6}$$

#### 2.2.2. Locomotive Network and Railway Power Conditioner

In the study of the low-frequency stability of locomotive network coupling systems, two-level locomotives are often considered. The China railway high-speed 3 (CRH3) locomotive with four power units is the research object of this paper. In [32], the impedance expression of a single power unit is presented. After converting the impedance of the China Railway High-speed 3 (CRH3) locomotive to the RPC AC side, the expression is:

$$\mathcal{Z}\_{\rm CRH3}(\mathbf{s}) = \frac{k\_{\rm a}^{\prime} k\_{\rm b}^{\prime} \mathcal{L} l\_{\rm s2}}{4 l\_{\rm s2}} = \frac{k\_{\rm a}^{\prime} k\_{\rm b}^{\prime} \left[1 + \mathcal{G}\_{\rm F} \mathcal{G}\_{\rm PWM} \mathcal{G}\_{\rm Ln} (\mathcal{G}\_{\rm CRH3u} \mathcal{G}\_{\rm switch} \mathcal{G}\_{\rm switch} + \mathcal{G}\_{\rm CRH3u})\right]}{4 \mathcal{G}\_{\rm Ln} (1 - \mathcal{G}\_{\rm delay} \mathcal{G}\_{\rm PWM})} \tag{7}$$

where the transformation ratio of the onboard transformer is *k*<sup>a</sup> and the transformation ratio of the RPC-side step-up transformer is *k*b.

At present, the equivalent circuit model is commonly used in low-frequency stability analysis of traction networks [31] and its expression corresponds to the RPC AC side:

$$Z\_s(s) = k\_b^2 \begin{bmatrix} R\_s + sL\_s & -\omega\_0 L\_s \\ \omega\_0 L\_s & R\_s + sL\_s \end{bmatrix} \tag{8}$$

The equivalent impedance of the locomotive network converted to the RPC AC side is recorded as *Z*load (*s*), where *Z*load (*s*) = *Z*<sup>s</sup> (*s*) + *Z*CRH3 (*s*).

RPC impedance modeling is conducted on the DC port of the inverter, that is, the input impedance of the RPC is studied. Considering that the back-to-back inverter in the RPC has a symmetrical structure and the impedance calculation and modeling method are similar, the β power supply arm is selected as the research object. For the emergency power supply scheme, the specific control method and topology of the inverter are shown in Figure 5.

**Figure 5.** RPC inverter main circuit topology and control.

In the above figure, *Z*load is the equivalent impedance of the locomotive network converted to the RPC AC side; *U*<sup>m</sup> is the amplitude of the modulated wave signal; and sin*ωt* is the sine wave frequency reference signal. Linearization is performed on the DC side and the variables are decomposed into a closed-loop impedance model in the dq coordinate system. According to the main circuit, the state equation under dq decoupling can be derived by:

$$\begin{cases} \begin{aligned} \mathcal{L}\_{\text{inv}} \frac{\text{di}\_{\text{Id}}}{\text{d}t} &= \mathcal{U}\_{\text{Id}} - \mathcal{U}\_{\text{0d}} + \omega \mathcal{L}\_{\text{inv}} \mathcal{I}\_{\text{Lq}} \\ \mathcal{L}\_{\text{inv}} \frac{\text{di}\_{\text{d}}}{\text{d}t} &= \mathcal{U}\_{\text{Iq}} - \mathcal{U}\_{\text{0q}} - \omega \mathcal{L}\_{\text{inv}} \mathcal{I}\_{\text{Ld}} \\ \mathcal{C}\_{\text{inv}} \frac{\text{du}\_{\text{0d}}}{\text{d}t} &= \mathcal{I}\_{\text{Ld}} - \mathcal{I}\_{\text{0d}} + \omega \mathcal{C}\_{\text{inv}} \mathcal{I}\_{\text{0q}} \\ \mathcal{C}\_{\text{inv}} \frac{\text{du}\_{\text{0q}}}{\text{d}t} &= \mathcal{I}\_{\text{Lq}} - \mathcal{I}\_{\text{0q}} - \omega \mathcal{C}\_{\text{inv}} \mathcal{I}\_{\text{0d}} \\ \mathcal{U}\_{\text{0d}} &= \mathcal{Z}\_{\text{load}} \cdot \mathcal{I}\_{\text{0d}} \\ \mathcal{U}\_{\text{0q}} &= \mathcal{Z}\_{\text{load}} \cdot \mathcal{I}\_{\text{0q}} \end{aligned} \tag{9}$$

Based on the basic derivation process of SPWM, the relationship between the output voltage *U*<sup>i</sup> before filtering and the modulation wave *U*<sup>m</sup> can be obtained using the switch average period method:

$$\begin{cases} \text{ } \mathcal{U}\_{\text{id}} = 0.5 \mathcal{U}\_{\text{dc}} \mathcal{U}\_{\text{md}} / \mathcal{U}\_{\text{tri}}\\ \text{ } \mathcal{U}\_{\text{iq}} = 0.5 \mathcal{U}\_{\text{dc}} \mathcal{U}\_{\text{mq}} / \mathcal{U}\_{\text{tri}} \end{cases} \tag{10}$$

In Equation (10), *U*tri is the amplitude of the triangular carrier wave, and the remaining symbols are the system variables under dq decoupling. In the dual-loop control strategy, the response speed of the inner loop is much faster than that of the outer loop. This article assumes that the response of the current inner loop is approximately the same within the average time of a single switching cycle. Due to the PI controller used in the voltage outer loop, the input–output relationship of the controller is:

$$\begin{cases} \left(k\_{\rm invp} + k\_{\rm invi}/s\right) \cdot \left(\mathcal{U}\_{\rm vd}^{\*} - \mathcal{U}\_{\rm 0d}\right) = \mathcal{U}\_{\rm md} \\ \left(k\_{\rm invp} + k\_{\rm invi}/s\right) \cdot \left(\mathcal{U}\_{\rm vq}^{\*} - \mathcal{U}\_{\rm 0q}\right) = \mathcal{U}\_{\rm mq} \end{cases} \tag{11}$$

where *U*\*vd and *U*\*vq are the product of the effective value of the voltage outer-loop output and the reference sine signal; *k*invp and *k*invi are proportional integral parameters of the voltage outer-loop PI controller.

According to the input power and output power of the RPC that are equal:

$$\mathcal{U}\_{\rm dc} I\_{\rm dc} = 1.5(\mathcal{U}\_{\rm 0d} I\_{\rm 0d} + \mathcal{U}\_{\rm 0q} I\_{\rm 0q}) \tag{12}$$

After processing the small signal in Equation (12), the following is obtained:

$$
\begin{bmatrix} I\_{\rm dc} \\ \mathcal{U}\_{\rm dc} \end{bmatrix}^{\rm T} \begin{bmatrix} \Delta u\_{\rm dc} \\ \Delta i\_{\rm dc} \end{bmatrix} = 1.5 \left( \begin{bmatrix} \mathcal{U}\_{\rm bd} \\ \mathcal{U}\_{\rm 0q} \end{bmatrix}^{\rm T} \begin{bmatrix} \Delta i\_{\rm bd} \\ \Delta i\_{\rm 0q} \end{bmatrix} + \begin{bmatrix} I\_{\rm bd} \\ I\_{\rm 0q} \end{bmatrix}^{\rm T} \begin{bmatrix} \Delta u\_{\rm bd} \\ \Delta u\_{\rm 0q} \end{bmatrix} \right) \tag{13}
$$

Equations (9) to (11) can be substituted into Equation (13) to obtain the input admittance on the DC side of the β arm inverter:

$$Y\_{\rm inv} = \frac{\Delta i\_{\rm dc}}{\Delta u\_{\rm dc}} = \frac{3}{2lI\_{\rm dc}} \left( \begin{bmatrix} \iota I\_{0\rm d} \\ \iota I\_{0\rm q} \end{bmatrix}^{\rm T} \left( G\_{\rm Zload} + G\_{\rm Zload}^{\rm T} \right) \begin{bmatrix} \frac{\Delta u\_{\rm bd}}{\Delta u\_{\rm dc}} \\ \frac{\Delta u\_{\rm dc}}{\Delta u\_{\rm dc}} \end{bmatrix} \right) - \frac{I\_{\rm dc}}{lI\_{\rm dc}} \tag{14}$$

As output impedance is the reciprocal of admittance, it can be inferred that *Z*inv = 1/*Y*inv and the variables in Equation (14) are as follows:

$$\begin{cases} \begin{bmatrix} \frac{\partial \mathbf{I\_{kd\_d}}}{\Delta \mathbf{I\_{kd\_d}}}\\ \frac{\partial \mathbf{u\_{0d}}}{\partial \mathbf{u\_{dc}}} \end{bmatrix} = 0.5[G\_{\text{Linv}}(\mathbf{G\_{\text{Linv}}} + G\_{\text{Zload}}) + E - 0.5lI\_{\text{dc}}G\_{\text{T}}]^{-1} \begin{bmatrix} \mathbf{U\_{md}}\\ \mathbf{U\_{mq}} \end{bmatrix} \\\\ \text{G}\_{\text{T}} = \begin{bmatrix} s & -\omega\\ \omega & s \end{bmatrix}^{-1} \begin{bmatrix} -(sK\_{\text{inv}\text{p}} + K\_{\text{inv}\text{i}}) & \omega K\_{\text{inv}\text{p}}\\ -\omega K\_{\text{inv}\text{p}} & -(sK\_{\text{inv}\text{p}} + K\_{\text{inv}\text{i}}) \end{bmatrix} \\\ \text{G}\_{\text{Linv}} = \begin{bmatrix} sL\_{\text{inv}} & -\omega L\_{\text{inv}}\\ \omega L\_{\text{inv}} & sL\_{\text{inv}} \end{bmatrix} \\\ \text{G}\_{\text{Cinv}} = \begin{bmatrix} sC\_{\text{inv}} & -\omega C\_{\text{inv}}\\ \omega C\_{\text{inv}} & sC\_{\text{inv}} \end{bmatrix} \\\ \text{G}\_{\text{Zload}} = \begin{bmatrix} Z\_{\text{load}} & 0\\ 0 & Z\_{\text{load}} \end{bmatrix}^{-1} \end{cases} \tag{15}$$

2.2.3. Verification of Established Impedance Model

To verify the established one-dimensional impedance model, the disturbance injection method was used for the simulation measurements [33]. By injecting disturbance signals with specific frequencies during the steady-state state of the system and using the Fourier transform to process these signals, the impedances of the three ports of the "PV–battery locomotive network" can be simultaneously measured at specific frequencies, and it should be noted that the impedances of the three ports of the "PV-battery locomotive network" refer to the output impedances *Z*PV towards the PV boost circuit port, *Z*Bat towards the battery DC/DC circuit port, and *Z*inv towards the inverter DC port. We will validate whether the theoretical deduced values of these impedances in the frequency domain match the simulation results.

Validation results are shown in Figure 6, the first row shows the corresponding impedance function's Bode plot, and the second row is the corresponding phase plot, with the blue curves representing the theoretical values of the model and the red circles representing the actual measured values. It can be seen that the modeling of the PV and battery system's output impedance and RPC's input impedance is consistent with the simulation test results.

**Figure 6.** Three-port impedance verification.

*2.3. Low-Frequency Stability Analysis of Proposed Emergency Power Supply System* 2.3.1. Stability Judgment Based on Generalized Nyquist Criterion

According to the mathematical model established earlier, the impedance equivalent circuit of the PV and battery locomotive traction is shown in Figure 7a, where *Z*<sup>0</sup> is equal to the parallel value of *Z*PV and *Z*Bat, and *Z*<sup>g</sup> is the impedance of the traditional locomotive network system source side. Due to the PI controller parameters of the converters in PV and battery systems that affect the output impedance of the system, *Z*0, *Z*PV, and *Z*Bat are different from *Z*<sup>g</sup> because the value of *Z*<sup>g</sup> depends on factors such as the material of the transmission line, its length, and the capacity of transformers, which are difficult to adjust in practical operation. In comparison, the values of *Z*0, *Z*PV, and *Z*Bat are related to the parameters of the PI controller, and adjusting the control parameters of the PI controller is less costly as well as more simple and fast. The introduction of a controllable equivalent impedance may exacerbate the LFO of the system or prevent system instability caused by impedance mismatch; this characteristic makes it particularly important to explore the impact of the control parameters and number of PV and battery parallel connections on the system's low-frequency stability. This section uses stability judgment methods based on impedance ratio criteria to conduct the research [27,28].

Firstly, it must be ensured that there are no RHP poles for *Z*0, *Z*PV, *Z*Bat, and *Y*inv and that each subsystem can operate independently and stably in the simulation. Under different traction modes, applying the one-dimensional impedance model established earlier, the Nyquisit curve based on the impedance ratio of *Z*0/*Z*inv, *Z*PV/*Z*inv, and *Z*Bat/*Z*inv can be obtained. Taking *Z*0/*Z*inv as an example, the expression for the system's impedance ratio can be represented by a transfer function of *G*(*s*) = *Z*0/*Z*inv. The stability of the system can be measured by plotting the Nyquist curve of *G*(*s*) and determining whether it encircles (−1, *j*0). Additionally, the variations in parameter values based on the impedance model discussed earlier can affect *Z*0, thus, modifying the Nyquist curve of *G*(*s*). Therefore, the parameters of each variable may have some effect on the stability of the system. Meanwhile, due to the lack of a clear concept of the stability margin in the impedance ratio criterion for MIMO equivalent systems, and to showcase the analysis process concisely and clearly, the critical amplitude margin (CAM) is defined to quantitatively analyze the system stability in order to reveal sensitive parameters and their influence laws. If CAM ≥ 0, the system is

stable; if CAM < 0, the system is unstable. The Nyquist plot of the impedance ratio transfer function for this system may change due to the effect of different parameters, leading to two different scenarios, as shown in Figure 7b.

**Figure 7.** (**a**) System impedance equivalent circuit; (**b**) CAM value scenario diagram.

Scenario 1: When the curve surrounds (−1, j0), the distance from the intersection point of the real axis of the curve to (−1, j0) is expressed as η1, CAM = −η1;

Scenario 2: When the curve does not surround (−1, j0), the corresponding distance is expressed as η2, CAM = η2.

In summary, changes in system variables can affect the Nyquist curve of *G*(*s*), thus, impacting the stability of the system. Such an impact can be translated into changes in CAM values to quantify the effect of a particular variable on the stability of the system. To further reveal the main controller parameters that affect the system stability, this article defines the sensitivity of the parameters as:

$$\varepsilon = \frac{\text{CAM}\_{\text{n}}}{(k\_{\text{s}} - k\_{\text{y}})/k\_{\text{y}}} \times 100\% \tag{16}$$

where *k*<sup>y</sup> is a certain original parameter in Table 1, *k*<sup>s</sup> is the value of the critical instability parameter with the smallest change from *k*y, (*k*<sup>s</sup> − *k*y)/*k*<sup>y</sup> is the multiple of the minimum change in the original parameters when the system is in critical instability, and CAMn is the critical amplitude margin corresponding to the original parameter *k*<sup>y</sup> under different operating conditions. The introduction of CAMn is aimed at correcting the sensitivity under different traction conditions (PV traction, battery traction, PV, and battery co-traction). This expression can be simply understood as the sensitivity to the electrical instability of the nearest value of the original parameter.


**Table 1.** Original simulation system parameters.

2.3.2. Influence of DC Converter and RPC Parameters on the Low-Frequency Stability of the Proposed System

Based on the method proposed in Section 2.3.1, Figure 8a shows the electrical sensitivity bar chart of the control parameters of the DC converter and RPC inverter under three different traction modes when the ratio between the PV module, battery module, and CRH3 locomotive is 1:1:1. Due to the sensitivity of the inner-loop control parameters of the system under this scheme being less than 1%, the Nyquist curve hardly changes with the changes in the parameters, and the corresponding test shows that the system remains

stable even after a large number of changes in the inner-loop control parameters; therefore, it is not included in the research object. The bar chart shows that the proportional gain of the outer ring is more sensitive to the system stability and the combined traction of the PV and battery will reduce the parameter sensitivity.

**Figure 8.** (**a**) Parameter sensitivity of the RPC and DC converters; (**b**) CAM change law of the control parameters.

The corresponding parameter influence law is shown in Figure 8b. Within the range of −0.9 to 2 times the original parameter value, the proportional gain of the outer ring of each converter exhibits monotonic, identical laws under different operating conditions, and the sensitivity (i.e., slope) is relatively high. However, the sensitivity of the integral gain of the outer ring is low, and the impact of the integral gain on the stability of the system exhibits diversification and non-monotonic characteristics under different traction conditions. Therefore, the outer-loop proportional gain of the system can be more easily used to control or prevent system impedance mismatch instability. The specific parameter adjustment criterion and case verification are expounded in Section 2.3.3.

In addition, LFO in locomotive network systems always occurs when the Nyquist curve of the impedance ratio (*Z*g/*Z*CRH3) passes through (−1, j0). To investigate whether the occurrence of LFO under the topology and control scheme proposed in this paper is consistent with the theoretical analysis of locomotive network systems when the locomotive is separately tractioned by the PV or battery, the impedance ratio expression for the PV traction locomotive is *Z*PV/*Z*inv, and the expression for battery traction locomotive is *Z*Bat/*Z*inv, where *Z*PV, *Z*Bat, and *Z*inv come from Equations (5), (6), and (14), and these impedance ratio expressions are used to plot Nyquist curves under different traction conditions. Then adjust the high sensitivity parameters *k*PVup and *k*Bup so that the Nyquist curve just surrounds (−1, j0), thereby obtaining the Nyquist curve shown in Figure 9 and the test waveform shown in Figure 10. The results show that when the Nyquist curve of the impedance ratio for the system described in this paper precisely intersects the (−1, j0) point, LFO at around 9.6 Hz also occurs. The principle and testing results are consistent under other parameters in Figure 8 but due to the large number of parameters, they are not listed in detail.

**Figure 9.** Nyquist curve before and after changes in high-sensitivity parameters: (**a**) working conditions of PV locomotive traction; (**b**) working conditions of battery locomotive traction.

**Figure 10.** (**a**) Experimental waveform of the critical instability of photovoltaic locomotive traction; (**b**) experimental waveform of the critical instability of battery pack locomotive traction.
