Preprint
Article

Network-Independent Grid Synchronous Stability Boundary and Spontaneous Synchronization

Altmetrics

Downloads

266

Views

92

Comments

0

Yu Yuan  *

This version is not peer-reviewed

Submitted:

11 April 2024

Posted:

11 April 2024

You are already at the latest version

Alerts
Abstract
The synchronization of complex networks is believed to be closely related to network topology. How to study synchronization of real network without clear topology and parameters? In this Letter, a synchronized stability boundary equation is derived that is system- and disturbance-independent and applicable to arbitrarily coupled grids. The results imply that spontaneous synchronization in a network may be network independent. These conclusions provide new research paths for network synchronization and analyze the synchronization stability of grids in a unified way.
Keywords: 
Subject: Physical Sciences  -   Other

Introduction

The study of synchronization began with Huygens. With the increase in the number of studies on collective behavior in complex systems, spontaneous synchronization in coupled systems has attracted public attention [1]. Synchronization in complex systems is widespread in the real world [1,2], e.g., the aggregation of flocks of birds and schools of fish, the flashing of fireflies, and the synchronization of generators.
Currently, it is widely believed that network topology is closely related to spontaneous synchronization [2]. To analyze these spontaneous synchronization behaviors, the interactions between individuals are first simplified into a coupled complex network, and then studied using the knowledge of network synchronization. However, this is actually difficult to achieve in practice [3,4]. Realistic interactions are often invisible, which may lead to incorrect network topologies. The complexity and nonlinearity of real networks can also make topological information incomplete [5]. Overly large networks make clear topology modeling very difficult. Therefore, the appropriate solution is to construct a network-independent synchronization analysis path that identifies collective synchronization and which individuals are not in this collective solely by the behavior of each individual in the system. In this study, this approach is experimentally verified on a realistic grid.
Synchronization is a prerequisite for the normal operation of a power grid [6]. Large power systems are complex coupled systems where uncertainties [7] and nonlinearities coexist. To analyze the stability of these systems, researchers have developed numerous insightful discriminatory methods [1,8,9,10,11,12], including finding stability region boundary [13,14,15] and describing the synchronization of generators using the spontaneous synchronization conditions of complex networks [2].
Finding the stability boundary of a system is an important fundamental problem [14,16,17]. A stable boundary is the union of unstable equilibrium points [18]. When an operating point is outside the boundary, the corresponding system is desynchronized [15]. A synchronous stability boundary, which is a core concept of grid stability, is closely related to many issues [19,20,21]. Therefore, studying synchronous stability boundary has a significant impact on the development of power systems. For decades, scientists and engineers have searched for an analytical equation to describe the boundary [13,14,22]. On the other hand, the spontaneous synchronization of complex systems has been used to explain the synchronized operation of generators in interconnected grids [23]. Therefore, many studies on power system stability are based on the knowledge of synchronization conditions on complex networks [2,23,24]. However, this scheme is sometimes considered to oversimplify real systems [9,23], while arguing that self-organization does not exist in the grid case [9].
Here, an equation is derived and visualized to describe the stability boundary of a power system in a unified way. This approach solves the problem of multi-swing stability discrimination of multiple generators in the grid. This equation proves that the synchronization stability boundary is independent of the network and disturbances. In addition, evidence that spontaneous synchronization occurs near the boundary and manifests itself as a specific spatiotemporal structure was found in the study. This is the first time evidence of self-organization has been found in a power system. This indicates that the synchronization conditions on a network are also independent of the network. In this context, for potential scenarios where it is difficult to construct network models, this research develops a new method for synchronization on complex networks. This approach can provide analysis for the synchronization of power grids in a unified way and be instructive for other disciplines.

Stability Boundary

