Multi‐harmonic oscillation and stability analysis of double‐input buck/buck‐boost inverter

Funding information National Natural Science Foundation of China, Grant/Award Number: 51577141; Natural Science Basic Research Plan in Shaanxi Province of China, Grant/Award Number: 2021JM Abstract In this study, multi-harmonic oscillation behaviour and stability problem in double-input buck/buck–boost inverter are investigated both theoretically and experimentally. First, the observer-pattern model is obtained by the local orthogonal and autonomous transformation in order to eliminate the time-variance effect from both fundamental component and the ‘hidden’ second-harmonic in the double-input cascade system. Based on the proposed model, bifurcation analysis method is employed to reveal the underlying mechanism of the multi-harmonic oscillation behaviour and identify the dominant harmonics. It is shown that the occurrence of Hopf bifurcation results in the double-adding harmonic oscillation, and the interaction between the double-adding harmonic components and the ‘hidden’ second harmonic component is also responsible for the multi-harmonic oscillation behaviour. Meanwhile, the occurrence of the multi-harmonic oscillations leads to the decrease of the system efficiency. For stability enhancement, the harmonic-stability sensitivity is presented to evaluate the influence of the key circuit parameters on the system stability quantitatively. In addition, the stability boundaries of the double-input inverter for these key parameters are provided to guide the system optimal design. Finally, the effectiveness of the above analysis is verified by those experimental results.


INTRODUCTION
For some remote communities and isolated residences far away from grid distribution lines, the standalone hybrid energy system with multiple energy sources is a more environmentfriendly and economical candidate than the grid-connected power system [1][2][3]. In this hybrid energy systems, multiple energy sources including solar energy, wind energy and even storage energy are connected to a common AC load or the grid via a number of inverters [4,5]. Due to some considerable advantages such as the integrated topology, centralised control and low-manufacturing size, multi-input inverters are capable of replacing several single-input inverters by combining different energy sources in the DC end instead of the AC end [6][7][8].
Since the multi-input inverters can enhance the system reliability and reduce the cost, they have attracted more and more attention recently [9][10][11][12]. Different from single-input single-output (SISO) inverters, multi-input inverters need the appropriate power management This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. IET Power Electronics published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology strategy to balance the power from different sources and regulate the output voltage, and significantly, they are often designed to operate under multiple topology structures and multiple operating modes [4][5][6][7][8][9][10][11][12][13]. As a consequence, stronger nonlinearity and more significant coupling effect become the salient characteristics of multi-input inverters, and they may exhibit complex oscillation in the form of multiple harmonics. These non-linear behaviours may degrade the system performance such as power quality, conversion efficiency and EMI. To avoid the occurrence of this phenomenon, it is essential to investigate the inherent mechanism of the multi-harmonic oscillation and identify the dominant harmonics.
Actually, harmonic stability is one of the main concerns in power converters, especially in inverters [14][15][16][17][18]. Much research effort has been devoted to investigating the harmonic oscillation and stability enhancement of inverter system in different applications. For such non-linear time-varying system as inverters, the non-linear analysis methods are powerful to analyse non-linear behaviours and system stability [19,20]. In [21], the Lyapunov-based large-signal stability study of inverters in microgrids is introduced. In [22], the Poincaré map is used to analyse the stability of the inverter based on Filippov's theory. However, these non-linear methods have computational complexity and are not directly applicable to the design of the inverters. In practice, such linear techniques as impedance-based method and eigenvalue method are widely used in stability analysis of inverters for engineers [23,24]. Since the seminal work of Middlebrook [25], the impedance-based method has become a very effective tool to analyse the stability of cascade systems in the frequency domain. In [26], an impedance-based small-signal stability criterion of the grid-connected system was proposed, which can assess the system stability based on the source-load impedance ratio. In [27], the impedance-based method is also used to analyse the influence of coupling effect and phaselocked loop on the stability of the inverter. Nevertheless, in the impedance-based method, the observability of some states is limited because of the black-box modelling of converters. Moreover, the impedance-based method cannot satisfy the demand in assessing the influence of the parameters on the harmonic stability quantitatively. Unlike the impedance-based method, the eigenvalue method based on state-space model provides a global view on the system dynamics in time-domain and enables to do further analysis for the effect of parameters. In [28], the eigenvalue trajectories of Jacobian matrix are calculated to assess the harmonic instability of the power conversion system. In [29], the non-linear dynamics and stability of a DC-AC inverter is investigated by using eigenvalue method. As a matter of fact, those previous contributions are beneficial to the judgement of the system stability, but little attention is paid to the inherent mechanism of the harmonic oscillations. In addition, the eigenvalue method widely used in stability analysis of SISO inverters cannot be directly applied to multi-harmonic instability analysis of multi-input cascade inverters. The main reason is the time-variance effect from the 'hidden' second harmonic in the source stage of multi-input inverters, which leads to much trouble in obtaining the time-invariant state-space model. Thus, it is necessary to develop effective modelling and analysis method to reveal the occurrence mechanism of the multiharmonic instability problem and assess the influence of the parameters.
In this study, a double-input buck/buck-boost inverter will be taken as an example to investigate the multi-harmonic oscillation and stability problem by using the proposed observerpattern model and bifurcation analysis. By the local orthogonal and autonomous transformation, the time-variance effect from both fundamental component and second-harmonic one is eliminated and the observer-pattern model is obtained. Then both bifurcation and loss analyses are performed to reveal the underlying mechanism of the multi-harmonic oscillation and identify the dominant harmonics. Particularly, the harmonic-stability sensitivity is introduced to evaluate the influence of the parameters on the system stability quantitatively. On this basis, the key parameters are identified and the stability boundaries are provided to guide the system optimal design. The rest of this study is organised as follows. The operating principles of the double-input buck/buck-boost inverter will be introduced in Section 2. The observer-pattern model will be developed in Section 3. In Section 4, the multi-harmonic oscillation phenomenon will be investigated by numerical simulations, bifurcation and loss analyses. In Section 5, the harmonicstability sensitivity will be presented to evaluate the influence of some parameters on the system stability. Finally, these theoretical results will be verified by circuit experiments.

