Modelling and stability assessment of the MMC‐HVDC energy interconnected system with the cyber delay of communication network

National Natural Science Foundation of China, Grant/Award Number: U1766210; National Natural Science Foundation of China, Grant/Award Number: 51377117; National Natural Science Foundation of China, Grant/Award Number: 51625702 Abstract Modular multilevel converters (MMCs) are favoured in high voltage direct current (HVDC) systems because of huge energy transmission capacity and few harmonics. However, the complicated control and the intricate communication network introduce cyber delay that affects the system dynamics. The communication structure is discussed and delay components are studied. The authors derive the MMC‐HVDC system model considering the cyber delay effect. Padé approximation is utilised to transform the transcendental characteristic equation of the delay model into an algebraic equation so that the dynamics can be observed through the eigenvalues. The comparison between the established model and the simulation in Simulink validates the proposed model's accuracy. The exact delay margin and the root cause of losing the stability are also discussed. The results show that the model is able to obtain the delay margin with a tiny error caused by the omitted high‐order harmonics and the cyber delay destabilises the system by lowering the controller performances.

modulate with the reference signals. The communication network only delivers the reference values, while the centralised modulation requires all switching signals shared in the network.
Previous studies on modelling the MMC-HVDC system [1,[14][15][16][17][18] have built some practical models. For example, [1,14] considered both the converters' internal dynamics and the control dynamics, and then derived the small-signal model. The authors of [15,16] derived the bilinear MMC model. With the proposed bilinear model, [15] predicted the behaviour of the MMC states one step ahead, and [16] designed a controller controlling both the external and internal dynamics. However, the models mentioned above didn't consider cyber delay. Reference [19] investigated the integration process of the electric vehicles and electrical grid, which revealed the distorted stability space and triggering mechanism of instability during the interaction. The authors of [20] proposed an improved small-signal modelling approach based on the characteristic equation of the PLL-based converter-dominated ac microgrids. The authors of [21] took the computational and switching delays into account and obtained the VSC-HVDC impedance model, then applied the passivity theory and the Nyquist stability criterion to study the DC side stability issue. The authors of [22] studied the generation mechanism of the delay and deduced the transfer function model of the MMC system including the delay term. With the obtained transfer function, [22] deduced the delay's effect and came up with the optimal design methods of the system controller parameters. Reference [4] pointed out that the delay is related to the sampling frequency and the switching frequency. Furthermore, to find out the proper frequency range, [4] built a kind of MMC delay model in the state space, analysed the model by using Padé approximation and achieved the quantitative stability assessments for MMCs. However, both [22] and [4] only considered a simple MMC model and only tested the model in the singlephase situation.

| Contribution and paper organisation
To address all kinds of problems mentioned above, the authors take both the control and converter dynamics and obtain the MMC-HVDC system state-space model. Based on the communication structure's discussion, cyber delay components are modelled. Integrating the cyber delay with the state-space model, the delay model is derived and then analysed by Padé approximation and its eigenvalues. The proposed model is able to reflect system dynamics under cyber delay effect compared to the models in [1,[14][15][16][17][18], and is more suitable for threephase, HVDC situation than the models in [4] and [22]. The proposed model is validated by the simulation in Simulink and the results show that the delay margin can be estimated through this model without time-consuming simulation.
The authors propose the MMC-HVDC system model that can reflect the dynamics under cyber delay's effect accurately and help to analyse the system. The communication network issue is discussed and EtherCAT in the ring topology is proved to be a wise choice. Delay's interaction with the converters is analysed and verification shows that the omitted high-order harmonics (≥3) can bring in a tiny stability margin error of no more than 3%. The results also show that cyber delay destabilises the system by lowering the control performances and unsteady mode, and in our case has an angular frequency of 237 rad/s approximately.
The remainder of the e authors' work is organised as follows. In Section 2, the whole system's communication structure is discussed and three components of the cyber delay are introduced and calculated. On account of the harmonic dynamics and the controllers' dynamics, Section 3 establishes the state-space model with varied harmonics and the delay model is derived by integrating the cyber delay with the state-space model. After Padé approximation, the simplified delay model is obtained. With the help of the Simulink simulation, proposed models are verified and some issues like cyber delay influence, delay margin and destabilization mechanism are studied in Section 4. Finally, section 5 draws the conclusions.

| Communication structure
To control the whole MMC-HVDC system, information sharing among the central controllers and all SMs is necessary. This makes the communication structure complicated and becomes a challenge for the control hardware [6]. If linking the central controllers to all SMs directly, that is, adopting the star topology (see Figure 1), synchronisation may be guaranteed and the communication speed is high but it will be impracticable and cumbersome to assemble because there will be a lot of long wires (interfaces, cables and fibre optics) that connect the central controllers to some SM [6]. Besides, due to the inherent attribute of the star topology, the corresponding communication hub of the central controller may be required to possess N or 2N communication ports. The authors of [23] proposed a scheme achieving communication using CAN protocol. This scheme makes the master module (the central controllers) and the slave modules (the submodules) to only communicate through the CAN bus, which is more like the tree topology. This scheme is low cost and communication ports requirement is less compared to the star topology. However, it can't satisfy the requirement to control a high number of SMs like MMC-HVDC [13,21]. For the sake of construction and engineering requirements, ring topology [6][7][8][9][10] is applied a lot as a kind of topology. It has a slower speed compared to the star topology but requires fewer wires and fewer ports, and it can function well even though there exist plenty of nodes. Among all the communication protocols fit for ring topology, EtherCAT is favoured because of its high speed [7], low synchronization jitter [8] and fault tolerance [9].
EtherCAT is an open-source protocol based on Ethernet which allows full-duplex communication and uses the classical master/slave configuration [24]. It has a high speed of 100 Mbps and is appropriate for almost all kinds of topology including star, tree, ring, etc. The distributed clocks CAO ET AL.
-87 synchronisation mechanism guarantees all of the nodes' clocks are almost synchronous with little error [8], which ensures the structure stability. In the communication ring, when some slave node or its cable is broken, the master node is still able to communicate with all other nodes by means of full-duplex communication (see Figure 2), which shows its high reliability and robustness.
The working procedure of the control system with the EtherCAT protocol is shown in Figure 3. Necessary signals in the AC system like voltages or currents are sampled into the central controllers. Then the central controllers calculate the voltage reference U ref based on the sampled data and reference values. [2] realised the voltage balancing in the distributed modulation with the designed proportional controller, which means it's feasible and be able to lighten the central controllers' burden to implement modulation in every sub-module instead of the central controllers. Thus the central controllers send the voltage references directly through the EtherCAT protocol and let LCUs modulate. By sending the references, the error will be much smaller than sending the on/off signals, which is better than the centralised control [2]. Besides, the distributed control also lightens the computation burden and reduces the information shared inside the system compared to the centralised control [13]. Figure 4 shows the EtherCAT frame signal transmission procedure. The central controllers send an EtherCAT frame signal (see Figure 5) containing voltage reference U ref for every SM through a ring whose nodes are the central controllers (master) and SMs' LCUs (slaves), which is shown in Figure 6. The LCU receives the corresponding segment, then updates the corresponding data segment with corresponding SM's information and resends this updated EtherCAT frame signal to the next LCU until it comes back to the central controllers [10].
The LCU finishes modulation and gets the switching signal 0/1, then transmits it to its corresponding SM. Thus the control of the MMC-HVDC system is accomplished through our communication structure. However, as can be seen in Figure 3, cyber delay exists in many stages in the system [4,6,20], which will affect the system's normal operation. So it's obligatory to study the cyber delay's character and its effect on the system.

| Cyber delay in the MMC-HVDC system
According to Figure 3, three kinds of dominating cyber delays are considered: cyber delay of sampling and data processing τ 1 , Usually sampling itself and controllers' calculation don't cost much time, but U ref won't be updated until one sampling period T s goes by. So it is assumed that where f s is the sampling frequency.

| Cyber delay of signal transmission τ 2
In the ring topology, it can be calculated by [6] where D is the data payload in bytes, T byte is the time to transmit one byte, the forwarding delay of every SM is T f w ∈ ½0:7μs; 1μs�. From Figure 5, it is known that there are 8 bytes data and 1 byte checksum signal for every LCU. Considering the frame signal has 1 byte control signal in the beginning and 1 byte checksum signal in the end, and one arm has N SMs, the following equation is obtained: The right side's second term of (1) is introduced by the data packaging process, and it's easy to verify that it's negligible compared to the right side's second term of (1). Considering the bandwidth of EtherCAT is 100 Mbps, The components of τ 2 can be interpreted as a data transmission process (DT byte ) with the addition of all nodes' forwarding delivery (NT f w ) as can be seen in Figure 7 (the right side's second term of (1) is omitted).

