A Novel Fast Path Planning Approach for Mobile Devices using Hybrid Quantum Ant Colony Optimization Algorithm

Mayukh Sarkar1, Jitesh Pradhan1, Anil Kumar Singh2, Hathiram Nenavath3

1Department of Computer Science and Engineering, National Institute of Technology Jamshedpur, India

2Department of Computer Science and Engineering, Motilal Nehru National Institute of Technology Allahabad, India

3Department of Electronics and Communication Engineering, National Institute of Technology Jamshedpur, India

With IoT systems' increasing scale and complexity, maintenance of a large number of nodes using stationary devices is becoming increasingly difficult. Hence, mobile devices are being employed that can traverse through a set of target locations and provide the necessary services. In order to reduce energy consumption and time requirements, the devices are required to traverse following a Hamiltonian path. This problem can be formulated as a Travelling Salesman Problem (TSP), an NP-hard problem. Moreover, in emergency services, the devices must traverse in real-time, demanding speedy path planning from the TSP instance. Among the well-known optimization techniques for solving the TSP problem, Ant Colony Optimization has a good stronghold in providing good approximate solutions. Moreover, ACO not only provides near-optimal solutions for TSP instances but can also output optimal or near-optimal solutions for many other demanding hard optimization problems. However, to have a fast solution, the next node selection, which needs to consider all the neighbors for each selection, becomes a bottleneck in the path formation step. Moreover, classical computers are constrained to generate only pseudorandom numbers. Both these problems can be solved using quantum computing techniques, i.e., the next node can be selected with proper randomization, respecting the provided set of probabilities in just a single execution and single measurement of a quantum circuit. Simulation results of the proposed Hybrid Quantum Ant Colony Optimization algorithm on several TSP instances have shown promising results, thus expecting the proposed work to be important in implementing real-time path planning in quantum-enabled mobile devices.

Mobile IoT Devices, Path Planning, Travelling Salesman Problem, Ant Colony Optimization, Quantum Computing.

§ INTRODUCTION

As the scale and complexity of the IoT systems are increasing, the maintenance of a large number of stationary sensors at target locations is becoming costlier. Hence, mobile devices such as Unmanned Ground Vehicles and Unmanned Aerial Vehicles are being deployed. Rather than being stationary, these devices can traverse through target locations and provide the necessary services <cit.>. However, providing the services at target locations in real time requires extremely fast planning of possible paths. Moreover, the devices need to traverse as little distance as possible while visiting each target location. This problem can be formulated as a Travelling Salesman Problem (TSP), a well-known NP-hard problem.

The intractability of the Travelling Salesman Problem involves the applications of various optimization techniques that can provide near-optimal routes in a graph of target nodes. The most notable techniques among them include Ant Colony Optimization (ACO). Colorni et al. proposed the possibility of optimization following the natural phenomena occurring in ant colonies in 1991 <cit.>. Dorigo and Gambardella successfully utilized this concept in optimizing the famous Travelling Salesman Problem in 1997 <cit.>. Since then, ACO has been utilized with great success in solving multitudes of problems that range from basic optimization up to as diverse as analog design <cit.>, reversible circuit synthesis <cit.>, etc.

Though ACO has been quite successful in obtaining near-optimum and optimum solutions for many optimization problems, including TSP, its execution needs to be sped up for real-time solutions such as path planning in mobile devices. One of the major bottlenecks in its speed, especially with a large number of nodes, is in selecting the next node. Any TSP-solving algorithm eventually needs to find the complete Hamiltonian cycle. However, while finding a path, selecting the next node among many possible neighbors needs to be executed fast to obtain the complete path as fast as possible. In the classical implementation of the ACO, finding the next neighbor involves calculating the cumulative probabilities of all neighbors and then selecting one based on a pseudorandom number. Quantum hardware can provide substantial speed-up in such scenarios by processing all probabilities together by a single execution of a quantum circuit, thanks to its inherent parallel nature obtained from superposition. Though the speed-up will be linear, with a large number of nodes and a sufficiently large number of ants, the overall accumulation of the speed-up starts to show its effects. This also has the advantage of selecting the next node based on pure randomization rather than using pseudorandomization as obtained from the classical counterpart. This paper attempts to achieve the same by proposing a novel hybrid classical-quantum technique of ACO.

Quantum computing started its journey from the proposal of a quantum equivalent of the Turing machine by Paul Benioff <cit.><cit.>. In 1982, Richard Feynman postulated that simulating quantum physics on classical hardware would be impossible and suggested using quantum mechanical hardware for more efficient computation <cit.>. The real power of quantum computing technique came into view when Peter Shor, in 1994, proposed quantum algorithms for prime factorization and discrete logarithm problems, hence showing the possibility of breaking RSA and Diffie-Hellman encryption techniques <cit.>. In 1996, Lov Grover showed the possibility of speed-up of search in an unstructured database using quantum system <cit.>. These results immediately brought quantum computing to the forefront of futuristic computational technique research. The declaration of quantum supremacy by Google using their Sycamore processor <cit.>, and the possibility of universal quantum logic operations with high average fidelity <cit.> has raised the curtain of those possibilities which have been only theoretical concept for long.

The current paper attempts to harness this power of quantum computation to speed up the well-known classical algorithm of ACO and employ it in possible real-time path planning applications as required in mobile devices. The rest of the paper is organized as follows. Section <ref> provides the necessary overview of quantum computing and the ant colony optimization algorithm. Several attempts have been made to provide a quantum version of ACO in literature. Section <ref> provides a brief overview of previous attempts and the novelties of the current paper. Section <ref> demonstrates the complete Hybrid Quantum Ant Colony Optimization algorithm as proposed in this paper, step-by-step, with an example of an instance of TSP. Section <ref> shows the performance of the proposed algorithm on several instances of the TSP obtained from a well-known collection of TSP instances, namely TSPLIB <cit.>. The result shows that the proposed algorithm can provide optimal or near-optimal solutions to TSP problems with sufficient speed on real quantum hardware. Section <ref> finally concludes the paper.

§ BACKGROUND INFORMATION

In this section, required background information on quantum computing and ant colony optimization is being demonstrated.

§.§ Quantum Computing

In classical computing, the basic unit of information, i.e., a bit, can stay in one of the two states, either 0 or 1, at any particular time. But in a quantum computer, a qubit can stay as a superposition of the two basis states $\ket{0}$ and $\ket{1}$, where state $\ket{0}$ corresponds to classical state 0 and state $\ket{1}$ corresponds to classical state 1. Mathematically, quantum state of a qubit can be represented as a vector $\ket{\psi} = \begin{bmatrix} \alpha \\ \beta \end{bmatrix} \in  \mathbb{C}^2 $, where $\alpha$ and $\beta$ are amplitudes associated with the states $\ket{0}$ and $\ket{1}$ respectively. Upon measurement, the quantum state collapses to one of the basis states and outputs the corresponding classical state. As for example, for the above-mentioned quantum state, upon measurement, the quantum state will collapse to state $\ket{0}$ and will output classical bit $0$ with probability $|\alpha|^2$, and it will collapse to state $\ket{1}$ and will output classical bit $1$ with probability $|\beta|^2$. Thus, the basis state $\ket{0}$ is associated with the vector $\ket{0} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}$ and basis state $\ket{1}$ is associated with the vector $\ket{1} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}$. An arbitrary quantum state $\ket{\psi}$ can thus be represented as follows.

\begin{equation*}
\ket{\psi} = \alpha \ket{0} + \beta \ket{1}
\end{equation*}

Moreover, as the measurement of a qubit should always result in either 0 or 1, it must follow that $|\alpha|^2 + |\beta|^2 = 1$.

On a multiqubit system, the qubits involved in the system, remain entangled to each other. As an example, on a two-qubit system, the state $\ket{\psi}$ becomes the superposition of all possible two-bit basis states $\ket{00}$, $\ket{01}$, $\ket{10}$ and $\ket{11}$, and can be written as follows.

\begin{equation*}
\ket{\psi} = \alpha_{00}\ket{00} + \alpha_{01}\ket{01} + \alpha_{10}\ket{10} + \alpha_{11}\ket{11}
\end{equation*}

, where $|\alpha_{00}|^2 + |\alpha_{01}|^2 + |\alpha_{10}|^2 + |\alpha_{11}|^2 = 1$. On measurement of both qubits, the classical result will be either $00$, $01$, $10$ or $11$ with the probabilities equal to the square of their corresponding amplitudes. Upon measuring a single qubit, as for example the first qubit, it will output $0$ with probability $|\alpha_{00}|^2 + |\alpha_{01}|^2$ and will output $1$ with probability $|\alpha_{10}|^2 + |\alpha_{11}|^2$. In case the first qubit measures to $0$, the complete quantum system collapses to the state $\frac{\alpha_{00}}{\sqrt{|\alpha_{00}|^2 + |\alpha_{01}|^2}} \ket{00} + \frac{\alpha_{01}}{\sqrt{|\alpha_{00}|^2 + |\alpha_{01}|^2}} \ket{01}$. The similar post-measurement state can be obtained in case the first qubit measures to $1$.

The state of a qubit gets modified by quantum gates, and a linear sequence of quantum gates eventually forms a quantum circuit. There cannot be any fanout or loop in a quantum circuit. The state of a quantum system can be manipulated using an adequately designed sequence of gates to harness the inherent parallelism obtained from the superposition property. A quantum gate must be a unitary matrix, and conversely, any unitary matrix is well-suited for a quantum gate. Only a small set of quantum gates is majorly used in quantum architectures, especially the universal set of quantum gates. There are uncountably many possible sets of universal quantum gates. As an example, a Clifford gate set $\lbrace H, S, CNOT \rbrace$, along with a non-Clifford gate such as $T$ gate, is a well-known universal gate set. Now, any single qubit gate $U$ can be represented as a sequence of three rotational gates up to a global phase factor, i.e., for any unitary $U$, there exist four real numbers $\alpha$, $\beta$, $\gamma$, and $\delta$, such that

\begin{equation*}
U = e^{i\alpha}R_z(\beta)R_y(\gamma)R_z(\delta)
\end{equation*}

As all gates in the mentioned gate set are single qubit gates except $CNOT$ gate, it can be said that any functionality can be implemented using two rotation gates and $CNOT$ gate. A complete overview of the basics of quantum computing, along with the proof of the above representation, can be found in <cit.>.

In the proposed work, a quantum circuit will be utilized to produce a quantum state provided the probabilities of nodes to be selected, and then will be measured to obtain the next node. As the probabilities involved in the ACO are real, the circuit will need to generate only those quantum states that are in real vector space. This allows us to restrict ourselves to only $R_y$ and $CNOT$ gates. The quantum selector proposed in <cit.> has been utilized here for the next node selection and is built up of only controlled-$R_y$ and $CNOT$ gates. An overview of the proposed quantum selector has been given in subsection <ref>.

§.§ Ant Colony Optimization

Ant Colony Optimization (ACO) is a well-known evolutionary algorithm that can obtain optimal or near-optimal solutions to many demanding optimization problems. This technique was first proposed in <cit.>. In this technique, multiple ants move through some path in a graph in order to obtain optimal solutions. Each ant determines its path as follows.

* If there can be multiple possible start nodes, the ant can start from any node randomly, or otherwise, it will start from the root. Then, it follows a path from the starting node as long as it does not reach some solution or as long as it does not reach the leaf node, depending on the problem statement. While sitting at a node, the ant always tries to find the optimum next node probabilistically. This probability of the subsequent nodes depends on two criteria: firstly, the likeliness of the path, which generally is some function of the path's cost. In a minimization problem, the likeliness varies inversely with the cost; in a maximization problem, the likeliness varies proportionately with the cost. The second selection criterion represents the experience of previous ants in finding paths for optimal or near-optimal solutions. For the first ant, this factor is initialized to the same value for all edges. This factor is commonly known as pheromone. Now, let from the node $a$, the selection of the next node $b$ involves the likeliness $w_{ab}$ and the pheromone $\eta_{ab}$. Then, the probability of the node will be calculated as in equation <ref>.