Regarding the synchronization, there is one different assumption not previously considered [21]: the system requires the synchronous power Pu to maintain synchronous stability, and the synchronous power is provided by the coupling power PΔu within the system. That is, when PΔuPu, the system is synchronized, otherwise it is out of step (see supplementary material for details).
The amplitude (voltage in this case), whose dynamics have been neglected for simplicity [9,23], actually plays a decisive role in synchronization stability. Counterintuitively, there is no direct relationship to phase. The boundary equation is
u K = u L u L u K = 2 cos δ K L 1 , u K > u L 0 , δ K L 0 , π 3 .
uK, uL are the voltage per unit of the Kth and Lth meta-generator(see supplementary material for details), respectively. δKL = |δKδL| ≥ 0 and δK = 2arctan(ωK) is the angle of rotation rate of the Kth meta-generator. ωK is the rotor speed per unit of the Kth meta-generator, which is a normalization quantity commonly used in power systems studies. Note that δK is not a phase and L = K + 1. By the definition of synchronization, δ1 = ⋯ = δK = δL = ⋯ = δn.
Eq.(1) is the analytical equation for the synchronous stability boundary. Specifically, when δKL > 0, the set where |uL| = |uK| is the isolated stability domain [13]. This may explain why symmetric networks have better synchronization capabilities [23].
The current analysis for synchronous stability on a network requires the creation of a network coupling matrix and then the master stability equation [25]. This approach is mainly applicable to a case where the network topology is complete. In contrast to previous reports [2,13,14], Eq.(1) is independent of the network topology, system parameters, perturbations and the number of subsystems. This shows that the stability boundary is independent of these factors. This provides the basis for the applicability of Eq. (1) to different grids. Additionally, this case demonstrates that the synchronization of complex networks can be analyzed even without a coupling matrix. As shown in Figure 1 and Figure S4, these conclusions are proved by numerical experiments on different standard arithmetic models.
In Figure 1, a New England test system (10-gen) was used. A three-phase short-circuit ground fault occurred at node 18.
Figure 1(b) shows that the synchronization stability boundary can distinguish between synchronized and desynchronized states and provides specific information about synchronized and desynchronized individuals or groups, respectively. This result shows that the system is critically stable at Δt = 0.154s. For n-1 operating points of n generators, the system is considered out of global synchronization if any of the operating points (uK, uL, δKL) is outside the boundary and far from the other operating points (σΔt(δi) >> σΔtε(δi), 0 < ε < Δt, where Δt is the fault duration) [see Figure 1(b) and Figure 2(a)]. This result suggests that the partial synchronization phenomenon [26,27,28] in which two states coexist is caused by the loss of synchronization stability between these individuals or groups. The different synchronization patterns can be easily identified in Figure 1(b), where m operating points outside the boundary indicate that the n meta-generators are sequentially divided into m+1 coherent groups. This approach can be directly applied to discriminate coherent generators [29].
Preprints 103679 g002
The results shown in Figure 1(b), (c) fit well with the results of the numerical experiments [see Figure 1(d), (e)]. This shows the ability of the boundary to discriminate the synchronization stability of the grid, and demonstrates that Eq. (1) is valid.
Eq. (1) has a variant as follows:
δ K L c r = arccos 1 u K u L 2 u K .
When |uK|, |uL| are sufficiently close [32], the stability margin of the Kth and Lth meta-generator angle rate difference δ K L c r also tends to 0. This indicates that the system may already be in a critical state during normal operation. In this case, even if the difference in the values between δK and δL is small, it takes only a very small perturbation to make the system asynchronous [33,34]. As |uK| − |uL| increase, the stability margin δ K L c r increases. This explains why the system has better synchronization stability when it is highly asymmetric [6,35]. This may also indicate the need to enhance the low voltage ride-through capability of generators in the future [36].
Similar results are obtained for the 3-generator test system (3-gen) via the same steps (see Figure S5). The New England test system and 3-gen are two completely different network systems, but both apply the same boundary equation. This demonstrates that the synchronization stability boundary is independent of the network topology. As a result, Eq. (1) allows a unified analysis of the synchronization stability of the grid and provides a new understanding of synchronization.
Synchronization is a pervasive topic in other disciplines such as social networks, biological systems, brain neural networks and other physical systems [1]. In many cases, it is impractical to model a network with a clear topology. Here, a synchronization-stabilizing boundary equation for real networked systems, independent of network topology, is developed. This helps to open up new avenues for synchronization studies. For example, studying the response of operating points to network and parameter changes, and to compare the post-response results with the boundary to discriminate synchronous behaviour and identify important nodes. Boundary equations in other fields may have different forms.

Spontaneous Synchronisation