| Cyber delay of modulation τ 3
Taking Phase-Shifted Pulse Width Modulation (PS-PWM) for example, [4] claimed that delay caused by modulation was half a switching period T sw , which means where f sw is the switching frequency.

| Summarised cyber delay τ
Considering several kinds of cyber delays above, it is possible to summarise them and treat the sum of them as the lumped cyber delay τ ¼ τ 1 þ τ 2 þ τ 3 . In order to explain the relationship between cyber delay and its influence factors visually, it is assumed that f s ¼ f sw ¼ f ,and τ ¼ τðf ; NÞ surface graph is made in Figure 8. It can be seen that the more f is and the less N is, the less cyber delay is. Furthermore, cyber delay changes rapidly with changing f when f is quite small, which should be noticed.

| MMC-HVDC SYSTEM DELAY MODEL
Here, the state-space model considering both the dynamics of the control system and converters will be derived first. Then on the basis of the state-space model, the cyber delay is integrated into the model and Padé approximation is used. Finally, the delay model which is convenient to analyse mathematically is obtained.

| State-space model
A typical MMC-HVDC system is shown in Figure 9 Taking phase A for example, the capacitor voltage dynamics of the upper arm (u cp ) and lower arm (u cl ) can be decided by where C is the SM capacitance, i pa , and i la are currents of upper and lower arm, respectively, and S pa and S la are insertion indexes of the upper and lower arm respectively, which can be expressed as where U ca is the output of the PQ control and U cira is the output of the circulating current suppressing controller (CCSC) which will be introduced simply below.
Considering the upper arm for example and ignoring higher-order harmonics, the capacitor voltage and the arm current can be expressed as where � u c , u 1 c , u 2 c , u 3 c are 0, 1, 2, 3 order components of the capacitor voltage in averaging situation and i cira is the circulating current of phase A.
Speaking of the control dynamics, current control (PQ control), CCSC and phase-locked-loop (PLL) are considered. In the MMC-HVDC system, PQ control means regulating the transmitted active power P and reactive power Q with 2 double-loop PI (Proportional Integral) controllers, each one consists of 2 PI units.
where L eq ¼ L tr þ L arm 2 . Circulating current inside the converters can result in extra loss, so CCSC uses 2 PI controllers to depress the negative effect that circulating current causes.
where i cirdref and i cirqref usually equal to zero. From the above equations, it is known that PQ control and CCSC are designed based on d-axis and q-axis signals. So PLL is needed to track the phase of the AC system so that the original AC signals are transformed into the d-axis and q-axis signals.
. ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Equations (9)-(22) make up the control system model which indicates the system's control dynamics.
Besides, the AC system dynamics should also be taken into consideration. Taking advantage of the Thevenin theorem, the AC system is simplified as an ideal voltage source with resistance R ac and inductance L ac in series.
where U ed , and U eq are d-axis and q-axis component of the ideal voltage source value, respectively, while U AC is the amplitude of the AC side's phase voltage. By integrating all the dynamics mentioned above and algebraic simplification, the linearized state-space model is obtained as stated by  [14] deduced a very similar state-space model in detail and as the authors concentrate on the following delay model, detailed algebraic simplification is no longer exhibited. Instead, detailed algebraic simplification of the delay model is given in the appendix.