Operating principle
The double-input buck/buck-boost inverter consists of two standalone converters, that is, one front-end double-input buck/buck-boost DC-DC converter and one single-phase Hbridge inverter, as shown in Figure 1. In the source-stage, the two input ports can deliver power from different voltage sources E 1 and E 2 to the downstream inverter. Two input ports are defined as port 1 # and port 2 # , respectively, and the downstream inverter is defined as port 3 # . Switches S 1 and S 2 are connected to E 1 and E 2 , respectively. In order to obtain regulated output voltage, the voltage-mode control is applied to the switch S 1 . Meanwhile, the averaged current-mode control is applied to the power switch S 2 for the faster dynamic response. Thus, the PWM signals of the source-stage subsystem are generated by comparing control signal v ref (t ) and v con1 (t ) with a sawtooth ramp voltage v ramp (t ). The control signals v ref (t ) and v con1 (t ) are given as where K 1 , K 2 are the proportional coefficients, T 1 , T 2 are the integral time constants in the controllers, V ref1 is the reference voltage, and K s1 , K s2 are the sensor gains of the output voltage v o1 and inductor current i L1 , respectively. The sawtooth ramp signal v ramp (t ) is given by where T s1 is the clock cycle of source-stage subsystem, and V m1 is the amplitude of v ramp (t ).
In the load-stage converter, the single-phase H-bridge inverter adopts the voltage control strategy, where the SPWM signal is achieved by comparing the control signal v con2 (t ) with the independent triangular voltage v tri (t ). The expression for the where K s3 is the sensor gain of the output voltage v o2 , the sinusoidal reference voltage v ref2 (t ) is V ref2 sin( t )( = 2 f l ), K 3 and T 3 are the proportional coefficient and integral time constant of the PI 3 controller, respectively. The triangular voltage v tri (t ) can be given as where T s2 is the clock cycle of load-stage and V m2 is the amplitude of v tri (t ).
Supposing that the double-input inverter operates in continuous conduction mode, eight switch states are possible within one switching cycle. Thus, the dynamics of the power stage can be governed by the following exact switch equations: Configuration 1: S 1 , S 2 , S 3 and S 6 turn on, S 4 and S 5 turn off Configuration 2: S 1 , S 2 , S 4 and S 5 turn on, S 3 and S 6 turn off Configuration 3: S 1 , S 3 and S 6 turn on, S 2 , S 4 and S 5 turn off Configuration 4: S 1 , S 4 and S 5 turn on, S 2 , S 3 and S 6 turn off Configuration 5: S 2 , S 3 and S 6 turn on, S 1 , S 4 and S 5 turn off Configuration 6: S 2 , S 4 and S 5 turn on, S 1 , S 3 and S 6 turn off Configuration 7: S 3 and S 6 turn on, S 1 , S 2 , S 4 and S 5 turn off Configuration 8: S 4 and S 5 turn on, S 1 , S 2 , S 3 and S 6 turn off

