PINN[PINN]Physics-Informed Neural Network
NN[NN]Neural Network
ODE[ODE]Ordinary Differential Equation

The ability to accurately approximate trajectories of dynamical systems enables their analysis, prediction, and control. Neural network (NN)-based approximations have attracted significant interest due to fast evaluation with good accuracy over long integration time steps. In contrast to established numerical approximation schemes such as Runge-Kutta methods, the estimation of the error of the NN-based approximations proves to be difficult. In this work, we propose to use the NN's predictions in a high-order implicit Runge-Kutta (IRK) method. The residuals in the implicit system of equations can be related to the NN's prediction error, hence, we can provide an error estimate at several points along a trajectory. We find that this error estimate highly correlates with the NN's prediction error and that increasing the order of the IRK method improves this estimate. We demonstrate this estimation methodology for Physics-Informed Neural Network (PINNs) on the logistic equation as an illustrative example and then apply it to a four-state electric generator model that is regularly used in power system modelling.

Dynamical system, error estimation, physics-informed neural network, Runge-Kutta method

§ INTRODUCTION

Differential equations provide a powerful modelling tool to describe how systems change. They are ubiquitous in science and engineering and are used to formulate the behaviour of dynamical systems [Strogatz, 2015]. While instantaneous changes are simple to evaluate, there, usually, exists no analytical solution to these differential equations. To predict future states of the system, we revert to numerical integration methods, see for example [Hairer et al., 1993, Hairer and Wanner, 1996, Iserles, 2008] for detailed treatment of the topic. A key question is the estimation of the approximation error of an integration scheme—both for a single time step and multiple time steps. For established schemes such as RK methods, this estimation hinges around the local approximation error $\varepsilon$ which is closely related to the step size $h$. This relation is expressed as the order $p$ of a scheme $\varepsilon = \mathcal{O}\left(h^{p+1}\right)$ and it describes the error characteristic when $h$ is sufficiently small [Hairer et al., 1993]. To obtain absolute error estimates, an additional RK scheme can be constructed to form an embedded RK method [Hairer et al., 1993]. The resulting error estimates are crucial in choosing time step sizes that yield a balance between a sufficiently accurate approximation and the required computational cost.

For dynamical systems found in power systems, the computational burden of solving the associated differential equations is significant—[Konstantelos et al., 2017] provide an indication of the scales. Grid operators, who need to ensure system stability, rely heavily on these numerical integration methods to choose control actions, but the computational burden becomes a limiting factor. Therefore, the choice of the integration method and accurate error estimate are critical; this becomes even more relevant in light of the energy transition with additional dynamic phenomena [Hatziargyriou et al., 2021]. NN-based approaches such as in [Cui et al., 2023, Moya and Lin, 2023, Stiasny et al., 2023] could have upsides in terms of prediction speed, however, the question of estimating and controlling the approximation error poses a major challenge to their adoption. Motivated by this example, this work shall provide a novel error estimation procedure to understand the error characteristics of learned solutions for dynamical systems.

Neural network-based integration: With NN-based integration methods we learn an approximate solution function to a differential equation or a system of differential equations. This learned solution offers the upside of fast and explicit evaluations while providing sufficiently accurate results over large time steps, given appropriate training. These benefits come at the expense of the need to train the model in advance. As suggested in [Lagaris et al., 1998] the training procedure can be set up as a collocation method. The NN constitutes the candidate function and is fit to match the differential equations at collocation points. This methodology has been revived in [Raissi et al., 2018] under the term PINN.
A related approach, that can be seen as a generalisation of PINNs, is termed operator learning [Lu et al., 2021, Kovachki et al., 2023]; we will subsequently use NN or PINN to refer to these kind of approaches. As these methods do not adhere to error characteristics by construction—unlike RK schemes—understanding and estimating the prediction errors is an open research question.

Approaches for analysing the error of PINNs:
The trivial approach consists in comparing the PINN's prediction with a highly accurate integration scheme that can be regarded as ground truth. However, as this assessment might involve significant computational cost, our interest lies on approaches for estimating the error characteristics without knowing the ground truth. For the simplest version of PINN, we obtain a local approximation error $\mathcal{O}(1)$, meaning that the error is independent of the time step size $h$. [Lagaris et al., 1998] proposed to change the NN's architecture to ensure numerical consistency for small $h$; reducing the time step size should reduce the error. However, achieving high orders $p$ would require many additional function evaluations which reduces the attractiveness of PINN from the standpoint of speed. Furthermore, these changes do not necessarily help with controlling the error for large time steps. This is undesirable as PINN become an attractive alternative to RK schemes only for large time steps.

The following approaches focus on providing guaranteed error bounds to certify a PINN. [Hillebrecht and Unger, 2022] formulate error dynamics between the exact solution and the PINN approximation. These dynamics can be bounded using the Lipschitz constant of the dynamical system, but for non-linear systems this might not result in tight error bounds. Several works [Yarotsky, 2017, De Ryck et al., 2021, De Ryck and Mishra, 2022] aim to provide theoretical upper bounds on the approximation errors of NN based on the function class that they describe, however, these results do not provide practical error estimates. [Eiras et al., 2023] verify that the residual of the differential equations based on the learned approximation is bounded and they show empirically that this bound is related to the maximum approximation error.

Proposed approach: In this work, we combine PINN with the collocation method of Gauss-Legendre Runge-Kutta schemes, a special form of implicit RK schemes, to estimate the error, similar to embedded RK schemes. To this end, we evaluate the PINN at time instances (nodes) defined by the IRK scheme. For a high-order IRK scheme, the function values at the nodes coincide with intermediate values along the trajectory of interest. We then evaluate how closely these predictions match the implicit system of equations defined by the IRK scheme. Based on the mismatch, we can estimate the error at the nodes and hence along the trajectory. As a numerical example, we illustrate the methodology on the logistic equation. A power system generator model with four states serves as a second example to showcase the resulting error estimation on a relevant use case.

Paper structure:<Ref> briefly introduces the notation as well as IRK methods and PINNs as numerical integration methods. <Ref> explains the error estimation procedure. <Ref> illustrate the methodology on the example of the logistic equation and the power system generator respectively. <Ref> concludes and discusses next steps.

§ APPROXIMATING THE TRAJECTORY OF A DYNAMICAL SYSTEM

We assume an autonomous dynamical system governed by

\begin{align}
\ddt \xstate(t) = \fupdateof{\xstate(t)}\label{eq:f_update}
\end{align}

with state $\xstate{} \in \Rdim{n}$, time $t \in \Rdim{}$ and the update function $\fupdate: \Rdim{n} \mapsto \Rdim{n}$. Given an initial condition $\xstate_0 = \xstate{}(t_0)$, the trajectory $\xstate(t_0 + h)$ can be obtained by integration of <ref>

\begin{align}
\xstate(t_0 + h) = \xstate_0 + \int_{t_0}^{t_0 + h} \fupdateof{\xstate(t)} dt\label{eq:trajectory_definition}.
\end{align}

We refer to $h = t - t_0$ as the time step size. Given that $\fupdate$ is continuous and Lipschitz continuous, the existence and uniqueness of $\xstate(h)$ is ensured. Usually, no analytical form for <ref> exists, hence, we revert to constructing an approximation $\xstatehat(h)$ of the trajectory. There is a great variety of established methods, see for example [Hairer et al., 1993, Iserles, 2008], in this work we focus on implicit RK methods (<ref>). In <ref>, we present PINN as a NN-based integration method.