| Delay model and Padé approximation
The existing cyber delay between controllers and circuits results in a phenomenon such that the circuits' control command at time t is actually calculated based on the data at the time t − τ. So, cyber delay is supposed to be integrated into the model by replacing all the sampled data gðtÞ in the control dynamics (9)-(18) with gðt − τÞ(g ¼ i sd ; i sq ; i cird ; i cirq ; U sd ; U sq ). Thus Equation (27) Theoretically speaking, our delay model, Equation (28) is able to indicate the system's dynamics under the cyber delay effect. By solving the characteristic equation of (28), which is The roots are obtained and it can be said that the system is asymptotically stable if and only if all the roots have the negative real part [25]. However, delay introduces infinite eigenvalues [4,5], making dynamics hard to be analysed. Hence, we apply Padé approximation to simplify the delay's s-domain expression e −τs as a rational fraction. As the transcendental term e −τs is missing, the transcendental characteristic Equation (29) can be turned into an algebraic equation. The authors make the numerator l remain zero-order and choose the order of the denominator k with 2 criteria in [4] so that the derived equations are as simple as possible. According to [26,27], the following equations can be derived: Taking k = 2 for example and if we define then we can get 2g τ þ 2τ _ g τ þ τ 2 € g τ ¼ 2g. Again, if we define g z ¼ _ g τ , then a pair of differential equations is obtained: After replacing gðtÞ with g τ and combining (27) and (32) together, and after a series of complicated algebraic simplification (see appendix), the simplified delay model is obtained: The 32 eigenvalues are obtained and the system stability can be assessed by observing the eigenvalues' location. Thus it is possible to analyse the system's dynamics and stability easily (the values of k and model's order are only used for explanation, the appropriate order will be chosen in the next section).

