\subsection{Proof of Proposition~\ref{prop:eva-reduc}}

We include here a formal definition of extended VA. An extended VA is a tuple $\bigcal{A} = (Q, \Sigma, \bigcal{X}, q_0, \Delta, F)$ where $Q$ is a finite set of states, $\Sigma$ is the input alphabet, $\bigcal{X}$ is a finite set of variables, $\Delta$ is a set consisting of {\it letter transitions} $(p, a, q)$ where $a\in\Sigma$ and $p,q\in Q$ or {\it extended variable transitions} $(p, S, q)$ where $S\subseteq \{\varop{x} , \varcl{x}\mid x\in\bigcal{X}\}$ and $S$ is nonempty; and $F\subseteq Q$. A run $\rho$ over a document $d = a_1\ldots a_n$ is a sequence
$$
\rho = q_0  \xrightarrow{S_1} p_0 \xrightarrow{a_1}
q_1  \xrightarrow{S_2} p_1 \xrightarrow{a_2} \cdots
\xrightarrow{a_n} q_n  \xrightarrow{S_{n+1}} p_n
$$
where each $S_i$ is a possibly empty set of markers, $(p_i, a_{i+1}, q_{i+1})\in \Delta$, and if $S_i$ is nonempty, then $(q_i, S_{i+1}, p_i)\in \Delta$, and if it is, then $p_i = q_i$. We say that a run is accepting if $p_n\in F$. Also, we say that a run $\rho$ is {\it valid} if variables are opened and closed in a correct manner: every marker $\varop{x}$ and $\varcl{x}$ must appear at most once each among the sets $S_1\ldots S_{n+1}$, and if they do, they both appear, and if $\varop{x}\in S_i$ and $\varcl{x}\in S_j$ then it holds that $i\leq j$. For a valid run $\rho$, we define the mapping $\mu^\rho$ that maps $x$ to $\spanc{i}{j}$ iff $\varop{x}\in S_i$ and $\varcl{x}\in S_j$. Lastly, we define the spanner $\sem{A}$ by the function $\sem{\bigcal{A}}(d) = \{\mu^\rho\mid \rho \text{ is an accepting and valid run of $\bigcal{A}$ over $d$}\}$. Like in varable-set automata, we say that $\bigcal{A}$ is sequential if every accepting run of $\bigcal{A}$ is valid. Also, we say that $\bigcal{A}$ is unambiguous if for every document $d$ and $\mu\in\bigcal{A}(d)$, there is only one accepting run $\rho$ of $\bigcal{A}$ over $d$ such that $\mu = \mu^\rho$.