§.§ Implicit Runge-Kutta (IRK) methods

As described in [Iserles, 2008], <ref> constitutes a quadrature problem and it can be approximated by a weighted sum of the function values at  nodes

\begin{align}
\xstate(t_0 + h) \approx \xstate{}_0 + h \sum_{i=1}^{\nstages} b_i \, \fupdateof{\xstate(t_0 + c_i h)}.
\end{align}

The location of each node is determined by the parameter $c_i$ which usually ranges in $[0, 1]$ and it is associated with the positive weight $b_i$. As the quadrature requires $\xstate(t_0 + c_i h)$ at the nodes, their values are approximated by $\rkstage{i} \approx \xstate(t_0 + c_i h)$ and are subsequently referred to as the $i$-th stage of the RK scheme. By constructing a polynomial of order  that passes through the node values i, a set of conditions can be derived—we construct a collocation method. The conditions are described by the following system of implicit equations

\begin{align}\label{eq:irk_system}
\rkstage{i} = \xstate_0 + h \sum_{j=1}^{\nstages} a_{ij} \fupdateof{\rkstage{j}}.
\end{align}

By solving for i, we can then calculate

\begin{align}
\xstatehat(t_0 + h) = \xstate{}_0 + h \sum_{i=1}^{\nstages} b_i \,\fupdateof{\rkstage{i}}.
\end{align}

The number of stages  and the choice of $a_{ij}, b_i, c_i$ determine the accuracy of the scheme. It can be shown, that the choice of Gauss-Legendre polynomials for $b_i, c_i$ leads to a scheme of numerical order $2\nstages{}$. This implies that the approximation error scales according to $\mathcal{O}\left(h^{2\nstages{}}\right)$. This relationship allows IRK schemes to achieve highly accurate approximations for large time step sizes $h$ and moderate numbers of stages $s$. The approximation error $\errorcollocation{i} \in \Rdim{n}$

\begin{align}\label{eq:collocation_error}
\errorcollocation{i} := \rkstage{i} - \xstate(t_0 + c_i h)
\end{align}

of the stages that result from solving the collocation scheme scales as $\mathcal{O}\left(h^\nstages\right)$ [Butcher, 2009]. We will use this relation in the error estimation in <ref> as it provides insightful information at intermediate points along the trajectory.

§.§ Neural-Network-based approximation

In the previous section, the candidate solution for the collocation problem is a polynomial. Instead, we now choose a NN to approximate

\begin{align}\label{eq:NN_approximation}
\xstatehat(t_0 + h) = \xstate{}_0 + h \, \NN{}_{\theta}(h).
\end{align}

The modifications of multiplying $\NN{}_{\theta}(h)$ with $h$ and the addition of the initial condition $\xstate{}_0$ can already be found in [Lagaris et al., 1998]. They ensure that the initial condition is satisfied for $h=0$ and the local approximation error will scale as $\mathcal{O}\left(h\right)$ if $\NN{}_{\theta}(h)$ is bounded.

In this work, we use a simple feed-forward neural network $\NN{}_{\theta}: \Rdim{d_0} \rightarrow \Rdim{d_L}$ with $L$ layers of size $d_l$. It describes a concatenation of affine transformations $\mathcal{A}_l : \Rdim{d_{l-1}} \rightarrow \Rdim{d_l}: z \mapsto W_l z + b_l$ with non-linear element-wise activation functions $\sigma: \Rdim{} \rightarrow \Rdim{}$

\begin{align*}
\NN{}_{\theta}(h) =\mathcal{A}_L \circ \sigma \circ \mathcal{A}_{L-1} \circ \dots \circ \sigma \circ \mathcal{A}_{1}(h).
\end{align*}

The weight matrices $W_l$ and bias vectors $b_l$ of the $L$ layers constitute the trainable parameters $\theta$ and for the present case, the time step size forms the input $(d_0 = 1)$ and the $n$ states of the system form the output $(d_L = n)$. The training problem to determine $\theta$ becomes

\begin{align}\label{eq:loss_pinn}
\min_{\theta} \frac{1}{N} \sum_{j=1}^{N} \norm{\ddt \xstatehat\left(t_0 + h^{(j)}\right) - \fupdateof{\xstatehat\left(t_0 + h^{(j)}\right)}}{2}^2
\end{align}

as introduced in [Lagaris et al., 1998] and termed PINN by [Raissi et al., 2018]. This can be seen as a collocation method for finding $\theta$ by evaluating the norm in <ref> at $N$ collocation points at time step sizes $h^{(j)}$. The temporal derivative $\ddt\xstatehat(t_0 + h^{(j)})$ can be computed using automatic differentiation (AD) [Baydin et al., 2018], hence, this training procedure does not require simulated data points. We solve the optimisation problem in <ref> with a gradient descent method, namely the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS) [Liu and Nocedal, 1989].

§ ERROR ESTIMATES OF PINN PREDICTIONS WITH IRK METHODS

The following error estimation approach is closely linked to the works in [Hairer and Wanner, 1996, de Swart and Söderlind, 1997, Butcher, 2009], however, we do not focus on the error estimation at the end of the time step $|\xstatehat(t_0 + h) - \xstate(t_0 + h)|$. Instead, we aim for  error estimates at the  locations of the RK stages $t_0 + c_i h$ which are spread across the interval $[t_0, t_0 + h]$ as $c_i \in [0, 1]$.

To this end, we evaluate the approximation function $\xstatehat$ defined by the PINN at the node locations $t_0+ c_i h$ that correspond to an IRK method of our choice. These predictions act as approximations

\begin{align}\label{eq:stage_approximation}
\rkstagehat{i}=\xstatehat{}\left(t_0+ c_i h\right)
\end{align}

of the IRK stages i and will, generally speaking, lead to a mismatch $\rkstagedefect{i} \in \Rdim{n}$

\begin{align}\label{eq:stage_error}
\rkstagedefect{i} := \rkstagehat{i} - \rkstage{i},
\end{align}

that we will call the stage error. Using <ref> and $\errorcollocation{i} = \mathcal{O}\left(h^{\nstages}\right)$, we can formulate the desired error estimate for the NN-based function approximation

\begin{align}\label{eq:error_approximation_NN}
\xstatehat(t_0 + c_i h) - \xstate(t_0 + c_i h) = \rkstagehat{i} - (\rkstage{i} + \errorcollocation{i}) =  \rkstagedefect{i} + \mathcal{O}\left(h^{\nstages}\right)
\end{align}

The subsequent paragraph describes the estimation of i without explicitly solving for i.

Estimation of i:
Due to the stage error i, the corresponding system of equations in <ref> evaluated with i will also not be satisfied; a residual $\rkresidual{i} \in \Rdim{n}$ for the system of equations for the $i$-th stage remains

\begin{align}
\rkresidual{i} := \rkstagehat{i} - \xstate_0 - h \sum_{j=1}^{\nstages} a_{ij} \fupdateof{\rkstagehat{j}}.\label{eq:irk_nn_residual}
\end{align}