| Delay model verification
Based on the delay model (28), the time-domain simulation M code of the single-terminal MMC-HVDC system is built. The single-terminal MMC-HVDC system model is also built in Simulink with the purpose of validating the model (28). The system parameters are shown in Table 1.U AC , N, and U dc are determined from Zhouding station of Zhoushan multiterminal VSC-HVDC transmission project [28]. The authors consider 3 cases in this subsection and the following subsection 4.2:(a) ideal situation, which means the sampling frequency and the switching frequency can be as high as possible regardless of loss. For example, f s ¼ f sw ¼ 3 � 10 4 Hz, then τ ¼ 1 the power reference P ref is increased from 400 to 420 MW.
The comparison results are shown in Figures 10--12, which demonstrates the delay in model's correctness. From Figures 10-12, it can be observed that i sd ,i cird and u cp start to change at t ¼ 1:8s; because the equilibrium point changes with increasing P ref . When τ ¼ 400μs or τ ¼ 800μs; i sd ,i cird and u cp change to a new steady value after a period of oscillation, which means the system is stable at the new equilibrium point. However, they fail to converge and keep oscillating with an increasing oscillating amplitude when τ ¼ 1000μs; which means the system loses stability facing the cyber delay in the communication network.
From Figures 10-12, Two conclusions are drawn and they are:(1) The system is stable when τ ¼ 400μs or τ ¼ 800μs: But the system responses are slower when τ ¼ 800μs; which  indicates that delay deteriorates system dynamics.
(2) Given that the delay reduction may require a higher switching frequency which introduces more loss, it's important to achieve the trade-off between dynamics and loss by designing the frequency and communication network appropriately. A suitable communication structure and a high-speed protocol like EtherCAT can lighten the pressure of choosing the switching frequency by making communication delay as small as possible.

| Simplified delay model verification
In this subsection, the simplified delay model is also verified by the method applied above in Section 4.1. The time-domain simulation M code is built based on (33). Usually, Padé approximation is verified by the critical eigenvalues. It can be concluded that an n-order Padé approximation is appropriate if its critical eigenvalues are almost equal to the critical eigenvalues of the corresponding (n+1)-order Padé approximation when the delay is a step away from the steady state and a step away from the unsteady state [4]. After the verification, it is found out that 870μs is the delay margin when using third-order Padé approximation. If 870μs is regarded as 'a step away from the steady state' and 880μs as 'a step away from the unsteady state', the eigenvalues of all orders Padé approximation can be analysed. The result is shown in Table 2. It is observed that the differences between k = 3 and k = 4 of both the 'steady state eigenvalues' and the 'unsteady state eigenvalues' are smaller than 0.02, which indicates that k = 3 is appropriate. Besides, this Padé approximation is also validated with the help of the simulations. The comparison results are shown in Figures 13-15.
From Figures 13 to 15, it can be seen that compared to the Simulink simulation, the simplified delay model is also able to describe the system dynamics. So it is possible to analyse the system dynamics taking advantage of the simplified delay model's eigenvalues, which will be shown in the following Section 4.3.