To study desynchronization process, the trajectories of the operating points near the critical point are investigated and found to exhibit a special spatiotemporal structure.
As shown in Figure 2(a), the decrease in σ(δi) from the highest point indicated that near the critical point, the system spontaneously directed the velocity of the subsystem to the mean value. This can be interpreted as a spontaneous synchronization effect due to system coupling [8,37]. As the stability boundary is approached, the system self-organizes toward synchronous evolution [23]. σ(δi) is discontinuous on the boundary, and in addition, spontaneous synchronization causes σ(δi) to decrease whether the running point crosses the potential barrier inwardly or outwardly[see Figure 2(a)]. Since the phase transition point [the magenta point in Figure 2.(a)] does not coincide with the critical point, the different starting points of spontaneous synchronization in these two directions may lead to hysteresis. This phenomenon is similar to the hysteresis phenomenon of explosive synchronization [38,39].
Self-organizing behavior emerges from the interactions of these meta-generators, and its effects are reflected in the perturbation trajectories at the operating point (see Figure 2). When the conditions are correct (as Δt increases, δL undergoes a second-order phase transition), the trajectories emerge with a special spatiotemporal structure. Before reaching the boundary, the operating points were in "decelerating motion", as if they were crossing a potential barrier in Figure 2(b). The potential barrier and cycle were the results of self-organizing in the δKδL and uKuL planes, respectively [see Figure 2(b), (c)]. Although spontaneous synchronization has been used to understand the synchronous operation of generators [23], this structure has not been reported. Due to the presence of long-range correlation effects near the critical point, the perturbed trajectories of all the meta-generators take on this structure at the same time (see Figure 3 for detailed). This suggests that spontaneous synchronization causes long-range correlations among all individuals in the system. Moreover, the phase transition point is not a critical point for desynchronization.
Spontaneous synchronization is considered the onset of network synchronization [1]. Currently, spontaneous synchronization is considered closely related to network topology [40]. However, the experimental results show that spontaneous synchronization occurs on the synchronization stability boundary. (see Figure 2 and Figure S7 for detailed). For coupled network systems, this strong correlation indicates that the location where spontaneous synchronization occurs is determined in uKuLδKL by Eq.(1). This may indicate that the mechanism of spontaneous synchronization is not necessarily related to the network topology, i.e., synchronization of the network may be independent of the network. This will challenge the traditional perception of synchronization in networks. Additionally, this finding also implies that there may be a close relationship between the critical stability of the dynamic system and the self-organized behavior.
The behavior of the operating point near the boundary is very complex. For example, it does not always result in the formation of a potential barrier [In some cases, δL undergoes a first-order phase transition in Figs. S2.(c) and S6.(c)]. The reasons for this difference, or rather, the specific conditions for the formation of a potential barrier require further research.

Supplementary Materials

The following supporting information can be downloaded at the website of this paper posted on Preprints.org.

Data availability

All the data that support the findings of this study are available at Figshare (DOI:10.6084/m9.figshare.23585961).

Conflicts of Interest

The author declares no conflicts of interest.