We begin by replacing $\rkstagehat{i} = \rkstage{i} + \rkstagedefect{i}$ and expanding $\fupdateof{\rkstage{i} + \rkstagedefect{i}}$ in a Taylor series expansion

\begin{align}\label{eq:error_expansion}
\begin{split}
\rkresidual{i} &=  \rkstage{i} + \rkstagedefect{i} - \left(\xstate(t_0) + h \sum_{j=1}^{\nstages} a_{ij} \fupdateof{\rkstage{j} + \rkstagedefect{j}}\right)\\
&= \rkstage{i} + \rkstagedefect{i} - \left(\xstate(t_0) + h \sum_{j=1}^{\nstages} a_{ij} \left[\fupdateof{\rkstage{j}} + \jacobianfx{\rkstage{j}} \rkstagedefect{j} + \mathcal{O}\left({\rkstagedefect{j}}^2\right)\right]  \right)\\
&= \rkstagedefect{i} -  h \sum_{j=1}^{\nstages} a_{ij}  \jacobianfx{\rkstage{j}} \rkstagedefect{j} + \mathcal{O}\left({\rkstagedefect{j}}^2\right)
\end{split}
\end{align}

The last simplification follows from the definition of the IRK scheme evaluated on i as it must equal $0$. As <ref> depends on j, we perform a Taylor series expansion of j using $\rkstage{j} = \rkstagehat{j} - \rkstagedefect{j}$. The relation between i and i becomes

\begin{align}
\rkresidual{i} &= \rkstagedefect{i} -  h \sum_{j=1}^{\nstages} a_{ij}  \jacobianfx{\rkstagehat{j}} \rkstagedefect{j} + \mathcal{O}\left({\rkstagedefect{j}}^2\right).
\end{align}

By collecting $\rkstagedefectall = \begin{bmatrix}
\rkstagedefect{1} & \dots & \rkstagedefect{\nstages}
\end{bmatrix} \in \Rdim{n\nstages}$ and $\rkresidualall = \begin{bmatrix}
\rkresidual{1} & \dots & \rkresidual{\nstages}
\end{bmatrix} \in \Rdim{n\nstages}$ we obtain

\begin{align}
\rkresidualall &= \left( I - h J\right) \rkstagedefectall  + \mathcal{O}\left({\rkstagedefectall}^2\right)\end{align}

with the identity matrix $I\in \Rdim{n\nstages \times n\nstages}$ and the Jacobian-like matrix $J\in \Rdim{n\nstages \times n\nstages}$ that is constructed based on the Jacobians $\jacobianfx{\rkstagehat{i}} \in  \Rdim{n\times n}$ and the weights of the IRK scheme $a_{ij}$

\begin{align*}
J = \begin{bmatrix} a_{11}  \jacobianfx{\rkstagehat{1}} & \dots &  a_{1\nstages}  \jacobianfx{\rkstagehat{\nstages}} \\ \vdots & \ddots & \vdots \\ a_{\nstages 1} \jacobianfx{\rkstagehat{1}} & \dots &  a_{\nstages\nstages} \jacobianfx{\rkstagehat{\nstages}} \end{bmatrix}.
\end{align*}

Neglecting the higher order terms of $\rkstagedefectall$ and inverting $(I - hJ)$, we obtain

\begin{align}\label{eq:error_estimate_delta}
\rkstagedefectall{} &= \left( I - h J\right)^{-1} \rkresidualall{}
\end{align}

which can then be used in <ref> to estimate the prediction error of the NN.

§ LOGISTIC EQUATION - ILLUSTRATION WITH AN ANALYTICAL SOLUTION

This section demonstrates the methodology on the logistic equation

\begin{align*}
\ddt x = x ( 1 - x)\label{eq:logistic}
\end{align*}

with the analytical solution for an initial conditions $x(t_0) > 0$

\begin{align*}
x(t_0 + h) = \frac{e^h}{c_0 +  e^h}; \qquad c_0 = \frac{1}{x(t_0)} - 1.
\end{align*}

We use the simple but nonlinear logistic equation due to the ease of visualisation and the existence of an analytical solution. The arguments can be equally made for dynamical systems with multiple states as we demonstrate later on. To build intuition, we show in <ref> the prediction accuracy of different integration schemes and in <ref> we illustrate the interpretation of the IRK stages as approximate intermediate values $\xstate(t_0 + c_i h)$. <Ref> then presents the error estimation.

§.§ Approximating a trajectory of the logistic equation

We apply different integration schemes to the logistic equation, starting with an initial condition $x_0 = 0.01$. <Ref> shows the trace of the predictions when a single time step of size $h$ is applied. We first want to note that—despite the simplicity of the trajectory—an explicit RK scheme such as the well-known fourth-order RK4 scheme fails to predict large time steps. Hence, we revert to IRK schemes for predictions over long time steps. With three stages ($s=3$), we can still observe an error, for eight stages ($s=8$), the predicted values are indistinguishable from the analytical solution. Lastly, we plot the prediction of a trained PINN and it also closely matches the analytical solution over the entire range of the tested time step sizes.

[Trace of predictions]

[Prediction errors]

Predictions from different numerical integration schemes for a time step of size $h$.

Plotting the approximation error $|\hat{x} - x|$ against the time step size in <ref> provides further insights. By showing both axes on a logarithmic scale the order $p$ of the different schemes become visible. For the PINN, such a convergence properties is not observed, in fact, the error is significantly higher for $h =5$ than for $h=10$. This example illustrates, that while PINNs can predict large time steps accurately, reducing the time step size can lead to larger errors. This creates an error characteristic that is more difficult to analyse and control than for a RK scheme.

§.§ Visualising the stages of an IRK scheme

We now perform a similar analysis for the stage value $k^{(i)}$ of the IRK schemes. We plot their values against the time $c_i h$ that corresponds to the $i$-th node. For $s=3$, we can observe that the stages $k^{(2)}$ and $k^{(3)}$ in <ref> follow the analytical solution well at first, but then deviate. For the higher order IRK scheme $(\nstages = 8)$ the stages $k^{(3)}$ and $k^{(5)}$ coincide closely with the analytical solution, that is the true trajectory.

[Trace of stage values $k^{(i)}$]

[Stage approximation errors]

Analysis of the prediction accuracy of IRK stages for a time step of size $h$.

<Ref> underlines this observation as we plot the stage approximation error $\eta^{(i)}$ against the time step size $h$. This error follows the order $\mathcal{O}\left(h^{\nstages}\right)$, which illustrates that interpreting the stage values $k^{(i)}$ as approximations of the trajectory itself at intermediate points becomes very precise for high order schemes, even for large time step sizes $h$. The proposed error estimate exploits this fact.

§.§ Estimating the stage error

With this intuition in mind, we now apply the methodology presented in <ref>. The question we want to ask is how well the estimated error i correlates with the prediction error $\xstatehat(t_0 + c_i h) - \xstate(t_0 + c_i h)$. A strong correlation improves the usefulness of the estimate as it is more reliable. As the analytical solution exists, we compute the two errors for a time step of $h=5$ (black dots) and $h=10$ (red dots) and show the result in <ref>.

Correlation between the estimated stage error $\delta^{(i)}$ and the corresponding prediction error $|\hat{x}(t_0 + c_i h) - x(t_0 + c_i h)|$ for IRK schemes with 3, 8, and 32 stages. The black dots represent a time step size of $h=5$ and the red dots a time step size of $h=10$.

