\section*{Introduction}
Non-abelian anyons are interesting for a variety of reasons. For example, they naturally appear in quantum field theory descriptions of knot theory \cite{Witten:1988hf}, they are believed to play an important role in the fractional quantum Hall effect \cite{Moore:1991ks}, and they underly a topological form of quantum computation \cite{wang2010topological}. More recently, they have attracted attention as providing possible lessons for quantum gravity \cite{Rudelius:2020orz}.

In this paper we will be exclusively concerned with the appearance of non-abelian anyons in a particular type of $2+1$-dimensional topological quantum field theory (TQFT): discrete gauge theories \cite{Dijkgraaf:1989pz,Roche:1990hs}. These are gauge theories based on some discrete gauge group, $G$, along with a Dijkgraaf-Witten 3-cocycle, $\omega\in H^3(G,U(1))$ (when $\omega$ is cohomologically non-trivial, the theory is said to be \lq\lq twisted"). The basic degrees of freedom are anyonic line operators (i.e., operators supported on one-dimensional loci of spacetime that have non-trivial braiding with each other) of three general types:
\begin{enumerate}
\item {\it Wilson lines}, which carry electric charge labeled by a linear irreducible representation of $G$, $\pi$. These operators have trivial magnetic charge.
\item {\it Magnetic flux lines} carrying magnetic charge labeled by a conjugacy class, $[g]$, of an element $g\in G$ with $g\ne1$. These operators have trivial electric charge. Depending on the choice of $\omega$, such operators may or may not exist.
\item {\it Dyonic lines} (or simply dyons) carrying a magnetic charge labeled by a conjugacy class, $[g]$, of an element $g\in G$ with $g\ne1$ and an electric charge labeled by an, in general, projective representation of the centralizer of $g$, $N_g$. In the case of an untwisted gauge theory (i.e., $\omega=0\in H^3(G,U(1))$), the representation is linear. We will describe, in some detail, when this statement continues to hold for certain dyons in twisted theories. Dyons are the most generic type of anyons in discrete gauge theories.
\end{enumerate}
As a physical toy model, one can think of dyons as Aharonov-Bohm systems with charges bound to magnetic flux lines \cite{Preskill}.

Our first rather basic observation is that the line operators in discrete gauge theories naturally relate close cousins in group theory: representations (and their characters) to centralizers. Therefore, discrete gauge theory is a natural way to organize and unify ideas in the theory of finite groups.

In what follows, we will focus on the case of finite simple groups. Via group extensions, these are the basic building blocks of all finite groups. The celebrated classification of finite simple groups guarantees that any such group fits into the following categories:
\begin{enumerate}
\item Abelian groups of prime order
\item Alternating groups
\item Lie groups over finite fields
\item Twenty-six sporadic groups
\end{enumerate}

In spite of this complete classification, there are still many open problems involving these groups. Of particular interest to us is the following old conjecture:

\noindent
{\bf Conjecture (Arad-Herzog):} Consider a non-abelian finite simple group, $G$, and non-trivial elements $g,h\in G$. Then,
\begin{equation}
[g]\cdot [h]\ne [gh]~,
\end{equation}
where $[g]$, $[h]$, and $[gh]$ are conjugacy classes of $g, h,$ and $gh$ respectively \cite{AradHerzog}.

\noindent
More pithily, Arad and Herzog (AH) conjectured that in non-abelian finite simple groups, the product of non-trivial conjugacy classes cannot be a single conjugacy class.

As we will argue in section \ref{conjecture}, this conjecture has the following implication (which we then prove in section \ref{proofs}):

\noindent
{\bf Theorem 1:} In a (twisted or untwisted) $2+1$-dimensional discrete gauge theory with a non-abelian finite simple gauge group, the fusion of any two lines carrying non-trivial magnetic flux (as in Figure \ref{dyonfusion}) cannot have a unique fusion outcome.

In other words, Theorem 1 asserts we cannot have
\begin{equation}\label{fusionL}
\mathcal{L}_{([g],\pi^{\omega}_g)}\times \mathcal{L}_{([h],\pi_h^{\omega})}=\mathcal{L}_{([k],\pi_{k}^{\omega})}~,\ \ \ g,h\ne1~,
\end{equation}
where, generically, all lines (denoted by $\mathcal{L}$) are non-abelian dyons \footnote{Of course, we may allow for pure fluxes to appear in \eqref{fusionL}.}. We will think of this theorem as a first cousin of the AH conjecture.

\begin{figure}[h!]

\tikzset{every picture/.style={line width=0.75pt}}

\begin{tikzpicture}[x=0.75pt,y=0.75pt,yscale=-0.8,xscale=1]

\draw [color={rgb, 255:red, 0; green, 0; blue, 0 }  ,draw opacity=1 ]   (302.88,78.94) -- (331,121.28) ;
\draw [color={rgb, 255:red, 0; green, 0; blue, 0 }  ,draw opacity=1 ]   (359.68,80.44) -- (331,121.28) ;
\draw [color={rgb, 255:red, 0; green, 0; blue, 0 }  ,draw opacity=1 ]   (331,121.28) -- (330.57,168.94) ;

\draw (272,52.4) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{L}_{([ g] ,\pi ^{\omega }_{g})}$};
\draw (344,53.4) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{L}_{([ h] ,\pi ^{\omega }_{h})}$};
\draw (301,173.4) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{L}_{([k] ,\pi ^{\omega }_{k})}$};

\end{tikzpicture}
\caption{Fusion of dyons}
\label{dyonfusion}
\end{figure}