| Eigenvalue analysis
As demonstrated in Section 4.2, the simplified delay model is able to describe the system dynamics, so we study its Jacobian matrix A ″ and extract critical information from its eigenvalues. Figure 16a shows partial root locus of A ″ when τ increases from 400μs to 1000μs, it can be seen that some eigenvalues pass through the imaginary axis from the left side during the increment of τ. Figure 16b is the waveform of q-axis AC current. Before t ¼ 1:8s; τ ¼ 400μs and the system is stable, at t ¼ 1:8s; τ increases to 1000μs and the system loses stability.
Both of them indicate the delay margin is some value in ½400; 1000� μs. To find out the exact margin, eigenvalues are plotted in steps of every, and simulation is also tested and the results are shown in Figures 17-19.
In Figure 17, at τ ¼ 870 μs, the critical eigenvalues are -0:1361 � j237:6(blue dots in the blue circle) whose real part is negative but the critical eigenvalues pass through the imaginary axis and become 0:4831 � j237:3(red dots in the red circle) when τ increases to 880 μs:This means delay margin is 870μs and the oscillation angular frequency when losing stability is 237 rad/s. Figure 18 shows the waveforms obtained in the simulation when suffering a disturbance at t ¼ 1:8s:We can see that the system is able to oscillate back to the original steady state when τ ¼ 840μs from Figure 18a. the system loses stability when τ ¼ 850μs and τ ¼ 860μs, respectively. In summary, the delay margin in simulation can be regarded as 850μs:From Figure 19, it can be seen that the system takes 0.272 s to oscillate for 10 periods, which means the oscillation angular frequency when losing stability is ð2π � 10Þ=0:272 ≈ 231 rad=s: Note that the obtained error of 20μs between the model and the simulation may be caused by the high-order harmonics of converters' current and voltage. The simulation's circuit includes all harmonics caused by switches' turning on or off but this model only considers dominant harmonics. Current affects voltage's dynamics through capacitors' charging and discharging and voltage affects current's dynamics through inductors in return. Thus, all order harmonics of the converters' current and voltage are coupled tightly. A certain order harmonic disturbance may affect other orders' harmonics. According to the proposed delay model, delay exerts an influence on the converters via the insertion indexes in (5) (6). The insertion indexes containing delay is multiplied by the current in (8), and thus the voltage dynamics are relative with the delay. Also, the current dynamics are dependent on the voltage, so the current is also related with the delay. The more the harmonic components are considered, the more the terms will be in (5) (6) (7) (8), and the bigger influence the delay is able to exert.
To demonstrate the idea stated above, the reduced-order model is obtained by setting u 3 c ¼ 0 compulsively and the margins calculated by Padé approximation are shown in Table 3. The results indicate the margin increases With the higher-order voltage components, the delay inside the insertion indexes interacts with more terms and makes a bigger impact, tends the system to lose stability. The higher-order harmonics are smaller and coupled less tight, which means delay's interaction with harmonics will not be increased highly. So, it is rational to assume that the higher order harmonics will be able to bring 20μs margin decrement. Furthermore, this case also demonstrates that in such a strong-coupled system, suppressing unnecessary harmonics can increase the delay margin, which may be a method of improving the system's stability. Besides, it is also worth noting that the pair of critical eigenvalues seems not to be the introduced eigenvalues. By decreasing delay to an extremely tiny value (20μs), it can be seen that the pair of critical eigenvalues become -38:69 � j251:5(red dots in the red circle of Figure 20), which are approximately equal to a pair of original eigenvalues of the traditional state-space model (τ ¼ 0) -39:16 � j252:2(blue dots in the red circle of Figure 20). This phenomenon indicates that the root cause for the system to lose its stability is delay, which lowers the controllers' performance and induces the original eigenvalues to move rightwards, rather than introducing eigenvalues located at the right side of the imaginary axis directly.
In summary, this section validates the accuracy of the proposed delay model and the simplified delay model by the comparison with Simulink simulation. The delay margin, the unstable mode's frequency and the instability reason are studied, through the Padé approximation. The error of the proposed model is explained and verified by the reduced-order model test. The more specific conclusions will be presented in the next section.

| CONCLUSION
The communication structure of the MMC-HVDC system is discussed and the inherent control delay and communication delay are modelled with the help of the EtherCAT protocol. Integrated with the modelled delay, the MMC-HVDC system delay model is derived in the state space. Padé approximation is used to simplify the model and the system dynamics are analysed through the simplified model's critical eigenvalues behaviours. The