The subplots show the correlation for different number of stages . The dashed grey line indicates identity between the estimated and true predicted error. It is clearly visible that the larger  and the smaller time step size $h$ improve the correlation which can be expected as the effect of the stage error $\eta^{(i)}$ reduces (see <ref>).

To identify the number of stages that is required for a reliable estimate, we plot the estimated error $\delta^{(i)}$ against the associated time value $c_i h$ for different values of $h$, shown in <ref>. We observe that for an IRK scheme with few stages (left plot), nearby values of $c_i h$ have orders of magnitude different error estimates $\delta^{(i)}$. By increasing the number of stages, the error estimates become more consistent and approach the error characteristic based on the analytical solution.

Consistency of stage prediction for different orders of the IRK schemes. The points indicate the error estimates $\delta^{(i)}$ for a range of values of $h$. The exact error $|\xstatehat - \xstate|$ is shown in orange. For low order IRK schemes, the error estimates for a given time $c_i h$ clearly contradict each other. For a high-order IRK scheme all error estimates form a consistent estimate that is in line with the exact error characteristics.

If the analytical solution is unknown, checking the consistency of the error estimate along a trajectory can be used for validating the estimation accuracy.

§ ERROR ESTIMATION FOR A POWER SYSTEM GENERATOR MODEL

As motivated in <ref>, NN-based approximations might be able to accelerate the simulation of dynamical phenomena in power systems. The PINNSim approach [Stiasny et al., 2023] relies on learning the dynamics of single components, such as generators, and then resolving their interactions. Being able to estimate errors of the learned dynamics will be an important step for future development of such approaches. As an example for the dynamics of a component, we use a two-axis generator model which is modelled by a state vector $\xstate \in \Rdim{4}$ governed by