Non-linear averaged equations
Since there are more than one control variable in the doubleinput buck/buck-boost inverter, the topology sequences within one switching cycle are variable. Regardless of the topology sequence this system is running on, its exact switch equations can be presented by the following expression: where S 2 turn off and S 3 and S 6 turn on, S 4 and S 5 turn off 0 S 4 and S 5 turn on, S 3 and S 6 turn off .
After averaging the state space equations (Equation 14) in one switching cycle, the state space-averaged model of the power stage is obtained as where d 1 and d 2 are the duty ratios of switch S 1 and S 2 in the source-stage subsystem, d 3 is the duty ratio of switch S 3 and S 6 , and thus (1 − d 3 ) is the duty ratio of S 4 and S 5 in the load-stage subsystem.
Taking a derivative of each side of Equations (1) and (2) corresponding to t , the models of controller can be expressed as The unified expressions of the two PWM modules in the source-stage converter can be displayed as For the load-stage inverter, taking a derivative of each side of Equation (4) corresponding to t , one obtains The unified expression of the SPWM module in the loadstage inverter can be displayed as Taking Equations (15), (16), (17), (18), (19), (20) and (21) into account, the entire circuit system is described by the following state equations: From Equation (22), it is clear that the non-linear averaged equations of the double-input inverter is divided into the two parts: the first four equations are used to describe the dynamics of DC state variables with 'hidden' second-harmonic component in the front-end converter. The fifth to the seventh equations are used to describe the dynamics of AC state variables with fundamental frequency component in H-bridge inverter.

OBSERVER-PATTERN MODEL
Obviously, the non-linear averaged equations in Equation (22) are time-variant, which becomes one main barrier for system modelling and analysis. Thus, eliminating the effect of timevariance is very necessary to system modelling and analysis. To achieve the above goals, an observer-pattern model will be proposed as follows.
Supposing that Equation (22) is the averaged model of the double-input buck/buck-boost inverter under the rOi reference frame. Let x r be one arbitrary physical state variable in loadstage subsystem, then we will reconstruct one imaginary variable x i , which maintains 90 degree phase shift with reference to x r . According to Park orthogonal transformation, we obtain where Substituting Equation (23) into the fifth to seventh equations in Equation (22), then we obtain Pre-multiplying both sides of Equation (24) by T , Equation (24) can be expressed as ] can be rewritten as where For simplicity, Equation (26) can be rewritten as In actual cases, the reactive power of this cascade system is almost insignificant for the system analysis. Accordingly, the input instantaneous active power of the load-stage subsystem P (t ) is . Based on the instantaneous output power, the 'hidden' second-harmonic in the source stage can be exposed by the following model of the source-stage converter: Clearly, due to the time-variant term cos(2 t ), Equation (28) is still a non-autonomous periodic system. Since cos(2 t ) is one specific solution of the following differential equations: where the initial conditions are g 1 (0) = 1 and g 2 (0) = 0.
Based on Equations (27), (28) and (29), we can obtain the observer-pattern model of the double-input buck/buck-boost inverter as follows: For the sake of convenience,

NUMERICAL AND THEORETICAL ANALYSIS
In what follows, circuit simulations will be performed to present the multi-harmonic oscillation behaviours in the double-input buck/buck-boost inverter. Furthermore, the occurrence mechanism of multi-harmonic oscillation will be revealed in detail by using bifurcation analysis method.