\begin{equation}
\label{eqtn-aco-next-node-prob}
P(a \rightarrow b) = \frac{w_{ab}^\alpha \eta_{ab}^\beta}{\sum_{i \in N(a)}w_{ai}^\alpha \eta_{ai}^\beta}
\end{equation}

, where $\alpha$ and $\beta$ are two hyperparameters that needs to be tuned for better results, and sum in the denominator involves all neighbours of node $a$ ($N(a)$ represents the neighbour set of node $a$).

* Once the ant completes its tour, the edges' pheromones change along with the success of the ant following the rule as in equation <ref>.

\begin{equation} \label{eqtn-aco-pheromone-update}
\eta_{ab} = (1 - \rho)*\eta_{ab} + f(w)
\end{equation}

, where $\rho$ is the evaporation rate of the paths, and $f(w)$ is the predetermined function of the overall cost of the final solution found if the ant is successful. If the ant fails to reach any solution, $f(w)$ is taken as zero. If the ant is successful, $f(w)$ varies inversely with $w$ for minimization problems and proportionately for maximization problems. In simpler terms, if the ant fails, the pheromones will only be evaporated, and nothing will be deposited. Furthermore, if the ant succeeds, some amount of pheromone gets deposited depending on the likeliness of the solution found, i.e., the better the solution, the higher the amount of pheromone deposited. For a large number of ants, this eventually allows all ants to converge toward the optimal solution by depositing more pheromones over the optimal paths.

§ RELATED WORKS

Wang et al. <cit.><cit.> introduced the novel concept of implementing Ant Colony Optimization on a Quantum Computer. In this work, the authors have implemented discrete binary ACO by representing the pheromone, i.e., the probability of the selected path, by qubits. A similar approach of quantum ACO has been employed in implementations of collision detection in virtual environment <cit.>, privacy-preserving data mining algorithm <cit.>, optimization of evacuation path <cit.>, rule-based fuzzy classifier <cit.>. Wang et al., in the same paper <cit.>, have also proposed a quantum-inspired ACO employing the same approach, and has been utilized in link prediction <cit.>. A similar approach to quantum-inspired ACO has been employed in optimizing queries in distributed database systems <cit.>. The major problem of these works is the updation of the qubit by rotation gates, whose rotation angle is determined by a function of the probability amplitudes of the qubit being applied. This technique, thus, involves knowing the state of the qubit. However, on real quantum hardware, measuring a qubit would collapse its state, and hence, these algorithms cannot be implemented on a real quantum computer.

A different approach to quantum implementation of ACO has been performed by representing the ant positions by coordinates on Bloch sphere <cit.><cit.>. A similar idea has also been employed in various optimization problems, such as in path planning of automated guided vehicles <cit.>. However, the major problem with this approach is its lack of consistency with the original concept of Ant Colony Optimization.

Ghosh et al. <cit.> have proposed another novel approach to implementing ACO on a quantum computer named MNDAS. This paper employs a set of qubits to encode all paths, a separate set of qubits to represent pheromones, and three qubits as registers, requiring many qubits in a sufficiently large instance of graphs. Moreover, the proposed architecture in the paper employs a large number of SWAP gates, resulting in the tendency of introducing noise <cit.>. Thus, such architecture is not implementable in current Noisy Intermediate-Scale Quantum (NISQ) era quantum computers.

De Andoin et al. <cit.> have proposed an implementable ACO on a hybrid quantum system. The algorithm in this paper also does not follow the conventional approach of ACO and hence faces problems in constrained optimization problems. The paper proposes generating valid solutions from invalid ones in such cases, but that requires the knowledge of the complete set of valid solutions, which is impossible in most optimization problems.

In the current paper, the authors have proposed a novel implementation of ACO on a hybrid quantum computer that follows the exact steps of classical ACO. This gives the freedom of employing ACO in solving optimization problems, irrespective of being constrained or unconstrained, in the same manner as its classical counterparts. The proposed technique in this paper does not encode ants, pheromones, or paths using qubits like previous approaches in the literature. Instead, the arithmetic computations in the original ACO are assigned to a classical processor, and the quantum coprocessor makes the probabilistic selection of the next node. This design is novel because quantum computers are not known to provide sufficient speed improvements in arithmetic computations over classical ones. Moreover, in the present NISQ era, the fault-tolerance of classical computers is much more reliable than quantum computers. Along with possible speed-up in selecting the next node, employing a quantum computer in the subsequent node selection allows us to perform the selection with pure randomization rather than employing pseudorandom numbers as in classical computers.

§ PROPOSED IMPLEMENTATION OF ANT COLONY OPTIMIZATION ON HYBRID QUANTUM SYSTEM

This section demonstrates the proposed implementation of the Ant Colony Optimization (ACO) algorithm on a hybrid Quantum-Classical system. Firstly, all the steps involved in the ACO algorithm can be implemented on a quantum computer because any classical algorithm can be simulated on a quantum computer efficiently using Fredkin Gate <cit.>. Moreover, several implementations of classical logical and arithmetic operations on a quantum computer have been proposed in literature <cit.><cit.><cit.>. However, in the current NISQ era, quantum computers are not fault-tolerant enough, and they can lose coherence before executing for sufficiently long <cit.>. Thus, the algorithms employing a quantum-classical hybrid approach in solving problems rather than employing quantum hardware solely, such as Variational Quantum Eigensolver (VQE) <cit.> and Quantum Approximation Optimization Algorithm (QAOA) <cit.> have made immediate impact in this NISQ era of quantum computation.

The proposed algorithm in this paper employs the exact steps of ACO as its classical counterpart, except a suitable quantum circuit is performing the next node selection. Besides the possibility of implementing ACO on a futuristic NISQ respecting quantum system, the proposed Hybrid Quantum Ant Colony Optimization (HQACO) has the following advantages.

* Possibility of solving optimization problems following the same approach as on classical computers. There is no need for some special preprocessing of input data other than the ones necessary for applying conventional ACO.

* Employing a quantum computer in selecting the next node has two major advantages. In the classical counterpart, a pseudorandom number is generated, and the cumulative probability of all possible next nodes is compared one by one with the generated number. The quantum selector will parallel process all possible next nodes in a single execution and select the next node with pure randomization respecting the exact probabilities. Thus, the quantum selector will speed up the next node selection and allow us to implement the ACO perfectly by enacting pure randomization rather than using pseudorandom numbers generated by classical random number generators.

The current section elaborates on the proposed algorithm as follows. Firstly, subsection <ref> demonstrates the design of a quantum selector that, upon providing a set of probabilities adding up to 1, will select one of the elements respecting their corresponding probabilities. The proposed HQACO algorithm utilizes this selector to select the next node. The complete step-by-step details of HQACO, with an example of the Travelling Salesman Problem, have been demonstrated in subsection <ref>.

§.§ Design of the Quantum Selector

In <cit.>, the author has demonstrated an algorithm that given a set of probabilities $\lbrace p_i \rbrace$ such that $\sum p_i = 1$, will generate a parameterized quantum circuit generating the output state $\sum{\sqrt{p_i} \ket{i}}$. As an example, given set of probabilities $\lbrace 0.3, 0.2, 0.4, 0.1 \rbrace$, the algorithm will generate a two-qubit circuit producing the state $(\sqrt{0.3} \ket{00} + \sqrt{0.2} \ket{01} + \sqrt{0.4} \ket{10} + \sqrt{0.1} \ket{11})$. Observe that, upon measurement of the two qubits, the classical output 00 will appear with probability 0.3, 01 will appear with probability 0.2, 10 will appear with probability 0.4 and 11 will appear with probability 0.1. In this subsection, a brief overview of the algorithm is being provided. Detailed discussion with its correctness can be found in <cit.>.

* Generate spherical coordinates from the set of amplitudes employing algorithm <ref>. As for example, from the set of amplitudes $\lbrace \sqrt{0.3}, \sqrt{0.2}, \sqrt{0.4}, \sqrt{0.1} \rbrace$, the algorithm will generate the spherical coordinate angles $\lbrace r_1, r_2, r_3 \rbrace = \lbrace 1.982313, 2.013707, 0.927295 \rbrace$. It can be readily verified that

\begin{align*}
cos(\frac{r_1}{2}) &= \sqrt{0.3} \\
sin(\frac{r_1}{2})cos(\frac{r_2}{2}) &= \sqrt{0.2} \\
sin(\frac{r_1}{2})sin(\frac{r_2}{2})cos(\frac{r_3}{2}) &= \sqrt{0.4} \\
sin(\frac{r_1}{2})sin(\frac{r_2}{2})sin(\frac{r_3}{2}) &= \sqrt{0.1}
\end{align*}

Array of amplitudes $A = \lbrace a_1, a_2, \cdots, a_n \rbrace$ of length $n$
Array of spherical angles $R = \lbrace r_1, r_2, \cdots r_{n-1} \rbrace$ of length $n-1$
$i\leftarrow 1$ $n-1$
$a_i$ is $0$
Append 0 to $R$.

Fill the rest of $R$ with $\pi$.

Set $r_i \leftarrow 2 * arccos(a_i)$.

Divide all elements of $A$ by $sin(r_i / 2)$.

* For single-qubit circuit generating 2-dimensional vector, the circuit will be, trivially a single $R_y$ gate with the angle set to the spherical angle generated in step 1. Observe that, from a 2-dimensional amplitude vector $\ket{\psi} = \begin{bmatrix} a_1 \\ a_2 \end{bmatrix}$, step 1 will generate a single angle $r_1$ such that $\begin{bmatrix} a_1 \\ a_2 \end{bmatrix} = \begin{bmatrix} cos(\frac{r_1}{2}) \\ sin(\frac{r_1}{2}) \end{bmatrix}$. The circuit is as shown in figure <ref>.