\begin{align*}
\ddt \begin{bmatrix}
E'_q \\ E'_d \\ \delta \\ \Delta \omega
\end{bmatrix} &= \begin{bmatrix}
\frac{1}{T'_{do}}\left(-E'_q - \left(X_d - X'_d\right) I_d + E_{fd}\right)\\
\frac{1}{T'_{qo}}\left(-E'_d + \left(X_q - X'_q\right) I_q\right)\\
\omega_s \Delta \omega\\
\frac{1}{2H}\left(P_m - E'_d I_d - E'_q I_q - \left(X'_q - X'_d\right) I_d I_q - D \Delta \omega\right)
\end{bmatrix}\\
\begin{bmatrix}
\end{bmatrix} &= \begin{bmatrix}
R_s & -X_q\\X_d & R_s
\end{bmatrix}^{-1} \begin{bmatrix}
E'_d - V \sin{(\delta  - \theta)}\\
E'_q - V \cos{(\delta  - \theta)}
\end{bmatrix}.
\end{align*}

We refer to [Sauer and Pai, 1998] for a detailed description of the modelling and the used parameters. The following results use the parameters associated with generator 2 in <cit.>.

In <ref>, we learned a single trajectory. Here, we aim to predict a collection of trajectories starting from initial conditions $\xstate_0 \in \mathcal{X}_0$ within the set $\mathcal{X}_0$, hence, we adjust <ref> to consider $\xstate_0$ as an additional input to the PINN

\begin{align}
\xstatehat{}(t_0 + h, \xstate_0) = \xstate_0 + h \, \NN{}_{\theta}(h, \xstate_0).
\end{align}

As a consequence, the approximation error will depend on $h$ and $\xstate_0$. Therefore, we sample initial conditions from $\mathcal{X}_0$ and apply the proposed error estimation technique for each $\xstate_0$. As a result we obtain $s$ error estimates along the trajectory of every tested $\xstate_0$. To report the error, we compute a norm over the prediction error $\xstatehat{} - \xstate{}$ (evaluated at $c_i h$) and the estimated error $\rkstagedefect{i}$

\begin{align*}
\norm{\xstatehat{} - \xstate{}}{2} &= \sqrt{\sum_j^n \xi_j^2 \left( \hat{x}_j - x_j \right)^2}\\
\norm{\rkstagedefect{i}}{2} &= \sqrt{\sum_j^n \xi_j^2 {\delta_j^{(i)}}^2}.
\end{align*}

We include a scaling factor $\xi_j$ in the norm so that all state variations have similar magnitude, see Stiasny, 2023 for more details. This factor was also employed when evaluating <ref> during training the PINN model (3 layers with 64 neurons each). As for <ref>, we collect the estimated errors at all $s$ stages. <Ref> shows the resulting error estimates $\norm{\rkstagedefect{i}}{2}$ and the prediction errors evaluated with a differential equation solver with tolerance settings of $10^{-12}$, as no analytical solution is available. To represent the dependency of the prediction error on the different initial conditions, we compute quantiles over small intervals of $c_i h$ which corresponds to the location of the stages.

[Estimated error]

[Error evaluated with a solver]

Error distribution of the PINN predictions with respect to the time step size $h$ for a two-axis generator model. We compare the estimated errors i2 to the prediction error evaluated with a highly accurate solver. As the error depends also on the initial condition, the grey shadings correspond to the errors for 100%, 80%, and 50% of these trajectories.

<Ref> uses $s=8$ and shows a great agreement of the estimated error with the true prediction error apart from a few instances where the proposed error estimation scheme over-estimates the maximum error, for example at $c_i h = 0.01$ and $c_i h = 0.05$. By increasing the number of stages, the estimation quality can be further improved.

§ CONCLUSIONS

The estimation of prediction errors poses a major challenge for PINN. This work introduces a simulation-free and practical error estimation approach, which can remove the barrier for incorporating PINN in established numerical solvers.

Our methodology only requires the evaluation of an implicit Runge-Kutta (IRK) scheme (whose coefficients are usually pre-computed) and a matrix inversion, which leads to a simple and easy-to-implement estimation method, see <ref>. We demonstrate that the proposed approach can accurately estimate the prediction errors along a trajectory and is also applicable for dynamical systems with multiple states. We also show that it is straightforward to assess the quality of this estimate by analysing the consistency of the error estimates that stem from different time step sizes $h$.

Future work: From a theoretical perspective, we see the need to further investigate the influencing factors on the estimation accuracy. That is, first, the analysis of the impact of higher-order terms of $\rkstagedefect{i}$ which could become significant if the PINN's prediction is not very accurate. Furthermore, the impact of system characteristics such as stiffness and oscillatory behaviour on the required number of stages $s$ should be investigated.

With respect to practical applications, the proposed estimation methodology can be used in monitoring the training progress, assessing the learned model, or providing estimates during prediction tasks. We plan to explore these use cases to develop heuristics for this method. Furthermore, the methodology could inform an additional term in the loss function similar to [Raissi et al., 2018, Stiasny, 2023]. Lastly, the error estimates could also be employed at the evaluation stage to correct the PINN's prediction.

This work was supported by the ERC Project VeriPhIED, funded by the European Research Council, Grant Agreement No: 949899. The authors would like to thank Samuel Chevalier, Ignasi Ventura Nadal, and Petros Ellinas for their feedback on the manuscript.

[Baydin et al., 2018]
Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.
Automatic differentiation in machine learning: a survey.
Journal of Machine Learning Research, 180 (153):0 1–43, 2018.

[Butcher, 2009]
J. C. Butcher.
Practical Runge–Kutta methods for scientific computation.
The ANZIAM Journal, 500 (3):0 333–342, 2009.

[Cui et al., 2023]
Wenqi Cui, Weiwei Yang, and Baosen Zhang.
A Frequency Domain Approach to Predict Power System Transients.
IEEE Transactions on Power Systems, pages 1–13, 2023.

[De Ryck and Mishra, 2022]
Tim De Ryck and Siddhartha Mishra.
Generic bounds on the approximation error for physics-informed (and) operator learning.
In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in neural information processing systems, volume 35, pages 10945–10958. Curran Associates, Inc., 2022.

[De Ryck et al., 2021]
Tim De Ryck, Samuel Lanthaler, and Siddhartha Mishra.
On the approximation of functions by tanh neural networks.
Neural Networks, 143:0 732–750, 2021.

[de Swart and Söderlind, 1997]
Jacques J.B de Swart and Gustaf Söderlind.
On the construction of error estimators for implicit Runge-Kutta methods.
Journal of Computational and Applied Mathematics, 860 (2):0 347–358, 1997.

[Eiras et al., 2023]
Francisco Eiras, Adel Bibi, Rudy Bunel, Krishnamurthy Dj Dvijotham, Philip Torr, and M. Pawan Kumar.
Provably correct physics-informed neural networks.
arXiv preprint arXiv:2305.10157, 2023.

[Hairer and Wanner, 1996]
Ernst Hairer and Gerhard Wanner.
Solving Ordinary Differential Equations II, volume 14 of Springer Series in Computational Mathematics.
Springer Berlin Heidelberg, Berlin, Heidelberg, 1996.

[Hairer et al., 1993]
Ernst Hairer, Gerhard Wanner, and Syvert P. Nørsett.
Solving Ordinary Differential Equations I, volume 8 of Springer Series in Computational Mathematics.
Springer Berlin Heidelberg, Berlin, Heidelberg, 1993.

[Hatziargyriou et al., 2021]
Nikos Hatziargyriou, Jovica Milanovic, Claudia Rahmann, Venkataramana Ajjarapu, Claudio Canizares, Istvan Erlich, David Hill, Ian Hiskens, Innocent Kamwa, Bikash Pal, Pouyan Pourbeik, Juan Sanchez-Gasca, Aleksandar Stankovic, Thierry Van Cutsem, Vijay Vittal, and Costas Vournas.
Definition and Classification of Power System Stability – Revisited & Extended.
IEEE Transactions on Power Systems, 360 (4):0 3271–3281, 2021.

[Hillebrecht and Unger, 2022]
Birgit Hillebrecht and Benjamin Unger.
Certified machine learning: A posteriori error estimation for physics-informed neural networks.
In International Joint Conference on Neural Networks, IJCNN 2022, Padua, Italy, July 18-23, 2022, pages 1–8. IEEE, 2022.

[Iserles, 2008]
Arieh Iserles.
A First Course in the Numerical Analysis of Differential Equations.
Cambridge University Press, 2 edition, November 2008.

[Konstantelos et al., 2017]
Ioannis Konstantelos, Geoffroy Jamgotchian, Simon H. Tindemans, Philippe Duchesne, Stijn Cole, Christian Merckx, Goran Strbac, and Patrick Panciatici.
Implementation of a massively parallel dynamic security assessment platform for large-scale grids.
IEEE Transactions on Smart Grid, 80 (3):0 1417–1426, 2017.

[Kovachki et al., 2023]
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.
Neural operator: Learning maps between function spaces with applications to pdes.
Journal of Machine Learning Research, 240 (89):0 1–97, 2023.

[Lagaris et al., 1998]
I.E. Lagaris, A. Likas, and D.I. Fotiadis.
Artificial neural networks for solving ordinary and partial differential equations.
IEEE Transactions on Neural Networks, 90 (5):0 987–1000, 1998.

[Liu and Nocedal, 1989]
Dong C. Liu and Jorge Nocedal.
On the limited memory BFGS method for large scale optimization.
Mathematical Programming, 450 (1-3):0 503–528, 1989.

[Lu et al., 2021]
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.
Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators.
Nature Machine Intelligence, 30 (3):0 218–229, 2021.

[Moya and Lin, 2023]
Christian Moya and Guang Lin.
DAE-PINN: a physics-informed neural network model for simulating differential algebraic equations with application to power networks.
Neural Computing and Applications, 350 (5):0 3789–3804, 2023.

[Raissi et al., 2018]
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.
Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.
Journal of Computational Physics, 3780 (C), 2018.

[Sauer and Pai, 1998]
Peter W. Sauer and M. A. Pai.
Power system dynamics and stability.
Prentice Hall, Upper Saddle River, N.J, 1998.

[Stiasny, 2023]
Jochen Stiasny.
Physics-Informed Neural Networks for Power System Dynamics.
PhD thesis, Technical University of Denmark, 2023.

[Stiasny et al., 2023]
Jochen Stiasny, Baosen Zhang, and Spyros Chatzivasileiadis.
PINNSim: A simulator for power system dynamics based on physics-informed neural networks.
arXiv preprint arXiv:2303.10256, 2023.

[Strogatz, 2015]
Steven H. Strogatz.
Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.
Westview Press, a member of the Perseus Books Group, Boulder, CO, second edition edition, 2015.

[Yarotsky, 2017]
Dmitry Yarotsky.
Error bounds for approximations with deep ReLU networks.
Neural Networks, 94:0 103–114, 2017.


# Walsh Meets OAM in Holographic MIMO

Charles Vanwynsberghe†, Jiguang He†, Chongwen Huang‡, and Merouane Debbah†
†Technology Innovation Institute, 9639 Masdar City, Abu Dhabi, United Arab
Emirates
‡College of Information Science and Electronic Engineering, Zhejiang
University, Hangzhou 310027, China

###### Abstract

Holographic multiple-input multiple-output (MIMO) is deemed as a promising
technique beyond massive MIMO, unleashing near-field communications,
localization, and sensing in the next-generation wireless networks. Semi-
continuous surface with densely packed elements brings new opportunities for
increased spatial degrees of freedom (DoFs) and spectrum efficiency (SE) even
in the line-of-sight (LoS) condition. In this paper, we analyze holographic
MIMO performance with disk-shaped large intelligent surfaces (LISs) according
to different precoding designs. Beyond the well-known technique of orbital
angular momentum (OAM) of radio waves, we propose a new design based on polar
Walsh functions. Furthermore, we characterize the performance gap between the
proposed scheme and the optimal case with singular value decomposition (SVD)
alongside perfect channel state information (CSI) as well as other benchmark
schemes in terms of channel capacity. It is verified that the proposed scheme
marginally underperforms the OAM-based approach, while offering potential
perspectives for reducing implementation complexity and expenditure.

###### Index Terms:

Holographic MIMO, LIS, near field, OAM, Walsh function.

##  I Introduction

Holographic multiple-input multiple-output (MIMO) is considered as a cutting-
edge technique beyond massive MIMO [1, 2]. It unlocks the great potential of
radiative near-field broadband communications by exploiting the numerous
spatial degrees of freedom (DoFs) even under pure line-of-sight (LoS)
condition. Thanks to densely packed antennas usually with less than half-
wavelength inter-element spacing, highly directional pencil beams can be
designed for ultra-high data rate transmissions.

Due to the extreme densification of both transmit and receive elements, the
channel state information (CSI) acquisition faces tremendous challenges. The
reason lies in that the channel sparsity in the angular-domain representation
vanishes when the receiver (Rx) moves from the far field to the near filed of
the transmitter (Tx) [3]. In terms of precoding design, adopting the
conventional right singular vectors of the channel (encounters high
computational complexity), requires the availability of CSI information, and
demands hardware flexibility to perform operations with variable amplitudes
and phases. In the literature, the potential of communication modes and
orbital orbital angular momentum (OAM) multiplexing was intensively
investigated [4, 5] with different antenna structures, e.g., uniform circular
arrays (UCAs) and rectangular arrays. The OAM takes the advantage of
electromagnetic (EM) wave property and generates orthogonal waves to multiplex
various data streams at the transmitter to take full advantage of the
available spatial multiplexing gain. Recently, two different types of OAM
designs, i.e., focused and unfocused, were introduced in [6]. Interestingly,
no CSI is required111It should be noted that the Tx-Rx distance is needed to
generate focused OAM modes., and continuous phase-shifters are sufficient to
implement the precoder because the amplitude of OAM modes is constant.

Motivated by [6], we in this paper develop a design of communication modes
based on polar Walsh functions. Such modes possess attractive properties: i)
they are separable in radial and angular coordinates, ii) they are orthogonal,
iii) their amplitude is constant, and their phases only take two values. The
point iii) is interesting to substantially reduce the hardware complexity of
the precoder and combiner[7]. In addition, we compare the channel capacity
with perfect CSI and the two benchmark schemes from [6]. With a marginal
channel capacity degradation with respect to the OAM modes [6], the polar
Walsh functions provide an interesting trade-off between communication
performance and hardware complexity.