Now let us describe the reduction. Let $\bigcal{A} = (Q, \Sigma, \bigcal{X}, q_0, \Delta, F)$ be an unambiguous sequential extended VA. We build an annotated automaton $\bigcal{A}' = (Q', \Sigma\cup\{\texttt{\#}\}, \Omega, \Delta', q_0, F')$ as follows.  Let $Q' = Q \cup \{q^*\}$ for some fresh state $q^*$, and let $F' = \{q^*\}$. Let $\Omega = 2^{\{\varop{x} , \varcl{x}\mid x\in\bigcal{X}\}}$ and let $\Delta' = \{(p,a,q)\mid a\in \Sigma\text{ and }(p,a,q)\in\Delta\} \cup \{(p,(a, S), q)\mid (p,S,q'),(q',a,q)\in\Delta \text{ for some }q'\in Q\} \cup \{(p,\texttt{\#}, q^*)\mid p\in F\} \cup \{(p,(\texttt{\#}, S), q^*)\mid (p,S,q)\in\Delta\text{ for some }q\in F\}$.

To see the equivalence between $\bigcal{A}$ and $\bigcal{A}'$, let $d = a_1\ldots a_n$ be a document over $\Sigma$, and let $\rho$ be an accepting run of $\bigcal{A}$ over $d$:
$$
\rho = q_0  \xrightarrow{S_1} p_0 \xrightarrow{a_1}
q_1  \xrightarrow{S_2} p_1 \xrightarrow{a_2} \cdots
\xrightarrow{a_n} q_n  \xrightarrow{S_{n+1}} p_n
$$
We define $\rho'$ as the following sequence:
$$
\rho' = q_0  \xrightarrow{b_1}
q_1  \xrightarrow{b_2} \cdots
\xrightarrow{b_n} q_n  \xrightarrow{b_{n+1}} q^*
$$
where $b_i  = (a_i, S_i)$ if $S_i$ is nonempty, and $b_i = a_i$ if it is for each $i\in[1, n]$. We define $b_{n+1} = (\texttt{\#}, S_i)$ if $S_i$ is nonempty, and $b_{n+1} = \texttt{\#}$ if it is. Since $\bigcal{A}$ is sequential, $\rho$ is a valid run which defines a mapping $\mu^\rho\in\sem{\bigcal{A}}(d)$. We can straightforwardly check that $\mu^\rho$ is equivalent to $\textsf{ann}(\rho')$.
It can also be seen directly from the construction that $\rho'$ is a run from $\bigcal{A}'$ over $d\texttt{\#}$, and since $\rho'$ is uniquely defined from $\rho$, we conclude that for every document $d$ each $\mu\in\sem{\bigcal{A}}(d)$ is equivalent to some unique $\nu\in\sem{\bigcal{A}'}(d\texttt{\#})$.

To see the equivalence on the opposite way, consider an accepting run $\rho'$ of $\bigcal{A}'$ over $d\texttt{\#}$ that is written as above, where each $b_i$ might be either $a_i$ or a pair $(a_i, S)$. From the construction, it can be seen that if $b_i = a_i$, for $i\in[1, n]$, there exists a transition $(q_{i-1},a_i, q_i)\in\Delta$. If $b_i = (a_i, S)$, for $i\in[1, n]$, there exist transitions $(q_{i-1},S, q'),(q',a_i,q_i)\in\Delta$ for some $q'\in Q$. Also, if $b_{n+1} = \texttt{\#}$ that means that $q_n\in F$ and if $b_{n+1} = (\texttt{\#}, S)$ there exists $(q_n,S,q')\in\Delta$ for some $q'\in F$. We define $\rho$ as a run of $\bigcal{A}$ over $d$ based on these transitions, and we note first that this run is accepting and valid, and since $\bigcal{A}$ is unambiguous, $\rho$ must be uniquely defined. For example, when replacing $(q_{i-1}, (a_i, S), q_i)$, we know there exist transitions $(q_{i-1}, S, q')$ and $(q', a_i, q_i)$ in $\Delta$. We also know this $q'$ must be unique, otherwise we could define a different accepting run with that defines the same mapping. We see that $\textsf{ann}(\rho')$ is equivalent to $\mu^\rho$, so we conclude that each annotation $\nu\in\sem{\bigcal{A}}(d)$ is equivalent to some unique mapping $\mu\in\sem{\bigcal{A}}(d)$. We also conclude that $\bigcal{A}'$ is unambiguous, otherwise for two runs $\rho_1'$ and $\rho_2'$ that render the same annotation, and that differ in some $q_i$, we would be able to build runs $\rho_1$ and $\rho_2$ that also differ in some $q_i$ that define the same mapping.

\subsection{Proof of Theorem~\ref{theo-algo-succ}}

For this proof, we will use an extension of Shift-ECSs\xspace called {\it succinct Shift-ECSs\xspace}.
A succinct Shift-ECSs\xspace as a tuple $\tilde{\bigcal{D}} = (\bigcal{S}, V, \mathcal{l}\!, \mathcal{r}, \lambda, \bot, \epsilon)$ as a Shift-ECS\xspace, with the difference that the output set has been replaced by a SERS $\bigcal{S}$. The set of strings associated to a node $v\in \tilde{D}$ is defined as follows: If $\lambda(v) = r\in R$, where $R$ is the set of representations in $\bigcal{S}$, then $\sem{\tilde{D}}(v) = \bigcal{L}(r)\times\{1\}$. The rest of the sets $\sem{\tilde{D}}(v)$ for union, product, and shift nodes $v$ remain unchanged, and so are the notions of duplicate-free, $k$-bounded and $\varepsilon$-safe nodes, along with the operations for $\mathsf{add}$, $\prod$, $\mathsf{union}$ and $\mathsf{shift}$.

\begin{proposition}\label{prop:ds-succ}
	The operations $\prod$, $\mathsf{union}$ and $\mathsf{shift}$ over succinct Shift-ECSs\xspace{} extended with empty- and $\epsilon$-nodes take constant time, and the operation $\mathsf{add}$ with a representation $r$ as argument takes time $O(|r|)$. Furthermore, if we start from an empty succinct Shift-ECS\xspace{} $\bigcal{D}$ and apply $\mathsf{add}$, $\prod$, $\mathsf{union}$, and $\mathsf{shift}$ over $\epsilon$-safe nodes, the resulting node $v'$ is always an $\epsilon$-safe node, and the set $\sem{\bigcal{D}}(v)$ can be enumerated with output-linear delay without preprocessing for every node $v$.
\end{proposition}
\begin{proof}
	For this proof, we need to be a bit more specific regarding the  enumeration algorithm $\bigcal{E}$ of an SERS. Precisely, we assume that $\bigcal{E}$ has an associated function $\texttt{yield}$ and a constant $c$, and its procedure is to receive $r$ and then allow up to $|\bigcal{L}(r)|$ calls to $\texttt{yield}$. Each of these calls produces a different annotation $\oout\in \bigcal{L}(r)$ and takes time $c\cdot (|\oout|+1)$. Each of these annotations carries a flag ${\tt end}$ which is true if, and only if, it is the last output of the set. We also assume that a sequence of $|\bigcal{L}(r)|$ calls to $\texttt{yield}$ might happen again with the same time bounds after the last call set ${\tt end}$ to true.

	We start this proof by re-stating that every definition in the Shift-ECS\xspace data structure is unchanged in its succinct version. The only difference is the definition of $\sem{\tilde{\bigcal{D}}}(v)$ for a bottom node $v\in\tilde{\bigcal{D}}$.

	First, we shall prove a version of Proposition~\ref{prop:lindelay} in this model, namely, that if a succinct Shift-ECS\xspace $\tilde{\bigcal{D}}$ is duplicate-free and $k$-bounded, then the set $\sem{\tilde{\bigcal{D}}}(v)$ can be enumerated with output-linear delay for every $v\in \bigcal{D}$. We do this by adapting the proof of Proposition~\ref{prop:lindelay} to handle this new type of bottom node. Let $\bigcal{S} = (\bigcal{R}, \Omega, |\cdot |, \bigcal{L}, \bigcal{E})$ be the SERS associated to $\tilde{\bigcal{D}}$ and let $c'$ be the constant associated to $\bigcal{E}$. In a bottom node $v$, the subroutines ${\tt first}(v, s)$, ${\tt next}(v, s)$ and ${\tt curr}(v, s)$ are defined as follows: To call ${\tt first}(v, s)$, we assume that the procedure ${\tt yield}$ from $\bigcal{E}$ has either not been called over $r$, or that the last call over $r$ set its flag ${\tt end}$ to true.
	The procedure ${\tt first}(v, s)$ then calls ${\tt yield}$ over $r = \lambda(v)$  and obtains an annotation $\oout$ in time $c'(|\oout| + 1)$ and prints $(\oout, s + 1)$ as output. The procedure ${\tt next}(v, s)$ does the same, but it the opposite is assumed. That is, that the last call to ${\tt yield}$ over $r$ did not set ${\tt end}$ to true. Both of these calls store the last annotation $\oout$ locally. Then, the procedure ${\tt curr}(v, s)$ returns $(\oout, s)$ using the $\oout$ that was stored from the last call. Each of these procedures returns ${\tt end}$ if and only if the call to ${\tt yield}$ that was just made (or the previous one in the case of ${\tt curr}$) set its own flag ${\tt end}$ to true. It can be seen now that calling ${\tt first}(v, s)$ and then ${\tt next}(v, s)$ until one of these calls returns ${\tt end}$, then this prints the entire set  $\mathsf{sh}(\sem{\tilde{D}}(v), s) = \bigcal{L}(r)\times\{s + 1\}$ without repetitions. This was the inductive assumption that used to prove correctness of the enumeration algorithm, so it holds.

	To show the time bounds, we just bring attention to the fact that the only difference in the algorithm time-wise is the time spent printing each pair $(\oout, i)$, as now this takes  $c'(|\oout| + 1)$. However, this does not affect output-linear delay since, when taking all pairs that belong to an output, the accumulate time is proportional to the total size of it plus its length.

	The rest of the proof pertains the operations $\mathsf{add}$, $\prod$, $\mathsf{union}$ and $\mathsf{shift}$, but it is not hard to see that they can be kept unchanged and maintain all the conditions in the statement, with the sole exception of $\mathsf{add}(r)$, which now adds a representation $r$ to the structure and takes $|r|$ time. This concludes the proof.
\end{proof}

The algorithm that we give to prove the statement is exactly Algorithm~\ref{alg:evaluation} line-by-line, with the sole exception of lines 10-11, which now read: $\textbf{for each }(p,a,\oout, q) \in \Delta$: $M_a[p,q] \gets \mathsf{union}(M_a[p,q],\, \mathsf{add}(\oout))$. The rest of the proof follows from the reasoning of the proof of Theorem~\ref{theo:evaluation}, and by noticing that the definition of the set $\sem{\tilde{D}}(v) = \bigcal{L}(r) \times \{1\}$ satisfies what is expected for an index $M_a[p,q]$. That is, that for every partial run $\rho$ of $\bigcal{T}\!$ over $a$ that starts on $p$ and end on $q$, all annotations from the set $\textsf{ann}(\rho)$ are included. The fact that the operations are done duplicate-free follows from the unambiguous condition over $\bigcal{T}\!$ as well.

\subsection{Proof of Proposition~\ref{prop:va-reduc}}

The SERS that we will consider are Enumerable Compact Sets as they were presented in~\cite{MunozR22}, defined over the output set $\Omega = \{\varop{x} , \varcl{x}\mid x\in\bigcal{X}\}$, with the sole difference that the sets of outputs stored in a node no longer contain strings in $\Omega^*$, but subsets of $\Omega$ instead. This difference is merely technical and has no influence in the time bounds of an ECS, as long as it is guaranteed to be duplicate-free (called {\it unambiguous} in~\cite{MunozR22}).

First we assume that all transitions in $\bigcal{A}$ are reachable from $q_0$, and all of them reach an accepting state. We know that in the VA $\bigcal{A}$, the graph induced by the variable transitions is a DAG, otherwise it would not be sequential. We start with an empty ECS $\bigcal{D}$ and define a matrix $K[p,q]$, that first starts with the empty node in each index except the diagonal, which have the $\varepsilon$ node. The idea is that at the end of this algorithm, $\bigcal{L}_{\bigcal{D}}(K[p,q])$ contains all sets of variable markers that can be seen in a path from $p$ to $q$ that does not contain any letter transition. We build this matrix by iterating over the variable transitions in $\bigcal{A}$ following some topological order of the DAG, starting from a root. For each variable transition $(p,V,q)$, let $u\gets\mathsf{add}(V)$ and for each $p'\in Q$ that has a lower value in the topological order assign $K[p',q] \gets \mathsf{union}(K[p',q], \prod(K[p', p], \mathsf{add}(V)))$. The time of this algorithm is $|Q|\times|\bigcal{A}|$.

After doing this, we perform a construction analogous to the one done in Proposition~\ref{prop:eva-reduc}, in which we replace the extended variable transitions $(p,S,q)$ by the indices in which $K[p,q]$ is not the empty node. The proof follows from that result as well.

\subsection{Proof of Proposition~\ref{prop:lindelay}}

\newcommand{\mathsf{len}}{\mathsf{len}}

The algorithm that we present is a depth-first traversal of the DAG, done in a particular way to ensure that after retrieving some output the next one can be found in $O(k\cdot \ell)$ time, where $\ell$ is the size of the next output. First we will assume that the algorithm receives a duplicate-free and $k$-bounded $\bigcal{D}$ and a bottom or product node $v$, and then extend it to the case of a union or shift node.

This algorithm uses subroutines ${\tt first}(v, s)$, ${\tt next}(v, s)$ and ${\tt curr}(v, s)$ defined only for product and bottom nodes which print the first element in $\mathsf{sh}(\sem{\bigcal{D}}(v), s)$ (according to an order that is induced by the algorithm itself), the element that follows the previous one printed, and the same element as the one previously printed respectively. We use the word ``print'' to refer to the fact that these elements are given as a global output of the algorithm, not in the local subroutine. This means that the pairs $(\oout, i)$ are given as output as soon as some of these calls is done over a bottom node $v$ for which $\lambda(v) = \oout$. We also assume that each of these subroutines returns (locally) a flag ${\tt end}$ whenever the output printed is the last one of the set. We will inductively assume that calling ${\tt first}(v, s)$, and then ${\tt next}(v, s)$ repeatedly until some of these calls returns ${\tt end}$ prints the entire set $\mathsf{sh}(\sem{\bigcal{D}}(v), s)$ without repetitions.

In a bottom node $v$, the subroutines ${\tt first}(v, s)$ and ${\tt curr}(v, s)$ each print the element $(\lambda(v), s + 1)$ and return ${\tt end}$. We assume that ${\tt next}(v, s)$ is never called on a bottom node.

We now describe the subroutines for a product node $v$. First, we will assume there are no shift nodes that are descendants of $v$ in $\bigcal{D}$, and afterwards we will adapt the algorithm for shift nodes. Also, consider a procedure $St\gets {\tt trav}(u)$ defined over a union node $u$ which does a traversal over the left descendants of $u$ until it sees a product or bottom node. These nodes are added to the stack $St$. Also, consider a procedure ${\tt move}(St)$ which pops the last two nodes from $St$, and let $u'$ be the second node. Then, $St'\gets{\tt trav}(r(u'))$ is performed and the nodes in $St'$ are added to $St$ in the same order. We assume that ${\tt move}(St)$ cannot be applied if $St$ has only one node. An example of the result of these procedures can be seen in Figure~\ref{fig-enum-stacks}.

\begin{figure}[t]
	\centering
	\input{enumfig.tex}
	\caption{An example iteration of ${\tt trav}$ and ${\tt move}$. The sequences of nodes joined by dashed lines represent a stack $St$, where the first one was obtained after calling ${\tt trav}$ over the topmost union node, and the following five are obtained by repeated applications of ${\tt move}(St)$.}
	\label{fig-enum-stacks}
\end{figure}

Now, ${\tt first}(v, s)$ operates as follows: We perform $St_1 \gets {\tt trav}(\ell(u))$, and $St_2 \gets {\tt trav}(r(u))$, and let $v_1$ and $v_2$ be the nodes at the top of each stack. Then the first output is produced by calling ${\tt first}(v_1, s)$ and ${\tt first}(v_2, s)$. We define ${\tt next}(v, s)$ as follows: First, check if the last call to $v_2$ returned ${\tt end}$, and if it did not, call ${\tt curr}(v_1, s)$ and ${\tt next}(v_2, s)$. If it did return ${\tt end}$, then check if $St_2$ has more than one element, and if it does, we apply ${\tt move}(St_2)$, and then call ${\tt curr}(v_1, s)$ and ${\tt first}(v_2', s)$ where $v_2'$ is the new node at the top of $St_2$. If it has only one node, we check if the last call to $v_1$ returned ${\tt end}$, and if it did not, we call ${\tt next}(v_1, s)$. If it did, then we assume $St_1$ has more than one element, so we apply ${\tt move}(St_1)$ and call ${\tt first}(v_1', s)$, where $v_1'$ is the new node at the top of $St_1$. After a call is done over $v_1$ in any of these cases, we perform $St_2 \gets {\tt trav}(r(u))$, and then call ${\tt first}(v_2, s)$. In ${\tt next}(v)$, the flag ${\tt end}$ is returned when both calls to $v_1$ and $v_2$ return ${\tt end}$ and both stacks $St_1$ and $St_2$ have only one element. We now have to prove that calling ${\tt first}(v, s)$ and then ${\tt next}(v, s)$ repeatedly until one of these returns ${\tt end}$ effectively prints the entire set $\mathsf{sh}(\sem{\bigcal{D}}(v), s)$ without repetitions. First, let us introduce the concept of a {\it union-descendant} of a node $u$. We say that $u$ is an union-descendant of itself, and that its children also are union-descendants of it, but only if $u$ is a union node. Alternatively, we say that $u'$ is a union-descendant of $u$ if $u'$ can be reached starting from $u$ and following a path where all the intermediate nodes are union nodes. Let $w$ be an element in $\mathsf{sh}(\sem{\bigcal{D}}(v), s)$. Since $\bigcal{D}$ is duplicate-free, there exists a unique $w_1$ and $w_2$ such that $w_1\in\mathsf{sh}(\sem{\bigcal{D}}(\ell(v)), s)$ and $w_2\in\mathsf{sh}(\sem{\bigcal{D}}(r(v)), s)$ and $w = w_1\cdot w_2$. Again, since $
\bigcal{D}$ is duplicate-free, there exist a unique union-descendant $v_1$ of $\ell(v)$ such that $w_1\in\mathsf{sh}(\sem{\bigcal{D}}(v_1), s)$. The same can be said for $w_2$ and a node $v_2$ which is a union-descendant of $r(v)$. From the algorithm it can be seen that for each pair of output nodes $v_1'$ and $v_2'$ which are union-descendants of $r(v)$ and $\ell(v)$ respectively the following calls are made: first ${\tt first}(v_1, s)$ and ${\tt first}(v_2, s)$ is called; then, ${\tt curr}(v_1, s)$ and ${\tt next}(v_2, s)$ is called until ${\tt next}(v_2, s)$ (or ${\tt first}(v_2, s)$) returns ${\tt end}$; then, ${\tt next}(v_1, s)$ and ${\tt first}(v_2, s)$ is called; then ${\tt curr}(v_1, s)$ and ${\tt next}(v_2, s)$ again until ${\tt next}(v_2, s)$ returns ${\tt end}$; this is done until ${\tt curr}(v_1, s)$ (or ${\tt first}(v, s)$) returns ${\tt end}$, and ${\tt next}(v_2, s)$ (or ${\tt first}(v_2, s)$) also returns ${\tt end}$. This, altogether, means that for each $w_1'\in \mathsf{sh}(\sem{\bigcal{D}}(\ell(v)), s)$ and $w_2'\in \mathsf{sh}(\sem{\bigcal{D}}(r(v)), s)$ the string $w_1'\cdot w_2'$ is printed, as we assumed that the calls ${\tt first}$, ${\tt next}$ and ${\tt curr}$ work as specified. This proves that each output $w\in \mathsf{sh}(\sem{\bigcal{D}}(v), s)$ is printed, and this happens exactly once since $w_1$ and $w_2$ are uniquely defined, are found in exactly one pair of output-node union-descendants $v_1$ of $\ell(v)$ and $v_2$ of $r(v)$, and in both of these nodes every output is printed exactly once.

To account for shift nodes, we modify ${\tt trav}$ so it stores pairs $(u, s')$ in the stack, where $u$ is a union or output node, and $s'$ is the sum of all shift nodes seen in the path from $v$ to $u$. Then, whenever ${\tt first}$, ${\tt next}$ or ${\tt curr}$ is called, it is instead called over ${\tt first}(u, s + s')$, ${\tt next}(u, s + s')$ and ${\tt curr}(u, s + s')$. To do this, we keep a variable $s^*$, and whenever a shift node is seen in the traversal, its shift value is added to $s^*$. Then, when a union node is seen, it is added to the stack alongside the value $s^*$. It can be seen that this modification accounts for all the needed shifts in any of the calls to ${\tt first}$, ${\tt next}$ and ${\tt curr}$ on an union-descendant of $\ell(v)$ and $r(v)$, so we conclude that the set $\mathsf{sh}(\sem{\bigcal{D}}(\ell(v)), s)$ is printed correctly.

To show the time bounds of this procedure, first we note the fact that each pair $(\oout, i)$ is printed directly so we do not take this time into account, as it only adds constant time per element in the output. Call $\mathsf{len}(w)$ the length of $w$, as in, the number of pairs $(\oout, i)$ in it. Let $c$ be the time needed to do the constant-time operations utilized in the algorithm, such as pushing and popping from the stack, navigating to some child of a node, etc. We inductively assume that for any product or bottom node $u$, a call to ${\tt first}(u, s)$ to ${\tt next}(u, s)$ and ${\tt curr}(u, s)$ each take at most $2c\cdot (k+2)\cdot (\mathsf{len}(w)-1)$ time, where $w$ is the output produced. Note that this follows trivially for a bottom node, since the time taken is merely what is needed to print the single pair associated to it, which we already accounted for, and the extra time is 0.
For a product node $v$, the call ${\tt first}(v, s)$ takes time at most $c\cdot 2(k+1)$ to build the stacks $St_1$ and $St_2$, and then $2c\cdot (k+2)\cdot (\mathsf{len}(w_1^*) + \mathsf{len}(w_2^*) - 2)$ to call ${\tt first}(v_1, s')$ and ${\tt first}(v_2, s'')$, which print $w_1^*$ and $w_2^*$ respectively. Therefore, by using that $\mathsf{len}(w^*) = \mathsf{len}(w_1^*) + \mathsf{len}(w_2)$ the call to ${\tt first}(v)$ takes time at most $2c\cdot (k+2)\cdot (\mathsf{len}(w^*)-1)$. For the call ${\tt next}(v, s)$, we have to account for several cases, but we simply note that the case which uses the most operations besides the calls to $v_1$ and $v_2$ is the one where both previous calls to $v_1$ and $v_2$ returned ${\tt end}$, and $St_2$ has only one node. In this case ${\tt move}(St_1)$ has to be called, and $St_1\gets {\tt trav}(r(v))$ has to be called, which take at most time $c(k+2)$ and $c(k+1)$ respectively. Adding this time to the time needed to make a call to $v_1$ and $v_2$ gives a total time at most $2c\cdot (k+2)\cdot (\mathsf{len}(w^*)-1)$. This proves that the set $\mathsf{sh}(\sem{\bigcal{D}}(v), s)$ is enumerated with output-linear delay.

Finally, we refer to the case in which the enumeration is done over a union or shift node $v$. This case is handled in the same way the product-node case was. More precisely, we start by calling $St\gets {\tt trav}(v)$, and then proceed in the exact same way as described in the subroutines for the product node, by calling ${\tt first}(v_1, s)$, then ${\tt next}(v_1, s)$ util they return ${\tt end}$ and then applying ${\tt move}(St)$ and calling ${\tt first}(v_1, s)$ again. This adds constant time $c\cdot(k+2)$ to produce each output, which maintains the output-linear delay bounds we showed. This completes the proof.

\subsection{Proof of Theorem~\ref{theo:data-structure}}

The operations $\mathsf{add}$, $\prod$, and $\mathsf{shift}$ were defined in the body of the paper, so what is remaining is to go over the cases for $\mathsf{union}$ that were not defined. That is, the cases where at least one node among $v_3$ and $v_4$ is a shift node with an output node as its left child.

Without loss of generality, assume that $v_3$ has an output node as its left child. We extend $\bigcal{D}$ with nodes $v'$, $v''$ and $v'''$, where $\ell(v) = \lambda(v_3)$, $\lambda(v') = \cup$ and $\lambda(v''')  = \lambda(v_4) - \lambda(v_3)$, $\ell(v) = v''$, $\ell(v) = \ell(v_3)$, $r(v) = v'''$ and $\ell(v''') = \ell(v_3)$. This construction is illustrated in FIgure~\ref{fig-theo-ds}. It can be seen that $v'$ is a safe node since the left child of $v''$ is an output node, and $v'''$ has output depth at most 2, either if $\ell(v_4)$ is an output node, or if $\ell(v_4)$ is a union node whose left child is an output node. Since $\bigcal{D}$ was assumed to be 3-bounded, and this operation only adds nodes with output depth at most 2, the new version of $\bigcal{D}$ is 3-bounded as well.

The fact that the set $\sem{D}(v)$ can be enumerated with output-linear delay for every node $v$ in $\bigcal{D}$ follows from Proposition~\ref{prop:lindelay}, since the operations are restricted to keep $\bigcal{D}$ duplicate-free and they result in a new version of $\bigcal{D}$ that is 3-bounded.

\begin{figure}
	\centering
	\begin{tikzpicture}[every label/.style={black!60}, ->,>=stealth',roundnode/.style={circle,inner sep=1pt},squarednode/.style={rectangle,inner sep=2pt}, scale=0.85]
		\node[squarednode] (0) at (0, 0) {$\ell(v_3)$};
		\node[squarednode, draw, label=180:\footnotesize{$v_3\!:$}] (2) at (2, 1) {$k_1$};
		\node[squarednode] (3) at (4, 0) {$\ell(v_4)$};
		\node[squarednode, draw, label=180:\footnotesize{$v_4\!:$}] (5) at (6, 1) {$k_2$};
		\node[squarednode, label=180:\footnotesize{$v'''\!:$}] (7) at (5.5, 1.7) {$(k_2-k_1)$};
		\node[squarednode, label=180:\footnotesize{$v''\!:$}] (8) at (3, 2.6) {$\cup$};
		\node[squarednode, draw, label=180:\footnotesize{$v'\!:$}] (9) at (4.5, 3.6) {$k_1$};
		\draw[dashed] (2.south west) to[out=-135,in=45] (0);
		\draw[dashed] (5.south west) to[out=-135,in=45] (3);
		\draw[dashed] (7.south west) to[out=-135,in=90] (3);
		\draw (8) to[out=-45,in=135] (7.north west);
		\draw[dashed] (8) to[out=-135,in=90] (0);
		\draw[dashed] (9.south west) to[out=-135,in=90] (8);

	\end{tikzpicture}
	\caption{Gadget used in Theorem~\ref{theo:data-structure}.}
	\label{fig-theo-ds}
\end{figure}

\subsection{Proof of Theorem~\ref{theo:data-structure-eps}}

Extending these four operations to deal with empty nodes is straightforward, and it has been explained in the body, so we will focus on operations that deal with the $\varepsilon$-node.

First, we add the condition that for any Shift-ECS\xspace $\bigcal{D}$ the $\varepsilon$-node is either disconnected to the rest of the structure (parent-less), or if it has a parent, then this must be a union node, and $\varepsilon$ must be its left child. We call this the {\it $\varepsilon$ condition}. We will see that any node that might be the result of applying the operations $\mathsf{add}$, $\prod$, $\mathsf{union}$ and $\mathsf{shift}$ is equivalent to a node in a Shift-ECS\xspace $\bigcal{D}$ that satisfies this.

The notions of output-depth, duplicate-free, and $k$-boundedness are mostly unchanged in this setting: If $v$ is the $\varepsilon$-node, then $\mathsf{odepth}(v) = 0$, and if $v$ is a union node with the $\varepsilon$-node as its left child, $\mathsf{odepth}(v) = 1$, the rest of the nodes have no descendants that are the $\varepsilon$-node, so the definition of output depth remains the same; and the definitions of duplicate-free and $k$-bounded remain unchanged as well. By assuming the $\varepsilon$ condition, we can state the following claim which is a direct consequence of Proposition~\ref{prop:lindelay}:

\begin{claim}\label{prop:lindelayeps}
	Fix $k\in\mathbb{N}$. Let $\bigcal{C}_k$ be the class of all  duplicate-free and $k$-bounded Shift-ECSs\xspace{} that satisfy the $\varepsilon$ condition. Then one can solve the problem $\textsc{ShiftECSEnum}[\bigcal{C}_k]$ with output-linear delay and without preprocessing (i.e. constant preprocessing time).
\end{claim}

We now define the notion of {\it $\varepsilon$-safe} nodes. These are nodes $v$ which satisfy one of three conditions: (1) $v$ is the $\varepsilon$-node, (2) $v$ is safe and none of its descendants is the $\varepsilon$-node, and (3) $v$ is a union node, its left child is the $\varepsilon$-node, and its right child is a safe node for which none of its descendants is the $\varepsilon$-node.
To guide the constructions, we call these conditions 1, 2 and 3.
These conditions will allow us to use the constructions already given for $\prod$, $\mathsf{union}$, and $\mathsf{shift}$ over safe nodes, such as node $v$, in the case of condition 2, or $r(v)$, in the case of condition 3.
The rest of this proof will be extending these operations to work over $\varepsilon$-safe nodes, and check that they can be done in constant time, return an $\varepsilon$-safe node, maintain the $\varepsilon$ condition in the Shift-ECS\xspace, and satisfy the specifications of each operation (i.e., that for the resulting node $v'$ the set $\sem{D}(v')$ contains what it should).

The operation $\mathsf{add}$ is defined the same as in a Shift-ECS\xspace without $\varepsilon$. The operation $\mathsf{add}(\varepsilon)$ is not defined since the $\varepsilon$-node is always part of $\bigcal{D}$.

The $\prod$ case is somewhat more involved. First, assume that some node among $v_1$ and $v_2$ is the $\varepsilon$-node (condition 1). Without loss of generality, assume that $v_1$ is. If $v_2$ is the $\varepsilon$-node as well, we simply return $v_1$; and if $v_2$ is not the $\varepsilon$-node (conditions 2 or 3), we simply return $v_2$. The requirements of the construction follow trivially. Now, assume that both $v_1$ and $v_2$ satisfy condition 2. For this case we simply return $\prod(v_1, v_2)$ as it was defined for a Shift-ECS\xspace without $\varepsilon$. Lastly, assume at least one node among $v_1$ and $v_2$ satisfies condition 3 (and none of them are the $\varepsilon$-node). These cases are depicted in Figure~\ref{fig-prod-multi-gadget}, and are formally described as follows, using the operations $\mathsf{union}$ and $\prod$ as they were defined for a Shift-ECS\xspace without $\varepsilon$:
\begin{itemize}
	\item[(a)] $v_1$ satisfies condition 2 and $v_2$ satisfies condition 3. Let $v'' \gets \prod( v_1, r(v_2))$, and let $v' \gets \mathsf{union}( v_1, v'')$.
	\item[(b)] $v_1$ satisfies condition 3 and $v_2$ satisfies condition 2. Let $v'' \gets \prod( r(v_1), v_2)$, and let $v' \gets \mathsf{union}( v'', v_2)$.
	\item[(c)] Both $v_1$ and $v_2$ satisfy condition 3. Let $v'' \gets \prod( r(v_1), r(v_2))$, let $v^3 \gets \mathsf{union}( v'', r(v_2))$, and let $v^4 \gets \mathsf{union}( r(v_1), v^3)$. Lastly, let $\lambda(v') = \cup$, $\ell(v') = \varepsilon$, and $r(v') = v^4$.
\end{itemize}

It can be seen that these constructions work as expected. That is, that $\sem{\bigcal{D}'}(v) = \sem{\bigcal{D}}(v_1)\cdot\sem{\bigcal{D}}(v_2)$. Also, if $v_1$ and $v_2$ are $\varepsilon$-safe, then $v'$ is $\varepsilon$-safe as well. Also, each construction does a fixed number of steps, so they take constant time.

\begin{figure}[t]
	\centering
	\begin{tikzpicture}[->,>=stealth',roundnode/.style={circle,draw,inner sep=1.2pt},squarednode/.style={rectangle,inner sep=2pt}, scale=1]
		\node [squarednode] (c) at (1.12, -0.6) {$(a)$};
		\node [squarednode] (0) at (0.8, 2.75) {$\mathsf{union}$};
		\node [squarednode] (6) at (-0.05, 2.75) {$v'\gets$};
		\node [squarednode] (1) at (0, 1) {$v_1$};
		\node [squarednode] (2) at (1.6, 1.75) {$\prod$};
		\node [squarednode] (3) at (1.6, 1) {$v_2$};
		\node [squarednode] (4) at (0.8, 0) {$\varepsilon$};
		\node [squarednode] (5) at (2.4, 0) {$r(v_2)$};
		\draw[dashed] (0) to[out=-135,in=90] (1);
		\draw (0) to[out=-45,in=120] (2);
		\draw (2) to[out=-45,in=90] (5);
		\draw[dashed] (3) to[out=-135,in=60] (4);
		\draw (3) to[out=-45,in=120] (5);
		\draw[dashed] (2) to (1);
	\end{tikzpicture}
	\begin{tikzpicture}[->,>=stealth',roundnode/.style={circle,draw,inner sep=1.2pt},squarednode/.style={rectangle,inner sep=2pt}, scale=1]
		\node [squarednode] (c) at (1.6, -0.6) {$(b)$};
		\node [squarednode] (0) at (2.5, 2.75) {$\mathsf{union}$};
		\node [squarednode] (6) at (1.65, 2.75) {$v'\gets$};
		\node [squarednode] (1) at (3, 1) {$v_2$};
		\node [squarednode] (2) at (1.90, 1.75) {$\ \prod$};
		\node [squarednode] (3) at (0.85, 1) {$v_1$};
		\node [squarednode] (4) at (0.2, 0) {$\varepsilon$};
		\node [squarednode] (5) at (1.5, 0) {$r(v_1)$};
		\draw (0) to[out=-60,in=90] (1);
		\draw[dashed] (0) to[out=-135,in=70] (2);
		\draw[dashed] (2) to[out=-120,in=90] (5);
		\draw[dashed] (3) to[out=-135,in=60] (4);
		\draw (3) to[out=-45,in=120] (5);
		\draw (2) to (1);
	\end{tikzpicture}
	\begin{tikzpicture}[->,>=stealth',roundnode/.style={circle,draw,inner sep=1.1pt},squarednode/.style={rectangle,inner sep=2pt},scale=1.1]
		\node [squarednode] (c) at (1.6, -0.5) {$(c)$};
		\node [squarednode] (0) at (1.5, 2.75) {$\cup$};
		\node [squarednode] (11) at (0.88, 2.8) {$v'\gets$};
		\node [squarednode] (1) at (0.75, 2.25) {$\varepsilon$};
		\node [squarednode] (2) at (2.25, 2.25) {$\mathsf{union}$};
		\node [squarednode] (3) at (3, 1.75) {$\mathsf{union}$};
		\node [squarednode] (4) at (2.1, 1) {$\prod$};
		\node [squarednode] (5) at (3.2, 1) {$v_2$};
		\node [squarednode] (6) at (2.5, 0) {$\varepsilon$};
		\node [squarednode] (7) at (4, 0) {$r(v_2)$};
		\node [squarednode] (8) at (0.5, 1) {$v_1$};
		\node [squarednode] (9) at (0, 0) {$\varepsilon$};
		\node [squarednode] (10) at (1, 0) {$r(v_1)$};
		\draw[dashed] (0) to (1);
		\draw (0) to (2);
		\draw[dashed] (2) to[out=-135,in=70] (10);
		\draw[dashed] (4) to (10);
		\draw (4) to[out=-45,in=145] (7);
		\draw[dashed] (5) to[out=-135,in=60] (6);
		\draw (5) to[out=-45,in=120] (7);
		\draw (3) to[out=-45,in=90] (7);
		\draw[dashed] (3) to (4);
		\draw (2) to (3);
		\draw (8) to[out=-50,in=100] (10);
		\draw[dashed] (8) to[out=-125,in=80] (9);
	\end{tikzpicture}
	\caption{Gadgets for $\prod$ as defined for an Shift-ECS\xspace with the $\varepsilon$-node.}
	\label{fig-prod-multi-gadget}
\end{figure}

We now address the $\mathsf{union}$ case:
\begin{itemize}
	\item Assume $v_3$ is the $\varepsilon$-node (condition 1). (1) If $v_4$ also is the $\varepsilon$-node, then we simply return $v'\gets\varepsilon$. (2) If $v_3$ satisfies condition 2, let $v'$ be a union node such that $\ell(v') = \varepsilon$ and $r(v') = v_4$. It can be seen that $v'$ satisfies condition 3 and that its $\varepsilon$-safe. (3) If $v_4$ satisfies condition 3, we simply return $v'\gets v_4$.
	\item Assume $v_3$ satisfies condition 2. (1) If $v_4$ is the $\varepsilon$-node, we return $v'\gets \mathsf{union}( v_4, v_3)$ using the construction for the case (2) in the previous item. (2) If $v_4$ satisfies condition 2, then we return $v'\gets \mathsf{union}( v_3, v_4)$ as it was defined for Shift-ECSs\xspace without $\varepsilon$. (3) If $v_4$ satisfies condition 3, let $v'' \gets \mathsf{union}( v_3, r(v_4))$, and let $v'$ be a union node such that $\ell(v') = \varepsilon$ and $r(v) = v''$. Again, it can be seen that $v'$ is $\varepsilon$-safe.
	\item Assume $v_3$ satisfies condition 3. (1) If $v_4$ is the $\varepsilon$-node, we simply return $v'\gets\varepsilon$. (2) If $v_4$ satisfies condition 2, we follow the construction given in case (3) of the previous item. (3) If $v_4$ satisfies condition 3, let $v''\gets\mathsf{union}( r(v_3), r(v_4))$ and let $v'$ be a union node such that $\ell(v') = \varepsilon$ and $r(v) = v''$. It can be seen that $v'$ is $\varepsilon$-safe.
\end{itemize}

Clearly, in each of these cases we have that $\sem{\bigcal{D}}(v') = \sem{\bigcal{D}}(v_3)\cup \sem{\bigcal{D}}(v_4)$, and if $v_3$ and $v_4$ are $\varepsilon$-safe, then $v'$ is $\varepsilon$-safe as well. Plus, they all take constant time and in each case the new version of $\bigcal{D}$ satisfies the $\varepsilon$ condition.

The operation $\mathsf{shift}( v, k)$ is defined as follows: If $v$ is the $\varepsilon$-node, we return $v'\gets\varepsilon$. If $v$ satisfies condition 2, we return $v'\gets\mathsf{shift}( v, k)$ as it was defined for Shift-ECSs\xspace without $\varepsilon$. If $v$ satisfies condition 3, let $v''\gets\mathsf{shift}(r(v), k)$ and let $v'$ be a union node such that $\ell(v') =\varepsilon$ and $r(v) = v''$. It can be seen that $v'$ is $\varepsilon$-safe.

We conclude the proof by stating that the operations $\mathsf{add}$, $\prod$, $\mathsf{union}$ and $\mathsf{shift}$, when done over $\varepsilon$-safe nodes, can be done in constant time, return an $\varepsilon$-safe node, maintain the $\varepsilon$ condition in the Shift-ECS\xspace, and satisfy the specifications of each operation. Furthermore, by virtue of Claim~\ref{prop:lindelayeps}, for each node $v$ in $\bigcal{D}$ the set $\sem{D}(v)$ can be enumerated with output-linear delay.

\subsection{Proof of Lemma~\ref{lem:ecsext}}

Let $\bigcal{D}$ be the succinct Shift-ECS\xspace associated to the query data structure for $S$ and $\bigcal{A}$ and let $\tilde{R} = R' \setminus R$ be the set of new rules in $S'$. Consider the DAG that is induced by this set. We can go over these new rules $\tilde{A}\to BC$ in a bottom-up fashion, starting from the ones where $B,C\in N$, and define the nodes $M_{\tilde{A}}[p,q] = (M_{B}\otimes \mathsf{shift}(M_C,|\operatorname{doc}(B)|))[p,q]$ for each $p,q\in Q$. Then, we build the nodes $v_{\tilde{A}}$ for each new non-terminal $\tilde{A}$ as defined. The new nodes are created by application of the Shift-ECS\xspace operations $\mathsf{add}$, $\mathsf{union}$, $\prod$ and $\mathsf{shift}$ which create the succinct Shift-ECS\xspace $\bigcal{D}'$ that defines the query data structure for $S'$ and $\bigcal{A}$. The procedure takes $O(|Q|^3\cdot|\tilde{N}|)$ time.

\subsection{Proof of Theorem~\ref{theo:evaluation}}

As the time bounds are explained in the body, this proof will only regard correctness of Algorithm~\ref{alg:evaluation}.

The proof will follow after proving that for each $M_X$, for an $X\in\Sigma\cup N$ that is reachable from $S_0$, the set $\sem{\bigcal{D}}(M_X[p,q])$ contains all partial annotations of a run of $\bigcal{A}$ over $\operatorname{doc}(X)$ that starts on $p$ and ends on $q$. First let $X = a\in\Sigma$. This case can be quickly verified by inspecting the procedure {\sc Terminal}$(a)$, as for each index, the set $\sem{\bigcal{D}}(M_a[p,q])$ contains $\varepsilon$, and the annotations $\oout$ in the pair $(\oout, 1)$ whenever they correspond. Now, let $X$ be a non-terminal and let $R(X) = \alpha$. Let $\alpha'$ be a nonempty prefix of $\alpha$, and let $Y$ be its last symbol. Note that for the empty string, a run of size one is always valid, so we can safely consider that the set $M_X$, as declared in line 16 effectively stores nodes properly for the runs that start and end on the same node, and read the empty string. The proof will follow inductively by proving that $M_X$, at the iteration of line 18 that $Y$ was read in line 19, satisfies that $\sem{\bigcal{D}}(M_X[p,q])$ contains all partial annotations from a run of $\bigcal{A}$ over $\operatorname{doc}(\alpha Y)$ that starts on $p$ and ends on $q$; that is, after assuming that in line 24, the set $\sem{\bigcal{D}}(M_X[p,q])$ contained all partial annotations from a run of $\bigcal{A}$ over $\operatorname{doc}(\alpha)$ that starts on $p$ and ends on $q$. We have assumed that the nodes in $M_Y$ store all annotations properly, so the inductive hypothesis is proven from the definition of the operator $\otimes$, as any possible annotation that should be added can be described by two partial runs, one that starts on $p$, reads $\operatorname{doc}(\alpha)$ and ends on $q'$, and a second that starts on $q'$, reads $\operatorname{doc}(Y)$ and ends on $q$, so it is properly added to the $\sem{\bigcal{D}}(M_X[p,q])$ as the operator $\otimes$ does consider this $q'$. Plus, it can be seen that the positions that were obtained from $M_Y$ are properly shifted when added to $M_X$ by assuming the $\operatorname{len}_X$ variable is storing the correct length. From this reasoning it can be concluded that the set $\sem{\bigcal{D}}(v)$, as enumerated in line 7, contains exactly the annotations from $\sem{A}(S)$.

To see that the operators $\mathsf{add}$, $\mathsf{union}$, $\prod$ and $\mathsf{shift}$ are done following all assumptions, we just have to note that $\mathsf{union}$ and $\prod$ are always called by respecting the duplicate-free property. In the case of $\prod$, this is straight forward, since it is only done in line 25, and the positions on the left side (i.e., the ones taken from $M_X$) are all less or equal than the value $\operatorname{len}_X$ has at that point, so by shifting the values from $M_Y$ by this amount, the sets from either side do not contain any annotation in common. In the case of $\mathsf{union}$, one can visually check that in lines 11 and 13, this is done properly, since transitions do not appear more than once, and one can check that in line 25, if any of the calls to $\mathsf{union}$ done inside the $\otimes$ operation is done with duplicates, this implies that there is some annotation that appears in two distinct calls, since each of the $\mathsf{union}$ calls is done over a different index in the matrices $M_Y$ and $M_X$, which contradicts the assumption that $\bigcal{A}$ is unambiguous. From this reasoning we conclude that the operators $\mathsf{add}$, $\mathsf{union}$, $\prod$ and $\mathsf{shift}$ are done properly, and the statement of the theorem follows.
\subsection{Proof of Lemma~\ref{ra:det}}

We use the standard determinization scheme for regular automata. Namely, let $\bigcal{A} = (Q, \Sigma, \Omega, \Delta, q_0, F)$, and let $\bigcal{A}' = (Q', \Sigma, \Omega, \Delta', \{q_0\}, F')$ where $Q' = 2^Q$, $\Delta' = \{(S, a, T) \mid S, T\subseteq Q \text{ and $T$ is the set of states $q$ for which there exists a $p\in S$ such that }(p,a,q)\in \Delta  \} \cup
 \{(S, (a,\oout), T) \mid S, T\subseteq Q \text{ and $T$ is the set of states $q$ for which there exists a } p\in S \text{ such that }$ $(p,(a,\oout),q) \in \Delta \}$, and $F' = \{S\mid S\cap F\neq\emptyset\}$. It can be seen that $\bigcal{A}'$ is deterministic. To show that $\sem{\bigcal{A}}(d) = \sem{\bigcal{A}'}(d)$ for every document $d\in\Sigma^*$, we define the {\it language} of an annotated automaton $\bigcal{A}$ as
 \begin{multline*}
 	 L(\bigcal{A}) = \{\bar{d}\mid \bar{d} \text{ is a document over $\Sigma$ annotated with symbols in $\Omega$ that}\\ \text{is accepted when fed to $\bigcal{A}$ as if it was a regular automaton}\}
 \end{multline*}
It is straightforward to see that $L(\bigcal{A}) = L(\bigcal{A}')$, and from here it can be deduced that for each document $d\in\Sigma^*$ it holds that $\sem{\bigcal{A}}(d) = \sem{\bigcal{A}'}(d)$.

\section{Introduction}\label{sec:introduction}

	\input{intro}

	\section{Setting and main problem}\label{sec:setting}

	\input{setting}

	\section{Enumerable Compact Sets with Shifts}\label{sec:ecs}

	\input{ecsnew}

	\section{Evaluation of annotatated automata over SLP-compressed strings}\label{sec:evaluation}

	\input{eval}

	\section{Applications in Regular Spanners}\label{sec:spanners}

	\input{appl}

	\section{Constant delay-preserving Complex Document Editing}\label{sec:edits}

	\input{edits}

	\section{Future Work}\label{sec:conclusions}

	\input{concl}

	\newpage
	\bibliographystyle{abbrv}


\section{Introduction}
\IEEEPARstart{Q}{uantum} computing is anticipated to outperform classical algorithms in processing speed and impact various industry sectors that need complex computation~\cite{spl01,spl02,KWAK2022}. In contrast to classical computation, where the computation unit (\textit{i.e.}, bit) is deterministically predetermined to be either $0$ or $1$, the computation unit in quantum computing (\textit{i.e.}, qubit) ranges from $0$ to $1$ and this leads an exponential scale in computation capability~\cite{aimlab2022arvix,kong2020review}. It's because the entanglement of qubits enables the simultaneous representation of multiple states. This distinguishing characteristic of quantum computing allows pursuing several paths concurrently in a single qubit, which is physically impossible for classical algorithms and necessitates multiple passes, increasing complexity orders~\cite{aimlab2022icdcs}. Therefore, even in the current decade of noisy intermediate scale quantum computation (NISQ), quantum machine learning (QML) has acquired linear or sublinear complexity as opposed to the polynomial complexity of conventional ML. Since conventional ML depends heavily on massive data, which is extremely hard to analyze and process, QML has drawn attention as a practical solution to these challenges. Various research has been conducted to utilize the nature of quantum computing on QML fully. For instance, a classification task, which is one of the representative machine learning problems, can be solved by QML based classifier~\cite{blank2020quantum}. In addition, QML can be used with not only itself but classical neural networks~\cite{schuld2019quantum}. Previous research showed that quantum computing performs complex computations in Hilbert space more efficiently than classical computing~\cite{schuld2019quantum,aimlab2022dynn}.
\begin{figure}[t!]
\centering
\includegraphics[width=0.75\linewidth]{Feature.PNG}
\caption{
Extracted feature maps according to various train strategies.}\label{fig:Featuremap}
\end{figure}
However, there is still a challenging problem that QML faces,~\textit{i.e.}, barren plateaus. The barren plateaus are a notorious problem in QML that occurs when the number of qubits increases. The barren plateaus vanish the gradients of the QML, making it impossible to guarantee trainability~\cite{mcclean2018barren}. As a solution, research in~\cite{pesah2021absence} proved that utilizing a quantum convolutional neural network (QCNN) with proper initialization can reduce the barren plateaus. Recent research in~\cite{DBLP:journals/qmi/HendersonSPC20} designed the QCNN with filters to extract the features of input data like classical convolutional neural networks (CNN). However, the QCNN is harsh to scale up the number of features because each feature is extracted by measurement of each qubit. Inspired by it, this paper proposes a new version of QCNN, \textit{scalable quantum convolutional neural network} (\textit{sQCNN}), and a new training algorithm,~\textit{reverse fidelity-train} (\textit{RF-Train}), which utilizes the concept of fidelity,~\textit{i.e.}, the nature of quantum computing. Therefore, sQCNN can fully use the intrinsic features of the input data. This paper has the following contributions.
\begin{itemize}
    \item First of all, we propose a scalable QCNN architecture with quantum computing, \textit{i.e.}, sQCNN, to achieve the scalability of filters while avoiding barren plateaus by maintaining QCNN architecture.

    \item In addition, we propose an sQCNN training algorithm, which is named to \textit{RF-Train}, in order to extract the intrinsic features with finite filters.

    \item Lastly, we conduct data-intensive experiments to corroborate the superiority of sQCNN with RF-Train in MNIST and FMNIST datasets, widely used in the literature.
\end{itemize}

Fig.~\ref{fig:Featuremap} shows each feature map of untrained QCNN, Vanilla-Trained QCNN, and RF-trained sQCNN, respectively. We describe the classification performance and Euclidean distance of each feature from these models in Sec.~\ref{sec:experiment}.

\section{Preliminaries of sQCNN}

\subsection{Classical Convolutional Neural Network (CNN)}
A classical CNN is mainly composed of 3 procedures, \textit{i.e.}, convolution layers, pooling layers, and fully connected layers.

\BfPara{Convolution Layers} In convolution layers, input data are convolved by a set of filters. Each filter is designed to extract the intrinsic feature of the input data. The number of filters are able to be adjusted under the consideration of elapsed calculation time in each layer in this classical CNN. The output is called feature map of this convolution computation.

\BfPara{Pooling Layers}
In pooling layers, dimensionality reduction is conducted on the convolved data. This procedure is essential as it reduces the computation time on next convolution layer. Moreover, it gives to the CNN the ability to learn representation invariant to small translations.

\BfPara{Fully Connected Layers}
After the computational procedure of convolution layers and pooling layers, conducting fully connected layers on input data enables the model to get the probability of initial input to belong to the corresponding class.

\begin{figure*}[t!]
\centering
\includegraphics[width=0.85\linewidth]{MA.PNG}
\caption{Comparison between QCNN with Vanilla-Train and sQCNN with RF-Train.}\label{fig:Model ARchi}
\end{figure*}

\subsection{Quantum CNN (QCNN)}\label{sec:quanvolution}

QCNN is a new design of CNNs for multi-dimensional vectors using quantum circuits as convolutional filters~\cite{DBLP:journals/qmi/HendersonSPC20}. To fully leverage massively parallel computations on the superposition of quantum states with a limited number of qubits, QCNN is designed to mimic the role of CNN. To make QCNN operate as CNN, spatial information should be processed by the quantum circuit adequately. The quanvolution filter is the primary method to utilize spatial information. The quanvolutional filter consists of three components,~\textit{i.e.}, encoder, parameterized quantum circuit (PQC), and measurement. By designing the filter with trainable parameters, the quanvolutional filter can be regarded as an extension of the filter in classical CNN. The architecture of quanvolutional filter is illustrated in Fig.~\ref{fig:Model ARchi}.

\BfPara{Encoder} The encoder in quanvolutional filter encodes classical information $\mathbf{x}$ into the state information of qubits. There are many encoding strategies, \textit{i.e.}, basis encoding, amplitude encoding, and angle encoding. In angle encoding, the state is denoted as follows,
\begin{equation}
    |\psi_\mathbf{x}\rangle = \prod_{m=0}^{ \ceil{ \text{size}(\mathbf{x}) / n_q} }U(\theta_m)U(\mathbf{x}_{n_q m:n_q(m+1)})|0\rangle^{n_q},
\end{equation} where $U(\theta_m)$ is the rotation gate of the $m$-th set of re-uploaded input described well in previous research~\cite{perez2020data}, and $n_q$ is the number of qubits in PQC. The terms $\text{size}(\mathbf{x})$ and $\mathbf{x}_{n_q m:n_q(m+1)}$ denote the vector size of input vector $\mathbf{x}$ and the vector which is composed by the first $n_q \cdot m $ to $n_q \cdot (m+1)$ elements, respectively. In this paper, we encode the states of the filters using this angle encoding strategy.

\BfPara{Parameterized Quantum Circuit} As with the universal approximation theorem, there is always a PQC that can represent the optimal objective function within a small error~\cite{benedetti2019parameterized}. Accordingly, QCNN uses PQC with trainable parameters as a filter and proposes to train the filter according to each data-driven task. Using the \textit{chain rule} and \textit{parameter shift rule}, we can write the derivative of the loss function in terms of the output of PQC~\cite{wierichs2022general}.

\BfPara{Measurement}
The measurement procedure of QCNN is the same as the QML. The output state of each n-qubit filter after training on PQC, $|\psi_l\rangle$, can be measured by a set of projection matrices $\mathbf{M}_i$, which form an orthogonal set,~\textit{i.e.}, $\mathbf{M}_i\mathbf{M}_j=0$ where $i\neq j$. If the set satisfies the requirement $\sum_i \mathbf{M}_i =\mathbf{1}$,
then any observable $\mathbf{M}$ has the spectral decomposition $\mathbf{M}=\sum_i i\mathbf{M}_i$. The possible outcomes correspond to the eigenvalues $i$ of $\mathbf{M}$. In this paper, we use the measurement operator $\mathbf{M}=\mathbf{I}^{\otimes n-1}\otimes \mathbf{Z}$, where $\mathbf{I}$ denotes the identity matrix and $\mathbf{Z} = \begin{bmatrix}
1 &  0 \\
0 &  -1
\end{bmatrix}$. Then, the measurement can be denoted as
$\langle O \rangle =\langle \psi_{l} | \mathbf{I}^{\otimes n-1}\otimes \mathbf{Z} |\psi_l\rangle$.
By pooling and fully connecting the measurement $\langle O \rangle$ on the FCN, the QCNN can classify the image classes.

\section{Training of sQCNN}
\subsection{Architecture of sQCNN}
Fig.~\ref{fig:Model ARchi} illustrates the architectural difference between QCNN and sQCNN. In classical QCNN with Vanilla-Train, each pixel data is encoded in each qubit. After feature extraction, we can obtain a feature corresponding to the measurement value of each qubit of PQC. Note that, the number of features extracted in this process is equal to the number of qubits in PQC, which is the same as the size of the kernel. Suppose a single quanvolutional layer and a square filter with length $M$. In this case, the number of features (\textit{i.e.,} equal to the number of channels) extracted through the QCNN is fixed to $n_q = M^2$, where $n_q$ is the number of qubits in the filter. In contrast to QCNN, we design sQCNN to be able to increase the number of filters. We assume that the increase of the number of filters outperforms the increase of qubits in a single filter (\textit{i.e.,} we corroborate this assumption via experiments in Sec.~\ref{sec:experiment}). The sQCNN utilizes multi-filters, and this enables sQCNN to adjust the number of extracted features. The number of features is denoted as $n_{f}n_{q}$, where $n_{f}$ is the number of filters.

\begin{algorithm2e}[t]\label{alg: RF-train}
    \SetCustomAlgoRuledWidth{0.44\textwidth}
\caption{Reverse Fidelity Train (RF-Train)}\label{alg:Training algorithm}
\textbf{Initialization.} \texttt{sQCNN} parameters, $w$;

 \For{ $e = \{1,2,\dots, E\}$}
 {
     \For{ $(x, y) \in \zeta^k$}
     {
         \For{$l, l^{'} \in \{1,2,\dots,L-1\}$}
         {
            Get features with $l$-th and $l^{'}$-th filter\;
            Calculate $\mathcal{L_{RF}}$\;
            Calculate loss gradients\;

         }
         Calculate $\mathcal{L}_e^k \leftarrow \mathcal{L}_{total}$\;
         $\ensuremath{\bm{\theta}}^k_{e+1}\leftarrow \ensuremath{\bm{\theta}}^k_{e}-\eta_e\nabla_{\theta^k_{e}}\ensuremath{\mathcal{L}}^{k}_{e}$\;
     }
 }
\end{algorithm2e}

\subsection{Reverse Fidelity Train (RF-Train)}

As sQCNN can adjust the number of features and measure each PQC-based quanvolutional filter as a set of features, we aim to utilize the architectural advantage of sQCNN fully. Motivated by the fidelity in quantum computing theory, we propose an sQCNN training algorithm named reverse fidelity-train (RF-Train). The RF-Train contains an RF regularizer that adjusts the fidelity between the quanvolutional filters of sQCNN. Note that fidelity is a nature of quantum computing that is a similarity metric between two quantum states. Here, it is possible to measure fidelity because we scale up the sQCNN by increasing the number of filters, not the number of qubits in a single PQC. Suppose two filters have the output states $|\psi_{l} \rangle$ and $|\psi_{l'}\rangle$, respectively. The fidelity between the two quanvolutional filters is denoted as $\Phi(\rho_{l}, \rho_{l'}) = |\langle \psi_{l}|\psi_{l'}\rangle|^2$, where
$\rho_{l} = |\psi_{l} \rangle \langle \psi_{l}| \textrm{ and }
\rho_{l'} =|\psi_{l'} \rangle \langle \psi_{l'}|$. The increasing similarity between the two quanvoltutional filters drives the fidelity to converge to 1. On the other hand, the fidelity converges to 0 when the similarity between the two filters decreases, indicating that the $l$-th filter does not follow the $l'$-th filter.
We assume that a reduction in the fidelity between output states of the quanvolutional filters enables the extraction of various intrinsic features. We corroborate this assumption with numerical experiments in Sec.~\ref{sec:experiment}. We define the RF regularizer as,
\begin{equation}
    \mathcal{L_{RF}} = 1 - {1\over L(L-1)} \sum_{l=1}^L\nolimits \sum^L_{l'\neq l}\nolimits \Phi(\psi_{q_{l}}, \psi_{q_{l^{'}}}),
\end{equation}
where $L$ is the number of filters. With the RF regularizer, the training procedure of sQCNN is described in Algorithm 1. The parameters ($\mathbf{x}$, $y$) are denoted as the input data and label, respectively. We adopt cross-entropy as,
\begin{equation}
\mathcal{L_{CE}} = -{1\over C}\sum^C_{c=1}\nolimits \log{p(y_{pred} = y_c|\mathbf{x})},
\end{equation}
where $C$ represents the number of classes. $y_{pred}$ and $y_c$ are the predicted class and the actual class, respectively. Consequently, we design the total loss of sQCNN as,
\begin{equation}
    \mathcal{L}_{total} ={\frac{1}{|\zeta|}} \sum_{(\mathbf{x}, y) \in \zeta}\nolimits [\mathcal{L_{CE}} -\lambda \mathcal{L_{RF}}],
\end{equation}
where $\zeta$, and  $\lambda$ denote the minibatch and an RF regularizer parameter, respectively.
\section{Performance Evaluation}\label{sec:experiment}
\subsection{Experimental Setting}
To corroborate the performance of the sQCNN with RF-Train, we design the experiments as follows:
\begin{itemize}
    \item We investigate the top-1 accuracy of sQCNN with various RF-regularizer parameters $\lambda$, and QCNN with Vanilla-Train on both MNIST and FMNIST datasets.
    \item To corroborate the impact of the RF regularizer parameter $\lambda$, we investigate the Euclidean distance between the extracted features due to the RF-regularizer parameter $\lambda$.
    \item The scalability of sQCNN is proven. In QCNN, an increase in qubits caused a barren plateau, which degraded the overall model performance. In contrast to QCNN, in sQCNN, the increase in filters results in a performance improvement, despite using the same number of qubits in QCNN. Note that the sQCNN shows performance improvement with the increased number of filters.
\end{itemize}

\subsection{Experimental Results}
\BfPara{Performance of sQCNN}
Fig.~\ref{exp:main} and Fig.~\ref{exp:main2} represent top-1 accuracy of various models with two filters on the MNIST and FMNIST datasets, respectively. Analyzing from the architectural point, both Fig.~\ref{exp:main} and Fig.~\ref{exp:main2} represent that RF-Trained sQCNNs outperform classical Vanilla-Trained QCNN. From the results, we confirmed that even using sQCNN, \textit{i.e.}, even without RF-Train, it can achieve performance improvement in a classification task. Moreover, when the RF-regularizer parameter $\lambda$ increases, sQCNN with a finite number of filters ($n_f =2$) shows performance improvement. sQCNN with high RF-regularizer parameter ($\lambda= 0.5)$ achieves $16 \%$ higher top-1 accuracy than sQCNN which does not utilize RF-Train ($\lambda=0$). Here, we observe that by diversifying the filters with RF-Train, we can improve the performance of sQCNN.
\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{accuracy-eps-converted-to.pdf}
    \caption{Top-1 accuracy (MNIST).}
    \label{exp:main}
\end{figure}
\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{bar-eps-converted-to.pdf}
    \caption{Top-1 accuracy on FMNIST dataset according to the number of filters (qubits). sQCNN ($\lambda = 0, 0.1$ and $0.5$) with filters ($n_f = 2, 3$ and $4$) and QCNN with number of qubits ($n_q =8, 12$ and $16$). }
    \label{exp:main3}
\end{figure}

\BfPara{Impact of RF-Train}
Table~\ref{tab:distance} represents the impact of RF regularizer parameter $\lambda$ on the classification performance. We observe that RF-Train increases the Euclidean distance between extracted features and experimentally confirmed that this diversity between features improves the classification task performance of sQCNN.

\BfPara{Scalability of sQCNN}
Fig.~\ref{exp:main3} represents that sQCNN achieves scalability. In QCNN, the top-1 accuracy decreases significantly as the number of qubits in a filter increases. As the number of qubits increases from 8 to 16, the performance of QCNN drops about $30\%$. In contrast to QCNN, sQCNN shows a stable performance when the number of filters increases. From the result, we corroborate that this scalability of SQCNN can be a significant characteristic in several tasks that require a large number of filters as well as simple MNIST and FMNIST data with one input channel.

\section{Conclusions and Future Work}
This paper proposes a scalable QCNN (sQCNN) architecture and a novel training algorithm (RF-Train) that enables sQCNN to diversify the extracted features. To achieve scalability while avoiding the barren plateaus which occur when the number of qubits in the filter increases, we utilize multiple filters with a finite number of qubits. To extract various features with the filters and maximize the performance of sQCNN, motivated by the quantum theory, we design an RF regularizer using the concept of fidelity. With extensive experiments, we corroborate the diversity of features that are extracted by using RF-Train and the scalability of sQCNN.

As future research, analyzing the scalability and trainability of sQCNN in various environments could be interesting.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth]{fmnist_acc-eps-converted-to.pdf}
    \caption{Top-1 accuracy (FMNIST).}
    \label{exp:main2}
\end{figure}

\begin{table}[t]
\centering
\caption{Top-1 accuracy and Euclidean distance comparison }
\resizebox{\columnwidth}{!}{
\begin{tabular}{c||ccc}
    \toprule[1pt]
     &\multicolumn{3}{c}{\textbf{RF-Train}}\\
    \textbf{Metric} & $\lambda=0$ & $\lambda=0.1$   &  $\lambda=0.5$ \\\midrule
    Top-1 accuracy (\%) & $76$& $78$& $\mathbf{82}$ \\
    Euclidean distance ($\times 10^{-2}$) & $0.4$& $0.7$& $\mathbf{1.1}$ \\
\bottomrule[1pt]
\end{tabular}}\label{tab:distance}
\end{table}

\bibliographystyle{IEEEtran}


\section{Introduction} \label{sec:IN}

Load frequency control (LFC) is one of the most important control problems in power system operations. The objective of LFC is to maintain the frequency of each area in an interconnected power system by adjusting the output of generators with automatic generation control (AGC) regulator or excitation controller \cite{shashi}. The LFC has been studied in various research to cope with conventional generators \cite{wen}, distributed energy resources \cite{zhu}, and electric vehicles (EVs) \cite{falahati}. However, most research considers a centralized framework, implying that a centralized dispatch center controls all generators \cite{kundur}. As a modern power system relies more on distributed generation (DG) and requires resilience to cyberattacks, the existing centralized control paradigm with one-point-failure faces many challenges. Distributed or decentralized frameworks can be considerable frameworks to improve the stability and security of the power system. Consequently, there is a lot of research that considers distributed or decentralized LFC problems, assuming the limited information exchanges among the interconnected areas \cite{ting, zhao, wang, andrade}. Especially in the case of peer-to-peer (P2P) based LFC, generation control in one area is determined by only the information from the areas connected in the information-exchange graph \cite{ting}.

The general LFC problem can be represented as a linear quadratic regulator (LQR) problem, which minimizes the frequency deviations and other factors such as power outputs, power inflow between interconnected areas, and control efforts such as AGC control signal \cite[Ch.~2]{bevrani}. When the model is known, this problem can be easily solved by adopting the Algebraic Riccati equation (ARE) \cite{ku} or by applying gradient-based methods \cite{bu,li}. However, finding the feedback gain will become complicated as we consider the model uncertainties, the constraints on the optimization problem, or the structure of the feedback gain. This paper aims to solve the P2P-based LFC problem, while dealing with these three challenges at once by considering uncertainty from the environment, the risk constraints, and structured output feedback.

In particular, reinforcement learning (RL) will be applied to solve this problem. There are two advantages that RL has \cite{sutton}. First, RL is model-free learning, i.e., we do not need to consider the model parameters. Instead, we will generate the trajectory, observe the reward and update the controller toward increasing the total reward. Second, RL is data-driven learning, i.e., instead of generating a probability distribution for the uncertainty and performing the Monte-Carlo method, we can directly use the data gathered from the grid and use it to train the controller. Recent research works such as \cite{sayak_arxiv, eth_distributed_learning, mukh_tac, KB} look into various distributed aspects of RL.

The main contributions of the paper are three-fold. First, we formulate the constrained LQR problem for LFC considering the mean-variance risk constraint and output structured feedback. Second, an RL-based Stochastic Gradient Descent max-oracle (SGDmax) algorithm with a zero-order policy gradient (ZOPG) is proposed to solve the formulated problem and find the feedback controller. Here, the learning process is applied to the emulator grid which approximates the dynamics of the physical grid. Last, the robustness test on the obtained feedback controller is checked to ensure that the controller trained from the emulator grid has enough robustness that can be directly applied to the physical grid.

The rest of this paper is organized as follows. In Section~\ref{sec:PS}, problem setup with two models is introduced along with the mathematical formulation of the constrained LQR problem including the objective function, risk constraints, and linear dynamics. Section~\ref{sec:ALG} introduces the SGDmax algorithm with ZOPG to estimate the gradient and find a feedback gain by a gradient-descent method. Section~\ref{sec:SR} shows the effectiveness of the proposed algorithm by adopting it to a load frequency control (LFC) problem. Finally, the paper is wrapped up in Section~\ref{sec:CON}.

\section{Problem Description} \label{sec:PS}
The objective is to find an optimal feedback controller $K$ of the physical system, represented as Model~2 in Fig.~\ref{fig:setup}. {\KBB We also estimate the parameters of the system from the historical data and construct Model~1, which is the offline emulator based on the actual physical system.} Here, we formulate offline learning based on uncertainty data from the physical system. From the initial feedback gain $K^0$, we will observe the output $y$ with on-policy learning, meaning that the controller will be directly applied to the system with the uncertainty. The output result will be used to update $K$ by gradient-based update. By updating $K$ iteratively, the solution will approach the optimal feedback gain.

Note that the setup can be extended to the case with the nonlinear dynamics-based physical system in Model~2. Even though the system has nonlinear dynamics, we can construct Model~1 as a linear system by estimating the linear parameters and finding $K$ that can show the best result. The obtained $K$ may not be an optimal feedback gain due to the approximation from the nonlinear to linear dynamics, but it will show a decent result with much less computation burden by directly finding the feedback controller from the nonlinear formulation.

After we obtain the trained $K$ from Model~1, we adopt it to Model~2 to run on the physical system. Here, the robustness of the feedback gain is a key factor as the parameters between Model~1 and Model~2 may have errors due to the approximation. If the robustness is not high enough, the feedback controller will not work properly in Model~2 and thus needs a very precise approximation in Model~1. On the other hand, high robustness refers to stability in Model~2 even though there are some errors in approximations.
\begin{figure}[t]
	\centering
	\includegraphics[width=0.9\linewidth]{fig_setup.png}
	\caption{Problem setup with emulator (Model 1) and physical system (Model 2)}
	\label{fig:setup}
\end{figure}
To formulate the algorithm, we start with the assumption of the emulator system with linear dynamics as:
\begin{align}
    &\dot{x} = Ax + Bu + w, x(0) = x_0\\
    & y = Cx
\end{align}
where $x\in \mathbb{R}^n, y \in \mathbb{R}^p, u \in \mathbb{R}^m$. We want to find static output feedback control:
\begin{align}
    u = -Ky = -KCx
\end{align}
that minimizes
\begin{align}
	R_0(K) &= \lim_{T \to \infty}\frac{1}{T}\mathbb{E} \sum_{t=0}^{T-1} [y^\top_t Q y_t + u^\top_t R u_t]  \label{eq:opt}
\end{align}
considering the mean-variance risk constraint:
\begin{align}
	R_c(K) = \lim_{T \to \infty}\!\frac{1}{T} \mathbb{E} \sum_{t=0}^{T-1} \left(y_t^\top Q y_t - \mathbb{E}[y_t^\top Q y_t \vert h_t]\right)^2 \leq \delta  \nonumber
\end{align}
Here, we consider the structured $K$ which is within the set $\mathcal{K}$ defined as
\begin{align}
	&\mathcal{K} = \{K: K_{a,b} = 0 \;\text{if and only if}\; (a,b) \notin \mathcal{E}) \}
\end{align}
where the structure pattern $\mathcal{E}$ is specified by the edges of a given communication or information-exchange graph. In other words, $K$ will have a sparse structure with the element as zero if the corresponding edges cannot exchange information.

By following \cite{tsiamis}, we can reformulate the constraint as
\begin{align}
	R_c(K) \!=\!  \lim_{T \rightarrow \!\infty} \frac{1}{T} \mathbb{E} \sum_{t=0}^{T-1} \left(4y_t^\top QCWC^\intercal Q^\intercal y_t+4y_t^\top Q M_3^y \right)\!\leq\!\bar{\delta}\! \label{eq:objr2}
\end{align}
with $\bar{\delta} = \delta-m_4^y + 4 \text{tr}\{(WQ_c)^2 \}$ and the (weighted) noise statistics given as
\begin{align}
	\bar{w}&=\mathbb{E}[w_t],\\
	W&=\mathbb{E}[(w_t-\bar{w})(w_t-\bar{w})^\top],\\
	M_3^y &= \mathbb{E}[(w_t-\bar{w})(w_t-\bar{w})^\top Q_c (w_t-\bar{w})],\\
	m_4^y&=\mathbb{E}[(w_t-\bar{w})^\top Q_c (w_t-\bar{w}) - \text{tr}(WQ_c)]^2.
\end{align}
with $Q_c = C^\intercal Q C$. Note that \eqref{eq:objr2} no longer depends on the past trajectory $h_t$, and has a quadratic form which is the same form as the objective function \eqref{eq:opt}.

To solve the optimization problem with the objective function as \eqref{eq:opt} and the risk-constraint as \eqref{eq:objr2}, we will formulate the Lagrangian function as below:
\begin{align}
\mathcal{L}(K, \lambda) =&R_0(K)+ \lambda [R_c(K)-\delta]\nonumber\\
 =&\lim_{T \rightarrow \infty} \frac{1}{T}\mathbb{E}\sum_{t=0}^{T-1} [y^\top_t (Q+4\lambda QCWC^\top Q^\top) y_t  \nonumber\\
 &+ 4\lambda y_t^\top Q M_3^y + u^\top_t R u_t]-\lambda\overline{\delta} \label{eq:lag}
\end{align}
By assuming that the problem is feasible and thus the dual variable $\lambda$ is finite, we can formulate the dual function with the bounded set $\mathcal{Y}:= [0,\Lambda]$ as
\begin{align}
\max_{\lambda \in \mathcal{Y}} \mathcal{D}(\lambda) = \max_{\lambda \in \mathcal{Y}}  \min_{K \in \mathcal{K}} \mathcal{L} (K, \lambda) \label{eq:dual}
\end{align}
Here, instead of \eqref{eq:dual}, we will consider the minimax counter problem:
\begin{align}
\min_{K\in\mathcal{K}}\Phi(K)~~\textrm{where} ~~\Phi(K):= \max_{\lambda \in \mathcal{Y}} \mathcal{L}(K, \lambda) \label{eq:dual_c}
\end{align}
Though \eqref{eq:dual} and \eqref{eq:dual_c} are different optimization problems, solving \eqref{eq:dual_c} will give the stationary point (SP) of \eqref{eq:dual} since the KKT conditions for \eqref{eq:dual} are related to the SP of \eqref{eq:dual_c}. Note that we can easily find $\lambda$ that maximizes $\mathcal{L}(K,\lambda)$ by choosing $\lambda = 0$ if the constraint is satisfied and $\lambda = \Lambda$ otherwise, since $\mathcal{L}(K,\lambda)$ is a linear function of $\lambda$.\\

\section{SGDmax Algorithm with ZOPG} \label{sec:ALG}

To solve \eqref{eq:dual_c}, we will adopt the SGDmax algorithm with ZOPG method to estimate the gradient \cite{KB}. Algorithm~\ref{alg:ZOO} shows how to find the zero-order policy gradient to update $K$ \cite{ZOO}. Specifically, we will choose the perturbation $U \in \mathcal{S}_\mathcal{K}$, with $\mathcal{S}_\mathcal{K}$ is the set of matrices that has the same structure in $\mathcal{K}$ with its magnitude being 1. After adopting the perturbation $rU$ with smoothing radius $r$ to current feedback gain $K$, the optimal $\lambda$ that maximizes the Lagrangian function value is chosen. By adopting $\lambda$, we can compute the estimated zero-order gradient. Note that the gradient is computed as $U$ multiplied by a scalar value, meaning that the estimated gradient follows the structure $\mathcal{K}$.
\begin{algorithm}[t]
	\SetAlgoLined
	\caption{Zero-Order Policy Gradient (ZOPG)}
	\label{alg:ZOO}
	\DontPrintSemicolon
	{\bf Inputs:} smoothing radius $r$,

	the policy $K$ and its perturbation  $U \in \mathcal{S}_{\mathcal{K}}$, both of $n_\mathcal{K}$ non-zeros.\;
	Obtain $\lambda' \leftarrow \argmax_{\lambda \in \mathcal{Y}} \mathcal{L}(K+rU, \lambda)$;\;

Estimate the gradient $\hat{\nabla}_{\mathcal{K}} \mathcal{L} (K; U) = \frac{n_\mathcal{K}}{r} \mathcal{L}(K+rU, \lambda') U$.\;
	{\bf Return:} $\hat{\nabla}_{\mathcal{K}} \mathcal{L} (K; U)$.
\end{algorithm}
With the gradient from Algorithm~\ref{alg:ZOO}, we adopt the SGDmax algorithm introduced in Algorithm~\ref{alg:SGD} to find SP $K$ using gradient method. From the zero-order gradient for each $U_s$, the gradient is averaged and used to update $K^j$ in the iteration $j$ with gradient descent. Note that at the initialization, the set $\mathcal{G}^1$ and regarding parameters $L_0, \ell_0$ and $\rho_0$ are defined as
\begin{align}
&\mathcal{G}^1 := \{K\in\mathcal{K} |\Phi_{\mu_0}(K)\leq 10~\Phi_{\mu_0}(K^0)\}, \nonumber \\
&L_0\!:=\!\sup\limits_{K\in\mathcal{G}^0}{L_K}, \ell_0\!:=\!\sup\limits_{K\in\mathcal{G}^0}\ell_K. ~\textrm{and}~ \rho_0\!:=\!\inf\limits_{K \in \mathcal{G}^0} \min\{\beta_K, \psi_K\} \label{eq:const}
\end{align}
with Moreau envelope $\Phi_{\mu_0}(k)$ of $\Phi(K)$ defined by solving the minimization problem
\begin{align}
    \Phi_{\mu_0}(K) := \min_{K'\in\mathcal{K}} \Phi(K') + \frac{1}{2{\mu_0}}\Vert K'-K\Vert^2, ~~\forall K \in \mathcal{K}, \label{eq:mor}
\end{align}
where $\mu_0 = 1/(2\ell_0)$. Here, the sublevel set $\mathcal{G}^1$ ensures the Lipschitz and smoothness properties for any $K \in \mathcal{G}^1$.{\KB Note that smoothing radius $r$, stepsize $\eta$, and the minimum number of iteration $J$ are determined by following \cite[Theorem 2]{KB} which is refered here as Lemma~\ref{lem:proof}.

\begin{lemma}
With an initial $K^0 \in \mathcal{K}$ and a given $\epsilon > 0$, by setting the parameters as
	\begin{align}
		&r \leq \min\Big\{\rho_0, \frac{L_0 \sqrt{M}}{\ell_0} \Big\}, \; \eta \leq \frac{\epsilon^2}{\alpha\ell_0(L_0^2+\ell_0^2r^2/M)},\nonumber\\
		\textrm{and}~&J = \frac{2\sqrt{10\alpha}\Phi_{\mu_0}(K^0)}{\eta \epsilon^2}\label{eq:thm2}
	\end{align}
with $\mathcal{G}^1, L_0, \ell_0$ and $\rho_0$ following \eqref{eq:const}, Algorithm~\ref{alg:SGD} converges to $\epsilon$-SP $K_\epsilon$ that satisfies $\Vert \Phi_{\mu_0}(K_\epsilon)\Vert \leq \epsilon$, with probability of at least $(0.9-\frac{4}{\alpha}-\frac{4}{\sqrt{10\alpha}})$.
\label{lem:proof}
\end{lemma}
}
\begin{algorithm}[t]
	\SetAlgoLined
	\caption{Stochastic Gradient Descent with max-oracle (SGDmax)}
	\label{alg:SGD}
	\DontPrintSemicolon
	{\bf Inputs:} A feasible and stable policy $K^0$, upper bound $\Lambda$ for $\lambda$, threshold $\epsilon$, and number of ZOPG samples $M$.\;
		Determine $L_0, \ell_0$, and $\rho_0$ with the set $\mathcal{G}^1$ and compute $r$, $\eta$, and $J$; \;
	\For{$j = 0, 1, \ldots, J-1$}{
		\For{$s=1,\ldots,M$}{
		Sample the random $U_s \in \mathcal{S}_\mathcal{K}$; \;
		Use Algorithm~\ref{alg:ZOO} to return $\hat{\nabla}\mathcal{L}_\mathcal{K}(K^j;U_s)$.
		}
		Update $K^{j+1} \leftarrow K^{j}-\eta \left(\frac{1}{M} \sum_{s=1}^M \hat{\nabla} \mathcal{L}(K^j;U_s)\right)$.}
	{\bf Return:} the final iterate $K^{J}$.
\end{algorithm}

{\KB \noindent Proof Sketch: Similar to \cite{KB}, the sketch of the proof for Lemma~\ref{lem:proof} is as follows.} From the result of recent works \cite{bu, malik, harvard}, the Lipschitz and smoothness properties of $\mathcal{L}(K,\lambda)$ are still valid by \cite[Lemma~1]{KB}. {\KBB Note that the output feedback involves only the linear combination from the original states, thereby such functional properties are maintained.} Therefore, we can bound the expectation of $\Phi_{\mu_0}(K^{j+1})$ from the expectation of $\Phi_{\mu_0}(K^{j})$ as below;
\begin{align}
	\mathbb{E}&\left[\Phi_{\mu_0}(K^{j+1})\right]\leq \mathbb{E}\left[\Phi_{\mu_0}(K^{j})\right]\nonumber\\
	& - \frac{\eta}{4} \mathbb{E} \Vert \nabla \Phi_{\mu_0}(K^{j})\Vert^2 + \eta^2 \ell_0 \left(L_0^2+{\ell_0^2 r^2}/{M}\right) \label{eq:pbd2}
\end{align}
Here, we can claim the convergence to an SP by showing these two aspects. First, we will verify the SP condition, i.e. $\Vert \nabla \Phi_{\mu_0}(K) \Vert \leq \epsilon$ in the expectation sense as the iteration continues. Second, we will verify that the probability that the algorithm fails to converge to the SP is bounded.
First, we can verify the SP condition by summing up \eqref{eq:pbd2} over iterations $j=0, ..., J-1$,
\begin{align}
    &\frac{1}{J}\sum_{j=0}^{J-1}\mathbb{E}\left[\Vert \nabla\Phi_{\mu_0}(K^j)\Vert^2\right] \nonumber\\
    \leq&\frac{4\left[\Phi_{\mu_0}(K^0)\!-\!\mathbb{E}[\Phi_{\mu_0}(K^J)]\right]}{J\eta}+\!4\eta\ell_0(L_0^2+\ell_0^2r^2/M),
\end{align}
which gives the upper bound on order of $\epsilon^2$ by choosing the stepsize $\eta=O(\epsilon^2)$.
Second, the probability that $\{K^j\}$ exceeds $\mathcal{G}^1$ is represented as in \eqref{eq:prob},
\begin{align}
    &\mathbb{P}\left( \frac{1}{J} \sum_{j=0}^{J-1} \Vert \nabla\Phi_{\mu_0}(K^j)\Vert^2 \geq \epsilon^2 \right) \nonumber\\
    = & \mathbb{P}\left( \frac{1}{J} \sum_{j=0}^{J-1}\Vert \nabla\Phi_{\mu_0}(K^j)\Vert^2 \geq \epsilon^2, \tau > J \right) \nonumber\\
    + &\mathbb{P}\left( \frac{1}{J} \sum_{j=0}^{J-1}\Vert \nabla\Phi_{\mu_0}(K^j)\Vert^2 \geq \epsilon^2, \tau \leq J \right), \label{eq:prob}
\end{align}
where $\tau := \min\{j\geq 0: K^j \notin \mathcal{G}^1 \}$, meaning the iteration index that $\{K^j\}$ exceeds $\mathcal{G}^1$.
Here, we can bound the first and second terms on the right-hand side by following the Markov's inequality and the Doob's maximal inequality with the supermartingale property. The detailed proof can be found at \cite[Appendix~B]{KB}.

\section{Simulation Study} \label{sec:SR}

To show the effectiveness of the proposed method, we consider the LFC in the multi-area power system. Assume there are $N$ number of radially connected areas through tie lines depicted in Fig.~\ref{fig:network} in case $N=6$ as an example. The corresponding parameters are indicated in \cite[Table I]{DLQR}.
\begin{figure}[t]
	\centering
	\includegraphics[width=0.65\linewidth]{fig_network.png}
	\caption{Power grid with six interconnected areas.}
	\label{fig:network}
\end{figure}
Assume the information-exchange graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ follows the physical interconnection between areas in Fig.~\ref{fig:network}. Each node $i\in \mathcal{V}$ represents an area and an edge $(i,j) \in \mathcal{E}$ between two areas indicates that they are connected through a tie-line.

For the aggregated generators in the area, the open-loop linearized dynamics of the $i$-th area are depicted in Fig.~\ref{fig:dynamics} \cite{bevrani, KB}.
\begin{figure}[t]
	\centering
	\includegraphics[width=\linewidth]{fig_dynamics.png}
	\caption{Block representation of load frequency control in $i$-th area \cite{KB}}
	\label{fig:dynamics}
\end{figure}
Notation $\Delta$ implies the deviation from the steady-state operation. When there is a frequency deviation $\Delta f_i$ happens in $i$-th area, the primary frequency control action $\Delta P_{f,i}$ is made proportional to $\Delta f_i$, as $\Delta P_{f,i} = -(1/R_i) \Delta f_i$. As a secondary control, automatic generation control (AGC) signal $\Delta P_{C,i}$ is determined based on the control design. The total control signal $\Delta u_{tot, i}$ combines primary and secondary control of generators, and eventually construct the signal total power output of generators denoted as $\Delta P_{G,i}$. Third, the unknown power load deviation $\Delta P_{L,i}$ and the total power inflow $\Delta P_{tie,i}$ are added to $\Delta P_{G,i}$ and affects the frequency deviation $\Delta f_i$ of the next time slot. Note that $\Delta P_{tie,i}$ is the summation of power inflow from the interconnected areas that are proportional to the difference in frequency as
\begin{align}
	\Delta P_{tie,i} = \sum_{j\in \mathcal{N}_i} K_{tie,i,j} (\Delta f_i - \Delta f_j)\Delta t,
\end{align}
{\KB where $j\in\mathcal{N}_i$ indicates the set of areas that are connected to $i$-th area, $\Delta t$ means the discrete time interval, and $K_{tie}$ refers to a synchronization coefficient between area $i$ and $j$.}
In addition to the dynamics indicated in Fig.~\ref{fig:dynamics}, we adopt the integral action of Area Control Error (ACE) $z_i$ defined as $z_i = \beta_i \Delta f_i + \Delta P_{tie,i}$, with bias factor $\beta_i$ often being defined as $D_i + 1/R_i$ \cite{DLQR_en}.

The state, control variable, and the unknown random disturbance of $i$-th area as $x_i = [\Delta f_i, \Delta P_{G,i}, \Delta P_{tie,i}, \int z_i]^\top$, $u_t = \Delta P_{C,i}$ and $w_i = \Delta P_{L,i}$ respectively. With an assumption that all areas are identical, the aggregated dynamics of the network can be represented as linear dynamics after removing the subscript $i$:
\begin{align}
	\dot{\tilde{x}} = (I_N \otimes A_1 + \mathcal{L} \otimes A_2)\tilde{x} + (I_N \otimes B_u) \tilde{u} + (I_N \otimes B_w)\tilde{w} \label{eq:agg_ss}
\end{align}

\noindent where $\tilde{x}\!=\!\text{col}(x_1,\cdots,x_N), \tilde{u}\!=\!\text{col}(\Delta P_{C,1},\cdots\!,\Delta P_{C,N}), \tilde{w} = \text{col}(\Delta P_{L,1},\cdots,\Delta P_{L,N})$ and the matrices $A_1, A_2, B_u$ and $B_w$ follow \cite[Section~V]{KB}.

Lastly, the output feedback $\tilde{y}$ is computed as $\tilde{y}=C \tilde{x}$, with $C$ indicating the output that can be observed. Here, since the frequency deviation $\Delta f_i$ is hard to track when considering the aggregated generators in an area, we will observe only the sum of total power output of generators and total power inflow, i.e. $\Delta P_{G,i}$ and $\Delta P_{tie,i}$. Later, we show an experiment when we include the area frequencies in the feedback. Also, we experiment with the design without using the conventional ACE inputs, however, indirectly, the tie-line power flow captures their impact. \color{black} Accordingly, $C$ represents the output that can be observed based on the information-exchange graph. {\KB In Fig.~\ref{fig:network}, the corresponding matrix $C$ will be a $6 \times 24$ matrix, with the elements $C_{i,4j-2}$ and $C_{i,4j-1}$ being 1 if the areas $i$ and $j$ are connected.
Here, each column represents the provided information from the neighborhood areas. For example, area~2 can use the information from area~3 and itself, which corresponds to column 2 that has $C_{2,6}, C_{2,7}, C_{2,9}$ and $C_{2,10}$ as 1.} Here, we adopt the linear control policy $\tilde{u} = -K \tilde{y} = -KC\tilde{x}$, with $K$ following the sparse structure $\mathcal{K}$. It implies that $K_{i,j} = 0$ when node $i$ and $j$ are not connected in the information-exchange graph.

\subsection{Training Result with Model~1}
To find the optimal feedback gain $K$, Algorithm~\ref{alg:SGD} is adopted with a stepsize $\eta = 10^{-4}$, smoothing radius $r=0.1$, and $M=100$ samples for ZOPG. Figure~\ref{fig:obj} shows that LQR objective value converges to a steady-state with a sufficient number of updates.
In addition, Fig~\ref{fig:element} depicts the trajectory of selected elements of $K$ during the updates. Even though there is some fluctuation due to the estimated gradient, all values converge as iteration continues, meaning that $K$ reaches the SP of the objective function.
\begin{figure}[t]
	\centering
	\includegraphics[width=0.95\linewidth]{fig_obj_value.png}
	\caption{LQR objective trajectory.}
	\label{fig:obj}
\end{figure}
\begin{figure}[t]
	\centering
	\includegraphics[width=0.95\linewidth]{fig_element.png}
	\caption{Updates of selected elements of $K$ during learning.}
	\label{fig:element}
\end{figure}
\subsection{Testing Result with Model~2}
With the feedback gain $K$ obtained from the training, test has been made by creating a scenario in a 20-second window. Fig~\ref{fig:load} shows one of the scenario on the changes in the load in each area. To show the effectiveness of the control, each area has one large random changes during a 20-second time period. For example, area~3 has huge increase in the load around 3 seconds according to Fig~\ref{fig:load}.
\begin{figure}[t]
	\centering
	\includegraphics[width=0.95\linewidth]{fig_load.png}
	\caption{Load deviations of areas}
	\label{fig:load}
\end{figure}
{\KB Fig~\ref{fig:freq}
shows the frequency deviation
of Area~2 and Area~3 with the obtained $K$ from the training. Note that Area~2 and Area~3 are interconnected according to Fig.~\ref{fig:network} with output values representing the sum of total power output of generator and total power inflow. As seen in Fig~\ref{fig:freq}, the frequency changes dramatically when there is a huge change in load as indicated in Fig.~\ref{fig:load}. Frequency deviation is stabled to zero in a short time period meaning that the feedback gain $K$ can successfully control the frequency from the uncertain load changes. Here, the connected area is also affected since they are connected to each other and thus affects the control decision based on the structure of $K$.} {\KBB Moreover, the performance in Fig~\ref{fig:freq} is similar to  Fig.~\ref{fig:freq2}, which includes the frequency deviation in the feedback. It shows that our proposed output feedback works quite well even when we do not have the information on frequency, with the influence of $\Delta P_{G,i}$ and $\Delta P_{tie,i}$ captured through peer-to-peer information exchange.}
\begin{figure}[t]
	\centering
	\includegraphics[width=0.95\linewidth]{fig_freq2.png}
	\caption{Frequency deviation of Area~2 and Area~3}
	\label{fig:freq}
\end{figure}
\begin{figure}[t]
	\centering
	\includegraphics[width=0.95\linewidth]{fig_freq2_f.png}
	\caption{Frequency deviation of Area~2 and Area~3 when including frequency deviation in the feedback}
	\label{fig:freq2}
\end{figure}

\subsection{Robustness Test Result}
As discussed, the parameters of the dynamics in Model~1 may have errors compared to Model~2 due to the estimation or the approximation from the nonlinear dynamics to linear dynamics. Therefore, checking out the robustness of the obtained feedback gain is crucial when adopting the obtained feedback gain to the physical system.

To test the robustness, we observe the frequency deviation of Area~3 with the changes in the parameters by 10\%, 15\% and 20\% of their actual values, respectively. Fig~\ref{fig:robust} shows the result of four cases. As seen in the figure, the controller successfully controls the frequency deviation until 15\% error, but the fluctuation increases in 20\% error case and thus takes much more time to reach the stable frequency. The result implies that even when the estimated parameters in Model~1 have errors with Model~2, the feedback gain from Model~1 can still perform a good frequency control in the physical system in around 15\% of the error.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.95\linewidth]{fig_robust.png}
	\caption{Frequency deviation of Area~3 for robustness test with Model 1}
	\label{fig:robust}
\end{figure}

\section{Conclusions} \label{sec:CON}
The paper presented a learning-based method to solve the P2P-based LFC problem, while considering three challenges including model uncertainty, the mean-variance risk constraint in load uncertainty, and output structured feedback. To train a feedback controller based on the data from a physical grid, an emulator grid that approximates the dynamics of the physical grid was formulated and the learning is performed in the emulator grid. Here, starting from the constrained LQR problem with output structured feedback, we reformulated the risk constraint to have the same form as the LQR objective function and construct the dual problem as a minimax problem. By adopting SGDmax algorithm with an estimated ZOPG, we obtained the trained output feedback that we can adapt to a physical grid with a robustness test. Numerical tests with radially interconnected areas in a power system showed that the feedback controller from learning can successfully control the frequency of the areas while satisfying the risk constraints. In addition, robustness test indicated that in the current setting of the model the obtained feedback controller can be adapted to a physical system with 20\% robustness, meaning that the controller obtained from the approximated emulator grid can be directly adopted to the physical system. {\KBB Future work includes uncertainty from distributed energy resources such as wind/solar or EVs, and adopts the proposed method to real power system models.}