\[
\begin{array}{c}
\Qcircuit @C=1em @R=1.6em {
\lstick{\ket{0}} & \gate{R_y(r_1)} & \qw & \rstick{\ket{\psi}}
\end{array}
\]

Single-qubit arbitrary state generator

* For two-qubit circuit generating 4-dimensional statevector $\ket{\psi} = [a_1\ a_2\ a_3\ a_4]^T$, first the three parameters $\lbrace r_1, r_2, r_3 \rbrace$ will be generated following algorithm <ref>. Then the circuit genrating $\ket{\psi}$ will be as in figure <ref>.

\[
\begin{array}{c}
\Qcircuit @C=1em @R=1.6em {
\lstick{} & \gate{R_y(r_1)} & \ctrl{1} & \gate{R_y(\pi+r_3)} & \qw \\
\lstick{} & \qw & \gate{R_y(-r_2)} & \ctrl{-1} & \qw\gategroup{1}{5}{2}{5}{2.0em}{\}} & \rstick{\ket{\psi}}
\inputgroupv{1}{2}{.8em}{.8em}{\ket{0}}
\end{array}
\]

Two-qubit arbitrary state generator

* For any n-qubit circuit with $n > 2$ generating $2^n$-dimensional statevector $\ket{\psi}$, the $2^n-1$ dimensional vector of spherical angles will be generated following algorithm <ref>. The spherical angle vector will be provided as the $2^n-1$ parameters to the state generator circuit as shown in figure <ref>. In the circuit, $G_1$ is a recursive (n-1)-qubit arbitrary state generator circuit with first $(2^{n-1} - 1)$ parameters $\lbrace r_1, r_2, \cdots, r_{2^{n-1}-1} \rbrace$, the multicontrolled $R_y$ gate has all the above (n-1) lines as controls and employs $(2^{n-1})^{th}$ parameter. Next, there is a group of (n-1) CNOT gates, each having control on $n^{th}$ line and target to each of the above (n-1) lines respectively. Finally, there is another recursive (n-1)-qubit arbitrary state generator circuit $G_2$ with last $(2^{n-1} - 1)$ parameters $\lbrace r_{2^{n-1}+1}, \cdots, r_{2^{n}-1} \rbrace$. Each gate on $G_2$ has an additional control from the $n^{th}$ line.

\[
\begin{array}{c}
\Qcircuit @C=1em @R=1.6em {
\lstick{} & \multigate{3}{G_1} & \ctrl{4} & \qw    & \qw & \qw          & \targ  & \multigate{3}{G_2} & \qw \\
\lstick{} & \ghost{G_1}        & \ctrl{3} & \qw    & \qw & \targ        & \qw    & \ghost{G_2}        & \qw \\
\cdots    & \nghost{G(_1}      & \cdots   & & \cdots                & \nghost{G_2}   &    & \cdots & \\
\lstick{} & \ghost{G_1}        & \ctrl{1} & \targ  & \qw  & \qw        & \qw    & \ghost{G_2}        & \qw \\
\lstick{} & \qw                & \gate{R_y(r_{2^{n-1}})} & \ctrl{-1} & \qw & \ctrl{-3} & \ctrl{-4} & \ctrl{-1} & \qw
\end{array}
\]

Multiqubit arbitrary state generator

Employing the stated algorithm, a 3-qubit arbitrary state generator circuit generating statevector $\lbrace a_1, a_2, \cdots, a_8 \rbrace$ can be designed as in figure <ref>, where the parameters generated from the input statevector, employing algorithm <ref> is $\lbrace r_1, r_2, \cdots, r_7 \rbrace$. The circuit generated using the stated algorithm can generate statevectors with only real amplitudes, which will be sufficient for the proposed HQACO algorithm, as the probabilites generated in the intermediate stages are always real numbers.

\[
\begin{array}{c}
\Qcircuit @C=1em @R=1.6em {
\lstick{} & \gate{R_y(r_1)} & \ctrl{1} & \gate{R_y(\pi+r_3)}  & \ctrl{2} & \qw & \targ & \gate{R_y(r_5)} & \ctrl{1} & \gate{R_y(\pi+r_7)} & \qw \\
\lstick{} & \qw & \gate{R_y(-r_2)} & \ctrl{-1} & \ctrl{1} & \targ & \qw & \qw & \gate{R_y(-r_6)}& \ctrl{-1} & \qw \\
\lstick{} & \qw & \qw & \qw & \gate{R_y(r_4)} & \ctrl{-1} & \ctrl{-2} & \ctrl{-2} & \ctrl{-1} & \ctrl{-2} & \qw
\end{array}
\]

Example of 3-qubit arbitrary state generator

§.§ HQACO Algorithm

In this subsection, the quantum selector circuit designed in subsection <ref> will be employed to perform Ant Colony Optimization on a hybrid quantum computer. In demonstrating the algorithm, an instance of the Travelling Salesman Problem on a small graph, as shown in figure <ref>, will be employed.

An Example Graph for HQACO

* A parametrized quantum selector with $q = max(\lceil log_2(v) \rceil, \lceil log_2(d) \rceil)$ qubits needs to be prepared with the structure as mentioned in subsection <ref>, where $v$ is the number of nodes, and $d$ is the maximum degree (outdegree for directed graph) of the graph. The given example graph has seven vertices, and the maximum degree is 4. Thus, it is sufficient to have a 3-qubit quantum selector, as depicted in figure <ref>. If the number of possible candidates, and hence, the number of probabilities to choose from, is less than $2^q$ (8, in the example), append 0 at the end of the amplitude array to make it of length $2^q$.

* The graph information is stored in a classical register $C_1$. Each edge has to be initialized with a pheromone value, which will be kept in another classical register $C_2$. Initially, all edges receive the same pheromone value. The initial pheromone value has been taken as $0.5$ for all the edges in the given example.

* Suitable values of the hyperparameters $\alpha$, $\beta$, and pheromone evaporation rate $\rho$ are chosen and kept in classical register $C_3$. In the given example, let us choose the values $\alpha=0.4$, $\beta=0.6$ and $\rho=0.01$ respectively. This will allow the selection to be more directed by pheromone values rather than edge costs.

* Perform the following sequence of operations for each ant.

* The start node for the ant is to be selected with equal probability. Set the parameters of the quantum selector so that the probability of each node becomes $\frac{1}{v}$. In the example with $7$ nodes, the $7$ parameters of the quantum selector will be the angles obtained using algorithm <ref> from the set of amplitudes $\lbrace \sqrt{\frac{1}{7}}, \sqrt{\frac{1}{7}}, \sqrt{\frac{1}{7}}, \sqrt{\frac{1}{7}}, \sqrt{\frac{1}{7}}, \sqrt{\frac{1}{7}}, \sqrt{\frac{1}{7}}, 0 \rbrace$. After calculating the parameters, the classical processor will set them to the quantum selector circuit for further processing.

* Quantum registers are initialized with $\ket{0}$. The quantum selector will produce the desired state. All the qubits of the selector will be measured and stored in a classical register $C_3$. The classical register holds the start node index with equal probabilities. Suppose the initial node has been measured as $4$ in the example considered.

*  The classical processor will read the node from the register $C_3$ and calculate the probabilities of the next transitions. To calculate this, the classical processor will store all unvisited neighbors and consider the inverse of the corresponding edge costs as the corresponding weights. The probability of a particular incident edge $e$ is then evaluated as $\frac{w(e)^\alpha ph(e)^\beta}{\sum_{i=1}^d w(e)^\alpha ph(e)^\beta}$, where $w(e)$ and $ph(e)$ are the edge weights and current pheromone values respectively, and $d$ is the degree of the current node. As for example, with start node $4$ and all pheromone values $0.5$, the probability of going to neighbor nodes $\lbrace 2, 3, 5, 6 \rbrace$ are $\lbrace 0.32625, 0.19402, 0.24725, 0.23247 \rbrace$, as the corresponding edge costs are $\lbrace \frac{1}{3}, \frac{1}{11}, \frac{1}{6}, \frac{1}{7} \rbrace$ respectively. Observe that the cost has been taken as the inverse of the edge weight, as it is a minimization problem.

*  Evaluate the angles using algorithm <ref> from the zero-appended amplitude set, where each amplitude is the square root of the corresponding probabilities. Thus, in the continuing example, the classical processor will set the parameters of the quantum selector after evaluating them from the amplitude set $\lbrace 0.57118, 0.44048, 0.49724, 0.48215, 0, 0, 0, 0 \rbrace$.

*  Reset the qubits in the quantum selector to $\ket{0}$, and allow them to pass through the quantum selector with set parameters. At the end, the qubits will be measured to obtain the next node index. For example, if the measurement result turns out to be $010$, the node at index 2 (0-indexed), i.e., node $5$, will be selected as the next node. Set the index in classical register $C_3$.

* Continue the steps <ref> to <ref>, as long as no unvisited next node can be found. After completion of the tour, reduce the pheromone of each edge by pheromone evaporation rate $\rho$.

* Two events can occur at this point: firstly, the ant has visited all nodes, in which case the ant is said to be successful, or secondly, the ant is stuck without any unvisited neighbor, but all nodes have not been visited. The tour is a failure in the second scenario, and the result will be discarded. In the first scenario, increase each edge's pheromone on the successful path by $\frac{1}{c}$, where $c$ is the total cost of the path.

* If a Hamiltonian cycle is to be obtained, another condition needs to be added for a path to be successful, and that is, the last node on the path must have the starting node as its neighbor.

In the given example, the above step has been performed with 100 ants, and the minimum Hamiltonian cycle $\lbrace 6, 5, 3, 1, 0, 2, 4, 6 \rbrace$ has been obtained with cost of 63.

§.§ Expansion to parallel ant systems

The proposed HQACO algorithm can be expanded naturally to parallel ant systems, where multiple ants will start in parallel rather than sequentially at each iteration. On a classical computer, this is achieved by multithreads, which will be executed on multiprocessors or multicores, whichever is available. On a quantum computer, the same can be achieved by placing multiple quantum co-processors, each co-processor being selected by each ant in parallel. Thus, in such a system, each ant will have a separate unit of classical-quantum hybrid processor.

§ EXPERIMENTAL RESULTS

This section demonstrates the experimental results of the proposed HQACO algorithm on several instances of the Travelling Salesman Problem obtained from TSPLIB <cit.>, a well-known collection of TSP instances for experimentations. The experiment has been performed using Python 3.11 with Qiskit 0.43 in a locally installed Linux environment. For larger instances, the Qiskit provider for Amazon Braket has been used and implemented on Amazon AWS. The table <ref> demonstrates the experimental results for 5 TSP instances, with the minimum path found and the overall estimation error. In each instance, 1000 ants have been used with $\alpha=0.4$, $\beta=0.6$ and $\rho=0.01$. Ten shots have been used during each measurement, and the output with the highest appearance count has been selected, as a large number of shots would make the output nearly deterministic by driving the output towards the highest probability.

of Nodes
2*Lower Bound
In Specification <cit.>
2c|Minimum Cycle Found using HQACO
in Estimation
2*Number of
Qubits used

Cost     Path

2*gr17.tsp     2*17     2*2085     2*2130     [12, 16, 6, 7, 5, 13, 14, 2, 10,      2*2.16%     5 for node selection

9, 1, 4, 8, 11, 15, 3, 0, 12]        4 for path selection

2*burma14.tsp     2*14     2*3323     2*3547     [12, 6, 11, 5, 4, 3, 2, 13,      2*6.74%     2*4

1, 9, 10, 8, 7, 0, 12]

2*gr21.tsp     2*21     2*2707     2*3216     [2, 8, 4, 15, 11, 3, 0, 5, 7, 6, 18,      2*18.8%    2*5

19, 10, 16, 9, 17, 12, 13, 14, 1, 20, 2]

2*bayg29.tsp     2*29     2*1610     2*1922     [23, 12, 15, 6, 22, 7, 27, 0, 5, 11, 8, 2, 28, 25, 4,      2*19.38%     2*5

20, 1, 19, 9, 3, 14, 17, 13, 21, 16, 18, 24, 10, 26, 23]

2*bays29.tsp     2*29     2*2020     2*2599     [6, 24, 18, 15, 14, 10, 13, 21, 16, 17, 3, 9, 19, 4,      2*28.66%     2*5

8, 11, 5, 27, 0, 23, 7, 22, 25, 28, 2, 1, 20, 12, 26, 6]

Experimental Results

In the first instance with 17 nodes, for each ant, a 5-qubit circuit is required to find the starting node, whereas a 4-qubit circuit is sufficient to find the path as each node has a degree of 16. In the simulation experiment, two separate circuits have been used for the same. The same can be followed and suggested in real quantum hardware rather than using a single 5-qubit circuit. Though it will increase the build cost, using a 5-qubit circuit, even where a 4-qubit circuit is sufficient, will increase the execution time and probability of noise in output.

The results have been sorted based on errors in estimation, and it can be observed that the error percentage increases with the number of qubits in the circuit. This output is reasonably expected, as the larger instances require more ants to converge, and the first 1000 ants in each case have shown apparent convergence towards a lower value. In the simulation, though larger circuits require exponential time to execute, they are expected to run extremely fast on real quantum hardware, as quantum simulation on classical hardware is exponential.
Hence, on real quantum hardware, a much larger number of ants can be utilized to obtain better results.

The above discussion shows that on real quantum hardware, the proposed HQACO algorithm can obtain near-optimum results of TSP instances quite fast, which is necessary for real-time path planning of mobile agents.

§ CONCLUSIONS

This paper proposes a novel implementation of the Ant Colony Optimization algorithm on a hybrid quantum computer. This algorithm alleviates two essential problems in implementing the same on a classical computer. Firstly, during path formation, each ant selects the next node probabilistically, which involves calculations of the cumulative probability of each node, making it linear in time. Secondly, the random number generated during node selection is not a pure random number, as classical computers can generate only pseudorandom numbers. Current work performs all arithmetic calculations, such as probability calculations, on a classical computer. Quantum computer in the NISQ era is not well-known enough to perform speedier arithmetic calculation reliably. Thus, the quantum phenomena have been applied at exactly that step where it would shine, i.e., the selection of the next node. A single execution of the quantum selector circuit, followed by measurement, will select the next node, respecting the assigned probabilities with pure randomness. Experimental results have shown that the speedy quantum circuit can reliably solve TSP instances with a suitable number of ants. Thus, this work is expected to be important in those areas where fast and accurate execution of ACO, such as solving TSP instances in real time for path planning of mobile devices.

§ ACKNOWLEDGEMENTS

This work would not have been possible without the financial support of MeitY Quantum Computing Applications Lab (QCAL). The authors are also thankful to Amazon AWS for providing the excellent environment for Quantum Computing research needed to perform the experiments using the Qiskit provider for Amazon Braket.


# System-Level Performance and Robustness of the Grid-Forming Hybrid Angle
Control

Ali Tayyebi∗†, Alan Magdaleno†, Denis Vettoretti∗, Meng Chen‡,
Eduardo Prieto-Araujo⋆, Adolfo Anta∗, and Florian Dörfler† ∗Electric Energy
Systems, Austrian Institute of Technology, Vienna, Austria.
†Automatic Control Laboratory (IfA), Swiss Federal Institute of Technology
(ETH), Zürich, Switzerland.
‡AAU Energy, Aalborg University, Aalborg, Denmark
⋆CITCEA, Universitat Politècnica de Catalunya (UPC), Barcelona, Spain.
Corresponding author’s e-mail address<EMAIL_ADDRESS>
###### Abstract

This paper investigates the implementation and application of the multi-
variable grid-forming hybrid angle control (HAC) for high-power converters in
transmission grids. We explore the system-level performance and robustness of
HAC concept in contrast to other grid-forming schemes i.e., power-frequency
droop and matching controls. Our findings suggests that, similar to the ac-
based droop control, HAC enhances the small-signal frequency stability in low-
inertia power grids, and akin to the dc-based matching control, HAC exhibits
robustness when accounting for the practical limits of the converter systems.
Thus, HAC combines the aforementioned complementary advantageous. Furthermore,
we show how retuning certain control parameters of the grid-forming controls
improves the frequency performance. Last, as separate contributions, we
introduce an alternative control augmentation that enhances the robustness and
provide theoretical guidelines on extending the stability certificates of HAC
to multi-converter systems.

##  I Introduction

The ambitious targets that are set to globally reduce the carbon footprint
require revolutionizing the foundations of legacy power systems. In other
words, the synchronous machines (SM)-based energy generation from fossil fuels
must be replaced with green and sustainable energy sources. The majority of
clean energy sources interface the power grids via dc-ac power converters.
Although, in contrast to the bulk SMs, power converters are fast, modular, and
highly controllable, they are subject to volatile energy resources and lack
the necessary robustness and maturity to ensure adequate and reliable power
delivery. It is envisioned that the advanced control architectures for power
converters can possibly address the aforementioned concerns [1, 2, 3, 4].

The so-called _grid-following_ converter controls that exploit an explicit
synchronizing mechanism are widely utilized [5, 6, 7]. However, the grid-
following converters exhibit robustness and stability issues in the converter-
dominated grids that are highlighted by a significant reduction of rotational
inertia i.e., _low-inertia_ grids; see [8, 9, 10] among others. Subsequently,
the concept of _ grid-forming converter (GFC)_ is introduced that provides
fast and robust frequency and voltage regulation to address the stability
challenges associated with the low-inertia systems [11, 12].

Several grid-forming control techniques have been proposed in recent years.
Restricting the focus to converter frequency definition under these controls,
one can highlight fundamental differences. The vast majority of grid-forming
techniques e.g., droop control, virtual synchronous machine (VSM), and
oscillator-based schemes shape the converter frequency based on the ac
quantities such as current, voltage, and power flows e.g., [13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23]. Recently, a class of controllers have been
proposed that define the converter frequency in proportion to a
linear/nonlinear dc voltage feedback [24, 25, 26, 27, 28, 29, 30]. Most
recently, the emerging multi-variable control trend is to combine dc and ac
information in designing the GFC frequency [31, 32, 33, 34, 35].

In this work, we focus on grid-forming hybrid angle control (HAC) that blends
linear dc and nonlinear ac feedback for defining converter frequency [31, 32].
We provide details on the implementation of HAC, highlight its performance in
contrast to other grid-forming controls via various simulation case studies,
and propose a complementary control augmentation for enhancing the robustness.
In addition, a high-level stability analysis is presented that provides
guidelines on extending the stability certificates of HAC (as in [32]) to an
interconnected system of converters. Finally, as a separate contribution in
its own right, we show that an appropriate tuning of converter control
parameters can globally enhance frequency stability across the low-inertia
power system.

The remainder of this paper is structured as it follows: Section II presents
our modeling approach for low-inertia power grids, Section III describes three
different grid-forming controls and highlights their dominant features.
Furthermore, it provides practical details on the implementation of HAC.
Section IV presents several comparative case studies. Last, the Appendix
introduces a variant of HAC with enhanced robustness and includes a high-level
stability investigation for multi-converter systems under HAC.

Figure 1: DC-AC power converter model schematic with ac quantities represented
in stationary $\alpha\beta$-coordinates [11, 36].

##  II Modeling Description

In this paper, we consider a transmission power grid that is fed by a mixture
of SM and converter-based generation units. For the sake of completeness, we
briefly review the SM, network, and converter modeling; the reader is referred
to [11] for further elaborations.

###  II-A Synchronous Machine

In this work, we consider a detailed nonlinear SM model that incorporates
second-order mechanical dynamics (associated with angle and frequency) and
sixth-order electrical dynamics (associated with the stator, damper, and field
winding fluxes). Furthermore, the SM is combined with a first-order dynamic
turbine model. The reader is referred to [11, Section II-B][37] for a detailed
description of the SM modeling. It is noteworthy that we include the standard
SM control mechanisms such as automatic voltage regulator (AVR), power system
stabilizer (PSS), and turbine governor control; see [11] for details.
Furthermore, the governor gain in [11, Equation 5] is exploited to realize
coordinated load-sharing with other generation units across the system; see
[11, Appendix] for a tuning criteria recommendation.

###  II-B Power Network

The generation units are considered to interface the transmission lines via
identical medium-to-high voltage (MV/HV) transformers. We consider dynamic
transmission lines represented with standard $\pi$-section models [37, 11].
Furthermore, the network loads are modeled by constant impedances. We
emphasize that the load models are voltage-dependent and restrict our focus to
the active power loads that primarily influence the frequency dynamics in
transmission grids. Note that due to the fast timescales of the GFCs and
potential adverse interactions with line dynamics, the quasi-steady-state
algebraic network model assumption is not valid [19, Subsection III-A], hence
it is necessary to consider dynamic line models for low-inertia power systems
[8, 38].

###  II-C High-Power DC-AC Converter

The power converter model in this work includes an explicit first-order
representation of the dc energy source dynamics, dc-link capacitance, power-
preserving (i.e., lossless) average model of the two-level dc-ac conversion
stage, and the LC output filter. Furthermore, we consider the current limits
of the dc energy source that in practice corresponds to the current limits of
the PV, battery, or wind generator systems. It is worth mentioning that in
this work the dc source model represents an aggregation of several energy
sources that supply a modular high-power converter; see Section IV for further
details. Figure 1 present the schematic of the converter model; see [11,
Subsection II-A] for further details and the differential equations
corresponding to the Figure 1.

The converter model illustrated in Figure 1 provides three control degrees of
freedom, namely: 1) the dc source current reference $i_{\text{{dc,r}}}$, 2)
the modulation signal magnitude i.e., $\mu=\lVert
m_{\alpha\beta}(\mu,\theta_{\text{{c}}})\rVert$, and the modulation signal
angle i.e., $\theta_{\text{c}}=\angle
m_{\alpha\beta}(\mu,\theta_{\text{{c}}})$. In Section III, we elaborate on the
definition of these control inputs under three different grid-forming control
strategies.

##  III Grid-Forming Frequency Controls

The vast literature on grid-forming control schemes can be distinctly assigned
into three categories:

  1. 1.

ac-based control techniques that define the converter frequency based on ac
measurements e.g., see [13, 14, 15, 16, 17, 19, 20, 18, 21, 22, 23],

  2. 2.

dc-based model-matching inspired techniques that relate the converter
frequency to the dc measurements e.g., see [24, 25, 26, 27, 28, 29, 30], and

  3. 3.

multi-variable hybrid control structures that exploit both dc and ac
measurements for a grid-forming frequency synthesis e.g., see [31, 32, 33, 34,
35].

In what follows, we review a candidate from each control category and
highlight their dominant features.

###  III-A AC-Based Power-Frequency Droop Control

The baseline power-frequency droop control that is inspired by the frequency
droop behavior of the SM, serves as a powerful yet simple control solution for
grid-forming converter applications [18, 21, 22, 23]. The power-frequency
droop control in its simplest form is described by

$\dot{\theta}_{\text{{c}}}=\omega_{\text{{c}}}:=\omega_{0}+d_{p-\omega}(p_{\text{{r}}}-p),$
(1)

where $\theta_{\text{{c}}}$ denotes the phase angle of the converter
modulation signal $m_{\alpha\beta}(\mu,\theta_{\text{{c}}})$ in Figure 1,
$\omega_{\text{{c}}}$ denotes the converter angular frequency, $\omega_{0}$
denotes the nominal system frequency, $d_{p-\omega}$ is the power-frequency
droop gain, $p_{\text{{r}}}$ is the control reference for the active power
flowing out of the converter’s ac-side terminal that is defined by
$p:=v^{\top}_{\alpha\beta}i_{\alpha\beta}$. The ac-based nature of the droop
control can be observed in (1) that highlights the frequency dependency on the
ac active power feedback. It is worth mentioning that low-pass filtering the
active power feedback in (1) is a common practice to safeguard the converter
control against measurement imperfections and high-frequency harmonics [39,
Introduction]. Section IV presents a case study that unveils the strong
influence of such low-pass filtering in reshaping the system-level post-
contingency frequency evolution.

######  Remark 1.

(Enhanced small-signal frequency stability)
Recent explorations [11, 9, 7] uncover that the ac-based grid-forming control
strategies e.g., droop control, VSM, and dispatchable virtual oscillator
control (dVOC) clearly improve the small-signal frequency stability (in terms
of the rate of change of frequency (RoCoF) and nadir performance metrics) of
the low-inertia system compared to the all-SMs conventional system. This is
primarily underpinned by the fast response timescale of the GFCs that enables
fast frequency regulation. Furthermore, it has been observed that the ac-based
grid-forming schemes exhibit slightly better frequency response compared to
the dc-based counterparts e.g., the matching control (presented in the next
Subsection); see [11, Subsection IV-C] for details. This behavior is due to
the fact that following a network contingency/disturbance, the ac feedback
(e.g., $p$ in (1)) quickly reflects the grid conditions in the converter
frequency definition and enables fast frequency regulation.

###  III-B DC-Based Model-Matching Control

Inspired by the structural similarities of the SM and converter dynamical
models, a family of model-matching control techniques have been proposed [24,
25, 26, 27, 28, 29, 30]. These control structures establish strong duality
between the converter and SM dynamical structure; see [29] for a detailed
derivation. Matching control in its original form is described by [24]

$\dot{\theta}_{\text{c}}=\omega_{\text{{c}}}:=\eta v_{\text{{dc}}},$ (2)

where $\eta:=\omega_{0}/v_{\text{{dc,r}}}$ in which $v_{\text{{dc,r}}}$
denotes the reference for the dc voltage $v_{\text{{dc}}}$. Note that (2)
highlights the dc-based nature of the matching strategy. Similar to droop
control, it might be necessary to low-pass filter the dc feedback in (2) to
safeguard the angle dynamics against potential dc-link voltage ripples.

######  Remark 2.

(Enhanced robustness)
Recall that the SMs can be operated with either flexible or constant
mechanical input torque, and regardless of this degree of freedom, the SMs
achieve robust synchronization with the power grid (thanks to their inherent
self-synchronizing feature [40, 26, 24, 29]). Under the strong converter-SM
duality induced by matching control, the converter can be operated with either
constant or flexible input dc current that is injected by the dc energy source
in Figure 1. This property of the matching control comes strongly into the
picture when considering the safety current constraints of the dc energy
source; see the saturation function in Figure 1. Recent works [11, 41, 33]
formally and numerically demonstrate that the matching controls (unlike the
ac-based schemes) preserve the closed-loop stability under active dc current
constraint. This behavior can be perceived as a built-in mode-switching
feature that simultaneously changes the converter operation mode from voltage
(i.e., grid-forming) to current (i.e., grid-following) source when the dc
constraint is activated; e.g., see [11, Figure 15]. To our knowledge this is
the most dominant feature of the dc-based controls i.e., superior robustness.

###  III-C Hybrid Angle Control

The hybrid ac/dc grid-forming control architecture appears as a natural
extension and the promising solution for combining the aforementioned
complementary benefits of the ac/dc-based control schemes [31, 32, 33, 34];
see Remarks 1 and 2. In what follows, we restrict our focus to the recently
proposed HAC that blends linear dc voltage and nonlinear ac angle error terms
for defining the converter frequency [31, 32]. The grid-forming HAC is
described by

$\dot{\theta}_{\text{{c}}}=\omega_{\text{{c}}}:=\omega_{0}+\underbracket{\gamma_{\text{{dc}}}\left(v_{\text{{dc}}}-v_{\text{{dc,r}}}\right)}_{\text{{dc
model-
matching}}}-\underbracket{\gamma_{\text{{ac}}}\sin\left(\dfrac{\delta-\delta_{\text{{r}}}}{2}\right)}_{\text{{ac
angle synchronization}}},$ (3)

where $\gamma_{\text{{dc}}}$ and $\gamma_{\text{{ac}}}$ respectively denote
the dc and ac control gains, and $\delta:=\angle
v_{\text{s},\alpha\beta}-\angle
v_{\alpha\beta}=\theta_{\text{{c}}}-\theta_{v}$ denotes the phase angle
difference (i.e., the relative angle) of the voltage before the filter and the
output voltage in Figure 1, and $\delta_{\text{{r}}}$ denotes the control
reference for $\delta$. Note that the dc term in (3) is identical to the
matching control variant proposed in [25] and the nonlinear ac term in (3)
resembles the Kuramoto-like angle synchronizing term associated with the
classic droop control (1) under certain assumptions [32, Remark 3][21]. The
practical implementation of (3) will be explored in the next section. The
reader is referred to [31, 32] for details on the design and properties of the
HAC. However, the following presents a brief summary:

  1. 1.

HAC (3) incorporates an inherent dc-ac power-balancing behavior i.e., assuming
$\omega_{\text{{c}}}\to\omega_{0}$ then
$\gamma_{\text{{dc}}}\left(v_{\text{{dc}}}-v_{\text{{dc,r}}}\right)-\gamma_{\text{{ac}}}\sin\left({\delta-\delta_{\text{{r}}}}/{2}\right)\to
0$. For instance, this means that if $v_{\text{{dc}}}>v_{\text{{dc,r}}}$ then
$\delta>\delta_{\text{{r}}}$ that allows for increased power injection into
the ac-side and subsequently stabilizes the dc voltage.

  2. 2.

The dc gain $\gamma_{\text{{dc}}}$ predominantly reinforces the frequency
dependency on the dc dynamics. It is expected that a nonzero dc gain enables
HAC to exhibit the robustness of dc-based controls (see Remark 2 for details
and Section IV for a numerical justification of this hypothesis.)

  3. 3.

The ac gain $\gamma_{\text{{ac}}}$ strongly influences the timescale
associated with the ac-side power flows. Therefore, it resembles the influence
of the droop gain on the performance of the droop control.

###  III-D DC and AC Voltage Control Schemes

The grid-forming droop, matching and hybrid angle controls (1)-(3) define the
phase angle that enters the converter modulation signal
$m_{\alpha\beta}(\mu,\theta_{\text{c}})$ in Figure 1. It remains to close the
loop by assigning the remaining control inputs i.e., the dc energy source
reference current $i_{\text{{dc,r}}}$ and modulation signal magnitude $\mu$.
For the sake of fairness in the forthcoming comparative investigation, we
consider identical complementary dc and ac voltage controls (that respectively
define $i_{\text{{dc,r}}}$ and $\mu$) for the aforementioned grid-forming
strategies.

Concerning the dc voltage control, we adopt the scheme that is proposed in
[11] (see [29] for a similar approach) i.e.,

$\displaystyle{i}_{{\text{{dc,r}}}}:=\underbracket{\kappa_{{\text{{dc}}}}\left(v_{{\text{{dc,r}}}}-{v_{{\text{{dc}}}}}\right)}_{\text{proportional
control}}+\underbracket{\frac{p_{\text{{r}}}}{v_{{\text{{dc,r}}}}}+\left(G_{{\text{{dc}}}}{v_{{\text{{dc}}}}}+\frac{v_{{\text{{dc}}}}i_{\text{{x}}}-{{p}}}{v_{{\text{{dc,r}}}}}\right)}_{\text{power
injection and loss feedforward}}$ (4)

where $k_{\text{{dc}}}$ denotes the proportional dc voltage control gain and
$i_{\text{{x}}}:=m^{\top}_{\alpha\beta}i_{\alpha\beta}$ denotes the net dc
current injection to the ac-side; see Figure 1. Note that the power injection
and loss compensation terms111The compensation scheme in (4) is conventionally
known as feedforward control, however due to the presence of the state-
dependent quantities (such as $v_{\text{{dc}}}$, $i_{\text{{x}}}$, and $p$) it
technically represents an algebraic state feedback. in (4) are not necessary
but improve the dynamic response and power set-point tracking [11, 29].

Next, inspired by the AVR mechanism of the SMs, we augment the grid-forming
strategies with a proportional-integral (PI) ac voltage magnitude control [11,
42] that is (with slight abuse of the notation) given by

$\mu:=\kappa_{\text{p}}\left(v_{\text{{r}}}-\lVert
v_{\alpha\beta}(t)\rVert\right)+\kappa_{\text{i}}\int_{0}^{t}\left(v_{\text{{r}}}-\lVert
v_{\alpha\beta}(s)\rVert\right)\text{d}s,$ (5)

where $\kappa_{\text{{p}}}$ and $\kappa_{\text{{i}}}$ respectively denote the
proportional and integral ac voltage control gains, and $v_{\text{{r}}}$
denotes the reference ac voltage amplitude. We remark that the robust ac
voltage regulation under (5) is favorable for a GFC (particularly in
islanded/microgrid applications). However, one can alternatively consider a
classic reactive power-voltage droop control as in [8]. Subsequently, the
converter modulation signal in Figure 1 is defined based the angle in (1)-(3)
and the magnitude prescribed by (5). Last, in case studies presented in the
Section IV we identically tune the dc and ac voltage controls (4) and (5) for
all grid-forming strategies.

###  III-E Hybrid Angle Control Implementation

In this subsection, we elaborate on the implementation of the grid-forming
HAC. Before embarking upon these discussions, we remark that the
implementation of droop control (1), matching control (2), and voltage
controls (4) and (5) is previously addressed; the reader is referred to [11,
43, 29, 18, 25] for details.

Concerning the implementation of HAC (3), previous works [31, 32] establish
the theoretical foundations. However, the implementation presented in [32, 31]
is not straightforward, thus in what follows, we make the reported ideas
clear.

To begin with, the dc component of the HAC (3) is easily constructed based on
the dc voltage measurement in Figure 1; see the the dc feedback control in
Figure 2. Prior to describing the implementation of the angle synchronizing
term in (3), note that $\delta=\theta_{\text{{c}}}-\theta_{v}$, thus one has
to implicitly/explicitly derive the angle information from the output ac
voltage measurement. Assume that the ac voltage control (5) is sufficiently
fast such that $\lVert v_{\alpha\beta}\rVert\approx v_{\text{{r}}}$ (i.e.,
$\lVert v_{\text{{abc}}}\rVert\approx v_{\text{{r}}}$ that follows from the
magnitude-preserving Clarke transformation [36, 32]).

As it is illustrated in Figure 2,
$\bar{v}_{\text{{abc}}}:=v_{\text{{abc}}}/v_{\text{{r}}}$ denotes the
normalized three-phase ac output voltage (note that $\bar{v}_{\text{{abc}}}$
represents a unity phasor in polar coordinates that rotates with angle
$\theta_{v}$). Next, we transform this quantity to a dq-coordinates that is
aligned with the converter modulation angle $\theta_{\text{{c}}}$. Thus, the
image of $\bar{v}_{\text{{abc}}}$ in the dq-rectangular coordinates (i.e., the
unit vector rotating with the angle
$\theta_{v}-\theta_{\text{{c}}}{=-\delta}$) implicitly contains the relative
angle information i.e.,

${\left(\bar{v}_{\text{{d}}},\bar{v}_{\text{{q}}}\right)=\left(\cos(\theta_{v}-\theta_{\text{{c}}}),\sin(\theta_{v}-\theta_{\text{{c}}})\right)=\left(\cos\delta,-\sin\delta\right).}$

Similar to low-pass filtering the active power feedback for the droop control
(1), we apply a first-order low-pass filter to
$\left(\cos\delta,-\sin\delta\right)$ where $\omega_{\text{{f}}}$ denotes the
cutoff frequency. Subsequently, the filter output represents an approximation
of the $\left(\cos\delta,-\sin\delta\right)$ that is
$\left(\widetilde{\cos\delta},-\widetilde{\sin\delta}\right)$222Note that the
LPF can by alternatively applied to $v_{\text{{abc}}}$ in Figure 2.. On the
other hand, we process a given relative angle reference $\delta_{\text{{r}}}$
that is consistent with desired power flows by the trigonometric functions
that results in
$\left(\cos\delta_{\text{{r}}},-\sin\delta_{\text{{r}}}\right)$ (see [32]
which shows how $\delta_{\text{{r}}}$ relates to the power and voltage set-
points). Finally, by exploiting the angle difference trigonometric identities
and the half-angle sine formula [32, Lemma 1 and Proposition 5] an
approximation of the angle term in (3) is obtained

$-\sin\left(\dfrac{\delta-\delta_{\text{{r}}}}{2}\right)\approx\dfrac{\sin\delta_{\text{r}}\widetilde{\cos\delta}-\widetilde{\sin\delta}\cos\delta_{\text{r}}}{\sqrt{2\left(1+\widetilde{\cos\delta}\cos\delta_{\text{r}}+\widetilde{\sin\delta}\sin\delta_{\text{r}}\right)}}.$
(6)

Respective multiplication of (6) and the dc voltage error with the gains
$\gamma_{\text{{ac}}}$ and $\gamma_{\text{{dc}}}$ provides the necessary
ingredients for implementing HAC (3). The reader is referred to [32, Section
V] for further discussions.

Figure 2: Block diagram of the HAC (3) based on the dc voltage and the output
ac voltage measurements in Figure 1.

##  IV Numerical Case Studies

In this section, we employ the standard IEEE 9-bus test system as described in
[8, 11, 43]; see Figure 3 and Table I for the model and control parameters.
This system represents a transmission grid in which the SMs are rated at
$100~{}[\mathrm{MVA}]$. Thus, in order to study the system-level influence of
the GFCs, it is ideally desired to consider a roughly similar rating for the
GFCs as that of the SMs. Therefore, we employ the dynamic aggregation
technique proposed in [44] that provides appropriate scaling laws for deriving
the aggregated parameters of a high-power converter system based on the
parameters of a smaller module.To this end, we consider a GFC module rated at
$500~{}[\mathrm{kVA}]$ that corresponds to a commercially available system;
see [11, Table 1] for the parameters. Next, we envision a parallel connection
of $200\times 500~{}[\mathrm{kVA}]$ units that is represented with a
$100~{}[\mathrm{MVA}]$ aggregated model with the same structure as in Figure
1. The reader is referred to [11, Figure 2 and Remark 1][9, 44] for further
details and similar model aggregation approaches. In what follows, the system-
level performance of the grid-forming hybrid angle control is explored in
various case studies. Note that the employed simulation model in
MATLAB/Simulink environment is publicly available [43].

Table I: Case study model and control parameters [43].

IEEE 9-bus test system base values
---
$S_{\text{b}}$ |  $100$ MVA | $v_{\text{b}}$ |  $230$ kV | $\omega_{\text{b}}$ |  $50$ Hz
synchronous machine (SM)
$S_{\text{r}}$ |  $100$ MVA | $v_{\text{r}}$ |  $13.8$ kV | $D$ | $0$
$H$ |  $3.7$ s | $d_{p}$ | $1\%$ | $\tau_{\text{g}}$ |  $5$ s
single converter module
$S_{\text{r}}$ |  $500$ kVA | $G_{{\text{{dc}}}},C_{{\text{{dc}}}}$ |  $0.83,0.008$ $\Omega^{-1}$,F | $v_{{\text{{dc,r}}}},v_{\text{ll-rms,r}}$ |  $2.44,1$ kV
$R$ |  $0.001$ $\Omega$ | $L$ |  $200$ $\mu$H | $C$ |  $300$ $\mu$F
$n$ | $200$ | $\tau_{{\text{{dc}}}}$ |  $50$ ms | $i^{\text{dc}}_{\max}$ |  $1.2$ pu
dc and ac voltage controls
$k_{{\text{{dc}}}}$ | $1.6\times 10^{3}$ | $k_{\text{p}}$ | $0.001$ | $k_{\text{i}}$ | $0.5$
droop control, matching control, and HAC
$\omega_{0}$ | $\omega_{\text{{b}}}$ | $d_{p-\omega}$ | $1\%$ | $\eta$ | $\omega_{0}/v_{\text{{dc,r}}}$
$\gamma_{\text{{dc}}}$ | $0.01\eta$ | $\gamma_{\text{{ac}}}$ | $205$ | $\delta_{\text{{r}}}$ | $0.0238$

Figure 3: IEEE 9-bus test system including SMs and GFCs; depending on the
choice of the generation technology at nodes 1, 2, and 3 this network model
represents: 1) a conventional system (i.e., all-SMs), 2) a SM-dominated system
(2SMs-1GFC), 3) a low-inertia system (i.e., 1SM-2GFCs), and 4) a purely
converter-based system (i.e., all-GFCs).