Notation: A bold lowercase letter ${\mathbf{a}}$ denotes the column vector, a
bold capital letter ${\mathbf{A}}$ denotes the matrix, $(\cdot)^{\mathsf{H}}$
and $(\cdot)^{*}$ denote the Hermitian transpose and conjugate, respectively,
and $\|\cdot\|_{2}$ is the Euclidean norm.

##  II System Model

We leverage OAM-like multiplexing for data transmissions in the LIS-assisted
communication network, shown in Fig. 1, which is capable of processing the
signals in the EM domain. The Tx and Rx are equipped with semi-continuous
antenna apertures, centered at origin and point $(0,0,d)$ with areas
$|\mathcal{S}_{t}|=\pi r^{2}_{t}$ and $|\mathcal{S}_{r}|=\pi r^{2}_{r}$,
respectively, where
$\mathcal{S}_{t}=\\{(x,y,0)\in\mathbb{R}^{3}:x^{2}+y^{2}\leq r_{t}^{2}\\}$,
$\mathcal{S}_{r}=\\{(x,y,d)\in\mathbb{R}^{3}:x^{2}+y^{2}\leq r_{r}^{2}\\}$,
and $r_{t}$ and $r_{r}$ are the radii of the Tx and Rx apertures. These two
antennas are thus parallel to the $x\text{-}y$ plane, centered onto the $z$
axis and separated by a distance $d$.

###  II-A Continuous Propagation Model

The generated wave $\phi({\mathbf{t}}_{0})$ at a certain Tx point
${\mathbf{t}}_{0}\in\mathcal{S}_{t}\subset\mathbb{R}^{3}$ and the received
wave $\psi({\mathbf{r}})$ at a certain Rx point
${\mathbf{r}}\in\mathcal{S}_{r}\subset\mathbb{R}^{3}$ satisfy [8, 6]

$\psi({\mathbf{r}})=\int_{{\mathbf{t}}_{0}\in\mathcal{S}_{t}}G({\mathbf{r}},{\mathbf{t}}_{0})\phi({\mathbf{t}}_{0})\mathrm{d}{\mathbf{t}}_{0},$
(1)

where $G({\mathbf{r}},{\mathbf{t}}_{0})$ denotes the Green’s function. Since
we focus on near-field propagation regime, $G({\mathbf{r}},{\mathbf{t}}_{0})$
characterizes the wave propagation from the point ${\mathbf{t}}_{0}$ to the
point ${\mathbf{r}}$ in terms of phase change and attenuation. Therefore,
$G({\mathbf{r}},{\mathbf{t}}_{0})$ is written in the form of

$G({\mathbf{r}},{\mathbf{t}}_{0})=\frac{\exp(-\jmath\kappa\|{\mathbf{r}}-{\mathbf{t}}_{0}\|_{2})}{4\pi\|{\mathbf{r}}-{\mathbf{t}}_{0}\|_{2}},$
(2)

where $\kappa=2\pi/\lambda$ with $\lambda$ being the wavelength and
$\jmath=\sqrt{-1}$.

Similar to singular value decomposition (SVD), we formulate a series of
eigenproblems as

$\varepsilon\phi({\mathbf{t}})=\int_{{\mathbf{t}}_{0}\in\mathcal{S}_{t}}K({\mathbf{t}},{\mathbf{t}}_{0})\phi({\mathbf{t}}_{0})\mathrm{d}{\mathbf{t}}_{0},$
(3)

where the kernel $K({\mathbf{t}},{\mathbf{t}}_{0})$ is defined as

$K({\mathbf{t}},{\mathbf{t}}_{0})=\int_{{\mathbf{r}}_{0}\in\mathcal{S}_{r}}G^{*}({\mathbf{r}}_{0},{\mathbf{t}})G({\mathbf{r}}_{0},{\mathbf{t}}_{0})\mathrm{d}{\mathbf{r}}_{0},$
(4)

and $\varepsilon$ is the eigenvalue associated with the eigenfunction
$\phi(\cdot)$. Finding the solutions of this continuous eigenproblem is not
straightforward. It is known from the literature that the circular prolate
spheroidal wave functions (CPSWFs) are analytical solutions to the
eigenproblem in Eq. (3). By resorting to CPSWFs, the DoFs of the channel
between the transmitter and receiver can be calculated numerically [8].

Figure 1: System model for OAM-like multiplexing based LIS communications.

It is well-known that the optimal communication modes are provided by the
solutions of the eigenproblem in Eq. (3) [9]. Their numerical values can be
obtained by discretizing the surfaces $\mathcal{S}_{t}$ and $\mathcal{S}_{r}$,
reformulating the Eq. (3) under a matrix form and calculating the SVD [5].
Although optimal, benefiting from these modes is challenging because it needs
to handle different continuous values of amplitudes and phases in the
precoder. Moreover, the exact knowledge of the CSI (i.e. Tx and Rx antenna
rotations and relative distance) is required to obtain the eigenfuntions. The
next section proposes alternative precoding schemes to tackle these
limitations.

##  III Alternative designs of transmission modes

We present the different options of hardware-friendly transmission modes for
the disk-shaped Tx LIS, with minimal requirements on the CSI. Following that
logic, we introduce two families of orthogonal modes:

  * •

OAM-based modes which depend on continuous phases only with amplitude being
constant;

  * •