Multi-harmonic oscillation phenomena
In actual application, the control parameters always have a significant impact on the system performance. Thus, the multiharmonic oscillation will be investigated with the variation of the integral time constant T 3 in the double-input inverter and other parameters values are listed in Table 1.
According to the circuit principle of double-input buck/buck-boost inverter, we can obtain the fact that the source-stage double-input buck/buck-boost DC-DC converter includes switching frequency of the sawtooth wave v ramp , whereas the load-stage H-bridge inverter contains the switching frequency of the triangular wave v tri and the line frequency of the reference voltage v ref2 . Note that the source-stage state variables contain the second-harmonic frequency, which is twice of the line frequency. Therefore, it can be concluded that there are switching frequency component, line frequency component and second-harmonic frequency component in the double-input buck/buck-boost inverter when the system is stable. Figure 2 shows the time waveforms and FFT of the variables v o1 and v o2 when T 3 = 0.0001. As illustrated in Figure 2(a) and (c), the output voltage v o1 is sinusoidal waves with the DC offset, whereas the output voltage v o2 is normal sinusoidal waves without the DC offset. In addition, the FFT of v o1 displays that the hidden second-harmonic frequency is 100 Hz in Figure 2(b), which is twice of the line frequency f l in Figure 2(d). Obviously, the double-input buck/buck-boost converter system operates stably when T 3 = 0.0001.
As T 3 decreases to 0.000083, the appearance of slow-scale unstable behaviour is observed in Figure 3. Different from the sinusoidal wave in Figure 2(a) and (c), an obvious oscillation happens in the waveforms shown in Figure 3(a) and (c). In addition, we can find that in Figure 3(b), three new harmonics, that is, f a1 = 536 Hz, f a2 = 436 Hz and f a3 = 336 Hz occur in the FFT of v o1 , which means that a multi-harmonic oscillation happens. Likewise, another three new harmonics, that is, f b1 = 586 Hz, f b2 = 486 Hz and f b3 = 386 Hz appear in the FFT of v o2 , as shown in Figure 3(d). In brief, this unstable phenomenon is referred to as multi-harmonic oscillation. Surprisingly, comparing Figure 3(b) and (d), it can be found that the difference between the three harmonics both in the sourcestage converter and in the load-stage inverter is just 100 Hz, that is, the second-harmonic frequency 2 f l , whereas the difference between the three harmonics in the source-stage converter and ones in the load-stage inverter is just 50 Hz, that is, the fundamental (line) frequency f l , respectively. It indicates that the multi-harmonic oscillation is close related to the original frequency components including the second-harmonic frequency and fundamental frequency.

Bifurcation analysis
From the above results, it is essential to uncover the underlying mechanism of the multi-harmonic oscillation in the doubleinput inverter. Bifurcation analysis is a powerful tool to capture the non-linear dynamics and assess the harmonic stability. In the following, bifurcation analysis will be applied to analyse the multi-harmonic oscillation observed above. According to F (X ) = 0, the equilibrium point can be calculated as , v e con2d , v e con2q , g e 1 , g e 2 ] T . The Jacobian matrix of the proposed observer-pattern model in Equation (31) can be written as follows: To find the eigenvalues , the following polynomial equation will be solved: Table 2 shows the eigenvalues with the variation of T 3 , it is obvious that the real parts of 11,12 are always zero. This reason is that 11,12 are originated from the constructed differential equations dg 1 dt = g 2 and dg 2 dt = −4 2 g 1 . Therefore, 11,12 are independent of the circuit parameters and have no effect on the harmonic instability. Note that the imaginary parts of 11,12 unmask the second-harmonic component 2 f l . In addition, when T 3 begins to decrease from 0.00098, the real parts of two pairs of eigenvalues 1,2 gradually increase together with 3,4 . In this case, 5,6 , 7,8 and 9,10 still stay on the left of complex plane. Figure 4 depicts the traces of 1,2 and 3,4 . It is evident that as T 3 decreases to 0.00009, 1,2 and 3,4 move through the imaginary axis simultaneously. This implies that a Hopf bifurcation behaviour happens in the double-input inverter and the double-adding harmonic oscillation with two different harmonics occur in the system. Furthermore, by the imaginary parts of 1,2 and 3,4 , the oscillation frequencies will be calculated as 537 and 437 Hz. Namely, the two dominant harmonics come to birth, which are consistent with the two harmonics f a1 and f a2 in Figure 3(b) except f a3 . Now the following question may be posed: How does the harmonic f a3 appear on earth? In fact, since the second-harmonic component 2 f l exists naturally in source-stage converter, the interaction between the secondharmonic component 2 f l and the above new harmonic f a2 should be responsible for the harmonic f a3 . Similarly, due to the interaction between three harmonics and the fundamental frequency, f b1 = f a1 + f l , f b2 = f a2 + f l and f b3 = f a3 + f l appear in load-stage inverter, as shown in Figure 3(d). Thus, the harmonics f a1 and f a2 are the dominant harmonics amongst these harmonics. Figure 5 demonstrates the phase portraits before and after the Hopf bifurcation, which inherently exposes the underlying mechanism of the multi-harmonic oscillation in the doubleinput inverter.