###  IV-A Purely Converter-Based System

In this scenario, we remove the SMs from the grid model shown in Figure 3 and
consider three identical GFCs at nodes 1-3. We implement the HAC strategy (as
in Figure 2) for all the GFCs and identically tune the controllers (resulting
in equal load-sharing). Figure 6 illustrates the frequency and active power
evolution at the generation nodes following a $0.75~{}[\mathrm{pu}]$ load-
disturbance at node 7. Note that we do not implement an explicit frequency
measurement mechanism and rather observe the internal frequency of the GFCs
i.e., (3).

Figure 6 firstly verifies the grid-forming nature of HAC technique i.e., the
autonomous operation of converters without relying on a reference frequency
provided by an external source (as in the case of grid-following devices).
Second, note that all GFCs reach the post-contingency equilibrium in
approximately $200~{}[\mathrm{ms}]$. This fast response timescale is crucial
for a satisfactory grid-forming performance [1]. Observe that the units that
are electrically closer to the disturbance location i.e., GFCs at node 2 and 3
react on a slightly faster timescale. Last, Figure 6 confirms the equal load-
sharing of the GFCs. See [11, 9] for similar observations concerning the
behavior of other grid-forming schemes in an all-GFCs network.

Figure 4: Frequency response of the all-GFCs IEEE 9-bus system configuration
under HAC strategy following a load disturbance (top), the active power time-
evolution associated with the GFCs at node 1, 2, and 3 (bottom).