Polar Walsh functions which depend on binary phases (i.e. $0$ and $\pi$) with
amplitude being constant.

We define the polar coordinates on the Tx LIS as
$(\rho,\theta)=(\sqrt{x^{2}+y^{2}},\arctan\frac{y}{x})$.

###  III-A OAM-based Transmission Modes

In [6], two different types of OAM transmissions, i.e., focused and unfocused,
are introduced and studied in terms of DoFs. The unfocused OAM modes are
defined for $n\in\mathbb{N}$ as:

$\phi_{n}(\rho,\theta)=\begin{cases}\dfrac{1}{\sqrt{\pi}r_{t}}e^{-\jmath\frac{n-1}{2}\theta},&\;\;\text{for
$n$ even},\\\
\dfrac{1}{\sqrt{\pi}r_{t}}e^{\jmath\frac{n}{2}\theta},&\;\;\text{for $n$
odd.}\end{cases}$ (5)

This expression is a straightforward extension to the disk case of the OAM
modes for UCAs [4]. Unfocused OAM modes are interesting because they do not
need to know the distance $d$ between the Tx and Rx antennas.

The scenario in which $d$ is known is in favor of focused OAM modes, whose
expression at the $n$-th state is:

$\phi_{n}(\rho,\theta)=\begin{cases}\dfrac{1}{\sqrt{\pi}r_{t}}e^{\jmath\kappa
d(1+\frac{\rho^{2}}{2d^{2}})}e^{-\jmath\frac{n-1}{2}\theta},&\;\;\text{for $n$
even},\\\ \dfrac{1}{\sqrt{\pi}r_{t}}e^{\jmath\kappa
d(1+\frac{\rho^{2}}{2d^{2}})}e^{\jmath\frac{n}{2}\theta},&\;\;\text{for $n$
odd.}\end{cases}$ (6)

Their expressions are based on an approximation of the exact CPSWFs that
enables constant amplitude and preserves orthogonality [10]. To be specific,
the phase varies also along the radial coordinate, so that the focal point is
located at a distance $d$ from the Tx LIS. Both types of OAM modes are
illustrated in Fig. 2 for $0\leq n\leq 3$.

In this paper, we aim at augmenting the hardware-friendly simplification, and
introduce a binary quantization of the modes. We will see that both radial and
angular coordinates are separate functions of the original linear Walsh
functions, and overall take only binary values $\pm\frac{1}{\sqrt{\pi}r_{t}}$.
A cost-efficient and low-complexity implementation can be expected
accordingly.

###  III-B Our Contribution: Transmission with Polar Walsh Functions

We present the Walsh functions in polar domain with separation in radial and
angular coordinates [11]. The orthogonality of a series of polar Walsh
functions is written as

$\displaystyle\int_{0}^{r_{t}}\int_{0}^{2\pi}\phi^{\mathcal{W}}_{m,n}(\rho,\theta)\phi^{\mathcal{W}}_{m^{\prime},n^{\prime}}(\rho,\theta)\rho\,\mathrm{d}\rho\,\mathrm{d}\theta$
$\displaystyle=\delta_{mm^{\prime}}\delta_{nn^{\prime}},$ (7)
$\displaystyle\;\forall m,m^{\prime},n^{\prime},n^{\prime},$

where $\delta_{mm^{\prime}}$ and $\delta_{nn^{\prime}}$ are Kronecker delta
functions. The index $m$ (resp $n$) of $\phi^{\mathcal{W}}_{m,n}$ stands for
the communication mode state according to the radial (resp. angular)
dimension. It can be shown from [11] that, for some positive integers $\mu$
and $\nu$, $\\{\phi^{\mathcal{W}}_{m,n}\\}$ forms an orthonormal basis for
$0\leq m\leq 2^{\mu}$ and $0\leq n\leq 2^{\nu}$ with:

$\phi^{\mathcal{W}}_{m,n}(\rho,\theta)=\dfrac{1}{\sqrt{\pi}r_{t}}.\prod\limits_{k=0}^{\mu-1}\mathrm{sgn}\bigg{[}\cos\Big{(}2^{k}\pi\left(\frac{\rho}{r_{t}}\right)^{2}m_{k}\Big{)}\bigg{]}.\\\
\prod\limits_{k=0}^{\nu-1}\mathrm{sgn}\bigg{[}\cos\Big{(}2^{k}\frac{\theta}{2}n_{k}\Big{)}\bigg{]},$
(8)

where $\mathrm{sgn}(\cdot)$ returns the sign of the argument,
$m=\sum\limits_{k=1}^{\mu-1}m_{k}2^{k}$ with $m_{k}$ being the binary bits (0
or 1) of the binary numeral for $m$, ditto for
$n=\sum\limits_{k=1}^{\nu-1}n_{k}2^{k}$. The polar Walsh function are
illustrated in Fig. 2 for $\mu=\nu=2$.

Among the complete set of functions, the radial Walsh functions are invariant
with respect to $\theta$, and generated by choosing $\nu=0$. Likewise, the
angular Walsh functions are invariant with respect to $\rho$, and generated by
choosing $\mu=0$.

To the best of our knowledge, these three sets of functions are the only ones
that (a) are defined on the disk, (b) are described with the minimal number of
values – two only – and (c) form an orthonormal basis. Finally it is also
independent from geometric parameters of the configuration. We emphasize that
polar Walsh functions should not be strictly referred to as OAM modes, because
they do not physically produce helical polarization of the EM waves.

Figure 2: Phase plots of three orthogonal functions for holographic disk
antennas: OAM unfocused, OAM focused ($r_{t}=5\lambda$, $d=20\lambda$), and
polar Walsh functions. Angular Walsh functions and radial Walsh functions are
a subset of the complete polar Walsh functions.

##  IV Capacity Formulation

Literature related to OAM-based holographic MIMO focuses on numerical analysis
of the DoFs. However, inter-mode interference may occur after propagation to
the Rx when transmission modes are not the eigenfunctions from Eq. (3). Here,
we focus on the capacity metric with different transmission schemes.

Since the LIS aperture is continuous, we need to perform spatial sampling on
the two-dimensional (2D) disk as to make the calculation of the channel
capacity feasible. We assume that $N_{t}$ and $N_{r}$ points uniformly samples
the Tx and Rx respectively, with a half-wavelength step size both $x$-wise and
$y$-wise222It is mentioned in [8] that half-wavelength sampling is sufficient
for equivalently modeling the discretized antenna aperture as its continuous
counterpart without loss of accuracy by following Nyquist space sampling
theorem.. Thus, we can formulate the Tx-Rx propagation channel with a matrix
${\mathbf{H}}\in\mathbb{C}^{N_{r}\times N_{t}}$, whose entries are calculated
based on the positions of the Tx and Rx elements and the model Eq. (2).
Therefore, the channel capacity $C$ in [bit/s/Hz] is written as

$C=\log_{2}\det\bigg{(}{\mathbf{I}}_{N_{r}}+\frac{P_{t}}{\sigma^{2}}{\mathbf{H}}{\mathbf{F}}\mathbf{Q}{\mathbf{F}}^{\mathsf{H}}{\mathbf{H}}^{\mathsf{H}}\bigg{)},$
(9)