So far, we have avoided discussing the fusion of Wilson lines. However, in light of \eqref{fusionL}, it is interesting to ask if we can fuse non-abelian Wilson lines ${\mathcal W}_{\pi}$ and ${\mathcal W}_{\pi'}$ (as in Figure \ref{Wilsonfusion}) to obtain a unique outcome
\begin{equation}\label{abcW}
{\mathcal W}_{\pi}\times{\mathcal W}_{\pi'}={\mathcal W}_{\pi''}~.
\end{equation}

\begin{figure}[h!]

\tikzset{every picture/.style={line width=0.75pt}}

\begin{tikzpicture}[x=0.75pt,y=0.75pt,yscale=-0.8,xscale=1]

\draw [color={rgb, 255:red, 45; green, 166; blue, 175 }  ,draw opacity=1 ]   (308.76,81.69) -- (336.89,124.03) ;
\draw [color={rgb, 255:red, 45; green, 166; blue, 175 }  ,draw opacity=1 ]   (365.57,83.19) -- (336.89,124.03) ;
\draw [color={rgb, 255:red, 45; green, 166; blue, 175 }  ,draw opacity=1 ]   (336.89,124.03) -- (336.46,171.69) ;

\draw (295.94,58.85) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{W}_{\pi }$};
\draw (322.94,175.85) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{W}_{\pi''}$};
\draw (354.94,59.85) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{W}_{\pi ^{'}}$};

\end{tikzpicture}

\caption{Fusion of Wilson lines}
\label{Wilsonfusion}
\end{figure}

As we will briefly explain in section \ref{fusionM}, \eqref{abcW} is equivalent to demanding that, at the level of group theory
\begin{equation}\label{charF}
\chi_{\pi}\cdot\chi_{\pi'}=\chi_{\pi''}~,
\end{equation}
where $\chi_{\pi},\chi_{\pi'}$, and $\chi_{\pi''}$ are, respectively, the characters of irreducible linear representations, $\pi, \pi'$, and $\pi''$, of $G$ with dimension greater than 1. Although it might seem strange that \eqref{charF} is possible (especially if one thinks of taking products of irreducible representations in $SU(N)$), it turns out that products of irreducible representations of finite simple groups can be irreducible \cite{zisser1993irreducible}.

The corresponding (twisted or untwisted) discrete gauge theory then has a product of Wilson lines as in \eqref{abcW}. One simple example of this phenomenon in theories with a non-abelian simple gauge group involves the fusion of a Wilson line carrying charge in the $8$-dimensional representation of $A_9$ with a Wilson line carrying charge in either of the $21$-dimensional representations. Intriguingly, the discrete gauge theories based on finite simple groups are prime \cite{naidu2009fusion}, so they do not consist of separate TQFTs with trivial mutual braiding. Therefore, \eqref{abcW} corresponds to some other structural properties of the $A_9$ discrete gauge theory. We will discuss these properties more generally in an upcoming work \cite{toAppear}.

Therefore, we learn that a version of the AH conjecture for characters alone cannot hold. However, our physical discussion above suggests studying one more type of fusion with a unique outcome (Figure \ref{Wilsonfluxfusion})
\begin{equation}\label{abcWd}
{\mathcal W}_{\pi}\times\mathcal{L}_{([g],\pi_g^{\omega})}=\mathcal{L}_{([h],\pi^{\omega}_h)}~, \ \ \ g\ne1~,
\end{equation}
where ${\mathcal W}_{\pi}$ is a non-abelian Wilson line, and the remaining anyons are non-abelian dyons. As a slightly simpler fusion, we may study the following fusion with a unique outcome
\begin{equation}\label{abcWmu}
{\mathcal W}_{\pi}\times\mu_{[g]}=\mathcal{L}_{([h],\pi_h^{\omega})}~, \ \ \ g\ne1~,
\end{equation}
where we have replaced the dyon on the left-hand side of \eqref{abcWd} with a non-abelian flux line. Here we have implicitly assumed that the flux line also exists in the theory (depending on the twist, this assumption may or may not hold).

\begin{figure}[h!]

\tikzset{every picture/.style={line width=0.75pt}}

\begin{tikzpicture}[x=0.75pt,y=0.75pt,yscale=-0.8,xscale=1]

\draw [color={rgb, 255:red, 45; green, 166; blue, 175 }  ,draw opacity=1 ]   (298.7,90.44) -- (326.82,132.78) ;
\draw [color={rgb, 255:red, 208; green, 2; blue, 27 }  ,draw opacity=1 ]   (355.5,91.94) -- (326.82,132.78) ;
\draw [color={rgb, 255:red, 0; green, 0; blue, 0 }  ,draw opacity=1 ]   (326.82,132.78) -- (326.39,180.44) ;

\draw (283.94,68.85) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{W}_{\pi }$};
\draw (345,70.4) node [anchor=north west][inner sep=0.75pt]    {$\mu _{[ g]}$};
\draw (300,183.4) node [anchor=north west][inner sep=0.75pt]    {$\mathcal{L}_{([ h] ,\pi _{h}^{\omega})}$};

\end{tikzpicture}

\caption{Fusion of a Wilson line with a magnetic flux line}
\label{Wilsonfluxfusion}
\end{figure}

This observation brings us to our second cousin of the AH conjecture:

\noindent
{\bf Theorem 2:} In any (twisted or untwisted) discrete gauge theory based on a non-abelian finite simple group, $G$, fusion of the types in \eqref{abcWd} and \eqref{abcWmu} is forbidden.

\noindent
{\bf Intuition:} One heuristic intuition behind this theorem is the following. As a consequence of theorem 1, theorem 2 implies that in discrete gauge theories based on non-abelian simple groups, the only allowed fusions with unique outcomes involving non-abelian anyons are those in \eqref{abcW}. Wilson lines have trivial braiding amongst themselves \footnote{Physically, this last statement is clear from the fact that Wilson lines do not carry magnetic flux. In the language of category theory, this statement follows from the well-known fact that Wilson lines form a symmetric fusion subcategory. In fact, it is a Lagrangian subcategory isomorphic to the category of finite dimensional representations of $G$ over $\mathbb{C}$, ${\rm Rep}(G)$.}. Therefore, even though the fusion in \eqref{abcW} does not arise from a factorization of the TQFT into separate theories with trivial mutual braiding, the Wilson lines themselves have trivial mutual braiding.

Just as theorem 1 follows from the AH conjecture, so too theorem 2 follows from a more basic theorem on finite simple groups which we refer to as the third cousin of the AH conjecture:

\noindent
{\bf Theorem 3:} Consider any non-abelian finite simple group, $G$, any irreducible linear representation, $\pi$, of $G$ having dimension greater than one, and the centralizer, $N_g$, of any $g\ne1$. The restricted representation, $\pi|_{N_g}$, is reducible.

\noindent
We refer to theorems 1-3 as \lq\lq cousins" of the AH conjecture since they are all related by TQFT.

Note that the above discussion is not relevant for abelian simple groups (the type 1 finite simple groups in the above-described classification) since these groups do not have conjugacy classes of length larger than one or representations of dimension larger than one. In other words, their fusion rules are those of a discrete finite group. As a result, we focus on non-abelian finite simple groups.

\noindent
{\bf Duality:} It is also interesting to understand how our above picture is compatible with a type of electric/magnetic duality that often features in discrete gauge theories. For example, the $S_3$ discrete gauge theory has a duality that exchanges the Wilson line charged under the 2-dimensional representation with the line having flux in the 3-cycle conjugacy class \cite{beigi2011quantum,nikshych2014categorical}. More general examples have been discussed in \cite{nikshych2014categorical,naidu2007categorical,Hu:2020ghf}. Clearly, theorems 1, 2, and 3 can only be compatible with such dualities if the Wilson lines participating in \eqref{abcW} are not exchanged with lines carying non-abelian flux. In fact, no such dualities exist in theories based on non-abelian finite simple gauge groups ({\bf Proof:} apply theorem 5.8 of \cite{naidu2007categorical} noting that non-abelian simple groups have no non-trivial abelian normal subgroups). This fact is a non-trivial check of the above picture and is a check of the AH conjecture (this latter claim holds since, if theorem 1 were not true, then the AH conjecture would be false) \footnote{This discussion does not preclude a duality between a twisted discrete gauge theory with a finite simple gauge group and some other type of TQFT (although presumably the dualities should live inside the class of theories considered in \cite{Freed:2020qfy}).}.

The structure of this paper is as follows. In the next section, we give a physicist's introduction to the machinery behind the fusion rules that are relevant for proving the above theorems and deriving theorem 1 from the AH conjecture. In section \ref{fusionM}, we give a brief overview for the more mathematically inclined reader. We derive theorem 1 from the AH conjecture in section \ref{conjecture} before proceeding to proofs of theorems 1, 2, and 3 in section \ref{proofs}. Theorem 3 turns out to be equivalent to theorem 2 and has its own, purely group theoretical proof, that we describe in section \ref{proofs}. The appendices include various additional checks we performed on theorem 3 and therefore 2.

\section{Physics: fusion in discrete gauge theories}\label{fusionP}
One modern perspective on how to go from a group, $G$, to a $2+1$-dimensional discrete gauge theory is to start from a $G$-symmetry-protected topological phase ($G$-SPT) and gauge $G$ \cite{Barkeshli:2014cna}. At the same time, it may be useful to keep in mind that many of the results we will need in this section predate this perspective and follow from the classic work \cite{Dijkgraaf:1989pz}.

The starting point is a set of surface defects in one-to-one correspondence with the elements $g\in G$. For simplicity, we label these defects by group elements as well. Fusion of these defects obeys the usual group multiplication law of $G$, so $g\times h=gh$. One may also consider deforming the associativity of defect fusion via a 3-cocycle
\begin{equation}\label{3cocycle}
\omega(g,h,k)\in H^3(G,U(1))~.
\end{equation}
The $H^3(G,U(1))$ cohomology group then labels the distinct $2+1$-dimensional $G$-SPTs.

Gauging $G$ corresponds to constructing conjugacy classes, $[g_i]$, for a set of representative $g_i\in G$ and pairing this data with an irreducible representation, $\pi^{\omega}_{g_i}$, of the centralizer of each $g_i$, $N_{g_i}$. These are, respectively, the magnetic and electric charges of the discrete gauge theory. The 3-cocycle in \eqref{3cocycle} is the Dijkgraaf-Witten twist (when $\omega=0$ in cohomology we have an untwisted gauge theory).

In this way, lines bounding the $G$-SPT surface operators are liberated and become anyons in the---depending on $\omega$---twisted or untwisted $G$ discrete gauge theory. These latter objects are given by the pair $([g], \pi^{\omega}_{g})$, where the square brackets around $g$ are there to emphasize that we are dealing with a conjugacy class (for any representative in $[g]$, the corresponding centralizers are isomorphic).

The question of whether the electric charge, $\pi^{\omega}_g$, is projective is determined by the reduction of $\omega$ to $N_g$
\begin{equation}
\label{2cocyclefrom3}
\eta_g(h,k):=\frac{\omega(g,h,k)\omega(h,k,g)}{\omega(h,g,k)}\in H^2(N_g,U(1))~,
\end{equation}
where $h,k\in N_g$. Indeed, this is the phase that appears in
\begin{equation}\label{etarep}
\pi^{\omega}_g(h)\pi^{\omega}_g(k)=\eta_g(h,k)\pi^{\omega}_g(hk)~.
\end{equation}
If $\eta_g$ is trivial in $H^2(N_g,U(1))$ the representation is linear \footnote{If $\eta_g(h,k)$ is a non-trivial 2-coboundary, the projective representations obtained will be in one-to-one correspondence with linear representations. These projective factors can be removed using a symmetry gauge transformation as detailed in \cite{Barkeshli:2014cna}. On the other hand, if $\eta_g(h,k)$ is non-trivial in cohomology, then the corresponding representations of $N_g$ must be higher dimensional.  To see this statement, suppose this were not the case. Then, solving \eqref{etarep} for $\eta_{g}(h,k)$ implies that $\eta_g$ is a 2-coboundary.}. For example, the group $PSL(2,4)$ has $\mathbb{Z}_3$ as the centralizer of elements in its length twenty conjugacy class. Since $H^2(\mathbb{Z}_3,U(1))=\mathbb{Z}_1$, the corresponding $\eta_g$ is trivial regardless of the choice of $\omega\in H^3(PSL(2,4),U(1))\simeq\mathbb{Z}_6\times\mathbb{Z}_{10}$. More generally, if $\omega$ is cohomologically non-trivial, then $\pi^{\omega}_g$ is typically projective.

In light of the discussion in the introduction, the most important thing for us to understand is the fusion of two anyons, $([g],\pi^{\omega}_g)$ and $([h],\pi^{\omega}_h)$. Intuitively, we have to fuse both the conjugacy classes as well as the representations that the anyons depend on. This involves identifying the conjugacy classes of the elements obtained by multiplying the elements in $[g]$ and $[h]$. Also, we have to consistently decompose the product $\pi^{\omega}_g \otimes \pi^{\omega}_h$ into irreducible representations of centralizers of $G$. The precise way to carry out these steps is given by \cite{Barkeshli:2014cna,Dijkgraaf:1989pz}
\begin{eqnarray}\label{genfusion}
N_{([g],\pi^{\omega}_g),([h], \pi^{\omega}_h)}^{([k],\pi^{\omega}_k)}&=&\sum_{(t,s)\in N_g \backslash G /N_h} m(\pi^{\omega}_k |_{N_{{}^tg} \cap N_{{}^sh} \cap N_k},\cr && {}^t\pi^{\omega}_{g}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes {}^s\pi^{\omega}_{h}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \cr&&\otimes \pi^{\omega}_{({}^tg,{}^sh,k)})~,
\end{eqnarray}
where ${}^t\pi^{\omega}_{g}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes {}^s\pi^{\omega}_{h}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes \pi^{\omega}_{({}^tg,{}^sh,k)}$ and $\pi^{\omega}_k |_{N_{{}^tg} \cap N_{{}^sh} \cap N_k}$ are (in general reducible) representations of $N_{{}^tg} \cap N_{{}^sh} \cap N_k$ (${}^t\pi^{\omega}_g$, ${}^s\pi^{\omega}_h$, and $\pi^{\omega}_k$ are representations of $N_{{}^tg}$, $N_{{}^sh}$, and $N_k$ which are then restricted to the intersection subgroup). Here we define ${}^tg:=t^{-1}gt$. The projectivity of the ${}^t\pi^{\omega}_{g}$, ${}^s\pi^{\omega}_{h}$, and $\pi^{\omega}_k$ representations is determined by the corresponding cohomology as in \eqref{2cocyclefrom3}. The representation $\pi^{\omega}_{({}^tg,{}^sh,k)}$ is one dimensional (it is a representation of the action of symmetries on the one-dimensional $V_{{}^tg{}^sh}^{k}$ fusion space in the $G$-SPT) and satisfies
\begin{eqnarray}
\pi^{\omega}_{({}^tg,{}^sh,k)}(\ell)\pi^{\omega}_{({}^tg,{}^sh,k)}(m)&=&{\eta_{k}(\ell,m)\over\eta_{{}^tg}(\ell,m)\eta_{{}^sh}(\ell,m)}\cr && \hspace{0.15cm} \cdot ~ \pi^{\omega}_{({}^tg,{}^sh,k)}(\ell m)~.
\end{eqnarray}
These projective factors guarantee that the two arguments of the $m(\cdot,\cdot)$ function can be meaningfully compared. Roughly speaking, this $m(\cdot,\cdot)$ function computes the inner products of the representations appearing as arguments (see \cite{Barkeshli:2014cna} for further details). Finally, let us note that the sum in \eqref{genfusion} is over the double coset, $N_g\backslash G/N_h$.

Another closely related quantity of interest is the modular data of a (twisted or untwisted) discrete gauge theory. It is given by \cite{Hu:2012wx}
\begin{eqnarray}
\label{smatrix}
S_{([g],\pi_{g}^{\omega}),([h],\pi_{h}^{\omega})}&=&\frac{1}{|G|}\sum_{\substack{k \in [g],\ \ell \in [h],\\ k \ell=\ell k}} \chi^{k}_{\pi_g^{\omega}}(\ell)^{*} \chi^{\ell}_{\pi_{h}^{\omega}}(k)^{*}~~, \nonumber\\ \theta_{([g], \pi_{g}^{\omega})}&=&\frac{\chi_{\pi^{\omega}_{g}}(g)}{\chi_{\pi^{\omega}_{g}}(e)}~,
\end{eqnarray}
where $\chi^{h}_{\pi_g^{\omega}}(\ell)$ is defined through the relation
\begin{equation}
\label{projchar}
\chi^{xgx^{-1}}_{\pi_g^{\omega}}(xhx^{-1}):= \frac{\eta_g(x^{-1},xhx^{-1})}{\eta_g(h,x^{-1})} \chi_{\pi^{\omega}_g}(h)~.
\end{equation}
 Here, $\theta$ is the topological spin, and $S$ is the modular $S$ matrix. From these definitions, one can check that the quantum dimensions of the anyons are
\begin{equation}
\label{qdim}
d_{([g],\pi^{\omega}_g)}={S_{([g],\pi_g^{\omega})([1],1)}\over S_{([1],1)([1],1)}}=|[g]|\cdot\text{deg }\pi^{\omega}_g~,
\end{equation}
where $|[g]|$ is the size of $[g]$, and $|\pi^{\omega}_g|$ is the dimension of $\pi^{\omega}_g$. Non-abelian anyons have $d_{([g],\pi^{\omega})}>1$ and necessarily satisfy
\begin{equation}
([g],\pi^{\omega}_g)\times([g^{-1}],(\pi^{\omega}_g)^*)=([1],1)+\cdots~,
\end{equation}
where the ellipses necessarily contain additional terms, $1$ is the trivial representation of $G$, and $(([g^{-1}],(\pi^{\omega}_g)^*)$ is the anyon conjugate to $([g],\pi^{\omega}_g)$.

As we will see in more detail when we prove theorems 1 and 2, anyons $([g],\pi_g^{\omega})$ and $([h],\pi_h^{\omega})$ that fuse to give a unique outcome satisfy the following condition with respect to the $S$ matrix
\begin{equation}
\label{Smatrixcond}
|S_{([g]\pi_g^{\omega}),([h],\pi_h^{\omega})}|=\frac{1}{|G|} d_{([g],\pi_g^{\omega})}d_{([h],\pi_h^{\omega})}~.
\end{equation}
Let us explore the consequences of this relation. To that end, using \eqref{qdim}, we have $d_{([g],\pi_g^{\omega})}d_{([h],\pi_h^{\omega})}=|[g]||[h]|\cdot\text{deg }\pi^{\omega}_g \cdot\text{deg }\pi^{\omega}_h$. Substituting in \eqref{Smatrixcond} and using \eqref{smatrix}, we have
\begin{eqnarray}
&\frac{1}{|G|}& |[g]||[h]|\cdot\text{deg }\pi^{\omega}_g \cdot\text{deg }\pi^{\omega}_h \nonumber \\
&=&\bigg |\frac{1}{|G|}\sum_{\substack{k \in [g],\ \ell \in [h],\\ k \ell=\ell k}} \chi^{k}_{\pi_g^{\omega}}(\ell)^{*} \chi^{\ell}_{\pi_{h}^{\omega}}(k)^{*} \bigg | \nonumber\\
&\leq & \frac{1}{|G|}\sum_{\substack{k \in [g],\ \ell \in [h],\\ k \ell=\ell k}} |\chi^{k}_{\pi_g^{\omega}}(\ell)| |\chi^{\ell}_{\pi_{h}^{\omega}}(k)|\nonumber \\
&\leq& \frac{|[g]||[h]|}{|G|}  \cdot\text{deg }\pi^{\omega}_g \cdot\text{deg }\pi^{\omega}_h
\end{eqnarray}
In the last inequality above, we have used \eqref{projchar} as well as the fact that projective characters satisfy $|\chi_{\pi^{\omega}_g}|\leq \text{deg } \pi_{g}^{\omega}$ \footnote{This statement is guaranteed as long as the projection factors defining the representations are roots of unity, which is satisfied in our case. Indeed, the 3-cocycle $\omega \in H^3(G,U(1))$ can be chosen to be valued in roots of unity without loss of generality.}. It is clear that \eqref{Smatrixcond} is satisfied if and only if the conjugacy classes $[g]$ and $[h]$ commute element-wise and the projective characters satisfy
\begin{equation}
\label{charconst}
|\chi_{\pi^{\omega}_g}(l)|= \text{deg } \pi_{g}^{\omega} \text{ and } |\chi_{\pi^{\omega}_h}(k)|= \text{deg } \pi_{h}^{\omega} ~
\end{equation}
$\forall ~ l\in [h], k \in [g]$. This result is a generalization of lemma 3.4 of \cite{naidu2009fusion}. As mentioned above, it will be crucial for our proofs of theorems 1 and 2.

In the language used in this section, we have that non-abelian Wilson lines, flux lines, and dyons correspond to
\begin{eqnarray}\label{dictionary}
{\mathcal W}_{\pi_1}&\leftrightarrow&([1],\pi_1)~, \ \ \ |\pi_1|>1~,\cr
\mu_{[g]}&\leftrightarrow&([g],1_g^{\epsilon})~, \ \ \ |[g]|>1~,\cr
\mathcal{L}_{([h],\pi_h^{\omega})}&\leftrightarrow&([h],\pi_h^{\omega})~,\ \ \ |[h]|\cdot|\pi_h^{\omega}|>1~.
\end{eqnarray}
We have dropped the $\omega$ superscript from $\pi_1$ in order to emphasize that the Wilson lines always transform under linear representations of $G$.
We attach the $\epsilon$ superscript on the trivial representation of the flux line because these objects only exist when the relevant $\eta_g$ in \eqref{2cocyclefrom3} is trivial in cohomology, and hence of the form $\eta_g(h,k)= \frac{\epsilon_g(h)\epsilon_g(k)}{\epsilon_g(h\cdot k)}$. Finally, $1_g^{\epsilon}$ is the irreducible projective representation of $N_g$ whose character is proportional to the trivial representation of $N_g$.

As a final comment, we note that, from the above modular data, it is easy to show that
\begin{equation}
\theta_{{\mathcal W}_{\pi}}=1~, \ \ \ {S_{{\mathcal W}_{\pi}{\mathcal W}_{\pi'}}\over S_{{\mathcal W}_1{\mathcal W}_{\pi'}}}=1~,
\end{equation}
where ${\mathcal W}_1=([1],1)$ is the trivial Wilson line. In other words, as alluded to in the introduction, the Wilson lines are all bosons and have trivial mutual braiding with each other (they have non-trivial braiding with other lines in the theory).

\section{Math: fusion in $\mathcal{Z}({\rm Vec}_{\omega}^G)$}\label{fusionM}
As a guide to the more mathematically inclined reader, we summarize certain aspects of the previous section in a more formal way. The reader who is only interested in the results or a more physical perspective on them is free to skip this section.

The mathematical framework underlying the physical discussion of \cite{Barkeshli:2014cna} used in the last section is the notion of a $G$-crossed braided category \cite{etingof2009fusion}. The gauging procedure corresponds to the mathematical notion of equivariantization.

In fact, to construct the (twisted or untwisted) discrete gauge theory we can use a simpler notion. Our starting point is the category of $G$-graded vector spaces, ${\rm Vec}^G_{\omega}$, with associator given by the $\mathbb{C}^{\times}$-valued 3-cocycle, $\omega$. The theory we obtain upon equivariantization is the modular tensor category (MTC) constructed by the process of taking the Drinfeld center \cite{bakalov2001lectures,ostrik2002module}. In particular, our gauge theory is just
\begin{equation}\nonumber
{\rm Twisted}\ G\ {\rm discrete\ gauge\ theory}\leftrightarrow\mathcal{Z}({\rm Vec}^G_{\omega})~.
\end{equation}
The various operators discussed in \eqref{dictionary} correspond to the simple objects of $\mathcal{Z}({\rm Vec}^G_{\omega})$ with categorical dimension larger than one. The simple objects corresponding to the trivial conjugacy class in $G$ (what we have called Wilson lines) have trivial topological spin, $\theta$, and are closed under fusion. This means they form a symmetric subcategory. In fact, as is well-known, these simple objects form a Lagrangian subcategory isomorphic to ${\rm Rep}(G)$, the category of finite dimensional representations of $G$ over $\mathbb{C}$. In particular, Wilson line fusion rules are those of the corresponding representation semiring. This observation explains the equivalence of \eqref{abcW} and \eqref{charF}.

\section{From fusion to theorem 1 and a relation between theorems 2 and 3}\label{conjecture}
Given the construction in section \ref{fusionP}, we will first explain why the AH conjecture implies that, in (twisted and untwisted) discrete gauge theories based on simple groups, the fusion of any two lines carrying magnetic flux must have more than one fusion outcome (i.e., theorem 1). After explaining this fact, we will explain the relation between theorems 2 and 3.

To understand the connection between (twisted and untwisted) discrete gauge theories and the AH conjecture, recall the fusion formula in \eqref{genfusion}. Since the arguments of the $m(\cdot,\cdot)$ function are representations of $N_{{}^tg} \cap N_{{}^sh} \cap N_k$, we can decompose them in terms of irreducible representations, $\pi^{\omega(i)}$, of this group
\begin{eqnarray}
&&{}^t\pi^{\omega}_{g}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes {}^s\pi^{\omega}_{h}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes \pi^{\omega}_{({}^tg,{}^sh,k)} \cr
&& \hspace{0.08cm} =\sum_i \alpha_i \pi^{\omega(i)} ~,\cr
&&\pi^{\omega}_k |_{N_{{}^tg} \cap N_{{}^sh} \cap N_k}=\sum_i \alpha^{'}_i \pi^{\omega(i)}~,
\end{eqnarray}
for some non-negative integers $\alpha_i, \alpha^{'}_i$. Then the definition of $m(\cdot,\cdot)$ in \cite{Barkeshli:2014cna} implies
\begin{eqnarray}
 &&m(\pi^{\omega}_k |_{N_{{}^tg} \cap N_{{}^sh} \cap N_k}, {}^t\pi^{\omega}_{g}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes {}^s\pi^{\omega}_{h}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \cr && \hspace{0.05cm} \otimes ~ \pi^{\omega}_{({}^tg,{}^sh,k)})= \sum_i \alpha_i \alpha^{'}_i~.
\end{eqnarray}
We know that $\pi^{\omega}_k$ is an irreducible representation of $N_k$. Also, $N_{{}^tg} \cap N_{{}^sh} \cap N_k$ is a subgroup of $N_k$. According to the Frobenius reciprocity theorem for projective representations of finite groups \cite{gregory} \footnote{We use this theorem in the twisted case; in the untwisted case we use the usual theorem for linear representations.}, we know that, given any irreducible representation, $\pi^{\omega(i)}$, of $N_{{}^tg} \cap N_{{}^sh} \cap N_k$, there is always an irreducible representation, $\pi^{\omega}_k$, of $N_k$ such that the decomposition of $\pi^{\omega}_k |_{N_{{}^tg} \cap N_{{}^sh} \cap N_k}$ into irreducible representations of $N_{{}^tg} \cap N_{{}^sh} \cap N_k$ contains $\pi^{\omega(i)}$. This reasoning shows that, given ${}^t\pi^{\omega}_{g}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes {}^s\pi^{\omega}_{h}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes \pi^{\omega}_{({}^tg,{}^sh,k)}$, there is always some irreducible representation, $\pi^{\omega}_k$, such that $ m(\pi^{\omega}_k |_{N_{{}^tg} \cap N_{{}^sh} \cap N_k}, {}^t\pi^{\omega}_{g}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes {}^s\pi^{\omega}_{h}|_{N_{{}^tg} \cap N_{{}^sh} \cap N_k} \otimes \pi^{\omega}_{({}^tg,{}^sh,k)})$ is non-zero. It follows that once we choose some conjugacy class, $[k]$, such that $[k] \in [g] \cdot [h]$, there is always some $\pi^{\omega}_k$ such that $N_{([g],\pi^{\omega}_g)([h],\pi^{\omega}_h)}^{([k],\pi^{\omega}_k)}\neq 0$. Here, $[g] \cdot[h]$ are the conjugacy classes obtained from taking a product of anyons with magnetic charges in $[g]$ and $[h]$.

Hence, in order to have a fusion rule of the type
\begin{equation}
\label{fusgen}
([g], \pi^{\omega}_g) \times ([h], \pi^{\omega}_h)= ([k], \pi^{\omega}_k)~,\ \ \ g,h\ne1~,
\end{equation}
where all magnetic fluxes on the LHS are non-trivial, we need the fusion of the orbits $[g] \cdot [h]$ to contain only a single orbit $[k]$ (note that $|[k]|$ need not be equal to $|[g]||[h]|$ \footnote{In the case of the fusion of pure fluxes, we do require $|[k]|=|[g]||[h]|$.}). Moreover, commutativity of the fusion rules requires $[k]=[h]\cdot [g]$. Hence, the double coset $N_g \backslash G /N_h$ should have only a single element. (Since the double coset is trivial, we will remove the $t,s$ dependence in the expressions below).
We also require that the decomposition of representations $\pi^{\omega}_k |_{N_g \cap N_h \cap N_k}$ and $\pi^{\omega}_{g}|_{N_g \cap N_h \cap N_k} \otimes \pi^{\omega}_{h}|_{N_g \cap N_h \cap N_k} \otimes \pi^{\omega}_{(g,h,k)}$ into irreps of $N_g \cap N_h \cap N_k$ to have only a single irrep (of multiplicity one) in common. That is, if
\begin{eqnarray}
\pi^{\omega}_{g}|_{N_{g} \cap N_{h} \cap N_k} \otimes \pi^{\omega}_{h}|_{N_{g} \cap N_{h} \cap N_k} \otimes \pi^{\omega}_{(g,h,k)}= \sum_i \alpha_i \pi^{\omega(i)} \nonumber\ \ \ \ \\
\pi^{\omega}_k |_{N_{g} \cap N_{h} \cap N_k}= \sum_i \alpha^{'}_i \pi^{\omega(i)}~,\ \
\end{eqnarray}
then there should be only one $i=i_0$ for which $\alpha_{i_0}= \alpha^{'}_{i_0}\ne0$. Furthermore, we require that $\alpha_{i_0}=1$.

So, in order to have a fusion of the type \eqref{fusgen}, we have two constraints:
\begin{enumerate}
\item $[g] \cdot [h]=[k]=[h]\cdot [g]$
\item $\exists! ~ \pi^{\omega}_k$ such that $m(\pi^{\omega}_k |_{N_{g} \cap N_{h} \cap N_k},\pi^{\omega}_{g}|_{N_{g} \cap N_{h} \cap N_k} \otimes \pi^{\omega}_{h}|_{N_{g} \cap N_{h} \cap N_k} \otimes \pi^{\omega}_{(g,h,k)})=1$
\end{enumerate}
The first constraint is on the conjugacy classes involved, and the second one is on the representations. The AH conjecture immediately implies that (1) is impossible for finite simple groups. Therefore, we see that
\begin{equation}\nonumber
{\rm AH\ conjecture}\Rightarrow{\rm no\ fusions\ as\ in\ \eqref{fusgen}\ for\ simple}\ G~.
\end{equation}
In particular, as claimed in the introduction, we see that
\begin{equation}\label{fusionLii}
\mathcal{L}_{([g],\pi^{\omega}_g)}\times \mathcal{L}_{([h],\pi^{\omega}_h)}\ne\mathcal{L}_{([k],\pi^{\omega}_{k})}~,
\end{equation}
where $\mathcal{L}_{([g],\pi^{\omega}_g)}=([g], \pi^{\omega}_g)$,  $\mathcal{L}_{([h],\pi^{\omega}_h)}=([h], \pi^{\omega}_h)$, and $\mathcal{L}_{([k],\pi^{\omega}_{k})}=([k], \pi^{\omega}_{k})$. So, in that language
\begin{equation}\nonumber
{\rm AH\ conjecture}\Rightarrow{\rm Theorem\ 1}~.
\end{equation}
Of course, this does not prove theorem 1 since the AH conjecture has not been proven. However, it is a non-trivial consistency check of the AH conjecture. We will prove theorem 1 in the next section.

Next, let us show how theorem 3 implies theorem 2. To understand this point, let us specialize the general fusion in \eqref{genfusion} to the product of a non-abelian Wilson line, ${\mathcal W}_{\pi_1}=([1],\pi_1)$, with a non-abelian flux line, $\mu_{[h]}=([h],1_h^{\epsilon})$. In order to have such a flux line in our theory we should, as discussed in section \ref{fusionP}, either consider an untwisted discrete gauge theory or else a theory in which $\omega$ is such that $\eta_h\in H^2(N_h,U(1))$ is cohomologically trivial.

To that end, we find
\begin{eqnarray}
N_{([1],\pi_1),([h], 1_h^{\epsilon})}^{([h],\pi^{\omega}_h)}&=&\sum_{(t,s)\in G \backslash G /N_h} m(\pi^{\omega}_h, {}^t\pi_{1}|_{N_h} \otimes {}^s1_{h}^{\epsilon} \cr&& \hspace{0.15cm}\otimes ~ \pi^{\omega}_{(1,h,h)}|_{N_h})~.
\end{eqnarray}
In this case, the double coset $G \backslash G/N_h$ is trivial. Hence, we have
\begin{equation}
N_{([1],\pi_1),([h], 1_h^{\epsilon})}^{([h],\pi^{\omega}_h)}= m(\pi_h, \pi_{1}|_{N_h} \otimes 1_{h}^{\epsilon}\otimes \pi^{\omega}_{(1,h,h)}|_{N_h} )~.
\end{equation}
In fact, the representation $\pi^{\omega}_{(1,h,h)}$ is trivial (this follows from the fixed nature of the $V_{1h}^h$ fusion space in the $G$-SPT \cite{Barkeshli:2014cna}). So the product of representations $\pi_{1}|_{N_h} \otimes 1_{h}^{\epsilon}\otimes \pi^{\omega}_{(1,h,h)}|_{N_h} $ is isomorphic to $\pi_{1} |_{N_h}\otimes 1_{h}^{\epsilon}$. Therefore, the expression above simplifies to
\begin{equation}\label{Nmrel}
N_{([1],\pi_1),([h], 1_h^{\epsilon})}^{([h],\pi^{\omega}_h)}= m(\pi^{\omega}_h, \pi_{1}|_{N_h} \otimes 1_h^{\epsilon})~.
\end{equation}
Note that $\pi_1$ is an irreducible representation of $G$. Its restriction to $N_h$ is in general reducible. So $m(\pi_h, \pi_1|_{N_h}\otimes 1_h^{\epsilon})$ gives the multiplicity of the irreducible representation, $\pi_h$, in the decomposition of the representation, $\pi_1|_{N_h}\otimes 1_h^{\epsilon}$, into irreducible representations of $N_h$. If $\pi_{1}|_{N_h}$ is irreducible, $m(\pi_h, \pi_1|_{N_h}\otimes 1_h^{\epsilon})= \delta_{\pi_h, \pi_{1}|_{N_h}\otimes 1_h^{\epsilon}}$. Hence, we have the following fusion rules
\begin{equation}
([1], \pi_1) \otimes ([h], 1_h)= ([h], \pi_1|_{N_h}\otimes 1_h^{\epsilon})~,
\end{equation}
if and only if $\pi_1|_{N_h}$ is an irreducible representation of $N_h$.

As a result, theorem 3 implies that we have more than one channel in the fusion
\begin{equation}\label{abcWmu2}
{\mathcal W}_{\pi_1}\times\mu_{[h]}=\mathcal{L}_{([h],\pi_h^{\omega})}+\cdots~.
\end{equation}
In fact, we may take the flux, $([h],1_h^{\epsilon})$, and replace it with a dyon, $([h], \pi^{\omega}_h)$. Note that, in some theories, such a dyon may exist while the flux line does not. We then find that the right-hand side of \eqref{Nmrel} becomes $m(\tilde\pi^{\omega}_h,\pi_1|_{N_h}\otimes \pi^{\omega}_h)$. Clearly, if the fusion in \eqref{abcWmu2} requires more terms on the right-hand side, so too will the fusion with the dyon replacing the flux. This is the content of theorem 2.

Similarly, by the logic of this section, if we satisfy theorem 2 for the untwisted discrete $G$ gauge theory, we then have that, for any irreducible linear representation, $\pi_1$, of $G$ having dimension greater than one, $\pi_1|_{N_h}$ is reducible. This is the content of theorem 3. In conclusion, we have
\begin{equation}\nonumber
{\rm Theorem\ 3}\Leftrightarrow{\rm Theorem\ 2}~.
\end{equation}

Let us also note that we have chosen $\pi_1$ to be an irreducible representation of $G$ with dimension $>1$ so that $([1],\pi_1)$ is non-abelian. Hence, for the above fusion rule to be consistent, $\pi_{1}|_{N_h}$ should be an irreducible representation of $N_h$ of the same dimension.

What remains is to prove at least one of theorems 2 or 3. In the next section we give independent proofs of theorems 2 and 3. Proceeding through theorem 2 first gives us a more TQFT-flavored proof. Proceeding through theorem 3 first gives us a more group theory-flavored proof. We then conclude the next section by proving theorem 1 as well.

\section{Proofs of the cousin theorems}\label{proofs}
From the discussion in the previous section, to prove theorems 2 and 3 we need only prove one of them. However, each route has its own merits, so we give independent proofs of each. We follow by proving theorem 1 (which is logically independent of the others).

Let us first prove theorem 2. To that end, suppose we have a fusion of the form given in \eqref{abcWmu}, which we reproduce below for ease of reference
\begin{equation}\label{abcWmu3}
{\mathcal W}_{\pi}\times\mu_{[g]}=\mathcal{L}_{([h],\pi_h^{\omega})}~, \ \ \ g\ne1~.
\end{equation}
In section \ref{conjecture}, we argued that, if such a fusion exists, the electric charge of the dyon on the right-hand side is given by a reduction of an irreducible representation of the gauge group $G$ (i.e., $\pi_h^{\omega}=\pi|_{N_g}\otimes 1_{h}^{\epsilon}$) and $h=g$. Next, we note that the $S$-matrix satisfies \cite{Kitaev:2006lla}
\begin{equation}
S_{{\mathcal W}_{\pi} \bar{\mu}_{[g^{-1}]}}={1\over|G|}{\theta_{\mathcal{L}_{([g],\pi_g^{\omega})}}\over\theta_{{\mathcal W}_{\pi}}\theta_{\mu_{[g]}}}d_{\mathcal{L}_{([g],\pi_g)}}={1\over|G|}{\theta_{\mathcal{L}_{([g],\pi_g^{\omega})}}\over\theta_{{\mathcal W}_{\pi}}\theta_{\mu_{[g]}}}d_{{\mathcal W}_{\pi}}d_{\mu_{[g]}}~,
\end{equation}
where $\bar\mu_{[g^{-1}]}$ is the conjugate of $\mu_{[g]}$. Therefore,
\begin{equation}\label{Smatrel}
|S_{{\mathcal W}_{\pi} \mu_{[g]}}|={1\over|G|}d_{{\mathcal W}_{\pi}}d_{\mu_{[g]}}~.
\end{equation}
Using \eqref{charconst}, we know that \eqref{Smatrel} implies $|\chi_{\pi}(g)|={\rm deg}\ \chi_{\pi}$, where $\chi_{\pi}$ is the character corresponding to the Wilson line's charge, and ${\rm deg}\ \chi_{\pi}=|\pi|>1$ is the dimension of $\pi$.

A standard argument in representation theory then implies that $\pi(g)=c\cdot\mathds{1}_{|\pi|}$, where $\mathds{1}_{|\pi|}$ is the $|\pi|\times|\pi|$ unit matrix, and $c$ is an $n^{\rm th}$ root of unity (the twist of the dyon). Next, choose some $k\in [G,g]:=\left<\ell g\ell^{-1}g^{-1}|\ell\in G\right>$. Clearly,
\begin{eqnarray}
\pi(k)&=&\pi(\ell g\ell^{-1}g^{-1})=\pi(\ell)\cdot c\cdot\mathds{1}_{|\pi|}\cdot\pi(\ell)^{-1}\cdot c^{-1}\cdot\mathds{1}_{|\pi|}\cr&=&\mathds{1}_{|\pi|}~.
\end{eqnarray}
Since $G$ is a simple group, we can choose $k\ne1$. As a result, $\pi$ is an unfaithful representation of $G$. Therefore, the kernel, ${\rm ker}(\pi)$, is a non trivial normal subgroup. Since $G$ is simple, we must have ${\rm ker}(\pi)=G$. But then, $\pi$ cannot be an irreducible representation. Note that we may repeat this proof verbatim by taking $\mathcal{L}_{([g],\pi_g^{\omega})}$ instead of the flux line. Therefore fusion of the form in \eqref{abcWd} is also forbidden. $\square$

By the discussion in section \ref{conjecture}, we have also proved theorem 3. Although this proof is mathematical, it also has a distinctly TQFT-flavor: notice the prominent role of the modular $S$ matrix (and also, to a lesser extent, the twists).

Alternatively, we may also give a direct group theoretical proof of theorem 3 (and therefore of theorem 2 via section \ref{conjecture}) as follows:

Since $G$ is a non-abelian simple group, its irreducible representations of dimension larger than one must be faithful (otherwise their kernels would be non-trivial normal subgroups). Now, consider some faithful non-abelian representation, $\pi$. Furthermore, take some $g\in G$ such that $g\ne1$ and consider the centralizer, $N_g$.

Suppose the restriction $\pi|_{N_g}$ is irreducible. Clearly $g$ is central in $N_g$. As a result, by Schur's lemma
\begin{equation}
\pi|_{N_g}(g)=c\cdot\mathds{1}_{|\pi|}~,
\end{equation}
where $c$ is an $n^{\rm th}$ root of unity. Since this is a restriction of a representation of $G$, we must also have in the parent group that
\begin{equation}
\pi(g)=c\cdot\mathds{1}_{|\pi|}~,
\end{equation}
and so it follows that
\begin{equation}\label{kereq}
\pi(hgh^{-1}g^{-1})=\mathds{1}_{|\pi|}~.
\end{equation}
Since the group is simple, $g\ne1$ cannot be in the (trivial) center of $G$. As a result, there exists $h$ such that $hgh^{-1}g^{-1}\ne1$. The result in \eqref{kereq} contradicts the fact that $\pi$ is faithful. $\square$

Let us now prove theorem 1. We reproduce the forbidden \eqref{fusionL} for ease of reference
\begin{equation}\label{fusionL2}
\mathcal{L}_{([g],\pi^{\omega}_g)}\times \mathcal{L}_{([h],\pi_h^{\omega})}=\mathcal{L}_{([k],\pi_{k}^{\omega})}~,\ \ \ g,h\ne1~,
\end{equation}
where, according to the discussion in the previous section, $[k]=[gh]$. Similarly to the case of theorem 2, we have that
\begin{eqnarray}
S_{\mathcal{L}_{([g],\pi^{\omega}_g)}\mathcal{L}_{([h^{-1}],(\pi_h^{\omega})^*)}}&=&{1\over|G|}{\theta_{\mathcal{L}_{([gh],\pi_{gh}^{\omega})}}\over\theta_{\mathcal{L}_{([g],\pi^{\omega}_g)}}\theta_{\mathcal{L}_{([h],\pi^{\omega}_h)}}}d_{\mathcal{L}_{([gh],\pi_{gh}^{\omega})}}\cr&=&{1\over|G|}{\theta_{\mathcal{L}_{([gh],\pi_{gh}^{\omega})}}\over\theta_{\mathcal{L}_{([g],\pi^{\omega}_g)}}\theta_{\mathcal{L}_{([h],\pi^{\omega}_h)}}}\cdot \cr&&\cdot d_{\mathcal{L}_{([g],\pi_{g}^{\omega})}} d_{\mathcal{L}_{([h],\pi_{h}^{\omega})}}~,
\end{eqnarray}
where $\mathcal{L}_{([h^{-1}],(\pi_h^{\omega})^*)}$ is the conjugate of $\mathcal{L}_{([h],\pi_h^{\omega})}$. Therefore,
\begin{equation}
|S_{\mathcal{L}_{([g],\pi^{\omega}_g)}\mathcal{L}_{([h],\pi_h^{\omega})}}|={1\over|G|} d_{\mathcal{L}_{([g],\pi_{g}^{\omega})}} d_{\mathcal{L}_{([h],\pi_{h}^{\omega})}}~.
\end{equation}
This last result allows us, as in the case of theorem 2, to use \eqref{charconst}. We then conclude that for any $\ell\in [g]$ and  $m\in [h]$, $\ell m=m\ell$ (i.e., that the two conjugacy classes $[h]$ and $[g]$ commute element-by-element).

Now, consider the product of conjugacy classes
\begin{equation}
[g]\cdot [g]=\sum_{[a]}N_{[g][g]}^{[a]}[a]~, \ \ \ N_{[g][g]}^{[a]}\in\mathbb{Z}_{\ge0}~.
\end{equation}
Clearly, we have that all elements on the left hand side commute with all elements of $[h]$. Therefore, the same is true of all elements in the conjugacy classes $[a]$. Now, consider taking pairwise products of all the $[a]$'s with themselves and with $[g]$ and so on. Eventually, we will come to a set of conjugacy classes closed under multiplication. This defines a normal subgroup $K\trianglelefteq G$ in which each element commutes with $[h]$. Since $G$ is simple, we must have that $K=G$. However, this means that $[h]$ commutes with all elements of the group and so we have a non-trivial center. This is a contradiction. $\square$

\section{Conclusion}
We have argued that $2+1$-dimensional discrete gauge theory is useful for putting conjectures and ideas involving finite simple groups into a broader context and unifying various relevant objects. Using this approach, we proved three theorems that TQFT relates to the AH conjecture.

In fact, we may also generalize the discussion in section \ref{conjecture} and show that the AH conjecture implies that, for any twisted or untwisted discrete gauge theory based on a non-abelian finite simple group, fusions of the form
\begin{equation}\label{fusionLfinal}
\mathcal{L}_{([g],\pi^{\omega}_g)}\times \mathcal{L}_{([h],\pi_h^{\omega})}=\sum_{\pi_{gh}^{\omega}} \mathcal{L}_{([gh],\pi_{gh}^{\omega})}~,\ \ \ g,h\ne1~,
\end{equation}
are not allowed.

Finally, we argued that the lack of electric-magnetic dualities involving discrete gauge theories with non-abelian finite simple groups is a consistency check of our picture above and of the AH conjecture itself.

One natural question that remains is to better understand to what extent ideas involving non-abelian anyons can be used to prove the AH conjecture (see \cite{moori2011products,guralnick2019conjugacy,beltran2019new} and references therein for interesting recent progress on the AH conjecture). Since discrete gauge theories feature in various physical systems, perhaps we can hope for a physics proof of this conjecture.

Theorem 3 is an example of the irreducible restriction problem for simple groups, in the special case of restriction of irreducible representations to centralizers. It will be interesting to explore its relationship with the Aschbacher-Scott program \cite{aschscott, kleschevtiep}.

Another interesting question is to understand to what degree fusion rules of the types we have been discussing constrain global properties of more general TQFTs. We will report on progress toward understanding this last question in \cite{toAppear}.

\acknowledgements{We are grateful for illuminating correspondence and discussions with D.~Aasen, M.~Fayers, S.~Ramgoolam, and I.~Runkel. M.~B. is funded by the Royal Society under the grant, \lq\lq New Constraints and Phenomena in Quantum Field Theory." M.~B. and R.~R. are funded by the Royal Society grant, \lq\lq New Aspects of Conformal and Topological Field Theories Across Dimensions." Our work is also partially supported by the STFC under the grant, \lq\lq String Theory, Gauge Theory and Duality."}

\newpage


\section{Introduction}
As a practically effective approach, model predictive control (MPC), or receding horizon control, has attracted notable attention in the field of industrial process control \cite{Karamanakos2014}.  To deal with optimal control problems, MPC can allow for industrial processes uncertainties and constraints much more straightforwardly than other methods \cite{Kothare1996,Mayne2000}. The core of all model-based predictive algorithms is to use ``open-loop optimal control" instead of ``closed-loop optimal control" within a moving horizon \cite{XI2013}. It brings a lot of robustness and reliability to allow the controller have the ability to recognize the control process. But, the acquisition of knowledge of a priori model seriously affects those performances of MPC. The practically measured data comes from the complicated processes and its utilization into MPC will greatly facilitate the design procedures, avoiding the need for initially accurate dynamic models \cite{Qin2012}.

Nowadays,  many research works focus on the data-driven predictive control, that is applicable on-line both to regulation and tracking control problems \cite{Xu2018,Yang2015,Shah2016,Wang2007,Hou2011,Lauri2010,Zhou2015,Aangenent2005,Ge2012,Ge2008,Kadali2003,zhou2019periodic}. More earlier works can refer to \cite{Skelton2018} and \cite{Aangenent2005}. A Markov data-based LQG control algorithm is suggested in \cite{Skelton2018} and data-based optimal control based on the system's Markov parameters is provided in \cite{Aangenent2005}. Both of them utilize the prior measurements to design the predictive control, along with its implementation on-line. Many results (see, \cite{Wang2007,Hou2011,Lauri2010}) still need some knowledge of the system information to fit the structure of the process model. Other works, like \cite{Li2013c}, focus on the polytopic uncertain systems with unmeasurable system states and a convex hull needs to be known. In \cite{Kadali2003,zhou2020synthesis}, data-driven subspace approach is introduced to design the predictive controller and in \cite{Gorges2017,Shah2016}, reinforcement learning approach is used to reduce the model-based dependence on predictive controller design procedures. We remark that all the aforementioned predictive control methods are designed based on the accurate model, as well as adequate uncertainty description of the linear or non-linear plant of the processes.

In this paper, the limitations of them are circumvented. We use the adequately measured data from the complicated industrial processes following the methods presented in \cite{Jiang2012,Jiang2017} and \cite{zhou2020secure,zhou2019secure}. To be precise, the data-driven learning technique of \cite{Jiang2012} will be employed to iteratively approximate the dynamical parameters, without requiring the \emph{a prior} knowledge of the system matrices. Then, the linear plant's version of \cite{Zhou2018} will be applied to predict the future trajectories, by following the continuous-time predictive control approach of \cite{Yang2015} but removing the assumption on partial knowledge of the system dynamics. Under this framework, the data-driven predictive control input can be generated on-line and can be used for the control of time-varying or nonlinear plants, since the algorithm is able to adapt to the actual dynamics by obtaining a linear model of the system at each sample.
The contributions of the paper are three-fold. (i) A data-driven approach is proposed to adaptively approximate the system matrices, without requiring the \emph{a prior} knowledge of the system. (ii) A continuous-time data-driven MPC approach is developed for the continuous-time linear system, using repeatedly the state and input formation on some fixed time intervals.
(iii) By implementing the proposed data-driven predictive control algorithm, both recursive feasibility of the optimization problem and closed-loop stability of the whole system are guaranteed.

The rest of this paper is organized as follows. In Section II, the problems are briefly formulated. In Section III, the data-driven predictive control approach with completely unknown dynamics algorithm is presented, and the feasibility and stability analysis are conducted. In Section IV, we apply the proposed approach to the optimal control problem of two continuous stirred tank reactor (CSTR). Conclusions are given in Section V.

\textbf{Notation}:  Through this note, $\mathbb{R}$ denotes the set of real numbers, $|\cdot |$ represents the Euclidean norm for a vector and the induced norm for a matrix. For a $a \in \mathbb{R}$, $\mathbb{R}_{\ge a}$ denotes the interval $[a , \infty)$ and $\mathbb{R}_{> a}$ denotes the interval $(a , \infty)$. For any vector $x$, $x^{\textrm{T}}$ denotes its transpose and ${\| x \|^2_P}$ represents ${{x^{\textrm{T}}}Px}$ for a real symmetrix and positive defined matrix $P$. For a real symmetric matrix $A \in \mathbb{R}^{n \times n}$, $A \succ 0$ or $A \succeq 0$ means that $A$ is positive definite or semi-positive definite, $\lambda_{M}(A)$ and $\lambda_{m}(A)$ denote the maximum and minimum eigenvalue of $A$, respectively. $\otimes$ indicates the Kronecker product operator, $\text{vec}(\cdot)$ denotes the vectorization operator and $\text{vec}^{-1}(\cdot)$ denotes the converse vectorization operator, i.e., $\text{vec}(A) = [a_1^{\textrm{T}}, \ldots, a_m^{\textrm{T}}]$ and $\text{vec}^{-1} [a_1^{\textrm{T}}, \ldots, a_m^{\textrm{T}}] = A$, where $a_i \in \mathbb{R}^n$ are the columns of a matrix $A \in \mathbb{R}^{n \times m}$. A continuous function $\alpha: \mathbb{R}_{\ge 0}\rightarrow \mathbb{R}_{\ge 0}$ is said to be a $\mathcal{K}$ function if it is strictly increasing, and $\alpha(s)>0$ for $s>0$ with $\alpha(0) = 0$. A continuous function $\alpha(\cdot)$ is said to be a $\mathcal{K}_{\infty}$ function if it is a $\mathcal{K}$ function, and $\alpha(s) \rightarrow \infty$ for $s \rightarrow 0$. For any piecewise continuous function $u(\cdot): \mathbb{R}_{\ge 0} \rightarrow \mathbb{R}^m$, $u^{[\textrm{r}]}$ denotes the $r$th order derivatives of $u(\cdot)$.

\section{Problem formulation and preliminaries}
Consider a continuous-time industrial process control system with the form
\begin{align}
{\dot x}(t) =&  A{x}(t) + B  u (t) , \quad x(t_0) = x_0, \label{eq:1}\\
y (t) = &  C x(t), \label{eq:2}
\end{align}
where $t \ge {t_0}$, ${x} \in \mathbb{R}^{n}$ is the measurable state; $ u\in \mathcal{U} \subset \mathbb{R}^{m}$ is the control input; $y \in \mathbb{R}^q$ is the controlled output fully available for feedback control design. Assume that $\mathcal{U} \subset \mathbb{R}^{m}$ is a nonempty compact convex set and contains the origin as its interior point. $A \in \mathbb{R}^{n \times n}$ and $B \in \mathbb{R}^{n \times m}$ are unknown system matrices with $(A,B)$ controllable, $(A, C)$ observable, satisfying $|A| \le A_M$, $|B| \le B_M$.

Some standard assumptions are made on \eqref{eq:1} and \eqref{eq:2}. Similar assumptions can be found in \cite{Jiang2012,Zhou2018,Yang2015} for solving (cooperative) output tracking problems.

\noindent\textbf{Assumption 1.} \label{as:3}
The input constraint has the box constrain form as $\underline{u}\le u\le\overline{u}$ with elementwise inequality and $\underline{u},\overline{u}$ the respective lower and upper bounds.

\noindent\textbf{Assumption 2.} \label{as:1}
	There exists a constant matrix $K_0$ such that $A-B K_0$ is a Hurwitz matrix with $ - K_0 x(t) \in \mathcal{U}$.

\noindent\textbf{Assumption 3.} \label{as:2}
	The input relative degree (IRD) of system \eqref{eq:1} is defined as $\rho$.

\noindent\textbf{Remark 1.} \label{re:1}
	In Assumption~\ref{as:3}, we refer the set $\mathcal{U}$ as a box constrain, which accurately describes nearly any set of standard mechanical actuators. Assumption~\ref{as:1} is made such that the initially feasibility can be achieved for the system \eqref{eq:1}. Assumption~\ref{as:2} imposes IRD for the system \eqref{eq:1}, which is used for simplify the solving of regulation problem \cite{Yang2015}.

For \eqref{eq:1} and \eqref{eq:2}, the output of the system should track the given reference $y_{\textrm{d}}(t), t\ge 0$. The tracking error can be given as $e(t) = y(t) - y_{\textrm{d}}(t), t\ge 0$. To that end, a sampled-data MPC, which is based on the repeated solution of an open-loop optimal control problem, is provided in this paper. At each time instant $t= t_k$, the state $y(t_k)$ is measured and then, the controller predicts the system behavior in the future over a prediction horizon $T$ by minimizing a certain objective cost function. The procedure is repeated at every sampling time instant $t_k$ for $k = 1, 2, \ldots$

For the system \eqref{eq:1} and \eqref{eq:2}, the cost function $J(t_k): = J(x(t_k), y_{\textrm{d}}(\cdot), \hat u_k(\cdot) )$, at time $t_k$, is defined as
\begin{align}\label{eq:3-1}
J(x(t_k), y_{\textrm{d}}(s), \hat u_k(s) ) =& \int_{t_k }^{ t_k+ T} L(x(t_k), y_{\textrm{d}}(s), \hat u_k(s)) \textrm{d} s \notag \\
& + F(y_{\textrm{d}}(t_k + T), \hat y(t_k + T))
\end{align}
where, $L(\cdot, \cdot, \cdot)$ and $F(\cdot, \cdot)$ denote the stage and terminal cost functions with the from
\begin{align*}
L(x(t_k), y_{\textrm{d}}(t), \hat u_k( t))  =  \|e (t)\|_Q^2  +\| \hat u_k (t) ) \|_R^2, \quad t \ge 0,
\end{align*}
where $Q  = Q^T \succ 0$ and $R = R^T \succeq 0$ are symmetric and sign definite weight matrices.

Then, at time $t_k$, the optimal control signal ${\hat u}_{k}^\star(s)$, $s\in [t_k, t_k + T]$, is obtained by solving the following finite-horizon optimal control problem as
\begin{subequations} \label{eq:4}
	\begin{align}
	&  \min\limits_{\hat u_k(s) \in \mathcal{U} } J(x(t_k), y_{\textrm{d}}(s), \hat u_k(s)  ) \notag \\
	&s.t. \quad {\dot {\hat x}}(s)  = \mathcal{H} \left(\hat x(s), \hat u_k(s) \right) {\Theta}   \label{eq:4A} \\
	& \qquad ~ \hat y (s) = C \hat x(s), {\hat x}(t_k) = x(t_k)   \label{eq:4B}   \\
	& \qquad~ {\hat u}(s) \in \mathcal{U}, \quad s\in [t_k, t_k +T ]   \label{eq:4C}
	\end{align}
\end{subequations}
where $\mathcal{H} (\cdot, \cdot): \mathbb{R}^n \times \mathbb{R}^m \rightarrow \mathbb{R}^{n \times ({n^2 + mn})} $ is defined as
\begin{align}\label{eq:5}
\mathcal{H}(x, u) = & \left[(x \otimes I_n)^{\textrm{T}} ~~  ( u \otimes I_n  )^{\textrm{T}} \right]
\end{align}
and $\Theta$ denotes the vector of the dynamical parameters for (1), defined as
\begin{align}\label{eq:6}
\Theta =\left[ \text{vec}(A)^{\textrm{T}}  ~~ \text{vec}(B)^{\textrm{T}} \right]^{\textrm{T}}   \in \mathbb{R}^{n^2 + nm}
\end{align}
where $\text{vec}(\cdot)$ denotes the vectorization operator, defined as for any $ P= [p_{ij}] \in \mathbb{R}^{n \times l} , l \in \{n, m\}$, $\text{vec} (P) = \hat P$, where
\begin{align}  \label{eq:7}
\hat P = [p_{11}, p_{21}, ..., p_{n1}, p_{12}, ..., p_{n2}, ..., p_{nl}]^{\textrm{T}}
\end{align}

If the controller has the exact knowledge of system matrices $A$ and $B$, then we can directly solve the finite-time optimal control problem (4) by \cite{Frank2012} to get the desired predictive control policy ${\hat u}_k^\star$. But, in this paper, we seek a data-driven approach to remove the model-based dependence of (6) on the predictive controller design. Thus, the primary objective of this paper is to explore an online learning-based approach to find the data-driven predictive control policy, without requiring any knowledge of the system dynamics, just by using  \emph{a prior}  data of inputs and outputs.

\section{Data-driven predictive control with completely unknown dynamics}
In this section, to facilitate the predictive controller design, the state ${\hat x}$ and the system parameters $\Theta$ of the predicted model (4a) are both estimated from the input-output measurements, using a simultaneous state and parameter estimator. Our proposed online-learning strategy does not rely on $\Theta$, i.e., either $A$ or $B$, which is totally data-driven approach.

\subsection{Dynamical Parameters Estimator}

To obtain the dynamical parameters $\Theta$ in (6), i.e., $A$ and $B$, we present our online learning strategy with the previous input-output data. By using more prior data for continuous-time system (1)  than the instants of solving (4), we try to seek a least-squares optimization solution to get the estimations of all the dynamical parameters $\Theta$.
To that end, at each time instant $t = t_k$, with enough prior date collected, we can find a minimal periodic time interval $\delta t$, such that for some $t_j< t_k$ of $j \in \{1, 2, ..., k-1\}$, there exists a constant integer $N_k$ satisfying $N_k \delta t = t_k - t_j$.

To begin with, for any time instant $t \le  t_k$, doing the integrals at the time interval $[t - \delta t, t ]$ along the trajectory of (1) leads to
\begin{align}\label{eq:8}
x(t) - x(t - \delta t)  = & A \int_{t - \delta t}^{t}  x \textrm{d} \tau +  B \int_{t - \delta t}^{t }  \hat u^\star  \textrm{d} \tau  \notag \\
= & \mathcal{H} \left( \int_{t - \delta t}^{t}  x \textrm{d} \tau, \int_{t - \delta t}^{t }  \hat u^\star  \textrm{d} \tau   \right) {\Theta}
\end{align}

By rearranging (8), we have the linear error system in the form of
\begin{align}\label{eq:9}
\mathcal{F}(t)  = \mathcal{G}(t) \Theta, ~~~ \forall t  \in \mathbb{R}_{[t_j, t_k]}
\end{align}
where the matrices $\mathcal{F} (\cdot): \mathbb{R}_{\ge 0} \rightarrow \mathbb{R}^n $ and $\mathcal{G} (\cdot): \mathbb{R}_{\ge 0} \rightarrow \mathbb{R}^{n \times (n^2 + m^2)} $ are defined,
\begin{align}
\mathcal{F}(t) = & \begin{cases}   x(t) - x(t - \delta t), \quad t \in [\delta t, \infty)  \\
0, \qquad \qquad \qquad \quad t<\delta t
\end{cases}   \label{eq:10} \\
\mathcal{G}(t) = & \left[(\mathcal{A}(t)  \otimes I_n)^\textrm{T} ~~ (\mathcal{B}(t)  \otimes I_n)^\textrm{T} \right]    \label{eq:11}
\end{align}
where $\otimes$ denotes the Kronecker product,  and the vectors $ \mathcal{A} (\cdot) : \mathbb{R}_{\ge 0} \rightarrow \mathbb{R}^n$ and $\mathcal{B} (\cdot): \mathbb{R}_{\ge 0} \rightarrow \mathbb{R}^n $ are defined,
\begin{align}
\mathcal{A}(t) = & \begin{cases}   \int_{t - \delta t}^{t } x(\tau ) \textrm{d} \tau, \quad t \in [\delta t, \infty)  \\
0, \qquad \qquad \quad \quad t<\delta t
\end{cases}    \label{eq:12}  \\
\mathcal{B}(t) = & \begin{cases}   \int_{t - \delta t}^{t } \hat u^\star(\tau )  \textrm{d} \tau, \quad  t \in [\delta t, \infty)  \\
0, \qquad \qquad \quad \quad  t<\delta t
\end{cases}    \label{eq:13}
\end{align}

By (12)-(13), we note that for any time instant $t  \le t_k$, by using the previous input-output measurements $x(\tau)$ and $u(\tau)$, $\tau \le t$, $\mathcal{F}(t)$ and $\mathcal{G}(t)$ in (10)-(11) are both available. Then, by taking periodic sampling $\delta t$ for (1) at each time $t = t_k - i \delta t $  with $ i \in \{0, 1, 2, ..., N_k\}$, such that $ t_j \le t \le t_k $, then we can get a set of composite data,
\begin{align}\label{eq:14}
\mathfrak{D}_k =   \bigcup \{\mathcal{F}_{i}, ~  \mathcal{G}_{i} \}_{i =0}^{N_k}
\end{align}
which satisfies
\begin{align}\label{eq:15}
\mathcal{F}_i = \mathcal{G}_{i} \hat \Theta + \epsilon , ~ \forall i \in \{1, 2, ..., N_k\}
\end{align}
where $\mathcal{F}_i = \mathcal{F} (t_k -i\delta t)$, $\mathcal{G}_i = \mathcal{G} (t_k -i\delta t)$, $\hat \Theta  = [ \hat A^T  ~~ \hat B^T ]^T$ is an estimate of $ \Theta$ in (6) and $\epsilon$ is the estimation error due to the data-driven approximation of $ \Theta$ by using $\mathfrak{D}_k$.

To bring $\epsilon$ to its minimum value, we have the following optimization problem (OP):
\begin{align}  \label{eq:16}
&  \min\limits_{\hat \Theta } ~~ \epsilon^\textrm{T} \epsilon  \\
&s.t. \quad (15) ~ \text{and} ~ (9)-(13)  \notag
\end{align}

Then, an adaptive least-square method is used here to give a solution to problem (16). But, before that, the full rank of dataset $\mathfrak{D}_k$ is defined as follows. To guarantee the existence of the solution of (16), we also give the full rank condition for dataset $\mathfrak{D}_k$.

\noindent\textbf{Definition 1. }\label{de:1}
	At each time instant $t_k$ for some integer $N_k$, the data stack $\mathfrak{D}_k$ is said to have full rank, if there exists an integer $N_0 >0$, such that for all $N_k \ge N_0$, we have the matrix $\mathfrak{A}_k$, defined as

	$\mathfrak{A}_k : = \sum_{i = 1}^{N_k} \mathcal{G}_{i} ^\textrm{T} \mathcal{G}_{i}  \in \mathbb{R}^{(n^2 + nm) \times (n^2 + nm) }$, satisfying
	\begin{align}\label{eq:17}
	\text{rank}(\mathfrak{A}_k ) = n^2 + nm
	\end{align}

Then, we have the following lemma  to give a sufficient condition to guarantee the full rank property of dataset $\mathfrak{D}_k$. The proof is similar to \cite{Ioannou1996}, and here it is omitted.

\noindent\textbf{Lemma 1.}\label{le:1}
	A data stack $\mathfrak{D}_k$ has full rank, if there exists a constant $\underline{d} > 0$, such that \begin{align}\label{eq:18}
	0 < \underline{d} < \gamma_{m}(\mathfrak{A}_k)
	\end{align}

Assume that the measurements of the inputs and the outputs are prior collected at an enough large number $N_k \gg n^2 + nm$ points of time $t_k - i \delta t$ with $i \in \{0, 1, ..., N_k\}$, which makes data stack $\mathfrak{D}_k$ has full rank by Lemma~\ref{le:1}. Then, at time $t = t_k$, one can use the data stack $\mathfrak{D}_k$ to evaluate the unknown dynamical parameters $\Theta$ in (6). By solving (16), the following learning-based update law is obtained
\begin{align}\label{eq:19}
\dot{\hat{\Theta}} =  \eta_{\theta} \sum\limits_{i =1}^{N_k} \mathcal{G}_i^\textrm{T} \left(\mathcal{F}_i - \mathcal{G}_i \hat{\Theta}  \right)
\end{align}
where $0 < \eta_{\theta} \in \mathbb{R}$ is a constant learning rate.

For ease of exposition, the original continuous-time system (1) is expressed in the form of,
\begin{align}\label{eq:20}
{\dot x}(t)  = \mathcal{H} \left(x(t), u(t) \right) {\Theta}
\end{align}
By considering the linear error system with parameters update law (19), it follows that
\begin{align}\label{eq:21}
{\dot x}(t)  = & \mathcal{H} \left( x(t) , u(t)  \right)  \hat  {\Theta} + w(x(t), u(t))   \notag \\
: = & \tilde {A} x(t) + \tilde B  u(t)  + w(t)
\end{align}
where
$ \tilde B = \text{vec}^{-1}(\hat B)$,  $ \tilde A = \text{vec}^{-1}(\hat A)$, $ \text{vec}^{-1}(\cdot)$ denotes the converse vectorization operator, that is, for any vector $\hat P \in \mathbb{R}^{n\times l} $ in (7), we have $  \text{vec}^{-1}(\hat P)= \tilde P = [p_{ij}] \in \mathbb{R}^{n\times l} $. And $w(t):  = w(x(t), u(t))$ is the continuous approximation error resulting from parameters uncertainty (19). If $\epsilon = 0$, then $\hat \Theta =\Theta$, it implies $\tilde A = A$ and $\tilde B = B$, thus we have $w(x, u) = 0$.

The following theorem analyses the property of the term $w$ with respect to (21).

\noindent\textbf{Theorem 1.}\label{th:1}
	The approximate error $w(t)$ in (20) is slowly time-varying, bounded and satisfies $\lim_{t \rightarrow \infty}  w(t) =0 $. \\
\emph{Proof. }
	By considering the closed-loop dynamics (20) and (21), we can refer to $w(t)$ as the unknown disturbance caused by the parameters uncertainty of (19). Then, $w(t)$ satisfies
	\begin{align}\label{eq:22}
	w(t) = \mathcal{H} \left(x(t), u(t) \right) \left( \hat{\Theta} -  \Theta  \right)
	\end{align}

	Letting $ \tilde{\Theta}  = \hat{\Theta} -  \Theta $ and bringing (19) to (22) lead to
	\begin{align}\label{eq:23}
	|\dot w(t)| \le  &| \mathcal{H} \left(x(t), u(t) \right) |  |\dot{ \tilde{\Theta}} | \notag \\
	\le & \eta_{\theta}   | \mathcal{H} \left(x(t), u(t) \right) |  | \sum\limits_{i =1}^{N_k} \mathcal{G}_i^\textrm{T}  \mathcal{G}_i |   | \dot{ \tilde{\Theta}} |  \notag \\
	= &  \eta_{\theta}   | \mathcal{H} \left(x(t), u(t) \right) |  | \mathfrak{A}_k|   | \dot{ \tilde{\Theta}} |
	\end{align}

	By Lemma 1, we have $\gamma_{m} (\mathfrak{A}_k) \le  | \mathfrak{A}_k| \le \gamma_{M} (\mathfrak{A}_k)$. Then, by (19), it implies $\dot {\tilde{\Theta} } = - \eta_{\theta} \mathfrak{A}_k \tilde{\Theta}$ and thus $\lim_{t \rightarrow \infty} \tilde{\Theta} (t) =0 $ and $\tilde{\Theta}  \le  e^{ - \eta_{\theta}  \gamma_{m} (\mathfrak{A}_k) }  \Theta(t_k - N_k \delta t)$. So, $ | \dot{ \tilde{\Theta}} | \le \eta_{\theta}  \gamma_{M} (\mathfrak{A}_k) e^{ - \eta_{\theta}  \gamma_{m} (\mathfrak{A}_k) }  \Theta(t_k - N_k \delta t): = C_{\Theta}$. Besides, for a fixed dataset $\mathfrak{D}_k$ in (14), we have $ | \mathcal{H} \left(x(t), u(t) \right)|$ in (23) bounded for some real constant $M \in \mathbb{R}$ at each time $t = t_k -i \delta t$. Thus, we have
	\begin{align}\label{eq:24}
	|\dot w(t)| \le  \eta_{\theta}   M C_{\Theta}  \gamma_{M} (\mathfrak{A}_k)
	\end{align}
	It implies that $w(t)$ is slowly time-varying and $\lim_{t \rightarrow \infty}  w(t) =0 $. By (23),(24), we also have $|w(t)| \le e^{ \eta_{\theta}   M C_{\Theta}  \gamma_{M} (\mathfrak{A}_k) } $, which means that $w(t)$ is bounded.

\noindent\textbf{Remark 2. }
	Note that in \eqref{eq:21}, the conventional receding horizon expression $A x(t) + Bx(t)$ depending on the unknown matrices $A, B$ is replaced by the term $ \mathcal{H} \left(\hat x(s), \hat u_k(s) \right) \tilde{\Theta}$, where $\tilde{\Theta}$ can be obtained by repeatedly learning from the states and inputs measurements. Furthermore, this learned results will not affects the convergence of the system by Lemma~\ref{th:1}. Therefore, \eqref{eq:19}  plays an important role in identifying the system dynamics from the \emph{a prior} data. As a result, the requirement of the system matrices in predicting the behavior of \eqref{as:2} can be replaced by the state and input information measured online.

\subsection{Receding-horizon Optimization}
To facilitate the data-driven predictive controller design for the system (21) with the dynamical parameters estimator (19), the receding-horizon predictive control problem of (4), at time instant $t = t_k$, can be reformulated,
\begin{subequations}  \label{eq:25}
	\begin{align}
	& {\hat u}_k^\star(s) = \arg \mathop {\min }\limits_{{{\hat u}}(t) \in \mathcal{U}} {J}({{x}}({t_k}),y_\textrm{d}(s), {{\hat u_k}}(s)) \notag \\
	&s.t. \quad {\dot {\hat x}}(s) = \tilde A {\hat x} (s) + \tilde B \hat u_k (s),  \label{eq:25A} \\
	& \qquad ~ \hat y (s) = C \hat x(s),  {\hat x}(t_k) = x(t_k),  \label{eq:25C}  \\
	& \qquad~ {\hat u}(s) \in \mathcal{U}, \quad s \in [t_k, t_k+ T].  \label{eq:25E}
	\end{align}
\end{subequations}

To solve the optimization problem (25),  under Assumption~\ref{as:2}, we define the decision variables as ${\bar u}_k (s)= [\hat u_k^\textrm{T}(s), (\hat u_k^{[1]})^\textrm{T}(s), \ldots, (\hat u_k^{[r]})^\textrm{T}(s)]$ for some  control order $r$ larger than $\rho\ge 1$. Note that the first term of ${\bar u}_k (\tau)$ is the to-be-optimized control input $\hat u_k (\tau)$ in (25). More generally, for the control law $\hat u_k (\tau) $ with a large enough control order $r$, we let $\hat u_k^{[l]} (\tau) =0 $ for any integer $l \ge r$.

Then, for the output prediction of optimization problem (25), by following (2), the future output $y(s) = y(t + \tau), t = t_k, k = 1, 2, \ldots$, in the moving horizon $\tau \in [0, T]$ is approximated by Taylor series expansion,
\begin{align}\label{eq:26}
y(t + \tau) = y(t) + \tau y^{[1]}(t) + \cdots + \frac{\tau^{r}}{ r!} y ^{[r] } (t) + O(\tau^{r})
\end{align}
where the $i$-th derivative of the output $y^{[i]}(t)$ with $i \in \{1, 2, ..., \rho, ..., r\}$ is obtained by
\begin{align}
{y^{[i]}} = & C{{\tilde A}^i}x + \sum\limits_{k = 0}^{i - 1} C {{\tilde A}^{i - 1 - k}}{w^{[k]}},i = 1,...,\rho  - 1  \label{eq:27}\\
{y^{[j]}} = & C{{\tilde A}^j}x + \sum\limits_{k = 0}^{j - \rho } C {{\tilde A}^{j - 1 - k}}\tilde B{u^{[k]}}  \notag \\
& + \sum\limits_{k = 0}^{j - 1} C {{\tilde A}^{j - 1 - k}}{w^{[k]}}, j = \rho, ..., r \label{eq:28}
\end{align}

By rewriting the output $y(t + \tau)$ in a compact form, it follows that
\begin{align}\label{eq:29}
y(t + \tau )  = \begin{bmatrix}  T_1 (\tau) & T_2 (\tau) \end{bmatrix}   \begin{bmatrix}
Y_1 \\  Y_2 \end{bmatrix}
\end{align}
where ${{T}}_1 (\tau) = \left[ 1, \tau, \ldots,  \frac{\tau^{\rho -1}}{(\rho - 1)!}  \right] $, $T_2(\tau) = \left[    \frac{\tau^{\rho}}{ \rho !}, \ldots,  \frac{\tau^{ r }}{ r!}  \right]$, ${Y}_1=\left[ y^\textrm{T}, (y ^{[1]})^\textrm{T} , \ldots,  (y ^{[\rho - 1 ]}) ^\textrm{T}  \right]^\textrm{T}$, ${Y}_2=\left[  (y^{[\rho]} ) ^\textrm{T}, (y^ {[\rho + 1 ]})^\textrm{T} , \ldots,  (y ^{[r]})^\textrm{T} \right]^\textrm{T}$, and
\begin{align}
Y_1 = & \mathcal{A}_1 x +  \mathcal{B}_1 \bar w   \label{eq:30} \\
Y_2 =&  \mathcal{A}_2 x +  \mathcal{B}_2 \bar w  +  \mathcal{B}_3 \bar u \label{eq:31}
\end{align}
where ${\bar w} (s)= \left[ w^\textrm{T}(s),  (w^{[1]})^\textrm{T}(s), \ldots,   (w^{[r]})^\textrm{T} (s) \right]^\textrm{T}$ and $\mathcal{A}_1, \mathcal{B}_1, \mathcal{A}_2, \mathcal{B}_2$ and $ \mathcal{B}_3$ are defined,
\begin{align}\label{eq:32}
\mathcal{A}_1 =& \begin{bmatrix}  C \\ C \tilde A   \\ \vdots \\  C \tilde A^{\rho - 1}  \end{bmatrix} ,    \mathcal{B}_1 =\begin{bmatrix}  C &  0 & \cdots & 0  \\ C \tilde A &  C & \cdots & 0  \\ \vdots & \vdots &  \ddots &  \vdots \\  C \tilde A^{ \rho -1 } & C \tilde A^{ \rho -2}  \tilde B & \cdots & C \end{bmatrix}  \notag \\
\mathcal{A}_2 = &  \begin{bmatrix}  C \tilde A^{\rho } \\ C\tilde A^{\rho +1 }    \\ \vdots \\  C \tilde A^{ r }  \end{bmatrix},
\mathcal{B}_2 \!=\! \begin{bmatrix}  C\tilde A^{\rho-1 }  & 0 \cdots & 0 \\ C\tilde A^{\rho }  &  C \tilde A^{\rho-1 }  & \cdots & 0  \\ \vdots & \vdots & \ddots &  \cdots \\  C \tilde A^{r  }  & C \tilde A^{ r -1}  & \cdots & C\tilde  A^{ \rho -1 }  \end{bmatrix}  \notag  \\
\mathcal{B}_3= & \begin{bmatrix}  C \tilde A^{\rho-1 } \tilde B & 0 &  \cdots & 0 \\ C\tilde A^{\rho }  \tilde B&  C \tilde A^{\rho-1 } \tilde B  & \cdots & 0  \\ \vdots & \vdots & \ddots & \cdots \\  C {\tilde A}^r  \tilde B & C {\tilde  A}^{r -1}  \tilde B & \cdots & C{ \tilde A}^{ \rho -1 } \tilde B  \end{bmatrix}
\end{align}

Then, for the reference signal $y_d(t + \tau)$, by Theorem 1, we have $w^{[k]} (\tau)\approx 0$ for $k=1, 2, \ldots, r$. Thus, $y_d(t + \tau)$ satisfies (29)-(31) with $Y_{1,d} = \left[ y_d^\textrm{T}, (y_d ^{[1]})^\textrm{T} , \ldots ,  (y_d ^{[\rho - 1 ]}) ^\textrm{T}  \right]$ and $ Y_{2, d} =\left[  (y_d^{[\rho]} ) ^\textrm{T}, (y_d^ {[\rho + 1 ]})^\textrm{T} , \ldots ,  (y _d^{[r]})^\textrm{T} \right]$.
But, along with the solution of optimization problem (25) without considering the parameters uncertainty from (21), it leads to
\begin{align}
& \hat y(t + \tau )  = \begin{bmatrix}  T_1 (\tau) & T_2 (\tau) \end{bmatrix}   \begin{bmatrix}
Y_1 \\  Y_2    \end{bmatrix}  \label{eq:33} \\
& Y_1 = \mathcal{A}_1 x, \quad Y_2 =  \mathcal{A}_2 x + \mathcal{B}_3 {\bar u}  \label{eq:34}\\
& \hat u(t + \tau) = T (\tau) \bar u  \label{eq:35}
\end{align}
where $T(\tau) = [T_1(\tau)~ T_2(\tau) ]$.

Then, by bring (34)-(35) into the given performance index (3), one can go to,
\begin{align}\label{eq:36}
J(t_k)  = & \int_{0}^{T}   \left[ \tilde Y_1^\textrm{T}, \tilde Y_2^\textrm{T} \right] \begin{bmatrix}
\Xi_1^\textrm{T} (\tau) \\  \Xi_2^\textrm{T}(\tau)  \end{bmatrix}  \left[ \Xi_1 (\tau) ,  \Xi_2 (\tau) \right]  \begin{bmatrix}
\tilde Y_1 \\ \tilde Y_2 \end{bmatrix}    \notag  \\
& + \int_{0}^{T}  \hat { u} ^\textrm{T}  T_3^\textrm{T} (\tau) R  T_3 (\tau)   \bar { u}   \textrm{d}\tau + F(\tilde Y_i (t_k + T))
\end{align}
where $\tilde Y_i =Y_i -Y_{i,d}$ for $i \in \{1, 2\}$, and $\Xi_i (\tau)= \sqrt{Q}  T_i(\tau) $ for $i \in \{1, 2\}$.
By defining ${\mathcal{T}}_{i, j} = \int_{0}^{T} \Xi_i^\textrm{T} \Xi_j \textrm{d}\tau$ with $i, j \in \{1, 2\}$ and ${\mathcal{T}} = \int_{0}^{T}  T^\textrm{T} R T \textrm{d}\tau$, we note that ${\mathcal{T}}_{1, 2} = {\mathcal{T}}_{2, 1}$. Thus, the performance index (37) can also be rewritten as
\begin{align}\label{eq:37}
J(t_k)  = & \tilde Y_1^\textrm{T} {\mathcal{T}}_{1,1} \tilde Y_1 + 2 \tilde{Y}_1^\textrm{T} {\mathcal{T}}_{1,2} \tilde {Y}_2  +  \tilde{Y}_2^\textrm{T} {\mathcal{T}}_{2,2} \tilde{Y}_2 \notag \\
& +   {\bar u} ^\textrm{T}  \mathcal{T}   {\bar u}   +  F(\tilde Y_i (t_k + T))
\end{align}

Due to equation (34),  taking partial derivative of $J(t_k)$ with respect to $ {\bar u}_k$  yields
\begin{align}\label{eq:38}
\frac{{\partial J}}{{\partial \bar u}} =& 2  \left(  (\frac{{\partial \tilde Y_2}}{{\partial \bar u}})^\textrm{T} {\mathcal{T}}_{2, 2}  \tilde Y_2  + {\mathcal{T}}_4 \right)  {\bar u}  + 2 (\frac{{\partial \tilde Y_2}} {{\partial \bar u}}  )^\textrm{T}  {\mathcal{T}}_{1,2}^\textrm{T}  \tilde{Y}_1     \notag  \\
= & 2  \left(   {\mathcal{B}}_3^\textrm{T}  {\mathcal{T}}_{2,2} {\mathcal{B}}_3 +  {\mathcal{T}}_4  \right)  {\bar u}  + 2 {\mathcal{B}}_3^\textrm{T}  {\mathcal{T}}_{1,2}^\textrm{T} \tilde{Y}_1  \notag   \\
& + 2   {\mathcal{B}}_3^\textrm{T}  {\mathcal{T}}_{2,2} \left( \mathcal{A}_2 x - Y_{2,d}\right)
\end{align}

By letting  ${\partial J}/{{\partial {\bar u}_k  }} = 0$, we can get the optimized predictive control law ${\bar u}_k^\star$ as
\begin{align}\label{eq:39}
{\bar u}_k^\star=& -  \left(   {\mathcal{B}}_3^\textrm{T}  {\mathcal{T}}_{2,2} {\mathcal{B}}_3 +  {\mathcal{T}}_4  \right) ^{-1}   {\mathcal{B}}_3^\textrm{T}  \left( {\mathcal{T}}_{2,2} \bar Y_2 - {\mathcal{T}}_{1,2}^\textrm{T} \tilde{Y}_1 \right)
\end{align}
where $\bar Y_2 =  Y_{2,d} -  \mathcal{A}_2 x $. Then, taking the first row of the optimized control law (39), the continuous-time predictive control law, applied to the plant, is given by
\begin{align}\label{eq:40}
\hat u_k^\star (t) = I_u   {\bar u}_k^\star
\end{align}
where $I_u = [1, 0, ..., 0]_{1 \times (r+1)}$.

\noindent\textbf{Remark 3.}
	Note that from (40), the existence of the optimized solution ${\bar u}_k^\star$ depends on the reversibility of matrix $\mathcal{M} = {\mathcal{B}}_3^\textrm{T}  {\mathcal{T}}_{2,2} {\mathcal{B}}_3 +  {\mathcal{T}}_4$, with ${\mathcal{B}}_3$ computed from (32), $\tilde A$ and $\tilde B$  calculated from (19).  Thus, before we implement the receding-horizon optimization, we first check the reversibility of matrix $\mathcal{M}$ by removing the repeated columns of dataset $\mathfrak{D}_{k}$, only left the distinct columns.

\subsection{Handing constraints}

To deal with optimal control problems, MPC can allow for industrial processes uncertainties and constraints much more straightforwardly than other methods {[}2{]}, {[}3{]}. Assumption~\ref{as:3} will allow us to use a specialized active-set algorithm which is more efficient and easier to implement. A box-constraint solver can be immediately generalized to any linear inequality constraints using slack variables. In the following, we formalize two classical ways to enforce the control limits.

\subsubsection{Saturating Functions}

A conventional attempt to enforce box constraints is to clamp the controls in the forward-pass. The element-wise clamping, or projection operator, is denoted by $sat\text{(\ensuremath{\cdot})}$,which is the input saturation function defined as
\begin{align}
sat(u) & =\left[sat(u_{1})\:sat(u_{2})\:\cdots\:sat(u_{m})\right]\label{eq:3}\\
sat(u_{i}(t)) & =\begin{cases}
u_{i}(t) & {\textstyle if}\,u_{i,min}<u_{i}(t)<u_{i,max}\\
u_{i,min} & {\displaystyle if}\,u_{i}(t)\le u_{i,min}\\
u_{i,max} & if\,u_{i}(t)\ge u_{i,max}
\end{cases}\nonumber
\end{align}
with $u=[u_{1}\:u_{2}\cdots u_{m}]$, and $u_{i,min}\le0$ and $u_{i,max}\ge0$
are the boundaries of $i$th control input of system (1).

It implies that it is tempting to simply replace the obtained control in the forward-pass with
\[
\hat{u}=sat(u^{*})
\]
However, the corresponding search direction may not be a descent direction anymore, harming convergence.

\subsubsection{Squashing Functions}

Another way to enforce box constraints is to introduce a sigmoidal squashing function $s(u)$ on the controls
\[
x_{i+1}=f(x_{i},s(u_{i}))
\]
where $s(\cdot)$ is an element-wise sigmoid with the vector limits
\[
\underset{u\to-\infty}{lim}s(u)=\underline{u},\quad\underset{u\to\infty}{lim}s(u)=\overline{u}
\]
For example, $s(u)=\frac{\overline{u}-\underline{u}}{2}tanh(u)+\frac{\overline{u}+\underline{u}}{2}$
is such a function. A cost term should be kept on the original $u$ and not only on the squashed $s(u)$. Otherwise it will reach very
high or low values and get stuck on the plateau. An intuition for the poor practical performance of squashing is given by the nonlinearity of the sigmoid. Since the backward pass uses a locally quadratic approximation of the dynamics, significant higher order terms will always have a detrimental effect on convergence.

\subsection{Our Proposed Algorithm}
Our proposed method is summarized as the following Algorithm $1$.

\smallskip
\begin{algorithm}
	\caption{Data-driven predictive control algorithm}\label{al:1}
	\KwData{current period $t$; a initially stable control $u^\star = -K_0 x$; prediction horizon $T$; terminal cost $F$;}
	\KwResult{Optimal MPC input $u^{\ast}_{t}$}
	$t_k = k \leftarrow 0$\;
	Collect data $\mathfrak{D}_k$ in (14) \;
	\While{}{Collect the data and form $\mathfrak{D}_k$\;
		\If{(18) is satisfied}{Generate the estimator (19) by using $\mathfrak{D}_k$\;
			Implement (25) to obtain the optimized control (40)\;
			Time evolves continuously with $t$\;}
		$k\leftarrow k+1$\;}
\end{algorithm}

\subsection{Performance Analysis}

Before proceeding further, we first introduce the following definition and lemma.

\noindent\textbf{Definition 2.}  For the system (21), given a compact set $\mathbb{E}$, with $\{0\}   \subset \mathbb{E}\subseteq \mathbb{R}^n$ and $\mathbb{E}$ being a robustly positively invariant set, if there exists a positive definite function $V(\cdot): \mathbb{R}^n \to \mathbb{R}_{\ge 0}$, such that,
	\begin{align}
	V( x) &\ge \alpha_1 (| x|), ~~V( x) \le \alpha_2 (| x|)  + c_1   \label{eq:41}\\
	V(\dot { x}) & -V( x)  \le - \alpha_3 (| x|) + \alpha_ 4(|w|) +c_2  \label{eq:42}
	\end{align}
	for all $t \in \mathbb{R}_{\ge 0} $ with $\alpha_1, \alpha_2, \alpha_3$ being $\mathcal{K}_{\infty}$ function, $\alpha_4$ being $\mathcal{K}$ function, and $c_1, c_2 \ge 0$. Then, the function $V(\cdot)$ is a regional input-to-state practical stability (ISpS)-type Lyapunov function in $\mathbb{E}$ for the system.

Based on Definition 2, we can have the following lemma, directly borrowed from \cite{Yang2015}.

\noindent\textbf{Lemma 2.}  Given a robust positively invariant set $\mathbb{E}$ for the system (21), if it admits an ISpS-type Lyapunov function $V(\cdot)$, then the system is regional ISpS in $\mathbb{E}$, and all the signals of the closed-loop system with the control input $u_{\mathbb{E}}$ are bounded, where  $u_{\mathbb{E}}(t)$ denotes the control such that the set $\mathbb{E}$ is an invariant region satisfying the constraints.

Note that by Assumption 1, the developed optimization problem of data-driven predictive control in (25) is initially feasible with $u_0 = -K_0 x$, then the global stability can be proved by using the Lemma 2.

\noindent\textbf{Theorem 2.}\label{th:2}
	Suppose the Assumption 1 and 2 hold for system with a robust positively invariant set $\mathbb{E}$, then the closed-loop control  plant (1) under the continuous-time MPC law (40), is globally asymptotically stable.

\emph{Proof.}
	The proof is composed by two parts, feasibility and convergence.

	\emph{Feasibility:} Consider any time $t_k$ such that the problem of (25) has a solution and the optimal input $\hat u_{k}^\star$ is implemented for time $[t_k, t_{k+1})$.  Assumed that at $t_{k+1}$, $\hat y(t_{k+1}) = y(t_{k+1})$. 	Therefore, the remaining piece of optimal input $\hat u_k^\star (s), \tau \in [t_{k+1}, t_{k} + T]$ satisfies the input constraints. Thus, we construct the control input as,
	\begin{align}\label{eq:43}
	\hat u_{k+1}(\tau) = \begin{cases}  \hat u_{k}^\star(\tau),  ~~ \tau \in [t_{k+1}, t_{k} + T]
	\\u_{\mathbb{E}}(\tau), ~~ \tau \in [t_{k} + T, t_{k+1} + T]
	\end{cases}
	\end{align}
	where $u_{\mathbb{E}}(\tau)$ makes the desired reference reached and the constraints satisfied. Thus, the predictive control problem is feasible at $t_{k+1}$. It implies that feasibility of the problem at $t_k$ implies the recursive feasibility at $t_{k+1}$.

	\emph{Convergence:} Let the optimal cost function at $t_k$ as the value function $V(x(t_k)) = J^\star(x(t_k), y_\textrm{d}(s), \hat u_k^\star(s) )$. If $V(x(t_k))$ is strictly decreasing, the tracking error $e$ will converge to the origin. To this end, we write the value function at $t_k$ as,
	\begin{align}\label{eq:44}
	V(t_k) =& \int_{t_k}^{t_k+ T} \left( \|e(\tau)\|_Q^2 +\| \hat u_{k}^\star(\tau) ) \|_R^2\right) \textrm{d} \tau  \notag \\
	& + F (y_d(t_k + T), \hat y(t_k + T))
	\end{align}

	Then, by applying $\hat u_{k+1}(t)$  in (44) to the system, beginning from $y(t_{k+1})$, one has
	\begin{align}\label{eq:45}
	J(t_{k+1}) =& \int_{t_{k+1}}^{t_{k+1}+ T} \left( \|e(\tau)\|_Q^2 +\| \hat u_{k+1}(\tau) ) \|_R^2\right) \textrm{d} \tau  \notag \\
	& + F (y_\textrm{d}(t_{k+1} + T), \hat y(t_{k+1} + T))
	\end{align}

	By substituting (45) in (46), it follows,
	\begin{align}\label{eq:46}
	J(t_{k+1}) =& V(t_k) - \int_{t_k}^{t_{k+1}} \left( \| e(\tau)\|_Q^2 +\| \hat u_k^\star(\tau) ) \|_R^2\right) \textrm{d} \tau  \notag \\
	& -   F (y_\textrm{d}(t_k + T), \hat y(t_k + T))       \notag \\
	& +   \int_{t_k + T}^{t_{k+1}+ T} \left( \|e (\tau)\|_Q^2 +\| \hat u_{k+1}(\tau) ) \|_R^2 \right) \textrm{d} \tau   \notag \\
	& + F (y_\textrm{d} (t_{k+1} + T), \hat y(t_{k+1} + T))
	\end{align}

	Note that $\mathbb{E}$ is a robust positively invariant set, it implies that for all $x \in \mathbb{E}$, we have
	\begin{align}\label{eq:47}
	(\partial F / \partial x) (Ax + Bu_{\mathbb{E}}) + L (x, y_\textrm{d}, u_{\mathbb{E}})  \le 0
	\end{align}
	where $L (x, y_\textrm{d}, u_{\mathbb{E}}) = \|e \|_Q^2 +\| u_{\mathbb{E}}  \|_R^2 $. By integrating the inequality (48) along the trajectory of $\dot {\hat x} = \tilde A \hat x + \tilde B u_{\mathbb{E}}$ and $\hat y = C\hat x$, we have $\int_{t_k + T}^{t_{k+1}+ T} \left( \|e (\tau)\|_Q^2 +\| u_{\mathbb{E}}(\tau) ) \|_R^2 \right) \textrm{d} \tau \le F (y_\textrm{d}(t_k + T), \hat y(t_k + T)) - F (y_\textrm{d} (t_{k+1} + T), \hat y(t_{k+1} + T)) $. Thus, by (47), one gets,
	\begin{align}\label{eq:48}
	V(t_k) - J(t_{k+1}) \! \le \! - \int_{t_k}^{t_{k+1}} \left( \| e(\tau)\|_Q^2 +\| \hat u_k^\star (\tau) ) \|_R^2\right) \textrm{d} \tau
	\end{align}

	Further, based on $V(t_{k+1}) = J^\star (t_{k+1}) $, we obtain,
	\begin{align}\label{eq:49}
	V(t_k) - V(t_{k+1}) \! \le \! - \int_{t_k}^{t_{k+1}} \left( \| e(\tau)\|_Q^2 +\| \hat u_k^\star(\tau) ) \|_R^2\right) \textrm{d} \tau
	\end{align}
	It implies that $V(t_{k+1})$ is strictly decreasing. Hence, the proof is complete.

\newcounter{MYtempeqncnt}
\begin{figure*}[!t]

	\normalsize

	\setcounter{MYtempeqncnt}{\value{equation}}

	\setcounter{equation}{49}
	\begin{align}
	\label{eqn_dbl_x}
	A =  \begin{bmatrix}
	-17.98 & -295.866 & 0    &   0 &   0 &    0\\
	0.0207  &  0.1889 &   0.0704  &       0      &   0      &   0\\
	0  &  0.3879  &  0.8000       &  0     &    0     &    0\\
	0.0977   & 0     &    0       &      -18.01 & -295.87 &  0\\
	0  &  0.0617   & 0      &        0.0131   & 0.0433  &  0.0589\\
	0  &  0   &      0      &        0    &     0.3787&   -0.622 \end{bmatrix}, \quad
	B   =  \begin{bmatrix}
	17.8996   &   -13.781  \\
	-0.0131   &  0.0101  \\
	0   &    0  \\
	17.8636& 17.8636  \\
	0.0082 &  0.0082 \\
	0 0 \end{bmatrix}

	\end{align}

	\setcounter{equation}{\value{MYtempeqncnt}}

	\hrulefill

	\vspace*{4pt}
\end{figure*}

\addtocounter{equation}{1}

\section{Application to two-CSTR process}\label{sec5}
Consider two continuous stirred tank reactor (CSTR) system with a full description in \cite{Cao2002,Wang2016}. The open-loop model is a six-state continuous model. The system matrices $A$ and $B$ are directly taken from \cite{Cao2002}, described in the form of (1), as you can see in (\ref{eqn_dbl_x}). The system output variables $y_1 = 362.995 x_2$ and  $y_2 = 362.995 x_4$, denoting the two tank outlet temperatures. The control problem is to maintain the two tank temperatures at desired values $y_d(t) =[y_{1d}(t)~y_{2d}(t)]^\textrm{T} $, where $y_{1d}(t) = 10$ when $0 \le t < 5$s and $y_{1d}(t) = 7$ when $ t \ge 5$s, $y_{2d}(t) = 10$ when $0 \le t < 5$s and $y_{2d}(t) = 4$ when $ t \ge 5$s. The constraints is,
\begin{align*}
\mathcal{U} = \{u =[u_1~ u_2]^\textrm{T}: |u_1| \le 80 , ~~|u_2| \le 70 \}
\end{align*}

In order to illustrate the efficiency of the proposed approach, the precise knowledge of $A$ and $B$ is not used in the design of the predictive controllers. Since the physical system is not stable, the initial stabilizing feedback gain is set as $K_0$,
\begin{align*}
K_0 = \begin{bmatrix}
\begin{smallmatrix}
-4.8949 & -3426.8 & -158.1712    &   -0.0320 &   -43.7963 &   -1.4675 \\
0.1  & 0  &   86.2934  &       1.1730      &   2.3886      &   104.8756
\end{smallmatrix}
\end{bmatrix}
\end{align*}

The weighting matrices $Q$ and $R$ are set to be $Q = diag([10~ 100~ 10~ 10~ 100~ 10])$ and $R = diag([1~ 1])$, respectively. In the simulation, the initial values are selected at the origin. The state and input information is collected over each interval of $0.01$s. When time arrives at $t =2$ s, all the inputs and outputs are repeatedly used to approximate the matrices $A$ and $B$ with $\eta_\theta = 0.85$. The predictive control also starts at $t = 2$ s with the prediction horizon $T = 1$ s. Since then, the control input is immediately updated by solving the problem (25), and the convergence of of $A_k: = \tilde{A}$ and $B_k: = \tilde{B}$ to their actual values is attained after 10 iterations. The procedure of solving (25) is repeated over a fixed interval of $0.1$s.
The convergence of $A_k$ and $B_k$ to their actual values is illustrated in Fig.1. The trajectories of the output variables and the flow rates are shown in Fig.2.
It can be seen that the data-driven predictive control algorithm can stabilize the system, without requiring the system matrices.

\begin{figure}
	\centering

	\includegraphics[width=3.6in]{parameters-eps-converted-to.pdf}\\
	\caption{Convergence of $A_k$ and $B_k$ to their actual values during the control process.}\label{fig:3}
\end{figure}
\begin{figure}
	\centering

	\includegraphics[width=3.6in]{performance-eps-converted-to.pdf}\\
	\caption{The trajectories of the output variables and the flow rates.}\label{fig:2}
\end{figure}

\section{Conclusion}\label{sec10}
In this paper, a data-driven predictive control approach for continuous-time linear system with completely unknown dynamics has been provided. This method solves the infinite-horizon optimal control problem, using the system inputs and outputs information collected online, without knowing the system matrices. The methodology developed in this paper may serve as a computational tool to study the finite-horizon adaptive optimal control of uncertain nonlinear systems. Some related work has appeared in \cite{Ge2008}, which was developed using neural networks, and also in our recent work \cite{Zhou2018}, which proposed a framework of distributed MPC to handle the asynchronous communication, using the \emph{a prior} information  associated with the interconnected neighbors to a distributed optimal design.

\bibliographystyle{IEEEtran}