Figure 5: RoCoF evolution of the SM at node 1 in 1SM-2GFCs system
configuration under different controls and with respect to variations in
network load disturbance (top), maximum frequency deviation (bottom).

Figure 6: The post-event frequency time-evolution of the SM at node 1 in
2SMs-1GFC system configuration under different controls for the GFC.

###  IV-B The Influence on the Frequency Performance Metrics

In what follows, we explore the influence of grid-forming architectures
(presented in Section III) on the frequency stability of low-inertia
configuration associated with the grid model shown in Figure 3 i.e., with SM
at node 1 and GFCs at nodes 2 and 3. More precisely, we 1) implement identical
droop controllers for both GFCs, 2) implement matching control for the GFCs,
and finally 3) consider HAC for the converters. Subsequently, for all the
1SM-2GFCs pairs under different controls, we apply five load-disturbances at
note 7 i.e., $\Delta p_{\ell,k}=15k~{}[\mathrm{MVA}]$ where $k=1,\ldots,5$.
Next, we observe the frequency dynamics via the angular frequency of the SM at
node 1; see [11, 9] for similar approaches. We evaluate the standard frequency
performance metrics, namely: 1) maximum frequency deviation, and 2) RoCoF that
are (with slight abuse of notation) defined by [45, 11]

$\displaystyle\lVert\Delta\omega_{\text{{SM}}}\rVert_{\infty}$
$\displaystyle:=\max\left(\omega_{0}-\omega_{\text{{SM}}}(t)\right),$ (7a)
$\displaystyle\big{|}\hat{{\dot{\omega}}}_{\text{{SM}}}\big{|}$
$\displaystyle\coloneqq\bigg{|}\frac{\omega_{\text{{SM}}}(t_{0}+\Delta
t)-\omega_{\text{{SM}}}(t_{0})}{\Delta t}\bigg{|},$ (7b)