where $\dfrac{P_{t}}{\sigma^{2}}$ denotes the received signal-to-noise ratio
(SNR), the $N_{t}\times N$ matrix ${\mathbf{F}}$ contains a collection of $N$
selected modes (similar to precoding in the conventional MIMO network), and
${\mathbf{Q}}$ is the $N\times N$ symmetric matrix that allocates the power to
each mode, such that $\operatorname{tr}({\mathbf{Q}})\leq 1$. In order to
analyze the capacity at constant received SNR, the total transmitted power is
adjusted to compensate the path loss, as

$P_{t}=\text{SNR}\cdot(4\pi d\sigma)^{2}.$ (10)

When the channel matrix ${\mathbf{H}}$ is known, the capacity is maximized
when ${\mathbf{F}}$ consists of the right singular vectors from
${\mathbf{H}}$, and when the power is distributed by water-filling [9]. In
that specific case $\mathbf{Q}$ forms a diagonal matrix. Another sub-optimal
yet reasonably efficient approach, named equal power allocation (EPA),
distributes the power on each mode equally [12, Appendix C] with

$\mathbf{Q}=\dfrac{{\mathbf{I}}_{N}}{N}.$ (11)

EPA is interesting in the general case of any choice of precoder
${\mathbf{F}}$, including OAM modes and Walsh functions. Indeed, finding the
matrix $\mathbf{Q}$ that maximizes $C$ under constrained transmission power
deals with the power allocation problem, which goes beyond the scope of our
study. First and foremost, we want to compare the possible candidates
described in Sec. III to design the precoder ${\mathbf{F}}$ independently from
the power allocation subject. Therefore, the next section analyzes the
capacity under EPA with different types of precoders ${\mathbf{F}}$.

##  V Numerical Results

In this section, we evaluate the channel capacity of the different
transmission modes described in Sec. III. The configuration is set with
$r_{t}=r_{r}=10\lambda$, $\lambda=1$ cm and the SNR equals $-20$ dB. For the
study, the following transmission schemes are compared under EPA, such that
${\mathbf{F}}$ contains:

  * •

the first $N=16$ right singular vectors of ${\mathbf{H}}$;

  * •

the first $N=16$ OAM unfocused modes ($0\leq n\leq 15$ in Eq. (5));

  * •

the first $N=16$ OAM focused modes ($0\leq n\leq 15$ in Eq. (6));

  * •

the first $N=16$ radial Walsh functions ($\mu=4$, $\nu=0$);

  * •

the first $N=16$ angular Walsh functions ($\mu=0$, $\nu=4$);

  * •

the first $N=16$ polar Walsh functions ($\mu=2$, $\nu=2$).

The comparison is also made with the optimal water-filling case.

Figure 3: Channel capacity comparison with disk-shaped LISs. The $x$-axis is
the ratio between the Tx-Rx distance and Rayleigh distance.

The capacity is plotted in Fig. 3 as a function of distance $d$ scaled to the
Rayleigh distance $d_{r}=\dfrac{r_{t}^{2}}{2\lambda}$. The vertical dotted
line discriminates between near field and far field regions. As expected, the
SVD modes with water-filling dominate all the other candidate transmission
schemes. The bell-shaped trend with a maximum in the near field (here $\approx
136$ bit/s/Hz) is a well known behavior of the continuous apertures when they
are closely spaced [8]. With EPA, SVD still outperforms all the other schemes,
since it benefits from the perfect CSI.

Except for the radial Walsh function that substantially has the lowest
performance, all the other precoding designs attain a capacity in a range
going from $56$ to $74$ bit/s/Hz. The focused OAM functions bring a slight
improvement over the unfocused ones thanks to the knowledge of $d$,
particularly in the near field. However, polar and angular Walsh functions
attain capacities at the same order of magnitude than OAM modes from [6], with
a marginal degradation in the near-field regime.

Interestingly, going from complex-valued (OAM) to binary-valued (Walsh) modes
has little impact on the communication performance, whereas it potentially
leads to significant hardware simplifications of the Tx precoder [7]. Note
that these trends are given for $N=16$, and choosing other small values (e.g.
$N<25$) does not compromise the stated conclusions. As potential extension of
the current study, the authors suggest to investigate a power allocation
method with the given the precoding constraint to further improve the
capacity. Investigating the actual benefits to the hardware complexity is also
of importance, since the required element densification in holographic MIMO
may lead to very large values of $N_{t}$.

##  VI Conclusion

In this paper, we have proposed a cost-efficient and hardware-friendly
holographic MIMO transmission scheme for disk apertures based on polar Walsh
functions. Rather than focusing on DoFs, we opt for a numerical analysis of
the capacity. Even though the proposed scheme exhibited slight performance
degradation compared to the state-of-the-art OAM schemes, our proposed
approach is more readily to be implemented/realized in the future LIS-assisted
network. In the future, we will focus on the post-processing of the received
signals and conduct rigorous transmit power allocation in order to achieve the
full potential of OAM based LIS communications.

## References

  * [1] S. Hu, F. Rusek, and O. Edfors, “Beyond massive MIMO: The potential of data transmission with large intelligent surfaces,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2746–2758, 2018.
  * [2] C. Huang, S. Hu, G. C. Alexandropoulos, A. Zappone, C. Yuen, R. Zhang, M. D. Renzo, and M. Debbah, “Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends,” _IEEE Wireless Commun._ , vol. 27, no. 5, pp. 118–125, 2020.
  * [3] M. Cui and L. Dai, “Channel estimation for extremely large-scale MIMO: Far-field or near-field?” _IEEE Trans. Commun._ , vol. 70, no. 4, pp. 2663–2677, 2022.
  * [4] O. Edfors and A. J. Johansson, “Is orbital angular momentum (OAM) based radio communication an unexploited area?” _IEEE Trans. Antennas Propag._ , vol. 60, no. 2, pp. 1126–1131, 2012.
  * [5] D. Dardari, “Communicating with large intelligent surfaces: Fundamental limits and models,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2526–2537, 2020.
  * [6] G. Torcolacci, N. Decarli, and D. Dardari, “OAM-based holographic MIMO using large intelligent surfaces,” in _proc. IEEE Global Communications Conference_ , 2022, pp. 651–655.
  * [7] Z. Wang, M. Li, Q. Liu, and A. L. Swindlehurst, “Hybrid precoder and combiner design with low-resolution phase shifters in mmwave mimo systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 12, no. 2, pp. 256–269, 2018.
  * [8] Z. Wu, M. Cui, Z. Zhang, and L. Dai, “Distance-aware precoding for near-field capacity improvement,” _arXiv_ , 2022.
  * [9] E. Telatar, “Capacity of multi-antenna Gaussian channels,” _European transactions on telecommunications_ , vol. 10, no. 6, pp. 585–595, 1999.
  * [10] J. Xu, “Degrees of freedom of OAM-based line-of-sight radio systems,” _IEEE Trans. Antennas Propag._ , vol. 65, no. 4, pp. 1996–2008, 2017.
  * [11] L. Hazra and P. Mukherjee, _Self-similarity in Walsh functions and in the farfield diffraction patterns of radial Walsh filters_. Springer, 2017.
  * [12] T. L. Marzetta, E. G. Larsson, and H. Yang, _Fundamentals of massive MIMO_. Cambridge University Press, 2016\.