Loss analysis
In what follows, the loss analysis will be used to assess the influence of the multi-harmonic oscillation on the system efficiency. Here the calculation of power losses is performed in a typical switch and the calculation method has been conducted as in [30]. The losses of the dual-input inverter can be divided into switching losses and conduction losses. The switching losses include the turn-on losses and the turn-off losses of the switches. Note that the switching losses of the diodes can be neglected because they are significantly lower than the other types of losses [31]. The turn-on energy loss of the switch S k in the nth switching cycle can be obtained where V n s k , I n s k ,on and t s k ,on are the off-state voltage on the switch S k , the current through S k after turning on S k and the turn-on transition time, respectively.
Likewise, the turn-off energy loss of S k in the nth switching cycle is where I n s k ,off and t s k ,off are the current through S k before turning off S k and the turn-off transition time, respectively.
The conduction losses mostly result from the on-state resistance of the switches (r s k ,on ), ESR of the inductors (r L k ) and the on-state voltage of the diodes (V n on,D k ). The conduction losses of the switch S k , the diode D k and the inductor L k can be calculated according to the following equations: where d n s k is the duty cycle of S k , T s is the switching cycle,I n D k is diode average current, and I n s k ,rms and I n L k ,rms are the root-mean-square (RMS) of the current through S k and L k , respectively. According to Equation (37), the total losses are P loss = P source,loss + P load,loss Then, the efficiency can be obtained by where p out is the output power. Figure 6 shows the system efficiency with the integral time constant T 3 varying under three switching frequencies f s . Obviously, the increase of f s will improve the efficiency. However, it can be seen that the efficiency will drop sharply under the three different f s when T 3 decreases to 0.00008. This indicates that no matter how the switching frequency is selected, the occurrence of these multi-harmonic oscillations causes the increase of the circuit loss, that is, the decrease of the system efficiency.

STABILITY ENHANCEMENT ANALYSIS
It is obviously that the variations of some parameters have a great influence on the system stability and other performance. Thus, it is necessary to evaluate the influence of the circuit parameters on the system stability, which can help power electronics engineers select proper circuit parameters so as to realise the stability enhancement. In what follows, the harmonic-stability sensitivity is proposed to reveal the relationships between the system stability and some parameters. On this basis, the key circuit parameters which have a major effect on the multi-harmonic oscillation are identified.
Since the real parts of the eigenvalues can show the system stability region, the harmonic-stability sensitivity is defined as where is the chosen circuit parameter and Re( i ) is the real part of the eigenvalue i . Re( i ) indicates Re( i )'s partial derivative with respect to . Figure 7 shows the harmonic-stability sensitivities of the double-input buck/buck-boost inverter. Since the eigenvalues 1,2 and 3,4 govern the dominant harmonics f a1 and f a2 respectively, it is important to identify the key parameters for 1,2 and 3,4 . We can see that for 1,2 and 3,4 , the parameters C 2 , R, K 3 and T 3 have the rather larger harmonic-stability sensitivities than other parameters. This implies that C 2 , R, K 3 and T 3 do the dominant effect on the multi-harmonic oscillation. Moreover, if K 3 , C 2 and R increase or T 3 decreases, 1,2 and 3,4 will move toward the imaginary axis, which indicates that the system will suffer from multi-harmonic oscillations more easily. Interestingly, 9,10 are highly sensitive to K 3 and T 3 , too. When K 3 increases or T 3 decreases, 9,10 will move toward left in the complex plane. However, for K 3 and T 3 , the harmonicstability sensitivity to 9,10 are much smaller than the ones to 1,2 and 3,4 . Thus, as K 3 increases or T 3 decreases, the system tends to lose stability more readily. Additionally, 5,6 and 7,8 are more sensitive to L 1 . The increase of L 1 will make 5,6 and 7,8 move toward the imaginary axis. Based on the above analysis, it can be found that the power-stage parameters C 2 , R and control parameters K 3 , T 3 are the key parameters for the multiharmonic oscillation in the double-input buck/buck-boost inverter.
Due to the intermittency of renewable energy, the input voltage may vary within a wide range. From Figure 7, the effect of the input voltage on the system performance can be assessed. It can be observed that 1,2 and 3,4 show much smaller harmonic-stability sensitivities to the input voltages E 1 and E 2 than other parameters. This indicates that E 1 and E 2 have negligible role on the multi-harmonic oscillation. For other eigenvalues, E 1 and E 2 is related with 7,8 and 5,6 , respectively. Figure 8 shows the traces of the eigenvalues with the change of E 1 and E 2 . As E 1 increases, the pair of the complex conjugate eigenvalues 7,8 move toward left in the complex plane. When E 2 increases, 5,6 leaves away from imaginary axis but 7,8 move toward the imaginary axis slowly. It is evident that within the change ranges of the two input voltage, the real parts of all the eigenvalues are negative, which means that the system is still stable.
For providing more design-oriented information for electrical engineers, some stability boundaries are presented to assess that in which way the system would suffer from multi-harmonic oscillation. Figure 9 illustrates the stability boundary under the key parameters T 3, C 2 and R. It can be found that the unstable region becomes larger with the load resistor R and capacitor C 2 increasing. It means that the system will be far away from the stability boundary when the load resistor R and capacitor C 2 increases. However, when T 3 increases, the buck/buckboost inverter is prone to reach the stable region, which is quite beneficial for system stability enhancement. Obviously, the trends of the boundaries in Figure 9 are in agreement with those harmonic-stability sensitivity analysis.