where $\omega_{\text{{SM}}}$ and $\Delta t$ respectively denote the SM
frequency and the RoCoF approximation time horizon. We consider $\Delta
t=150~{}[\mathrm{ms}]$ that allows to observe the influence of fast GFCs
dynamics on the system frequency; the reader is referred to [9] for a detailed
discussion on the choice of RoCoF window.

Figure 6 highlights the evolution of metrics (7) in 1SM-2GFCs system, under
different grid-forming controls, and with respect to the disturbance
variation. Concerning the RoCoF, the techniques incorporating an ac feedback
i.e., droop control (1) and HAC (3) result in better performance compared to
the purely dc-based matching control. This is underpinned by their fast
disturbance sensing and frequency modification feature. In contrast, matching
control (2) only reacts when the network disturbance is propagated to the
converter dc dynamics.

Figure 6 confirms that HAC inherits the advantage of the other ac-based
controls e.g., droop control; see Remark 1 and [11] for further details. On
the other hand, maximum frequency deviation performance is almost identical
when employing different controls. The reader is referred to [11] for a
detailed discussion on this aspect. Although Figure 6 uncovers a performance
variation that depends on the grid-forming controls in 1SM-2GFCs setup, all
controls outperform the all-SMs configuration by a clear margin; e.g., see
[11, Figures 11 and 12]. The positive influence of the grid-forming controls
on frequency stability is also reported in [9, 8].

###  IV-C Synchronous Machine Dominated System

In this test case, we consider a SM-dominated generation profile by including
identical SMs at nodes 1 and 3, and a GFC at node 2 in Figure 3 leading to a
2SMs-1GFC system. We consider the same load-disturbance as in Subsection IV-A
and enforce equal load-sharing for all units. Figure 6 shows the frequency of
the SM at node 1 when the GFC is under different controls. The network base
load and disturbance are chosen such that the desired post-disturbance GFC’s
active power does not result in exceeding the converter dc source current
limit; see Figure 1. However, due to the presence of slowly reacting SMs
(because of large turbine time constants), the GFC dominantly supports the
load during the transient that results in violating the dc limit. As it is
also reported in [41, 33, 11], under active dc current constraint, droop
control exhibits instability by depleting the dc-link energy. This aggressive
behavior is because droop control (1) is agnostic to the dc dynamics.
Subsequently, the converter instability propagates through the network an
destabilizes the SMs. This adverse interaction can be counteracted by
augmenting droop control with an ac current limiting mechanism e.g., threshold
virtual impedance [46]. However, in our experience the limiting performance is
fragile and depends on various factors e.g., the disturbance size and
location.

In contrast, as it can be observed in Figure 6, the controls that exploit dc
feedback in their angle dynamics e.g., matching control and HAC exhibit
robustness with respect to the dc current constraint. This behavior is
numerically and formally explored in [11, 41, 33]. In a nutshell, it is
underpinned by an inherent mode-switching behavior that transforms the grid-
forming control to a following one and allows injecting constant current while
dc source is saturated. The reader is referred to [11] for a detailed
discussion. Interestingly, HAC results in better frequency deviation
performance compared to matching control. This is because HAC provides global
stability certificates [32, 31] and includes a hybrid ac/dc structure that
strikes a balance between robustness and performance.

######  Remark 3.

(Adverse timescales interaction)
In case studies presented in [11, 41, 33], the instability of the purely ac-
based grid-forming techniques strongly depends on the disturbance magnitude.
However, Figure 6 provides an alternative insight. More precisely, the base
load and disturbance scenario specification does not derive the converter
beyond its limit in all-GFCs or 1SM-2GFCs systems. However, in SM-dominated
system the adverse interplay of the fast GFC and slow SMs results in excessive
transient current injection by the converter. Hence, aside the disturbance
characteristics, the penetration level of the converter-based generation is an
influencing factor for this instability mechanism.

###  IV-D The Influence of GFC parameters on Frequency Behavior

Several works e.g., [47, 48, 39] highlight the necessity and benefits of low-
pass filtering the ac measurements that are exploited in synthesizing
converter control. In the following, we show how an appropriate tuning of the
converters’ low-pass filters improves the frequency stability and can be
perceived as a decentralized frequency shaping strategy.

To begin with, we consider the 1SM-2GFCs configuration as in Subsection IV-B.
Next, we equip the GFCs at nodes 2 and 3 with HAC and consider different
cutoff frequencies for the underlying low-pass filters; see Figure 2. Figure 8
illustrates the frequency response of the SM with respect to the converters’
cutoff frequency variation. Note that this observation suggests that an
appropriate tuning of the LPFs can globally reshape the frequency response
across the system. This is due to the fact that the LPFs with an appropriate
tuning provide a certain amount of virtual inertia and thus enhance the SM
frequency response. The reader is referred to [48] for detailed discussion. We
remark that the positive influence of $\omega_{\text{{f}}}$-variation on
frequency stability must be compared against the enhancement that is achieved
via changing the converters’ controls. In particular, the maximum frequency
deviation metric is significantly reduced in Figure 8 as $\omega_{\text{{f}}}$
is reduced, in contrast, for a fixed $\omega_{\text{{f}}}$ the improvement due
to changing grid-forming strategies is minimal; see Figure 6. Last, our
numerical investigations confirms an almost identical behavior as when the
GFCs are controlled by droop control.

Figure 7: Frequency evolution of SM at node 1 when the GFCs are controlled by
HAC strategy with 5 different LPF cutoff frequency and without a LPF.

Figure 8: The PSS-free response of the SM at node 1 when the GFCs are
controlled by either HAC or droop control while considering two different LPF
cutoff frequencies. Figure 9: The simplified two-converter system
configuration where $i_{\text{{x}},j}:=m(\mu_{j},\theta_{j})^{\top}i_{\ell}$
and $v_{\text{{x}},j}:=m(\mu_{j},\theta_{j})v_{\text{{dc}},j}$ for $j=1$ and
$2$.

###  IV-E The PSS-Free Behavior of the Low-Inertia System

It has been recently uncovered that the PSSs of the SMs in a low-inertia
system might adversely interact with certain timescales of the converter
systems [8]. Moreover, it is highlighted that an adaptive PSS might be a
crucial element in ensuring system stability in transitioning to 100$\%$
converter-based generation [49, 9, 50]. However, the online/offline
modification of the well-established PSS control architectures might be an
expensive, challenging, and practically infeasible task. In this subsection,
we investigate if the LPF retuning strategy can possibly allow for a PSS-free
operation of the SMs in a low-inertia network configuration.

To this end, we consider the 1SM-2GFCs configuration and remove the PSS form
the SM model at node 1. Similar to previous cases, we separately consider
droop control and HAC for the GFCs. Next, we consider two different cutoff
frequencies for the underlying LPFs. Figure 8 depicts the frequency response
of the SM. It can be seen that under $\omega_{\text{{f}}}=5~{}[\mathrm{Hz}]$,
HAC results in slightly less oscillations compared to droop control. However,
when selecting $\omega_{\text{{f}}}=0.2~{}[\mathrm{Hz}]$, the SM frequency
nadir (that might potentially trip the low-frequency protection mechanisms) is
almost removed. Although, under $\omega_{\text{{f}}}=0.2~{}[\mathrm{Hz}]$ the
low-frequency oscillation are still present but the appropriate LPF tuning
reshapes the envelope on the frequency response and constrains the frequency
oscillations within the stability margin. Last, the enhanced response is
achieved regardless of the control strategies of the GFCs. We remark that this
insight is a result of a preliminary exploration and requires an in-depth
analysis.

##  V Conclusions and Outlook

In this paper, the system-level performance and robustness of grid-forming HAC
is explored. We verified that the multi-variable dc-ac HAC inherits the
enhanced performance and robustness of ac-based and dc-based controls,
respectively. Furthermore, we highlighted that how retuning the low-pass
filters for grid-forming controls can significantly enhance the frequency
performance across the system. Last, an alternative augmented controller is
introduced and a simplified stability analysis for a two-converter system is
presented. Our future work includes: 1) implementation of HAC for power
converters in more complex power grid models, 2) power hardware validation of
HAC concept, and 3) detailed stability analysis multi-converter system under
HAC.

###  -A Dynamic Inverse Droop Control Augmentation

The HAC implementation described in the Section III, explicitly relies on the
foreknowledge of the relative angle reference $\delta_{\text{{r}}}$. The
reader is referred to [32, Proposition 6] on deriving $\delta_{\text{{r}}}$
based on the power and voltage set-points. However, this techniques relies on
the system parameters and might exhibit robustness issues. In what follows, we
introduce a complementary high-level feedback control mechanism that achieves
the same objective with enhanced robustness.

To begin with we introduce an integral control that defines
$\delta_{\text{{r}}}$ in relation to the active power mismatch i.e.,

$\delta_{\text{{r}}}=\int_{0}^{t}\kappa_{p\delta}(p_{\text{{r}}}-p),$ (8)

where $\kappa_{p\delta}$ denotes the integrator gain. It is worth mentioning
that augmenting the HAC with (8) transforms the grid-forming HAC to a phase
locked loop (PLL)-free grid-following controller that achieves robust active
power reference tracking. This is due to the fact that the integrator (8)
disables the natural $p-\omega$ droop behavior of the HAC [32, Proposition 7].

Although robust reference tracking and disturbance rejection might be
desirable for certain applications, drooping behavior that enables load-
sharing between generation units in transmission grids is vital. Thus, we
combine (8) with an inverse $\omega-p$ droop control (reminiscent of the speed
droop control of the SM e.g., see [11, Equation 5]) that is

$p_{\text{{r}}}:=p^{\star}+d_{\omega-p}(\omega_{\text{{c}}}-\omega_{0}),$ (9)

where $p^{\star}$ denotes the power reference at nominal frequency,
$d_{\omega-p}$ is the inverse droop gain, and $\omega_{\text{{c}}}$ is the
internal feedback given by HAC (3). We emphasize that the combination of (3)
with the dynamic inverse droop control (8) and (9) is a heuristic approach,
and requires a separate detailed stability analysis as in [32]. However, our
numerical investigations confirms that the cascade control structure achieves
the aforementioned control specification. Furthermore, as in standard cascaded
control systems (e.g., [51]) the controls (3), (8), and (9) must be tuned in
harmony while respecting the required timescales separation between the nested
loops. Last, we remark that since the main focus of this work is to explore
the behavior of the standard HAC (3), the presented case studies in the
Section IV only incorporates (3).

###  -B Stability Analysis of the Interconnected Converters

The stability analysis of the HAC in [32, 31] is centered around a model
configuration that includes a single converter connected to an infinite bus or
a dynamic center of inertia (CoI) grid model. The converter-CoI model
abstractly represents the mixed converter-SM configurations of the IEEE 9-bus
system shown in Figure 3. However, the stability and synchronization analysis
of a multi-converter setup requires a separate investigation and is presented
as it follows.

To begin with, we consider a simplified two-converter model as in Figure 9;
see [30] for a similar model configuration. Without loss of generality we
assume the line current is flowing out of the converter 1 and into the
converter 2. It is noteworthy that the configuration in Figure 9 represents
two converters with RL output filters that are connected through a RL line
section. Thus, the filter elements and the line section are merged into a
single RL elements.

We equip both converters with HAC (3) and thus the overall dynamical model is
given by

$\displaystyle\dot{\theta}_{1}=\omega_{0}+\gamma_{\text{{dc}},1}(v_{\text{{dc}},1}-v_{\text{{dc,r}}1})-\gamma_{\text{{ac}},1}\sin\left(\dfrac{\delta-\delta_{r}}{2}\right),$
(10a)
$\displaystyle\dot{v}_{\text{{dc}},1}=C_{\text{{dc}},1}^{-1}\left(i_{\text{{dc}},1}-G_{\text{{dc}},1}v_{\text{{dc}},1}-m(\mu_{1},\theta_{1})^{\top}i_{\ell}\right),$
(10b)
$\displaystyle\dot{\theta}_{2}=\omega_{0}+\gamma_{\text{{dc}},2}(v_{\text{{dc}},2}-v_{\text{{dc,r}}2})+\gamma_{\text{{ac}},2}\sin\left(\dfrac{\delta-\delta_{r}}{2}\right),$
(10c)
$\displaystyle\dot{v}_{\text{{dc}},2}=C_{\text{{dc}},2}^{-1}\left(i_{\text{{dc}},2}-G_{\text{{dc}},2}v_{\text{{dc}},2}+m(\mu_{2},\theta_{2})^{\top}i_{\ell}\right),$
(10d)

where $\delta:=\theta_{1}-\theta_{2}$. In addition, the dc source current and
the modulation signal control input pairs are defined as

$\displaystyle i_{\text{{dc}},1}$
$\displaystyle:=-\kappa_{\text{{dc}},1}\left(v_{\text{{dc}},1}-v_{\text{{dc,r}}1}\right),$
(11a) $\displaystyle m(\mu_{1},\theta_{1})$
$\displaystyle:=\mu_{1}\left(\cos\theta_{1},\sin\theta_{1}\right)^{\top},$
(11b) $\displaystyle i_{\text{{dc}},2}$
$\displaystyle:=-\kappa_{\text{{dc}},2}\left(v_{\text{{dc}},2}-v_{\text{{dc,r}}2}\right),$
(11c) $\displaystyle m(\mu_{2},\theta_{2})$
$\displaystyle:=\mu_{2}\left(\cos\theta_{2},\sin\theta_{2}\right)^{\top}.$
(11d)

Furthermore, for simplicity of exposition we assume quasi-steady-state line
dynamics, thus the current in stationary $\alpha\beta$-coordinates is given by

$i_{\ell}:=R^{-1}\left(v_{\text{{dc}},1}m(\mu_{1},\theta_{1})-v_{\text{{dc}},2}m(\mu_{2},\theta_{2})\right).$
(12)

We subsequently 1) combine the absolute angle dynamics (10a) and (10c) into
relative angle dynamics, 2) replace $i_{\ell}$ in (10b) and (10d) with the
expression in (12), and 3) incorporate the controls (11) that results in the
overall closed-loop dynamics:

$\displaystyle\dot{v}_{\text{{dc}},1}=$ $\displaystyle-
C_{\text{{dc}},1}^{-1}\Big{(}\kappa_{\text{{dc}},1}\left(v_{\text{{dc}},1}-v_{\text{{dc,r}}1}\right)+\left(G_{\text{{dc,1}}}+R^{-1}\mu_{1}^{2}\right)v_{\text{{dc}},1}$
$\displaystyle-R^{-1}\mu_{1}\mu_{2}\cos(\delta)v_{\text{{dc}},2}\Big{)},$
(13a) $\displaystyle\dot{v}_{\text{{dc}},2}=$ $\displaystyle-
C_{\text{{dc}},2}^{-1}\Big{(}\kappa_{\text{{dc}},2}\left(v_{\text{{dc}},2}-v_{\text{{dc,r}}2}\right)+\left(G_{\text{{dc}},2}+R^{-1}\mu_{2}^{2}\right)v_{\text{{dc}},2}$
$\displaystyle-R^{-1}\mu_{1}\mu_{2}\cos(\delta)v_{\text{{dc}},1}\Big{)},$
(13b) $\displaystyle\dot{\delta}=$
$\displaystyle+\gamma_{\text{{dc}},1}(v_{\text{{dc}},1}-v_{\text{{dc,r}}1})-\gamma_{\text{{dc}},2}(v_{\text{{dc}},2}-v_{\text{{dc,r}}2})$
$\displaystyle-(\gamma_{\text{{ac}},1}+\gamma_{\text{{ac}},2})\sin\left(\dfrac{\delta-\delta_{r}}{2}\right).$
(13c)

Assume that there exists an equilibrium and it coincides with the desired
control references (similar to [32, Theorem 1]) i.e.,

$\left(\delta^{\star},v_{\text{{dc}},1}^{\star},v_{\text{{dc}},2}^{\star}\right):=\left(\delta_{\text{{r}}},v_{\text{{dc,r}}1},v_{\text{{dc,r}}1}\right).$
(14)

Next, we associate an energy function with the closed-loop dynamics (13), that
is given by

$\displaystyle\mathcal{V}:=$
$\displaystyle\frac{1}{2}\left(C_{\text{{dc}},1}\left(v_{\text{{dc}},1}-v_{\text{{dc,r}}1}\right)^{2}+C_{\text{{dc}},2}\left(v_{\text{{dc}},2}-v_{\text{{dc,r}}2}\right)^{2}\right)$
$\displaystyle+2\left(1-\cos\left(\dfrac{{\delta-\delta_{\text{{r}}}}}{2}\right)\right).$
(15)

Following the analysis recipe as in [32, Theorem 2], it is possible to show
that there exists positive lower bounds $\gamma_{\text{{ac}},\min}$ and
$\kappa_{\text{{dc}},\min}$ such that if
$\gamma_{\text{{ac}},1}+\gamma_{\text{{ac}},2}>\gamma_{\text{{ac}},\min}$ and
$\kappa_{\text{{dc}},1},\kappa_{\text{{dc}},2}>\kappa_{\text{{dc}},\min}$ then
the energy function $\mathcal{V}\to 0$ for (almost) all the solutions of (13).
Note that when $\mathcal{V}\to 0$ it yields that
$\left(\delta,v_{\text{{dc}},1},v_{\text{{dc}},2}\right)\to\left(\delta_{\text{{r}}},v_{\text{{dc,r}}1},v_{\text{{dc,r}}1}\right)$.
This result verifies the asymptotic convergence of converter dynamics to the
desired dc voltage and relative angle references.