References

  1. A. Koronovskii, O. I. Moskalenko, and A. E. Hramov, Synchronization in Complex Networks, Tech. Phys. Lett. 38, 924 (2012).
  2. Dörfler, F.; Chertkov, M.; Bullo, F. Synchronization in complex oscillator networks and smart grids. Proc. Natl. Acad. Sci. 2013, 110, 2005–2010. [CrossRef]
  3. Wu, K.; Hao, X.; Liu, J.; Liu, P.; Shen, F. Online Reconstruction of Complex Networks From Streaming Data. IEEE Trans. Cybern. 2020, 52, 5136–5147. [CrossRef]
  4. L. L. Linyuan and T. Zhou, Link Prediction in Complex Networks: A Survey, Phys. A Stat. Mech. Its Appl. 390, 1150 (2011). [CrossRef]
  5. Xu, Y.; Zhou, W.; Fang, J. Topology identification of the modified complex dynamical network with non-delayed and delayed coupling. Nonlinear Dyn. 2011, 68, 195–205. [CrossRef]
  6. Molnar, F.; Nishikawa, T.; Motter, A.E. Asymmetry underlies stability in power grids. Nat. Commun. 2021, 12, 1–9. [CrossRef]
  7. Martínez, I.; Messina, A.R.; Vittal, V. Normal Form Analysis of Complex System Models: A Structure-Preserving Approach. IEEE Trans. Power Syst. 2007, 22, 1908–1915. [CrossRef]
  8. Zhu, L.; Hill, D.J. Synchronization of Kuramoto Oscillators: A Regional Stability Framework. IEEE Trans. Autom. Control. 2020, 65, 5070–5082. [CrossRef]
  9. Casals, M.R.; Bologna, S.; Bompard, E.F.; D', G.; Agostino, N.; Ellens, W.; Pagani, G.A.; Scala, A.; Verma, T. Knowing power grids and understanding complexity science. Int. J. Crit. Infrastructures 2015, 11. [CrossRef]
  10. Gurrala, G.; Dimitrovski, A.; Pannala, S.; Simunovic, S.; Starke, M. Parareal in Time for Fast Power System Dynamic Simulations. IEEE Trans. Power Syst. 2015, 31, 1820–1830. [CrossRef]
  11. Gurrala, G.; Dinesha, D.L.; Dimitrovski, A.; Sreekanth, P.; Simunovic, S.; Starke, M. Large Multi-Machine Power System Simulations Using Multi-Stage Adomian Decomposition. IEEE Trans. Power Syst. 2017, 32, 3594–3606. [CrossRef]
  12. Wang, B.; Fang, B.; Wang, Y.; Liu, H.; Liu, Y. Power System Transient Stability Assessment Based on Big Data and the Core Vector Machine. IEEE Trans. Smart Grid 2016, 7, 2561–2570. [CrossRef]
  13. Yu, Y.; Liu, Y.; Qin, C.; Yang, T. Theory and Method of Power System Integrated Security Region Irrelevant to Operation States: An Introduction. Engineering 2020, 6, 754–777. [CrossRef]
  14. Yang, P.; Liu, F.; Wei, W.; Wang, Z. Approaching the Transient Stability Boundary of a Power System: Theory and Applications. IEEE Trans. Autom. Sci. Eng. 2022, 20, 2268–2279. [CrossRef]
  15. Al-Ammar, E.A.; El-Kady, M.A. Application of Operating Security Regions in Power Systems, IEEE PES Transm. Distrib. Conf. Expo. Smart Solut. a Chang. World (2010).
  16. Kundur, P.; Paserba, J.; Ajjarapu, V.; Andersson, G.; Bose, A.; Canizares, C.; Hatziargyriou, N.; Hill, D.; Stankovic, A.; Taylor, C.; et al. Definition and Classification of Power System Stability IEEE/CIGRE Joint Task Force on Stability Terms and Definitions. IEEE Trans. Power Syst. 2004, 19, 1387–1401. [CrossRef]
  17. B. Student Member and G. A. Senior Member, On the Nature of Unstable Equilibrium Points in Power Systems, IEEE Trans. Power Syst. 8, 738 (1993). [CrossRef]
  18. Chiang, H.-D.; Wu, F.; Varaiya, P. A BCU method for direct analysis of power system transient stability. IEEE Trans. Power Syst. 1994, 9, 1194–1208. [CrossRef]
  19. Shubhanga, K.N.; Kulkarni, A.M. Application of Structure Preserving Energy Margin Sensitivity to Determime the Effectiveness of Shunt and Serles FACTS Devices. IEEE Power Eng. Rev. 2002, 22, 57–57. [CrossRef]
  20. Bhui, P.; Senroy, N. Real Time Prediction and Control of Transient Stability Using Transient Energy Function. IEEE Trans. Power Syst. 2016, 32, 1–1. [CrossRef]
  21. Al Marhoon, H.H.; Leevongwat, I.; Rastgoufard, P. A Fast Search Algorithm for Critical Clearing Time for Power Systems Transient Stability Analysis, 2014 Clemson Univ. Power Syst. Conf. PSC 2014 (2014).
  22. Rimorov, D.; Wang, X.; Kamwa, I.; Joos, G. An Approach to Constructing Analytical Energy Function for Synchronous Generator Models With Subtransient Dynamics. IEEE Trans. Power Syst. 2018, 33, 5958–5967. [CrossRef]
  23. Motter, A.E.; Myers, S.A.; Anghel, M.; Nishikawa, T. Spontaneous synchrony in power-grid networks. Nat. Phys. 2013, 9, 191–197. [CrossRef]
  24. Li, X.; Wei, W.; Zheng, Z. Promoting synchrony of power grids by restructuring network topologies. Chaos: Interdiscip. J. Nonlinear Sci. 2023, 33. [CrossRef]
  25. Chen, G. Searching for Best Network Topologies with Optimal Synchronizability: A Brief Review. Ieee/caa J. Autom. Sin. 2022, 9, 573–577. [CrossRef]
  26. Y. Kuramoto and D. Battogtokh, Coexistence of Coherence and Incoherence in Nonlocally Coupled Phase Oscillators, Physics (College. Park. Md). 4, 380 (2002).
  27. Martens, E.A.; Thutupalli, S.; Fourrière, A.; Hallatschek, O. Chimera states in mechanical oscillator networks. Proc. Natl. Acad. Sci. 2013, 110, 10563–10567. [CrossRef]
  28. Panaggio, M.J.; Abrams, D.M. Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators. Nonlinearity 2015, 28, R67–R87. [CrossRef]
  29. Ding, L.; Gonzalez-Longatt, F.M.; Wall, P.; Terzija, V. Two-Step Spectral Clustering Controlled Islanding Algorithm. IEEE Trans. Power Syst. 2012, 28, 75–84. [CrossRef]
  30. Ajala, O.; Dominguez-Garcia, A.; Sauer, P.; Liberzon, D. A Second-Order Synchronous Machine Model for Multi-Swing Stability Analysis, 51st North Am. Power Symp. NAPS 2019 (2019).
  31. Alberto, L.; Bretas, N. Required damping to assure multiswing transient stability: the SMIB case. Int. J. Electr. Power Energy Syst. 2000, 22, 179–185. [CrossRef]
  32. Karatekin, C.Z.; Uçak, C. Sensitivity analysis based on transmission line susceptances for congestion management. Electr. Power Syst. Res. 2008, 78, 1485–1493. [CrossRef]
  33. Mei, S.; Ni, Y.; Wang, G.; Wu, S. A Study of Self-Organized Criticality of Power System Under Cascading Failures Based on AC-OPF With Voltage Stability Margin. IEEE Trans. Power Syst. 2008, 23, 1719–1726. [CrossRef]
  34. Dobson, I.; Carreras, B.; Lynch, V.; Newman, D. An Initial Model for Complex Dynamics in Electric Power System Blackouts, Proc. Hawaii Int. Conf. Syst. Sci. 51 (2001).
  35. T. Nishikawa and A. E. Motter, Symmetric States Requiring System Asymmetry, Phys. Rev. Lett. 117, 114101 (2016).
  36. Zhou, L.; Swain, A.; Ukil, A. Reinforcement Learning Controllers for Enhancement of Low Voltage Ride Through Capability in Hybrid Power Systems. IEEE Trans. Ind. Informatics 2019, 16, 5023–5031. [CrossRef]
  37. Dörfler, F.; Bullo, F. Synchronization in complex networks of phase oscillators: A survey. Automatica 2014, 50, 1539–1564. [CrossRef]
  38. Zou, Y.; Pereira, T.; Small, M.; Liu, Z.; Kurths, J. Basin of Attraction Determines Hysteresis in Explosive Synchronization. Phys. Rev. Lett. 2014, 112, 114102. [CrossRef]
  39. Chialvo, D.R. Emergent complex neural dynamics. Nat. Phys. 2010, 6, 744–750. [CrossRef]
  40. Fan, H.; Wang, Y.; Wang, X. Eigenvector-based analysis of cluster synchronization in general complex networks of coupled chaotic oscillators. Front. Phys. 2023, 18, 1–15. [CrossRef]