EXPERIMENTAL RESULTS
In order to verify the above numerical and theoretical results, an experimental prototype is built with the specifications listed in Table 1. The photograph of the prototype is depicted in Figure 10. In this experiment, all diodes are MUR860. The MOS-FET of the front-end converter is 2SK3569 and the H-bridge in the downstream inverter is PM50CLA060. TMS320F28335 is chosen as the digital controller .The output voltages v o1 and v o2 in the source stage and the load stage will be transmitted to the ADC module of the TMS320F28335 through the voltage conditioning circuit, which consists of the voltage sensor LV25-P and the operational amplifier OP07CP. The inductor current in the source stage will also be sent to the ADC module through the current conditioning circuit, which contains the current sensor HKY25. When T 3 is 0.00012, the time-domain waveforms and FFT results of v o1 and v o2 are shown in Figure 11. It can be seen that the output voltage v o1 is sinusoidal waves with the DC offset whereas the output voltage v o2 is normal sinusoidal waves. In Figure 11, the fundamental component f l = 50 Hz and secondharmonic component 2 f l = 100 Hz appear in the frequency spectrums of v o1 and v o2 , respectively. This indicates that the When T 3 decreases to 0.000083, the experimental waveforms of v o1 and v o2 are shown in Figure 12. Obviously, the multiharmonic oscillation phenomenon takes place in v o1 and v o2 . From the frequency spectrums, the harmonics f A1 , f A2 and f A3 in v o1 are mainly around 510, 410 and 310 Hz and the harmonics f B1 , f B2 and f B3 in v o2 are about 560, 460 and 360 Hz, respectively. These measured results and the mathematical relationship amongst these harmonics are in agreement with those results in Figure 3.
Keeping T 3 as 0.000083 and decreasing K 3 to 0.02, the experimental waveforms shown in Figure 13 are different from the ones in Figure 12. As mentioned in Section 5, when K 3 decreases, the system approaches towards the stability operation. It follows in Figure 13 that the magnitudes of harmonics are much small except f l and 2 f l , and they are negligible. This indicates that a smaller K 3 is beneficial to suppressing the new harmonics. Thus, the influence of the parameters on the system stability can be effectively predicted by those above analysis in Section 5.
When the load resistor R increases to 12 Ω and other parameters are the same as the ones in Figure 13, the multi-harmonic oscillation appears, as shown in Figure 14. It should be noted that the system will suffer from multi-harmonic oscillation more easily with the increase of R. It is in accordance with the stability boundaries shown in Figure 9.