## References

  * [1] F. Milano, F. Dörfler, G. Hug, D. J. Hill, and G. Verbič, “Foundations and challenges of low-inertia systems,” in _Power Systems Computation Conference (PSCC)_ , 2018.
  * [2] B. Kroposki, B. Johnson, Y. Zhang, V. Gevorgian, P. Denholm, B.-M. Hodge, and B. Hannegan, “Achieving a 100% renewable grid: Operating electric power systems with extremely high levels of variable renewable energy,” _IEEE Power and energy magazine_ , vol. 15, no. 2, pp. 61–73, 2017.
  * [3] R. Rosso, X. Wang, M. Liserre, X. Lu, and S. Engelken, “Grid-forming converters: Control approaches, grid-synchronization, and future trends—a review,” _IEEE Open Journal of Industry Applications_ , vol. 2, pp. 93–109, 2021.
  * [4] T. Ackermann, T. Prevost, V. Vittal, A. J. Roscoe, J. Matevosyan, and N. Miller, “Paving the way: A future without inertia is closer than you think,” _IEEE Power and Energy Magazine_ , vol. 15, no. 6, pp. 61–69, 2017\.
  * [5] Y. Li, Y. Gu, and T. C. Green, “Rethinking grid-forming and grid-following inverters: A duality theory,” 2021, Preprint available at https://arxiv.org/abs/2105.13094.
  * [6] D. Pattabiraman, R. H. Lasseter., and T. M. Jahns, “Comparison of grid following and grid forming control for a high inverter penetration power system,” in _IEEE Power Energy Society General Meeting (PESGM)_ , 2018.
  * [7] A. Tayyebi, F. Dörfler, F. Kupzog, Z. Miletic, and W. Hribernik, “Grid-forming converters – inevitability, control strategies and challenges in future grid applications,” in _Workshop on Microgrids and Local Energy Communities (CIRED)_ , 2018.
  * [8] U. Markovic, O. Stanojev, P. Aristidou, E. Vrettos, D. S. Callaway, and G. Hug, “Understanding small-signal stability of low-inertia systems,” _IEEE Trans. Power Syst._ , 2021.
  * [9] A. Crivellaro, A. Tayyebi, C. Gavriluta, D. Groß, A. Anta, F. Kupzog, and F. Dörfler, “Beyond low-inertia systems: massive integration of grid-forming power converters in transmission grids,” in _IEEE Power Energy Society General Meeting (PESGM)_ , 2020.
  * [10] D. Pattabiraman, R. H. Lasseter, and T. M. Jahns, “Impact of phase-locked loop control on the stability of a high inverter penetration power system,” in _IEEE Power Energy Society General Meeting (PESGM)_ , 2019.
  * [11] A. Tayyebi, D. Groß, A. Anta, F. Kupzog, and F. Dörfler, “Frequency stability of synchronous machines and grid-forming power converters,” _IEEE Trans. Emerg. Sel. Topics Power Electron._ , vol. 8, no. 2, pp. 1004–1018, 2020.
  * [12] P. Unruh, M. Nuschke, P. Strauß, and F. Welck, “Overview on grid-forming inverter control methods,” _Energies_ , vol. 13, no. 10, p. 2589, 2020.
  * [13] B. B. Johnson, M. Sinha, N. G. Ainsworth, F. Dörfler, and S. V. Dhople, “Synthesizing virtual oscillators to control islanded inverters,” _IEEE Trans. Power Electron._ , vol. 31, no. 8, pp. 6002–6015, 2015.
  * [14] M. Sinha, F. Dörfler, B. B. Johnson, and S. V. Dhople, “Uncovering droop control laws embedded within the nonlinear dynamics of van der pol oscillators,” _IEEE Trans. Control Netw. Syst._ , vol. 4, no. 2, pp. 347–358, 2017.
  * [15] M. Chen, D. Zhou, and F. Blaabjerg, “Active power oscillation damping based on acceleration control in paralleled virtual synchronous generators system,” _IEEE Transactions on Power Electronics_ , vol. 36, no. 8, pp. 9501–9510, 2021.
  * [16] Q. C. Zhong and G. Weiss, “Synchronverters: inverters that mimic synchronous generators,” _IEEE Trans. Ind. Electron._ , vol. 58, no. 4, pp. 1259–1267, 2011.
  * [17] P. Lorenzetti, Z. Kustanovich, S. Shivratri, and G. Weiss, “The equilibrium points and stability of grid-connected synchronverters,” 2021, Preprint available at https://arxiv.org/abs/2103.07952.
  * [18] M. Chandorkar, D. Divan, and R. Adapa, “Control of parallel connected inverters in standalone AC supply systems,” _IEEE Trans. Ind. Appl._ , vol. 29, no. 1, pp. 136–143, 1993.
  * [19] M. Colombino, D. Groß, J. Brouillon, and F. Dörfler, “Global phase and magnitude synchronization of coupled oscillators with application to the control of grid-forming power inverters,” _IEEE Trans. Autom. Control_ , vol. 64, no. 11, pp. 4496–4511, 2019.
  * [20] S. D’Arco, J. A. Suul, and O. B. Fosso, “A virtual synchronous machine implementation for distributed control of power converters in smart grids,” _Electric Power Systems Research_ , vol. 122, pp. 180–197, 2015.
  * [21] J. W. Simpson-Porco, F. Dörfler, and F. Bullo, “Synchronization and power sharing for droop-controlled inverters in islanded microgrids,” _Automatica_ , vol. 49, no. 9, pp. 2603–2611, 2013.
  * [22] K. Smith, S. Jafarpour, and F. Bullo, “Transient stability of droop-controlled inverter networks with operating constraints,” _IEEE Trans. Autom. Control_ , 2021.
  * [23] C. N. Rowe, T. J. Summers, R. E. Betz, D. J. Cornforth, and T. G. Moore, “Arctan power–frequency droop for improved microgrid stability,” _IEEE transactions on power electronics_ , vol. 28, no. 8, pp. 3747–3759, 2012.
  * [24] C. Arghir, T. Jouini, and F. Dörfler, “Grid-forming control for power converters based on matching of synchronous machines,” _Automatica_ , vol. 95, pp. 273–282, 2018.
  * [25] S. Curi, D. Groß, and F. Dörfler, “Control of low-inertia power grids: a model reduction approach,” in _IEEE Conference on Decision and Control (CDC)_ , 2017.
  * [26] I. Cvetkovic, D. Boroyevich, R. Burgos, C. Li, and P. Mattavelli, “Modeling and control of grid-connected voltage-source converters emulating isotropic and anisotropic synchronous machines,” in _IEEE Workshop on Control and Modeling for Power Electronics (COMPEL)_ , 2015.
  * [27] L. Huang, H. Xin, Z. Wang, K. Wu, H. Wang, J. Hu, and C. Lu, “A virtual synchronous control for voltage-source converters utilizing dynamics of dc-link capacitor to realize self-synchronization,” _IEEE Trans. Emerg. Sel. Topics Power Electron._ , vol. 5, no. 4, pp. 1565–1577, 2017.
  * [28] J. D. Watson and I. Lestas, “Control of interlinking converters in hybrid AC/DC grids: network stability and scalability,” _IEEE Trans. Power Syst._ , vol. 36, no. 1, pp. 769–780, 2020.
  * [29] C. Arghir and F. Dörfler, “The electronic realization of synchronous machines: Model matching, angle tracking, and energy shaping techniques,” _IEEE Trans. Power Electron._ , vol. 35, no. 4, pp. 4398–4410, 2020.
  * [30] T. Jouini and F. Dörfler, “Parametric local stability condition of a multi-converter system,” _IEEE Trans. Autom. Control_ , 2020.
  * [31] A. Tayyebi, A. Anta, and F. Dörfler, “Almost globally stable grid-forming hybrid angle control,” in _IEEE Conference on Decision and Control (CDC)_ , 2020.
  * [32] A. Tayyebi, A. Anta, and F. Dörfler, “Hybrid angle control and almost global stability of grid-forming power converters,” 2020, preprint.
  * [33] Y. Gao, H.-P. Ren, and J. Li, “Grid-forming converters control based on dc voltage feedback,” 2020, Preprint available at https://arxiv.org/abs/2009.05759.
  * [34] I. Subotić and D. Groß, “Power-balancing dual-port grid-forming power converter control for renewable integration and hybrid ac/dc power systems,” 2021, Preprint available at https://arxiv.org/abs/2106.10396.
  * [35] M. Chen, D. Zhou, A. Tayyebi, E. Prieto-Araujo, F. Dörfler, and F. Blaabjerg, “Generalized multivariable grid-forming control design for power converters,” 2021, Preprint available at https://arxiv.org/abs/2109.06982.
  * [36] A. Yazdani and R. Iravani, _Voltage-sourced converters in power systems: modeling, control, and applications_. John Wiley & Sons, 2010.
  * [37] P. Kundur, _Power system stability and control_. McGraw-hill, 1994.
  * [38] D. Groß, M. Colombino, J. Brouillon, and F. Dörfler, “The effect of transmission-line dynamics on grid-forming dispatchable virtual oscillator control,” _IEEE Trans. Control Netw. Syst._ , vol. 6, no. 3, pp. 1148–1160, 2019.
  * [39] B. Johnson, M. Rodriguez, M. Sinha, and S. Dhople, “Comparison of virtual oscillator and droop control,” in _IEEE Workshop on Control and Modeling for Power Electronics (COMPEL)_ , 2017.
  * [40] N. Barabanov, J. Schiffer, R. Ortega, and D. Efimov, “Conditions for almost global attractivity of a synchronous generator connected to an infinite bus,” _IEEE Trans. Autom. Control_ , vol. 62, pp. 4905–4916, 2017.
  * [41] S. Samanta and N. R. Chaudhuri, “On stability analysis of power grids with synchronous generators and grid-forming converters under dc-side current limitation,” 2021, Preprint available at https://arxiv.org/abs/2103.09966.
  * [42] A. Sajadi, R. W. Kenyon, M. Bossart, and B.-M. Hodge, “Dynamic interaction of grid-forming and grid-following inverters with synchronous generators in hybrid power plants,” in _2021 IEEE Kansas Power and Energy Conference (KPEC)_ , 2021.
  * [43] A.Tayyebi, D. Groß, and A. Anta, “GridFormingConverters: Implementation of grid-forming control techniques in IEEE 9-bus system,” Git repository, 2019, https://github.com/ATayebi/GridFormingConverters.
  * [44] V. Purba, B. B. Johnson, S. Jafarpour, F. Bullo, and S. V. Dhople, “Dynamic aggregation of grid-tied three-phase inverters,” _IEEE Trans. Power Syst._ , vol. 35, no. 2, pp. 1520–1530, 2020.
  * [45] B. K. Poolla, D. Groß, and F. Dörfler, “Placement and implementation of grid-forming and grid-following virtual inertia and fast frequency response,” _IEEE Trans. Power Syst._ , vol. 34, no. 4, pp. 3035–3046, 2019\.
  * [46] D. Groß and F. Dörfler, “Projected grid-forming control for current-limiting of power converters,” in _Conference on Communication, Control, and Computing (Allerton)_ , 2019.
  * [47] M. Raza, E. Prieto-Araujo, and O. Gomis-Bellmunt, “Small-signal stability analysis of offshore ac network having multiple vsc-hvdc systems,” _IEEE Transactions on Power Delivery_ , vol. 33, no. 2, pp. 830–839, 2018\.
  * [48] S. D’Arco, J. A. Suul, and O. B. Fosso, “A virtual synchronous machine implementation for distributed control of power converters in smart grids,” _Electric Power Systems Research_ , vol. 122, pp. 180–197, 2015.
  * [49] M. Amer, M. Ulrich, and F. Rolf, “Controller tuning in power systems using singular value optimization,” in _IFAC World Congress_ , 2020.
  * [50] R. Stanev and K. Nakov, “Power system stabilizers for inverter dominated future power systems,” in _International Symposium on Electrical Apparatus Technologies (SIELA)_ , 2020.
  * [51] I. Subotić, D. Groß, M. Colombino, and F. Dörfler, “A Lyapunov framework for nested dynamical systems on multiple time scales with application to converter-based power systems,” _IEEE Trans. Autom. Control_ , 2020.

  *[HAC]: hybrid angle control
  *[SM]: synchronous machines
  *[GFC]: grid-forming converter
  *[VSM]: virtual synchronous machine
  *[AVR]: automatic voltage regulator
  *[PSS]: power system stabilizer
  *[dVOC]: dispatchable virtual oscillator control
  *[RoCoF]: rate of change of frequency