Figure 1. Stability boundary. (a). Visualization of the boundary. The stability boundary (blue surface) and the planes 0 − uLδKL, uK − 0 − δKL and uKuL − 0 (pink surface) together enclose the stability domain. (b). Synchronization stabilization discrimination. The plane |uK| = |uL| is not shown. Δt is the fault duration. At Δt = 0.154s, all operating points are clustered at the boundary(cyan dots). At Δt = 0.155s, the three operating points are outside the boundary and away from the other points(magenta dots). dt) = 0.001s is a step length commonly used in power system studies. (c). Multiswing stabilization discrimination (Δt = 0.155s). The numbers 9, 10 denote the moment T9 = 9s, T10 = 10s, respectively. The figure shows the trajectory of (u1, u2, δ1,2). The operating point crosses the boundary outwardly, and δkl rapidly increases after the voltage changes. (d) and (e) are the results of the numerical experiment. The maximum value of δi is represented by the magenta curve and the minimum value is represented by the cyan curve. δminδiδmax. (d), δi of all meta-generators are close to each other and have approximately the same rate of change. (e), in time interval (9s, 10s), δmaxδmin becomes sharply larger.
Figure 1. Stability boundary. (a). Visualization of the boundary. The stability boundary (blue surface) and the planes 0 − uLδKL, uK − 0 − δKL and uKuL − 0 (pink surface) together enclose the stability domain. (b). Synchronization stabilization discrimination. The plane |uK| = |uL| is not shown. Δt is the fault duration. At Δt = 0.154s, all operating points are clustered at the boundary(cyan dots). At Δt = 0.155s, the three operating points are outside the boundary and away from the other points(magenta dots). dt) = 0.001s is a step length commonly used in power system studies. (c). Multiswing stabilization discrimination (Δt = 0.155s). The numbers 9, 10 denote the moment T9 = 9s, T10 = 10s, respectively. The figure shows the trajectory of (u1, u2, δ1,2). The operating point crosses the boundary outwardly, and δkl rapidly increases after the voltage changes. (d) and (e) are the results of the numerical experiment. The maximum value of δi is represented by the magenta curve and the minimum value is represented by the cyan curve. δminδiδmax. (d), δi of all meta-generators are close to each other and have approximately the same rate of change. (e), in time interval (9s, 10s), δmaxδmin becomes sharply larger.
Preprints 103679 g001
Figure 2. Unique spatiotemporal structure. The fault duration Δt increased from 0.140 s to 0.154 s. The arrow shows the direction of increase in Δt (the 18-node three-phase short circuit to ground fault). dt) = 0.001s. (a). σ(δi) denotes the standard deviation of δ, which started to decrease at 0.147 s and rose by 1300% at 0.155 s when the system became unstable. (b). In the δ1δ2 plane, the elliptical area marks the position of the potential barrier. From 0.147s onward the interval between operating points decreased in the direction of increasing Δt (shaded area). (c). In the u1u2 plane, the graph is a cycle structure that appears at 0.147 s.
Figure 2. Unique spatiotemporal structure. The fault duration Δt increased from 0.140 s to 0.154 s. The arrow shows the direction of increase in Δt (the 18-node three-phase short circuit to ground fault). dt) = 0.001s. (a). σ(δi) denotes the standard deviation of δ, which started to decrease at 0.147 s and rose by 1300% at 0.155 s when the system became unstable. (b). In the δ1δ2 plane, the elliptical area marks the position of the potential barrier. From 0.147s onward the interval between operating points decreased in the direction of increasing Δt (shaded area). (c). In the u1u2 plane, the graph is a cycle structure that appears at 0.147 s.
Preprints 103679 g003
Figure 3. Second-order phase transition, long-range correlation and the potential barrier. (a). Near the boundary, δi undergone a second-order phase transition at Δt = 0.147s. Amplified points were the phase transition points. The disturbed trajectories of all points were almost identical. (b)-(i). In the δKδL planes, all operating points crossed the potential barrier at Δt = 0.147s. The images are almost identical due to the effect of long-range correlation.
Figure 3. Second-order phase transition, long-range correlation and the potential barrier. (a). Near the boundary, δi undergone a second-order phase transition at Δt = 0.147s. Amplified points were the phase transition points. The disturbed trajectories of all points were almost identical. (b)-(i). In the δKδL planes, all operating points crossed the potential barrier at Δt = 0.147s. The images are almost identical due to the effect of long-range correlation.
Preprints 103679 g004
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated