# Efficient single-stage optimization of islands in finite-$\beta$ stellarator
equilibria

C. B. Smiet<EMAIL_ADDRESS>J. Loizu  E. Balkovic  Ecole
Polytechnique Federale de Lausanne, Lausanne, Switzerland  A. Baillod
Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA

###### Abstract

We present the first single-stage optimization of islands in finite-$\beta$
stellarator equilibria. Stellarator optimization is traditionally performed as
a two-stage process; in the first stage, an optimal equilibrium is calculated
which balances a set of competing constraints, and in the second stage a set
of coils is found that supports said equilibrium. Stage one is generally
performed using a representation for the equilibrium that assumes nestedness
of flux surfaces, even though this is not warranted and occasionally
undesired. The second stage optimization of coils is never perfect, and the
mismatch leads to worse performing equilibria, and further deteriorates if
additional constraints such as force minimization, coil torsion or port access
are included. The higher fidelity of single-stage optimization is especially
important for the optimization of islands as these are incredibly sensitive to
changes in the field. In this paper we demonstrate an optimization scheme
capable of optimizing islands in finite $\beta$ stellarator equilibria
directly from coils. We furthermore develop and demonstrate a method to reduce
the dimensionality of the single-stage optimization problem to that of the
first stage in the two-stage approach.

††preprint: AIP/123-QED

One of the challenges in stellarator design is the vastness of the design
space; encompassing families of entirely different magnetic configurations
including but not limited to (quasi)-isodynamic configurations Goodman  _et
al._ (2023), quasisymmetric configurations Rodriguez (2022); Rodriguez  _et
al._ (2022), and isodrastic configurations Burby  _et al._ (2023). Designing a
stellarator equilibrium can be conceptualized as an optimization in this
infinite-dimensional space of possibilities, which is a challenge in its own
right, but has proven successful in experiment Beidler  _et al._ (2021), and
yielded many interesting configurations Goodman  _et al._ (2023); Landreman
and Paul (2022).

The magnetic field supporting the plasma equilibrium needs to be generated by
a set of coils external to it, which in a two-stage process is designed a
posteriori. This is done by matching the coil solution to the plasma solution
on the plasma boundary, by minimizing the squared total magnetic flux through
this boundary. The field generated by a finite number of discrete coils can
closely approximate the required field, but the match will never be exact.
Reduction of the total coil length exponentially increases the error in this
match due to coil ripple, and destroys the perfect confinement characteristics
of the optimized equilibrium Wechsung  _et al._ (2022). Construction costs on
the other hand increase with the quantity of conductor needed, and are reduced
by having fewer total coils of shorter length.

The requirements of a reactor place constraints on the coils that can
potentially degrade equilibrium performance but can be mitigated through a
single-stage approach. Neutral beams require locations with increased spacing
between the coils, view ports and feedtroughs must be placed where they are
needed, and device maintenance considerations set aditional constraints on the
coil geometry Brown  _et al._ (2015). A better reactor can be designed if the
coils are directly included in the optimization, such that these trade-offs
can be weighed, and so that the unavoidable coil ripple is distributed in a
way to minimize the deterioration of performance.

Single stage optimization, where the coils are included during the equilibrium
optimization, can be achieved in several different ways, comprehensively
presented in Henneberg  _et al._ (2021a). The most precise, but also costly
method is the direct-from-coil optimization where the coils are directly
varied and the equilibrium supported by those coils is optimized. Recently
Jorge  _et al._ presented optimizations where a finite-$\beta$ fixed-boundary
equilibrium was optimized in parallel with the coils, achieving better
performing stellarator equilibria than equivalent two-stage optimizations
Jorge  _et al._ (2023). Another single-stage optimimzation method was reported
by Giuliani  _et al._ (2022) which utilizes adjoint methods to rapidly
optimize vacuum fields for quasisymmetry. The methods’ efficiency is attested
by the 300,000 (and counting) coil configurations made available in the open-
source QUASR database Giuliani (2024).

When optimizing a stellarator equilibrium it is important to consider magnetic
islands, which can reduce the confinement in the core. Many equilibrium
solvers, including the widely used VMEC Hirshman and Whitson (1983), are
unable to represent magnetic islands. This can be overcome by coupling a code
capable of resolving islands and chaotic fields to the optimization as
demonstrated by Landreman  _et al._ (2021a) for vacuum configurations
($\beta=0$). Islands are extremely sensitive to small changes in the magnetic
field as was demonstrated in W7-X, where a field perturbation as small as
$1:10^{-5}$ caused a discernible change in the island structure Pedersen  _et
al._ (2016).

Whereas magnetic islands are detrimental to confinement in the core, they can
be essential in the plasma edge, where a reactor furthermore requires a
controlled method of handling the heat exhausted from the plasma. This can be
achieved by an island divertor where a large island chain is intentionally
generated, or a nonresonant divertor Bader  _et al._ (2017) where the last
closed flux surface transitions into a chaotic field. To determine the points
where the plasma interfaces with the wall, and hence the diveror geometry, the
equilibrium representation needs to represent more than nested surfaces.
Furthermore, ensuring that the coils do not interfere with the divertor can
only be accomplished by a single-stage optimization that is capable of
resolving magnetic islands and magnetic chaos.

In this paper we present a first demonstration of a single-stage optimization
capable of addressing all the issues above: Optimizing a free-boundary
finite-$\beta$ equilibrium directly from and self-consistent with a set of
coils, and targeting properties of the islands in the optimization. We do this
by removing the islands from around a $\langle\beta\rangle=1.49\%$ rotating
ellipse equilibrium. We build upon our recent work Baillod  _et al._ (2022);
Baillod (2023); baillod2023thesis where magnetic islands in a finite-$\beta$
free-boundary equibrium were optimized without consideration of coils. In this
paper we extend this to direct-from-coil optimization of islands in
stellarator equilibria. We further demonstrate a dimensionality-reduction
scheme that reduces the complexity and cost of the fully single-stage
optimization to that of the direct free-boundary optimization.

##  I MRxMHD Equilibria coupled to coils

In this work we use the Multi-region Relaxed MagnetoHydroDynamics (MRxMHD)
Bruno and Laurence (1996) equilibrium formulation, solved by the Stepped
Pressure Equilibrium Code (SPEC) Hudson  _et al._ (2011) to describe
stellarator equilibria. An MRxMHD equilibrium consists of a number of discrete
and nested ideal interfaces which sustain a finite pressure jump, and a force-
free field in between. Unlike formulations such as VMEC Hirshman and Whitson
(1983), that enforce continuously nested flux surfaces in their
representations MRxMHD allows for solutions that contain magnetic islands and
magnetic field line chaos whilst still constraining pressure and current
profiles. This makes SPEC a powerful tool to predict equilibrium $\beta$
limits Baillod (2023) and study linear and nonlinear, ideal and resistive
instabilities Loizu and Bonfiglio (2023); Loizu  _et al._ (2020); Kumar  _et
al._ (2021), and helical states in reversed field configurations Dennis  _et
al._ (2014); Liu  _et al._ (2024).

We take as starting point the same equilibrium presented in Baillod  _et al._
(2022), namely a rotating-ellipse field in which the pressure has been
increased to $\langle\beta\rangle=1.49\%$, shown in figure 1 (a). This
equilibrium has 5 field periods, is stellarator-symmetric Dewar and Hudson
(1998), has 7 pressure interfaces, and is contained in a computational
boundary with major radius 10m, minor radius 1m and an elongation of 2.5. For
a full description we refer to Baillod  _et al._ (2022); Baillod (2023).

In free-boundary mode, SPEC requires as an input the normal component of an
externally produced magnetic field $(\mathbf{B}_{\rm ext}\cdot\mathbf{n})$ on
a computational boundary. SPEC uses a spectral representation, here with a
poloidal resolution $M=6$ and toroidal resoltuion $N=8$. This means all
doubly-periodic odd functions (and in particular, the normal component of the
field on the boundary $(\mathbf{B}_{\rm ext}\cdot\mathbf{n})(\theta,\phi)$ )
are represented using a double Fourier sine series of the form:

$(\mathbf{B}_{\rm
ext}\cdot\mathbf{n})(\theta,\phi)=\sum_{n=1}^{N}V_{0,n}\sin(-nN_{P}\phi)+\\\
\sum_{m=1}^{M}\sum_{n=-N}^{N}V_{m,n}\sin(m\theta-nN_{P}\phi)$ (1)

(even functions are described by the analogous cosine series). Here $\phi$ is
a toroidal angle, $\theta$ a poloidal angle, and $N_{P}$ the number of field
periods.

In Baillod  _et al._ (2022) an optimzization was performed using as degrees-
of-freedom the representation of the external field, which with the resolution
$M=6$, $N=8$ constitutes the 110 Fourier coefficients of equation 1. In a
single-stage optimization, these Fourier components are not degrees-of-
freedom, but they are determined by the coils. The degrees-of-freedom then
become the parameters determining the geometry of the coils, from which the
normal field to the boundary is calculated using the Biot-Savart law.

$(7/5)$$(8/5)$$(9/5)$(a)

\begin{overpic}[width=433.62pt,percent]{start_coils.png} \put(5.0,75.0){ (b)}
\put(40.0,9.0){\large$(\mathbf{B}_{\rm coils}\cdot\hat{\mathbf{n}})/B_{0}$}
\end{overpic}

Figure 1: Initial unoptimized configuration and coils supporting this
equilibrium. a: Poincaré section of the field, with the $(m/n)=(7/5)$, $(8/5)$
and $(9/5)$ islands labeled. The SPEC interfaces are from outside in:
Computational boundary (black), plasma boundary (red), 6 internal interfaces
(teal) and the magnetic axis (teal dot). b: coils consistent with the
equilibrium in (a), and the normal component of the coil magnetic field on the
computational boundary. The color scale is identical to figures 2-4 and is
normalized to $B_{0}$, the magnetic field strength on the magnetic axis.

We first generate an initial coil set that supports the initial equilibrium.
The geometry of the coils is parametrized by three Fourier series in their
Cartesian $x$, $y$, and $z$-coordinates, with the curve describing the $i$-th
coil given by:

$\displaystyle x^{i}(\tau)$
$\displaystyle=X^{i}_{c,0}+\sum_{k=1}^{K}\left[X^{i}_{c,k}\cos(k\tau)+X^{i}_{s,k}\sin(k\tau)\right],$
(2) $\displaystyle y^{i}(\tau)$
$\displaystyle=Y^{i}_{c,0}+\sum_{k=1}^{K}\left[Y^{i}_{c,k}\cos(k\tau)+Y^{i}_{s,k}\sin(k\tau)\right],$
$\displaystyle z^{i}(\tau)$
$\displaystyle=Z^{i}_{c,0}+\sum_{k=1}^{K}\left[Z^{i}_{c,k}\cos(k\tau)+Z^{i}_{s,k}\sin(k\tau)\right].$

We generate a total of five unique coils with Fourier order $K=6$, that are
copied according to stellarator symmetry and field periodicity to form a total
of 50 coils. This results in 195 Fourier modes needed to represent a coil set.

The coils are initialized as circles enclosing the computational boundary, and
the coils are then optimized such that $(\mathbf{B}_{\rm
coils}\cdot\hat{\mathbf{n}})$ matches the value of $(\mathbf{B}_{\rm
ext}\cdot\hat{\mathbf{n}})$ from the initial equilibrium on the computational
boundary using a FOCUS-like algorithm Zhu  _et al._ (2018a). This algorithm
consists of the gradient-based minimization of the quadratic flux, regularized
by a penalty on total coil length, i.e. the minimization of:

$f_{\rm FOCUS}=\int_{\mathcal{S}}\left(\left(\mathbf{B}_{\rm
coils}-\mathbf{B}_{\rm
ext}\right)\cdot\mathbf{n}\right)^{2}\mathrm{d}s+\mathcal{L}^{2}_{>L_{0}}$ (3)

where the integral calculates the total quadratic flux through the surface
$\mathcal{S}$(= computational boundary) of the difference between a target
field $\mathbf{B}_{\rm ext}$ (given by the SPEC solution) and $\mathbf{B}_{\rm
coils}$ produced by the coils. $\mathcal{L}_{>L_{0}}$ is a one-sided penalty
function that penalizes if the total coil length exceeds $L_{0}$, and that is
zero below $L_{0}$, with $L_{0}=1000$m.

An optimum is achieved with a r.m.s. difference between the coil field and the
target field of $\delta B/B_{0}\simeq 10^{-7}$ with $B_{0}$ the magnetic field
strength on the magnetic axis. This produces a SPEC equilibrium
indistinguishable from the starting equilibrium in Baillod  _et al._ (2022).
In fact figure 1 (a) shows the equilibrium calculated from the coils. These
initial coils, as well as the normal field they produce on the boundary are
shown in figure 1 (b). We note that the initial field was constructed such
that the external field on the computational boundary is exactly opposite to
the plasma field (and thus the computational boundary coincides with a
magnetic surface) and hence figure 1 (b) shows a non-zero amplitude of the
coil field.

##  II Optimizing the equilibrium

We define our optimization problem with the goal of removing the islands
labeled in figure 1. For this we use Cary and Hanson’s method Cary and Hanson
(1986) of reducing the Greene’s Residue Greene (1979, 1968). Greene’s residue,
$\mathcal{R}_{m,n}$, of an $(m,n)$ island chain (with $m$ the poloidal and $n$
the toroidal mode number) is a quantity that is calculated at the x- or
o-point of said island chain, by calculating the linearized mapping
$\mathsf{M}$ of the field line through the x- or o-point. The Residue is
defined as
$\mathcal{R}_{m,n}=\tfrac{1}{2}-\tfrac{1}{4}\textrm{Tr}(\mathsf{M})$. For
details see for example Greene (1968); Smiet  _et al._ (2019). Greene’s
residue is larger than zero for o-points, and smaller than zero for x-points,
and when the Residue is zero, the $(m,n)$ island chain has reduced to the
unbroken $m/n$ magnetic surface. The x- and o-points are fixed points of the
Poincaré map applied $m$ times, and these are found using a Newton method. We
target the residues of the $\mathcal{R}_{9/5}$, $\mathcal{R}_{8/5}$ and
$\mathcal{R}_{7/5}$ island chains, that corresponds to the labeled islands in
figure 1.

We perform the optimizations using the simsopt stellarator optimization code
Landreman  _et al._ (2021b) that has recently been coupled to SPEC. We define
our degrees-of-freedom for the optimization problem in two different ways: A
_brute-force_ optimization using all coil degrees-of-freedom, and a reduced-
dimensional optimization that we describe in section III. For the _brute-
force_ optimization, the degrees-of-freedom are the 195 Fourier coefficients
that describe the coils according to equation 2. The optimization target is a
least-squares functional with the following terms:

$f=\mathcal{R}_{9/5}^{2}+\mathcal{R}_{8/5}^{2}+\mathcal{R}_{7/5}^{2}+\\\
0.1(\mathcal{L}-L_{0})^{2}+0.1(\mathcal{C}-0.1\mathcal{C}_{0}^{2})$ (4)

where $\mathcal{L}$ is the total coil length, $L_{0}=1010$m, a little larger
than the starting coils, and $\mathcal{C}$ is the sum of the mean squared
curvatures of the coils for which $\mathcal{C}_{0}$ is its initial value.

The equilibrium is optimized using a trust-region-reflective least-squares
optimization with the Jacobian calculated by a one-directional finite-
difference evaluation of the components of equation 4. The step size for the
finite-difference evaluation is set to $10^{-4}$, and the trust region is
decreased for degrees-of-freedom corresponding with higher Fourier orders $k$
according to the relation $d_{trust}=1/20e^{0.3k}$.

The last two terms of equation 4 were found by trial-and-error, and were found
to be necessary for the optimization to converge and to avoid non-convergence
of the equilibrium solver when providing unsolvable fields. Inspection of the
steps taken in optimizations without these terms shows that the optimizer
attempts unfeasible configurations, with wild coils caused by high amplitude
high order Fourier modes. Since both coil length and coil complexity are
strongly increased by the higher order Fourier modes, their inclusion in the
functional helps taming the steps attempted by the optimizer. The order-
dependent trust region has a similar effect, and with these measures in place,
the optimization concludes successfully.

The optimization achieves an improved coil configuration in which the targeted
islands are greatly reduced in size. The Poincaré section of this field
(figure 2 (a) ) shows very thin remnants of the islands. The optimized coils
(2 (b), black curves) are only shifted by a small amount with respect to the
starting coils (2 (b) and 1 (b) white curves). Figure 2 (b) also shows the
computational boundary, colored by the normal component of the coil magnetic
field on it, using the same color scale as figure 1. The optimized coils
produce a stronger field on the computational boundary than the initial coils,
consistent with the observation that the computational boundary no longer
coincides with a magnetic surface. We also note that the optimization has
starkly reduced the magnetic stochasticity in the Poincaré plot, not only in
and around the targeted surfaces, but also in the core of the equilibrium.

\begin{overpic}[width=433.62pt,percent]{COIL_poincare.png} \put(5.0,55.0){
(a)} \end{overpic}

\begin{overpic}[width=433.62pt,percent]{COIL_coils.png} \put(5.0,55.0){ (b)};
\put(40.0,9.0){\large$(\mathbf{B}_{\rm coils}\cdot\hat{\mathbf{n}})/B_{0}$}
\end{overpic}

Figure 2: Optimized equilibrium and coils using the _brute-force_ optimization
scheme where the degrees-of-freedom are all elements of the coil
parametrization. a) Poincaré of the equilibrium achieved after optimization.
b) The coils that achieve this optimum (black) and the initial coils (white),
plotted around the computational boundary (blue-red) colored by the magnitude
of the coil field field normal to the surface. The color scale is identical to
figures 1,3 and 4.

The problems encountered in this _brute-force_ approach however hint at a
fundamental issue with single-stage optimization: The optimization problem has
195 degrees-of-freedom. The free-boundary equilibrium, however, depends only
on the field on the computational boundary, which is represented by equation 1
with finite Fourier resolution of 110 unique modes. This implies that the
optimization has an inherent null space, as there will be many linear
combinations of coil degrees-of-freedom that cannot be represented by the
limited Fourier resolution on the boundary. Not only the limited resolution is
at issue here; it is a known issue that the coil optimization problem is ill-
posed; there are variations of the coil geometry that do not affect the field
on the boundary Landreman (2017). This null space hinders the optimization
because it allows for the unfeasible steps taken by the optimizer: variations
that do not affect the equilibrium can be freely added to the step.

##  III Reducing the dimensionality of the optimization space

In order to mitigate the issues encountered in the _brute-force_ optimization,
we apply a transformation to the basis of the optimization space. The mapping
from the coil degrees-of-freedom to the magnetic field on the computational
boundary is a mapping from an $I=195$ to a $J=110$-dimensional space:

$g:\mathbb{R}^{I}\rightarrow\mathbb{R}^{J}$ (5)

$g$ can be linearized around the starting coil set defined by the state vector
$\mathbf{x}_{\rm c}\in\mathbb{R}^{I}$:

$g(\mathbf{x}_{\rm c}+\delta\mathbf{x})\simeq g(\mathbf{x}_{\rm
c})+\mathsf{\mathsf{J}}\delta\mathbf{x}$ (6)

where the Jacobian of the mapping $\mathsf{\mathsf{J}}_{ij}=\partial
g^{(j)}/\partial\mathbf{x}_{\rm c}^{(i)}$ (note that $(i)$ and $(j)$ are
dimension labels) is calculated using one-sided finite-differencing on the
evaluation of $g$ with step size $10^{-7}$.

The Jacobian $\mathsf{\mathsf{J}}$ is a non-square matrix that encodes the
effect of coil displacements on the field on the computational boundary. We
now perform singular value decomposition (SVD) on this Jacobian matrix:

$\mathsf{\mathsf{J}}=\mathsf{\mathsf{U}}\mathsf{\mathsf{\Sigma}}\mathsf{\mathsf{V}}^{*}.$
(7)

SVD decomposes $\mathsf{\mathsf{J}}$ into two square, real, and orthogonal
matrices $\mathsf{\mathsf{U}}$ ($J\times J$) and $\mathsf{\mathsf{V}}^{*}$
($I\times I$) and a diagonal matrix $\mathsf{\mathsf{\Sigma}}$ containing the
singular values. This decomposition has profound significance to the physical
problem at hand; rows of $\mathsf{\mathsf{V}}^{*}$ form an orthonormal basis
of ‘coil space’, and the $i$-th row of $\mathsf{\mathsf{V}}^{*}$ (right-
singular vector of $\mathsf{\mathsf{J}}$) is a linear combination of Fourier
modes (a displacement of the coils) that effects a change in the normal field
given by the $i$-th column of $\mathsf{\mathsf{U}}$, with a transfer
efficiency given by the singular value $\mathsf{\mathsf{\Sigma}}_{ii}$. The
SVD algorithm conveniently ranks the singular values from high to low, with
the highest singular-value corresponding to the right-singular vector
embodying the coil displacement that most efficiently effects change to the
field on the computational boundary. Because $I>J$, there are many right-
singular vectors that have zero singular value; i.e. these are the
displacements that do not influence the field on the computational boundary,
the null space of our optimization problem!

Zhu  _et al._ (2018b) used similar methods in analyzing coil sets at an
optimum; they used the singular values of a Hessian matrix to quantify coil
sensitivity and identify critical tolerances in construction Zhu  _et al._
(2018b, 2019).

We now redefine our optimization problem to use the non-zero right-singular
vectors as degrees-of-freedom. Further reduction by limiting to a subset of
the right-singular vectors could further speed up the computation and will be
considered in future work. Using all right-singular vectors we have reduced
the dimensionality of the problem back to the same dimensionality as needed to
represent a doubly-periodic function in our chosen Fourier resolution, i.e.
110 degrees-of-freedom. This is the minimal dimensionality for a first stage
optimization in a two-stage process, as it is the minimal dimensionality of a
fixed-boundary representation Henneberg  _et al._ (2021b) with similar
spectral complexity. Our reduced representation furthermore removes coil
displacements that do not affect the field on the boundary, i.e. it removes
the inherent null space in coil optimization Landreman (2017).

In this new optimization space we no longer need the regularizing penalties on
total coil length and coil complexity, and we optimize the stellarator
equilibrium with the least-squares functional:

$f=\mathcal{R}_{9/5}^{2}+\mathcal{R}_{8/5}^{2}+\mathcal{R}_{7/5}^{2},$ (8)

without the terms on the coils present in equation (4).

In the optimized equilibrium, the targeted island chains can no longer be seen
at all, as shown in figure 3 (a). The entire equilibrium is also shifted from
its original position, and the computational boundary is clearly no longer a
flux surface. The coils are not significantly displaced from their starting
positions, as shown in figure 3 (b). The normal component of the coil field on
the computational boundary, shown by the color in figure 3 (b) is increased
over the _brute-force_ optimization. The achieved optimum is an order of
magnitude more performant than the _brute-force_ approach, as can be seen from
the optimization results in table 1. From this table we can also see that the
coil complexity is increased by less than 2% and the total coil length is even
reduced, though these were not optimization targets.

The reduced-dimensional optimization space has significantly increased the
robustness of the calculation. The _brute-force_ optimization functional in
equation (4) was found through trial-and-error. Optimization functionals
without the stabilizing weights on coil complexity and length would not
converge, get stuck in local minima, or lead to crashes of the equilibrium
solver. The reduced-dimensional optimization has no need for stabilizing
terms.

\begin{overpic}[width=433.62pt,percent]{REDUCED_poincare.png} \put(5.0,55.0){
(a)} \end{overpic}

\begin{overpic}[width=433.62pt,percent]{REDUCED_coils.png} \put(5.0,55.0){
(b)} \put(40.0,9.0){\large$(\mathbf{B}_{\rm
coils}\cdot\hat{\mathbf{n}})/B_{0}$} \end{overpic}

Figure 3: Optimized equilibrium and coils using the dimensionality-reduced
optimization scheme. a: Poincaré section of the field. b: The coils that
achieve this optimum (black) and the initial coils (white), plotted around the
computational boundary (blue-red) colored by the magnitude of the field normal
to the surface. The color scale is identical to figures 1,2 and 4.

##  IV Comparison with two-stage optimization

We finally demonstrate the advantage of single-stage optimization by
performing a similar optimization in two stages. To increase the realism of
this example, we generate initial coils that are 40% shorter (coils closer to
the plasma than the coils in section I), with a total length of $600$m. For
the first stage of this optimization, we repeat the island optimization
procedure laid out in Baillod  _et al._ (2022). Our degrees-of-freedom are the
110 Fourier components describing the magnetic field normal to the
computational boundary, and the least-squares functional is identical to
equation (8). Similar to the results of Baillod  _et al._ (2022), this
optimization yields an equilibrium in which two new island chains appear with
mode numbers $(m,n)=(10,5)$ and $(6,5)$ respectively. A second optimization is
performed on this intermediate optimum with the least-squares functional

$f=\mathcal{R}_{10/5}^{2}+\mathcal{R}_{9/5}^{2}+\mathcal{R}_{8/5}^{2}+\mathcal{R}_{7/5}^{2}+\mathcal{R}_{6/5}^{2}.$
(9)

The optimization runs successfully and the Poincaré section of this field is
shown in figure 4 (a) (top half). The island size is reduced, but the optimum
achieved is not as performant as the single-stage or the reduced single stage
optimizations.

Table 1: Optimization Results Optimization | $\mathcal{R}_{9/5}$ | $\mathcal{R}_{8/5}$ | $\mathcal{R}_{7/5}$ | $\mathcal{L}$ | $\mathcal{C}$
---|---|---|---|---|---
initial configuration | 7.7$\times 10^{-2}$ | -2.11$\times 10^{-1}$ | 2.95$\times 10^{-1}$ | 1000 | 0.995
_brute-force_ optimization | 3.2$\times 10^{-3}$ | -1.78$\times 10^{-2}$ | 5.4 $\times 10^{-2}$ | 1010 | 0.847
reduced-dimensional optimization | -5.32$\times 10^{-4}$ | 7.78$\times 10^{-3}$ | -2.35$\times 10^{-3}$ | 997.7 | 1.02
first stage (shorter coils) | -4.23$\times 10^{-3}$ | -1.07$\times 10^{-2}$ | 7.77$\times 10^{-3}$ |  |
second stage (shorter coils) | 2.61$\times 10^{-2}$ | -5.81$\times 10^{-2}$ | 5.22$\times 10^{-2}$ | 600 | 1.7
single-stage (shorter coils) | 1.49$\times 10^{-2}$ | 8.97$\times 10^{-3}$ | 9.39$\times 10^{-3}$ | 600.2 | 1.72

This first-stage optimized equilibrium is then used as a target for a second-
stage coil optimization which seeks to match the field on the computational
boundary. We start from a set of circular coils, and we use the FOCUS-like
algorithm described in section I. To prevent the optimizer from getting stuck
in a local minimum, the optimization is performed four times; first with
$\mathcal{L}_{0}=600$, then the optimized coils are re-optimized with
$\mathcal{L}_{0}=610$, and then twice more with $\mathcal{L}_{0}=600$.

\begin{overpic}[width=433.62pt,percent]{TWOSTAGE_poincare.png} \put(5.0,55.0){
(a)} \put(5.0,52.0){ first stage} \put(5.0,5.0){ second stage}
\end{overpic}\begin{overpic}[width=433.62pt]{SHORTER_poincare.png}
\put(5.0,55.0){ (c)} \put(5.0,5.0){single-stage} \end{overpic}

\begin{overpic}[width=433.62pt]{TWOSTAGE_coils.png} \put(5.0,55.0){(b)}
\put(25.0,40.0){second stage} \put(65.0,40.0){target field}
\put(40.0,9.0){\large$(\mathbf{B}_{\rm
coils,ext}\cdot\hat{\mathbf{n}})/B_{0}$}
\end{overpic}\begin{overpic}[width=433.62pt]{SHORTER_coils.png}
\put(5.0,55.0){(d)} \put(40.0,9.0){\large$(\mathbf{B}_{\rm
coils}\cdot\hat{\mathbf{n}})/B_{0}$} \end{overpic}

Figure 4: Comparison of single-stage and two-stage optimization of islands in
stellarator equilibria. a: Poincaré section of the single-stage optimized
field (top half a), and of the second-stage field consistent with the coils
(bottom half a). The islands that were mostly removed in the first stage are
inadvertently re-introduced in the second stage. b) coils of the two-stage
optimized field (left side) and the resultant normal field to the
computational boundary, and the target thereof from the single-stage
optimization (right half). c) Poincaré section of the single-stage optimized
field. d: single-stage optimized coils. The color scale in c and d is
identical to figures 1-3.

The Poincaré section of the equilibrium after the second stage optimization is
shown in figure 4 (a) (bottom half). We see that the islands that were mostly
removed in the first stage optimization are re-introduced in the second stage.
The coils from the second stage are shown in figure 4 (b) (left half),
together with the normal field they produce on the computational boundary. The
right half of figure 4 (b) shows the target field, namely the result of the
first stage that needs to be matched by the second stage, but a close match is
not possible due to unavoidable coil ripple.

We compare this result to the equivalent single-stage optimization. First,
coils are generated for the initial equilibrium with a single optimization
with $\mathcal{L}_{0}=600$. We then use our dimensionality reduction scheme to
generate a 110-dimensional optimization space from this initial coil position.
The target functional for the single-stage optimization consists only of the
three initial island chain Residues (i.e. equation (8)). The optimization is
performed similarly to the other reduced optimization, and runs to completion
after two Jacobian evaluations.

The reduced optimization of the shorter coil set is successful, resulting in a
reduced island width visible in the Poincaré section (figure 4 (c) ). This
single-stage optimization yields better results than the two-stage process;
table 1 shows that the residues are roughly halved. The coils of the single-
stage optimization show significant ripple on the computational boundary, as
shown in figure 4 (d), comparable to the ripple of the two stage optimized
coils, yet the island residues are significantly decreased. We postulate that
this is due to the fact that the optimizer can distribute the coil ripple on
the computational boundary in a way that has minimal effect on the islands,
whereas when the second-stage targets an already optimized field, the mismatch
is distributed randomly, leading to a larger effect on the islands. We also
note that the Poincaré section of the single-stage optimized field shows
significantly reduced chaos, also in the core of the equilibrium near the
axis, hinting that the single-stage process changes the field in a more robust
way, and is able to find solutions that the two-stage process cannot.

##  V Conclusions and discussion

We have demonstrated a first of a kind single-stage optimization of a
finite-$\beta$ stellarator equilibrium with islands. We have furthermore shown
how the dimensionality of a single-stage optimization can be reduced to that
of the two-stage problem, and that this reduction scheme improves the
optimization performance by removing a null space.

All optimizations were performed on a cluster using a total of 288 cores per
run, split into 24 groups of 12 cores. The brute-force optimization took 4000
CPU-hours to complete, whereas the reduced-dimensional optimization took 3500
cpu-hours to terminate. The optimizations are terminated when the step size
becomes too small. The reduced-dimensional optimization evaluated
significantly more finite-difference Jacobians (the computationally most
costly step), and as a result reached a significantly better optimum
configuration. The total wall clock time was on the order of 10 hours.

This dimensionality reduction scheme has broader applicability than the
problem at hand, and we expect it to to enhance most single-stage
optimizations, as the number of degrees-of-freedom to describe a coil set are
generally larger than those needed to describe the plasma boundary. It can be
further improved in many ways which will be the subject of future research,
including but not limited to:

  1. 1.

Using only a subset of the singular vectors.

  2. 2.

Biasing the Jacobian calculation so the most relevant (to the problem at hand)
coil displacements have largest singular values.

  3. 3.

Projecting the reduced degrees-of-freedom into a subspace that conserves
desirable coil metrics such as total length or coil complexity.

This scheme is furthermore not limited to any representation of the coils. The
Fourier representation in (1) can be replaced with any other, for example
spline-based Lonigro and Zhu (2022); Takahashi (2024), coil representations.
Although splines have the advantages of having local support, less degeneracy
in the representation, and similar magnitude effect on the produced field for
each degree-of-freedom, coils represented in this way will still have the
inherent null space consisting of displacements that do not affect the field
at the plasma Landreman (2017); Landreman and Boozer (2016).

The dimensionality reduction scheme is most efficient when one can assume the
that the linearization embodied by the Jacobian $\mathsf{\mathsf{J}}$ remains
valid, i.e. the coils are only displaced a small amount. This makes the method
very well suited for the problem of removing islands for which only minor
adjustments to the coils are sufficient. If the coils need to be significantly
displaced, or if the coils are optimized starting from a simplified (’cold
start’) configuration, the linearization on which we performed the SVD can no
longer be assumed to be valid. In such a case, the mapping could be re-
calculated regularly. Since the calculation of the Jacobian
$\mathsf{\mathsf{J}}$ (from coil degrees-of-freedom to field-on-boundary) does
not involve free-boundary solves and can be evaluated in seconds, this can
still significantly speed up an optimization.

The capability to include islands in a single-stage finite-$\beta$ stellarator
optimization opens up new possibilities for stellarator design. We
demonstrated how this scheme can be used to remove unwanted islands, but
specific islands can now also be intentionally introduced. One example is the
optimization of either an island divertor or a nonresonant divertor Bader  _et
al._ (2017). The capability to resolve the magnetic structure at finite
$\beta$ directly from from the coils that support it, will make it possible to
enforce the necessary space for a divertor solution, and ensure its proper
functioning in future reactors.

We have demonstrated our single-stage optimization in a simplified geometry,
consisting of a rotating ellipse. In future work we will apply our techniques
to more strongly shaped and reactor-relevant stellarator geometries

###### Acknowledgements.

This work was supported by a grant from the Simons Foundation (1013657, JL).
This work has been carried out within the framework of the EUROfusion
Consortium, via the Euratom Research and Training Programme (Grant Agreement
No 101052200 — EUROfusion) and funded by the Swiss State Secretariat for
Education, Research and Innovation (SERI). Views and opinions expressed are
however those of the author(s) only and do not necessarily reflect those of
the European Union, the European Commission, or SERI. Neither the European
Union nor the European Commission nor SERI can be held responsible for them.

## References

  * Goodman  _et al._ (2023) A. Goodman, K. C. Mata, S. A. Henneberg, R. Jorge, M. Landreman, G. Plunk, H. Smith, R. Mackenbach, C. Beidler, and P. Helander, Constructing precisely quasi-isodynamic magnetic fields, Journal of Plasma Physics 89, 905890504 (2023).
  * Rodriguez (2022) E. Rodriguez,  _Quasisymmetry_ , Ph.D. thesis, Princeton University (2022).
  * Rodriguez  _et al._ (2022) E. Rodriguez, E. Paul, and A. Bhattacharjee, Measures of quasisymmetry for stellarators, Journal of Plasma Physics 88, 905880109 (2022).
  * Burby  _et al._ (2023) J. Burby, R. MacKay, and S. Naik, Isodrastic magnetic fields for suppressing transitions in guiding-centre motion, Nonlinearity 36, 5884 (2023).
  * Beidler  _et al._ (2021) C. Beidler, H. Smith, A. Alonso, T. Andreeva, J. Baldzuhn, M. Beurskens, M. Borchardt, S. Bozhenkov, K. J. Brunner, H. Damm, _et al._ , Demonstration of reduced neoclassical energy transport in wendelstein 7-x, Nature 596, 221 (2021).
  * Landreman and Paul (2022) M. Landreman and E. Paul, Magnetic fields with precise quasisymmetry for plasma confinement, Physical Review Letters 128, 035001 (2022).
  * Wechsung  _et al._ (2022) F. Wechsung, M. Landreman, A. Giuliani, A. Cerfon, and G. Stadler, Precise stellarator quasi-symmetry can be achieved with electromagnetic coils, Proceedings of the National Academy of Sciences 119, e2202084119 (2022).
  * Brown  _et al._ (2015) T. Brown, J. Breslau, D. Gates, N. Pomphrey, and A. Zolfaghari, Engineering optimization of stellarator coils lead to improvements in device maintenance, in  _2015 IEEE 26th Symposium on Fusion Engineering (SOFE)_ (IEEE, 2015) pp. 1–6.
  * Henneberg  _et al._ (2021a) S. A. Henneberg, S. R. Hudson, D. Pfefferlé, and P. Helander, Combined plasma–coil optimization algorithms, Journal of Plasma Physics 87, 905870226 (2021a).
  * Jorge  _et al._ (2023) R. Jorge, A. G. Goodman, M. Landreman, J. Rodrigues, and F. Wechsung, Single-stage stellarator optimization: Combining coils with fixed boundary equilibria, Plasma Physics and Controlled Fusion (2023).
  * Giuliani  _et al._ (2022) A. Giuliani, F. Wechsung, G. Stadler, A. Cerfon, and M. Landreman, Direct computation of magnetic surfaces in boozer coordinates and coil optimization for quasisymmetry, Journal of Plasma Physics 88, 905880401 (2022).
  * Giuliani (2024) A. Giuliani, Direct stellarator coil design using global optimization: application to a comprehensive exploration of quasi-axisymmetric devices, Journal of Plasma Physics 90, 905900303 (2024).
  * Hirshman and Whitson (1983) S. P. Hirshman and J. Whitson,  _Steepest descent moment method for three-dimensional magnetohydrodynamic equilibria_ , Tech. Rep. (Oak Ridge National Lab.(ORNL), Oak Ridge, TN (United States), 1983).
  * Landreman  _et al._ (2021a) M. Landreman, B. Medasani, and C. Zhu, Stellarator optimization for good magnetic surfaces at the same time as quasisymmetry, Physics of Plasmas 28 (2021a).
  * Pedersen  _et al._ (2016) T. S. Pedersen, M. Otte, S. Lazerson, P. Helander, S. Bozhenkov, C. Biedermann, T. Klinger, R. C. Wolf, and H.-S. Bosch, Confirmation of the topology of the wendelstein 7-x magnetic field to better than 1: 100,000, Nature communications 7, 13493 (2016).
  * Bader  _et al._ (2017) A. Bader, A. Boozer, C. Hegna, S. Lazerson, and J. Schmitt, Hsx as an example of a resilient non-resonant divertor, Physics of Plasmas 24 (2017).
  * Baillod  _et al._ (2022) A. Baillod, J. Loizu, J. Graves, and M. Landreman, Stellarator optimization for nested magnetic surfaces at finite $\beta$ and toroidal current, Physics of Plasmas 29 (2022).
  * Baillod (2023) A. Baillod,  _Equilibrium pressure limits in stellarators_ , Tech. Rep. (EPFL, 2023).
  * Bruno and Laurence (1996) O. P. Bruno and P. Laurence, Existence of three-dimensional toroidal mhd equilibria with nonconstant pressure, Communications on pure and applied mathematics 49, 717 (1996).
  * Hudson  _et al._ (2011) S. R. Hudson, R. Dewar, M. Hole, and M. McGann, Non-axisymmetric, multi-region relaxed magnetohydrodynamic equilibrium solutions, Plasma Physics and Controlled Fusion 54, 014005 (2011).
  * Loizu and Bonfiglio (2023) J. Loizu and D. Bonfiglio, Nonlinear saturation of resistive tearing modes in a cylindrical tokamak with and without solving the dynamics, arXiv preprint arXiv:2305.16785 (2023).
  * Loizu  _et al._ (2020) J. Loizu, Y.-M. Huang, S. Hudson, A. Baillod, A. Kumar, and Z. Qu, Direct prediction of nonlinear tearing mode saturation using a variational principle, Physics of Plasmas 27 (2020).
  * Kumar  _et al._ (2021) A. Kumar, Z. Qu, M. J. Hole, A. M. Wright, J. Loizu, S. R. Hudson, A. Baillod, R. L. Dewar, and N. M. Ferraro, Computation of linear mhd instabilities with the multi-region relaxed mhd energy principle, Plasma Physics and Controlled Fusion 63, 045006 (2021).
  * Dennis  _et al._ (2014) G. R. Dennis, S. R. Hudson, and M. J. Hole, Modeling the single-helical axis state in the reversed-field pinch, IEEE Transactions on Plasma Science 42, 2514 (2014).
  * Liu  _et al._ (2024) K. Liu, G. Yu, Y. Huang, W. Mao, Y. Xie, X. Nie, H. Li, T. Lan, J. Xie, W. Ding, _et al._ , Effects of magnetic helicity on 3d equilibira and self-organized states in ktx reversed field pinch, arXiv preprint arXiv:2401.14604 (2024).
  * Dewar and Hudson (1998) R. Dewar and S. Hudson, Stellarator symmetry, Physica D: Nonlinear Phenomena 112, 275 (1998).
  * Zhu  _et al._ (2018a) C. Zhu, S. R. Hudson, Y. Song, and Y. Wan, Designing stellarator coils by a modified newton method using focus, Plasma Physics and Controlled Fusion 60, 065008 (2018a).
  * Cary and Hanson (1986) J. R. Cary and J. D. Hanson, Stochasticity reduction, The Physics of fluids 29, 2464 (1986).
  * Greene (1979) J. M. Greene, A method for determining a stochastic transition, Journal of Mathematical Physics 20, 1183 (1979).
  * Greene (1968) J. M. Greene, Two-dimensional measure-preserving mappings, Journal of Mathematical Physics 9, 760 (1968).
  * Smiet  _et al._ (2019) C. Smiet, G. Kramer, and S. Hudson, Mapping the sawtooth, Plasma Physics and Controlled Fusion 62, 025007 (2019).
  * Landreman  _et al._ (2021b) M. Landreman, B. Medasani, F. Wechsung, A. Giuliani, R. Jorge, and C. Zhu, Simsopt: A flexible framework for stellarator optimization, Journal of Open Source Software 6, 3525 (2021b).
  * Landreman (2017) M. Landreman, An improved current potential method for fast computation of stellarator coil shapes, Nuclear Fusion 57, 046003 (2017).
  * Zhu  _et al._ (2018b) C. Zhu, S. R. Hudson, S. A. Lazerson, Y. Song, and Y. Wan, Hessian matrix approach for determining error field sensitivity to coil deviations, Plasma Physics and Controlled Fusion 60, 054016 (2018b).
  * Zhu  _et al._ (2019) C. Zhu, D. A. Gates, S. R. Hudson, H. Liu, Y. Xu, A. Shimizu, and S. Okamura, Identification of important error fields in stellarators using the hessian matrix method, Nuclear Fusion 59, 126007 (2019).
  * Henneberg  _et al._ (2021b) S. A. Henneberg, P. Helander, and M. Drevlak, Representing the boundary of stellarator plasmas, Journal of Plasma Physics 87, 905870503 (2021b).
  * Lonigro and Zhu (2022) N. Lonigro and C. Zhu, Stellarator coil design using cubic splines for improved access on the outboard side, Nuclear Fusion 62, 066009 (2022).
  * Takahashi (2024) T. Takahashi, Designing gradient coils with the shape derivative and the closed b-spline curves, Magnetic Resonance Imaging (2024).
  * Landreman and Boozer (2016) M. Landreman and A. H. Boozer, Efficient magnetic fields for supporting toroidal plasmas, Physics of Plasmas 23 (2016).


# Path-Based Conditions for Local Network Identifiability – Full Version

Antoine Legat and Julien M. Hendrickx A. Legat and J. M. Hendrickx are with
ICTEAM Institute, UCLouvain, Belgium. Work supported by the “RevealFlight” ARC
at UCLouvain, and by the Incentive Grant for Scientific Research (MIS)
“Learning from Pairwise Data” of the<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This work focuses on the generic identifiability of dynamical networks with
partial excitation and measurement: a set of nodes are interconnected by
transfer functions according to a known topology, some nodes are excited, some
are measured, and only a part of the transfer functions are known. Our goal is
to determine whether the unknown transfer functions can be generically
recovered based on the input-output data collected from the excited and
measured nodes.

We propose a decoupled version of generic identifiability that is necessary
for generic local identifiability and might be equivalent as no counter-
example to sufficiency has been found yet in systematic trials. This new
notion can be interpreted as the generic identifiability of a larger network,
obtained by duplicating the graph, exciting one copy and measuring the other
copy. We establish a necessary condition for decoupled identifiability in
terms of vertex-disjoint paths in the larger graph, and a sufficient one.

##  I INTRODUCTION

This paper addresses the identifiability of dynamical networks in which node
signals are connected by causal linear time-invariant transfer functions, and
can be excited and/or measured. Such networks can be modeled as directed
graphs where each edge carries a transfer function, and known excitations and
measurements are applied at certain nodes.

We consider the identifiability of a network matrix $G(q)$, where the network
is made up of $n$ node signals $w(t)$, external excitation signals $r(t)$,
measured nodes $y(t)$ and noise $v_{1}(t),v_{2}(t)$ related to each other by:

$\displaystyle\begin{split}w(t)&=G(q)w(t)+Br(t)+v_{1}(t)\\\
y(t)&=Cw(t)+v_{2}(t),\end{split}$ (1)

where matrices $B$ and $C$ are binary selections defining respectively the
$n_{B}$ excited nodes and $n_{C}$ measured nodes, forming sets $\mathcal{B}$
and $\mathcal{C}$ respectively. Matrix $B$ is full column rank and each column
contains one 1 and $n-1$ zeros. Matrix $C$ is full row rank and each row
contains one 1 and $n-1$ zeros. The nonzero entries of the network matrix
$G(q)$ define the network topology: some of them are known and collected in
$G^{0}(q)$, and the others are the unknowns to identify, collected in
$G^{\Delta}(q)$, such that $G(q)=G^{0}(q)+G^{\Delta}(q)$. The known edges are
collected in set $E^{0}$, the unknown ones in $E^{\Delta}$, and
$E=\\{E^{0},E^{\Delta}\\}$ is the set of all edges.

We assume that _the input-output relations between the excitations $r$ and
measurements $y$ have been identified_, and that the network topology is
known. From this knowledge, we aim at recovering the unknown entries of
$G(q)$.

The model (1) has recently been the object of a significant research effort.
If the whole network is to be recovered, the notion of network identifiability
is used [1]. If one is interested in identifying a single module, topological
conditions are derived in [2, 3]. In this paper, we do not consider the impact
of noise $v_{1},v_{2}$, but studying the influence of rank-reduced or
correlated noise under certain assumptions yields less conservative
identifiability conditions [4, 5, 6, 7].

An approach dual to ours is to assume that network dynamics are known, and aim
at identifying the topology from input/output data. This problem is referred
to as topology identification, and is addressed in e.g. [8, 9].

It turns out that the identifiability of the network, i.e. the ability to
recover a module or the whole network from the input-output relation, is a
generic notion: Either _almost all_ transfer matrices corresponding to a given
network structure are identifiable, in which case the structure is called
_generically identifiable_ , or none of them are. A number of works study
generic identifiability when all nodes are excited or measured, i.e. when $B$
or $C=I$ [10, 11]. Considering the graph of the network, path-based conditions
on the allocation of measurements/excitations in the case of full excitation/
measurement are derived in [12] /[13]. Reformulating these conditions by means
of disjoint trees in the graph, [14] develops a scalable algorithm to allocate
excitations/measurements in case of full measurement/excitation. In case of
full measurement, [15] derives path-based conditions for the generic
identifiability of a subset of modules, under the presence of noise.
Abstractions of dynamic networks yield conditions on nodes to measure for
identifiability of a target module [16].

The conditions of [12] apply for generic identifiability, i.e. identifiability
of almost all transfer matrices corresponding to a given network structure.
[17] extends the path-based conditions under full excitation to determine the
identifiability _for all_ (nonzero) transfer matrices corresponding to a given
structure, and [18] provides conditions for the outgoing edges of a node, and
the whole network under the same conditions.

In all these works, the common assumption is that all nodes are either
excited, or measured. In [19], this assumption is relaxed and generic
identifiability with partial excitation and measurement is addressed for
particular network topologies. Partial excitation and measurement is also
addressed by [20], which derives conditions for generic identifiability in
terms of disconnecting sets, akin to what was done in full
excitation/measurement [13].

In the general case of arbitrary topology, partial excitation and measurement,
[21] introduces the notion of local identifiability, i.e. only on a
neighborhood of $G(q)$. Local identifiability is a generic property, necessary
for generic identifiability and no counterexample to sufficiency is known to
the authors, i.e. no network which is locally identifiable but not globally
identifiable. [21] derives algebraic necessary and sufficient conditions for
generic local identifiability for both the whole network, and a single module.

The algebraic conditions of [21] allow rapidly testing local identifiability
for any given network, but we wish to find a combinatorial characterization
for generic identifiability, that is expressed purely in terms of graph-
theoretical properties, akin to what was done in the full excitation case e.g.
in [12]. Such characterization would in particular pave the way for optimizing
the selection of nodes to be excited and measured, akin to the work in [13] in
the full measurement case. [20] already provides conditions in terms of
disconnecting sets, but we believe that the vertex-disjoint paths conditions
of [12] can be extended to local identifiability under partial excitation and
measurement.

In this paper, we derive path-based local identifiability conditions from the
results of [21], in the general case of arbitrary topology, partial excitation
and measurement. We extend the results of [21] when some transfer functions
are known _a priori_. A more general notion of local identifiability is
introduced: _decoupled_ identifiability, necessary for local identifiability.
Interestingly, no counterexample to sufficiency is known to the authors,
despite extensive testing (code available at [22]). Then, necessary and
sufficient path-based conditions for decoupled identifiability are derived.
These conditions are given in terms of connected paths and vertex-disjoint
paths, and extend what one had in the full excitation case e.g. in [12].

Assumptions: We consider the problem modeled in (1). Consistently with
previous works, we assume that the network is well-posed, that is
$(I-G(q))^{-1}$ is proper and stable, and we assume that
$CT(q)B=C(I-G(q))^{-1}B$ has been identified exactly. We do not suppose having
access to any information related to the effect of the noise signals
$v_{1},v_{2}$. The additional information that could be gathered from this
knowledge in our context is left for further works.

Consistently with [21], we consider in this paper a single frequency $z$, so
that all transfer functions are modeled simply by a complex value, and the
matrices $G$ and $T(G)=(I-G)^{-1}$ are complex matrices rather than matrices
of transfer functions. Conceptually, our generic results directly extend to
the transfer function case: if one can recover a $G_{ij}(z)$ at a given
frequency $z$ for almost all $G$ consistent with a network, then one can also
recover it at all other frequencies, and hence recover the transfer function.
We intend to remove this simplification or to formalize this intuitive
argument in a further version of this work. In the remainder of this document,
we omit $(q)$ to lighten notations. Also, the proofs of this paper are
collected in the Appendix.

##  II LOCAL IDENTIFIABILITY

We start from the definition of identifiability, see e.g. [12], which we
extend to the case where some transfer functions are known ($G^{0}$), and some
are not ($G^{\Delta}$), as in [2]. In the remainder of this paper, we denote
$T(G)=(I-G)^{-1}$ and we sometimes drop the $(G)$ when there is no ambiguity.

###### Definition 1

The transfer function $G_{ij}$ is _identifiable_ at $G$ from excitations
$\mathcal{B}$ and measurements $\mathcal{C}$ if, for all network matrix
$\tilde{G}$ with same zero and known entries as $G$, there holds

$\displaystyle C\,T(\tilde{G})\,B=C\,T(G)\,B\Rightarrow\tilde{G}_{ij}=G_{ij}.$
(2)

The network matrix is identifiable at $G$ if each unknown transfer function
$G_{ij}$ is identifiable at $G$, i.e. if the left-hand side of (2) implies
$\tilde{G}^{\Delta}=G^{\Delta}$.

We remind a notion of identifiability amenable to linear analysis: local
identifiability, which corresponds to identifiability provided that
$\tilde{G}$ is sufficiently close to $G$. Again, we extend the definition of
[21] to the more general case where some transfer functions are already known.

###### Definition 2

The transfer function $G_{ij}$ is _locally identifiable_ at $G$ from
excitations $\mathcal{B}$ and measurements $\mathcal{C}$ if there exists
$\epsilon>0$ such that for any $\tilde{G}$ with same zero and known entries as
$G$ satisfying $||\tilde{G}-G||<\epsilon$, there holds

$\displaystyle C\,T(\tilde{G})\,B=C\,T(G)\,B\Rightarrow\tilde{G}_{ij}=G_{ij}.$
(3)

The network matrix is locally identifiable at $G$ if each unknown transfer
function $G_{ij}$ is locally identifiable at $G$, i.e. if the left-hand side
of (3) implies $\tilde{G}^{\Delta}=G^{\Delta}$.

As stressed in [21], local identifiability is a necessary condition for
identifiability. It is _a priori_ a weaker notion, yet no example of network
locally identifiable but not globally identifiable is known to the authors.
Moreover, one can show that if a network is locally identifiable, then it can
be recovered up to a discrete ambiguity, i.e. the set of unknown $G^{\Delta}$
corresponding to a measured $CTB$ would be discrete.

Genericity: Given a network and sets $\mathcal{B},$ $\mathcal{C}$ of excited
and measured nodes, we say that an edge is _generically_ (locally)
identifiable if it is (locally) identifiable at all $G$ consistent with the
graph and known transfer functions, except possibly those lying on a lower-
dimensional set [23, 24] (i.e. a set of dimension lower than $|E^{\Delta}|$).
In the remainder of this paper, we say that a property is _generic_ if it
either holds (i) for _almost all_ variables, i.e for all variables except
possibly those lying on a lower-dimensional set, or (ii) for no variable.

For example, take a polynomial $p$. The _nonzeroness_ of $p(x)$ is a generic
property of $x$: either (i) $p(x)\neq 0$ for all $x$ except its roots, or (ii)
$p$ is the zero polynomial, which returns zero for all $x$.

A handy consequence of this definition is the following: showing that a
generic property holds for one variable $x$ implies that it holds for almost
all $x$.

###  II-A Algebraic condition

Proposition II.1 below, adapted from [21], states that local identifiability
is a generic property which can be checked by computing the rank of the matrix
$K$:

$\displaystyle K(G)\triangleq\left(B^{\top}T^{\top}(G)\otimes
C\,T(G)\right)I_{G^{\Delta}},$ (4)

where $\otimes$ denotes the Kronecker product and the matrix
$I_{G^{\Delta}}\in\\{0,1\\}^{n^{2}\times|E|}$ selects only the columns of the
preceding $n_{B}n_{C}\times n^{2}$ matrix corresponding to unknown edges:

$\displaystyle
I_{G^{\Delta}}\triangleq\left[\operatorname{vec}(G^{\Delta}(\mathbf{e}_{1}))\,\cdots\,\operatorname{vec}(G^{\Delta}(\mathbf{e}_{|E^{\Delta}|}))\right],$
(5)

where $G^{\Delta}(\mathbf{v})$ is the matrix with same zero structure as
$G^{\Delta}$, its nonzero entries collected in vector $\mathbf{v}$, and
$\mathbf{e}_{e}$ is the standard basis vector filled with zeros except 1 at
$e$-th entry.

###### Proposition II.1

Exactly one of the two following holds:111Observe that this implies: $K(G)$
has full rank for almost all $G$ _if and only if_ the network is generically
locally identifiable, but the way Proposition II.1 is stated is stronger. The
same phrasing remark holds for Proposition III.1 and Corollary IV.1.

  1. (i)

$K(G)$ has full rank for almost all $G$ and the network is generically locally
identifiable;

  2. (ii)

$K(G)$ has full rank for no $G$ and the network is never locally identifiable.

Moreover, $K(G)$ has full rank if and only if the following implication holds
for all $\Delta$ with same zero entries as $G^{\Delta}$:

$\displaystyle C\,T(G)\,\Delta\,T(G)\,B=0\Rightarrow\Delta=0.$ (6)

The proof is given in [21] when all transfer functions are unknown, and can
straightforwardly be adapted when some transfer functions are already known.

In this paper, we study network identifiability, i.e. the identifiability of
_all_ unknown transfer functions $G^{\Delta}$. Algebraic conditions for the
identifiability of a single transfer function $G_{i,j}$ are derived in [21],
but so far we were unable to interpret them in terms of paths in the graph. If
the network is not identifiable, we do not have path-conditions to find out
which transfer functions are problematic.

##  III DECOUPLED IDENTIFIABILITY

Consider condition (6) of Proposition II.1. We know that $C\,T(G)\,B$ is the
closed-loop transfer matrix of the network with excitations $\mathcal{B}$ and
measurements $\mathcal{C}$. It motivates the introduction of a larger network,
whose closed-loop transfer matrix is given by $C\,T(G)\,\Delta\,T(G)\,B$. It
is built by duplicating the network, exciting the left copy, measuring the
right copy and adding the unknown transfer functions in the middle (from left
to right), see Fig. 1.

Figure 1: The decoupled network $\hat{G}$ related to $G$. The network is
duplicated, excitations are applied on the left copy, measurements taken on
the right, and unknowns in the middle.

Since the network is duplicated, we may allow the left and right copies to
have different unknown transfer functions, while keeping the same topology and
known transfer functions. We will see that relaxing the equality of left and
right copies leads to a simpler analysis, hence we consider this more general
notion of identifiability in this paper.

First, we introduce the _decoupled network_ , whose closed-loop transfer
matrix is $C\,T(G)\,\Delta\,T(G^{\prime})\,B$, see Fig. 1.

###### Definition 3

Consider a network of $n$ nodes with excitation matrix $B$, measurement matrix
$C$ and network matrix $G=G^{0}+G^{\Delta}$, where $G^{0}$ collects the known
transfer functions and $G^{\Delta}$ collects the unknown transfer functions.
Its _decoupled network_ is composed of $2n$ nodes:
$\\{1,\dots,n,1^{\prime},\dots,n^{\prime}\\}$. Its network matrix is defined
by

$\displaystyle\hat{G}(G,G^{\prime})\triangleq\begin{bmatrix}G&G^{\Delta}\\\
0&G^{\prime}\end{bmatrix},$

where $G^{\prime}$ has the same zero and known entries as $G$, and their
unknown entries are given fixed parameters. Transfer matrices $G$ and
$G^{\prime}$ are then fully known, while $G^{\Delta}$ contains the unknown
transfer functions. Excitations are applied on the left subgraph
($G^{\prime}$), and measurements on the right one ($G$), i.e. its excitation
and measurement matrices are

$\displaystyle\hat{B}\triangleq\begin{bmatrix}0&0\\\ 0&B\end{bmatrix},$
$\displaystyle\hat{C}\triangleq\begin{bmatrix}C&0\\\ 0&0\end{bmatrix}.$

An example of decoupled network is given in Fig. 2.

$1$$2$$3$$G_{21}\qquad$$G_{23}$$G_{32}$ (a)

$1^{\prime}$$2^{\prime}$$3^{\prime}$$1$$2$$3$$G_{21}^{\Delta}$$G_{32}^{\Delta}$
(b)

Figure 2: (a): A basic network example: unknown transfer functions are in
dashed blue, arrows from the left denote excitations and the arrow to the
right is a measurement. (b): Its decoupled network: The graph is duplicated,
excitations are applied on left copy and measurement taken on the right one.
The unknown transfer functions are in the middle, in dashed blue: they link
the excited and measured subgraphs.

We are now ready to introduce decoupled identifiability, which we will prove
to be a generic property, necessary for generic (local) identifiability. Take
condition (6) of Proposition II.1, and allow the two $T(G)$s to have different
unknown transfer functions. The problem is no longer quadratic in $T(G)$, but
linear in both $T(G^{\prime})$ and $T(G)$:

###### Definition 4

A network is _decoupled-identifiable_ at $(G,$ $G^{\prime})$, with $G$ and
$G^{\prime}$ sharing the same zero and known entries, if for all $\Delta$ with
same zero entries as $G^{\Delta}$, there holds:

$\displaystyle C\,T(G)\,\Delta\,T(G^{\prime})\,B=0\Rightarrow\Delta=0.$ (7)

Similarly to local identifiability, decoupled identifiability is a generic
property: either it holds for almost all $(G,G^{\prime})$, or for no
$(G,G^{\prime})$. It is proved in Proposition III.1 below, which relies on the
rank of the following matrix, defined analogously to $K$ in (4), but with
$G^{\prime}$ in the excitation part:

$\displaystyle\hat{K}(G,G^{\prime})\triangleq\left(B^{\top}T(G^{\prime})^{\top}\otimes
C\,T(G)\right)I_{G^{\Delta}},$

where $\otimes$ is the Kronecker product and $I_{G^{\Delta}}$ is defined in
(5).

###### Proposition III.1

Exactly one of the two following holds:

  1. (i)

$\hat{K}(G,G^{\prime})$ has full rank for almost all $(G,G^{\prime})$ and the
network is generically decoupled-identifiable;

  2. (ii)

$\hat{K}(G,G^{\prime})$ has full rank for no $(G,G^{\prime})$ and the network
is never decoupled-identifiable.

A first important proposition is that generic decoupled identifiability is
necessary for generic local identifiability:

###### Proposition III.2

If a network is generically locally identifiable, then it is generically
decoupled-identifiable.

Interestingly, no counterexample to sufficiency is known to the authors,
despite numerous tests: we have randomly generated $10^{6}$ networks with up
to $100$ nodes, and checked the generic local identifiability and generic
decoupled identifiability of each network by computing
$\operatorname{rank}\,K$ and $\operatorname{rank}\,\hat{K}$: for every
network, both matched. In other words, experiments seem to show that generic
decoupled identifiability is equivalent to generic local identifiability. Code
available at [22].

Moreover, since generic local identifiability is necessary for generic
identifiability [21], so is generic decoupled identifiability. Hence, the
necessary conditions derived for generic decoupled identifiability in this
paper also hold for generic identifiability, and for identifiability. Besides,
no example of network that is generically locally identifiable, but not
generically identifiable is known to the authors.

Definition 3 introduces the decoupled network, Definition 4 presents the
notion of decoupled identifiability. The proposition below unifies those two
notions.

###### Proposition III.3

The network $G$ is generically decoupled-identifiable if and only if its
decoupled network $\hat{G}(G,G^{\prime})$ is generically identifiable for
almost all $(G,G^{\prime})$.

The following example illustrates propositions III.1 and III.3.

###### Example III.1

For the network of Fig. 2 (a), one can check that $\hat{K}$ has generic rank
2. Since there are 2 unknown transfer functions, Proposition III.1 ensures
that we have generic decoupled identifiability. It can be interpreted on the
decoupled network, depicted on Fig. 2 (b), as the generic identifiability of
the unknown transfer functions $G^{\Delta}_{21}$ & $G^{\Delta}_{32}$.

##  IV PATH-BASED CONDITIONS
FOR DECOUPLED IDENTIFIABILITY

For ease of presentation, we consider the case where there are exactly
$n_{B}n_{C}$ unknown edges, i.e. as many as the number of (excitation $b$,
measurement $c$) pairs. If there are more unknown edges, then the network is
not identifiable since there are more unknowns than (in, out) data. The
situation with less unknown edges than $n_{B}n_{C}$ will be addressed in a
more complete version of this work. In this section, we drop the
$G,G^{\prime}$ arguments and denote $T(G^{\prime})$ by $T^{\prime}$. Also, we
refer to unknown transfer functions $G^{\Delta}_{i,j}$ as unknown edges
$\alpha$.

Since $|E^{\Delta}|=n_{B}n_{C}$, $\hat{K}$ is a square matrix, hence
$\operatorname{rank}\,\hat{K}=|E^{\Delta}|$ is equivalent to $\det\hat{K}\neq
0$. Proposition III.1 can then be rewritten in terms of determinant:

###### Corollary IV.1

Exactly one of the two following holds:

  1. (i)

$\det\hat{K}\neq 0$ generically and $G$ is generically decoupled-identifiable;

  2. (ii)

$\det\hat{K}$ is always zero and $G$ is never decoupled-identifiable.

In order to interpret $\det\hat{K}$, we develop the system of (7):

$\displaystyle\sum_{\alpha\in
E^{\Delta}}T^{\prime}_{\alpha,b}T_{c,\alpha}\,\Delta_{\alpha}=0\qquad\forall\
b\in\mathcal{B},\ c\in\mathcal{C}.$ (8)

Each equation of (8) represents a pair (excitation $b$, measurement $c$), and
corresponds to a row of $\hat{K}$. Each column of $\hat{K}$ matches an unknown
edge $\alpha$, and the entry of $\hat{K}$ corresponding to $(b,c)$ and
$\alpha$ is given by $T^{\prime}_{\alpha,b}T_{c,\alpha}$.222By abuse of
notation, $T^{\prime}_{\alpha,b}$ stands for the transfer function between
excitation $b$ and _start node of_ edge $\alpha$, and $T_{c,\alpha}$ denotes
the one between _end node of_ edge $\alpha$ and measurement $c$.

Besides, the determinant is expressed as the sum over all possible row-column
permutations by the Leibniz formula:

$\displaystyle\det\hat{K}=\sum_{\sigma\in
S}\operatorname{sgn}(\sigma)\prod_{\alpha\in
E^{\Delta}}T^{\prime}_{\alpha,\sigma_{B}(\alpha)}T_{\sigma_{C}(\alpha),\alpha},$
(9)

where each row-column permutation corresponds to a bijective assignation
$\sigma:E^{\Delta}\rightarrow\mathcal{B}\times\mathcal{C}$, composed of (Fig.
3 (a)):

  * •

an excitation assignation $\sigma_{B}:E^{\Delta}\rightarrow\mathcal{B}$ in
which $n_{C}$ unknown edges $\alpha$ are assigned to each excitation $b$,

  * •

a measurement assignation $\sigma_{C}:E^{\Delta}\rightarrow\mathcal{C}$ in
which $n_{B}$ unknown edges $\alpha$ are assigned to each measurement $c$.

Each bijective assignation $\sigma$ is composed of a $\sigma_{B}$ and a
$\sigma_{C}$, but not every $(\sigma_{B},\sigma_{C})$ pair gives a bijective
assignation $\sigma$, e.g. see Fig. 3 (b). We say that $\sigma_{B}$ and
$\sigma_{C}$ are compatible if they form a bijection. $S$ is the set composed
of all bijective assignations $\sigma$, and $\operatorname{sgn}(\sigma)$
equals $+1$ if the number of transpositions in assignation $\sigma$ is even,
and $-1$ otherwise. A transposition is the swap of two elements. Each $\sigma$
is obtained by combining a certain number of transpositions (although such
decomposition is not unique, the number of transpositions always has same
parity).

(a)

(b)

Figure 3: Each bijective assignation $\sigma$ is composed of a $\sigma_{B}$
and a $\sigma_{C}$: see (a), but not every $(\sigma_{B},\sigma_{C})$ pair
gives a bijective assignation $\sigma$. In (b), the orange and red edges are
assigned to the same (excitation, measurement) pair: $\sigma_{B}$ and
$\sigma_{C}$ are not compatible.

The graph-theoretical conditions of this section rely on _vertex-disjoint
paths_ : we say that a group of paths are vertex-disjoint if no two paths of
this group contain the same vertex [12], see Fig. 4. The following lemma links
vertex-disjoint paths in the graph with the generic rank of transfer matrices:

$\mathcal{I}$$\mathcal{J}$ Figure 4: One can route up to three vertex-disjoint
paths from $\mathcal{J}$ to $\mathcal{I}$ (in purple), i.e.
$\beta_{\mathcal{J}\rightarrow\mathcal{I}}=3$.

###### Lemma IV.1

[25, 12] Let $T_{\mathcal{I},\mathcal{J}}$ be the transfer matrix from nodes
of set $\mathcal{J}$ to nodes of set $\mathcal{I}$. We have

$\displaystyle\operatorname{generic\
rank}\,T_{\mathcal{I},\mathcal{J}}=\beta_{\mathcal{J}\rightarrow\mathcal{I}},$
(10)

where the generic rank of $T_{\mathcal{I},\mathcal{J}}$ is the rank for
_almost all_ $G=G^{0}+G^{\Delta}$, and
$\beta_{\mathcal{J}\rightarrow\mathcal{I}}$ stands for the maximum number of
vertex-disjoint paths that can be routed from $\mathcal{J}$ to $\mathcal{I}$.

Note that the rank is taken for almost all unknown transfer functions
$G^{\Delta}$ _and almost all_ known transfer functions $G^{0}$, otherwise (10)
would not hold for problematic values of $G^{0}$, as shown in Example 1 of
[15].

Combining Leibniz formula (9) with Lemma IV.1 yields the following basic
proposition. An assignation $\sigma$ is _connected_ if for each unknown edge
$\alpha$ there is a path from its assigned excitation $\sigma_{B}(\alpha)$ to
its assigned measurement $\sigma_{C}(\alpha)$, in which the unknown edge
$\alpha$ is included.

###### Proposition IV.1

A necessary and a sufficient condition:

  * •

If a network is generically decoupled-identifiable, then there is _at least_
one connected bijective assignation $\sigma$.

  * •

If there is _only_ one connected bijective assignation $\sigma$, then the
network is generically decoupled-identifiable.

Proposition IV.1 relies on (9), but this expression can be further developed,
and algebraic manipulations allow to derive a stronger condition in terms of
vertex-disjoint paths. The following lemma is the main building brick of this
stronger condition. In the following lemma, an assignation
$\sigma_{B}:E^{\Delta}\rightarrow\mathcal{B}$ (respectively
$\sigma_{C}:E^{\Delta}\rightarrow\mathcal{C})$ is _connected_ if there is a
path between each unknown edge $\alpha$ and its assigned excitation
$\sigma_{B}(\alpha)$ (resp. measurement $\sigma_{C}(\alpha)$).

###### Lemma IV.2

If a network is generically decoupled-identifiable, then there is at least one
assignation $\sigma_{B}$ (resp. $\sigma_{C}$) such that:

  1. (a)

$n_{C}$ (resp. $n_{B}$) unknown edges $\alpha$ are assigned to each excitation
$b$ (resp. measurement $c$)

  2. (b)

$\sigma_{B}$ (resp. $\sigma_{C}$) is connected

  3. (c)

for each excitation $b$ (resp. measurement $c$), there are $n_{C}$ (resp.
$n_{B}$) vertex-disjoint paths between the edges assigned to $b$ (resp. $c$)
and measurements $\mathcal{C}$ (resp. excitations $\mathcal{B}$).

If there is only one such assignation, then this condition is also sufficient
for generic decoupled identifiability.

Combining Lemma IV.2 for $\sigma_{B}$ and $\sigma_{C}$ yields the main result
of this paper.

###### Theorem IV.1

If a network is generically decoupled-identifiable, then there is at least one
assignation $\sigma$ such that:

  1. (a)

$n_{C}$ unknown edges $\alpha$ are assigned to each excitation $b$

  2. (b)

$n_{B}$ unknown edges $\alpha$ are assigned to each measure $c$

  3. (c)

$\sigma$ is connected

  4. (d)

for each excitation $b$, there are $n_{C}$ vertex-disjoint paths between the
edges assigned to $b$ and measurements $\mathcal{C}$.

  5. (e)

for each measurement $c$, there are $n_{B}$ vertex-disjoint paths between
edges assigned to $c$ and excitations $\mathcal{B}$.

If there is only one such assignation, then this condition is also sufficient
for generic decoupled identifiability.

##  V DISCUSSION

First, we highlight some subtle difference between the conditions of Theorem
IV.1 and those of Proposition IV.1:

  1. (i)

In contrast to Proposition IV.1, the assignation $\sigma$ of Theorem IV.1 is
not necessarily bijective: two edges assigned to the same excitation can be
assigned to the same measurement.

  2. (ii)

The assignations of condition a do not necessarily match the vertex-disjoint
paths of condition e.

  3. (iii)

The assignations of condition b do not necessarily match the vertex-disjoint
paths of condition d.

We believe that there could be a stronger version of Theorem IV.1, in which
$\sigma$ is bijective, as in Proposition IV.1. It would unify Proposition IV.1
with Theorem IV.1, and extend the vertex-disjoint path conditions of [12].

Besides, Proposition IV.1 and Theorem IV.1 provide necessary conditions for
generic decoupled identifiability and sufficient ones. As shown in Proposition
III.2, generic decoupled identifiability is necessary for generic local
identifiability, which was itself shown to be necessary for generic
identifiability in [21]. Therefore, the necessary conditions derived in
Section IV apply to (generic) (local) identifiability.

Regarding sufficiency, we remind that in the systematic tests we have
conducted (code at [22]), we have found no network which is generically
decoupled-identifiable but not generically locally identifiable. No
counterexample to the sufficiency of generic local identifiability for generic
identifiability has been found either yet [21], but we could not conduct
numerical tests for this one since we lack of a general criterion to check
generic identifiability. Consequently, one might hope that the sufficient
conditions of Section IV apply to generic (local) identifiability, but this
remains an open question.

##  VI CONCLUSION

This work was motivated by one main open question: determining path-based
conditions for generic local identifiability in networked systems.

The decoupled version of generic identifiability we have introduced allowed to
look at generic identifiability from a new angle, on a larger graph which
decouples excitations and measurements. In particular, we have derived
necessary conditions in terms of vertex-disjoint paths in the larger graph and
sufficient ones. The necessary conditions extend to generic (local)
identifiability, but whether the sufficient conditions extend as well remains
an open question.

We believe that the two identifiability conditions of this paper could be
merged into a unifying stronger condition. It would extend results of full
excitation or measurement.

A further open question would be to obtain graph theoretical conditions for
local identifiability of a subset of edges when not all edges can be
recovered, a question for which [21] gives an algebraic necessary and
sufficient condition.

## References

  * [1] H. H. Weerts, A. G. Dankers, and P. M. Van den Hof, “Identifiability in dynamic network identification,” IFAC-PapersOnLine, 2015.
  * [2] H. Weerts, P. M. Van den Hof, and A. Dankers, “Single module identifiability in linear dynamic networks,” in 2018 IEEE Conference on Decision and Control (CDC), pp. 4725–4730, IEEE, 2018.
  * [3] M. Gevers, A. S. Bazanella, and G. V. da Silva, “A practical method for the consistent identification of a module in a dynamical network,” IFAC-PapersOnLine, vol. 51, no. 15, pp. 862–867, 2018.
  * [4] H. H. Weerts, P. M. Van den Hof, and A. G. Dankers, “Prediction error identification of linear dynamic networks with rank-reduced noise,” Automatica, vol. 98, pp. 256–268, 2018.
  * [5] M. Gevers, A. S. Bazanella, and G. A. Pimentel, “Identifiability of dynamical networks with singular noise spectra,” IEEE Transactions on Automatic Control, vol. 64, no. 6, pp. 2473–2479, 2018.
  * [6] P. M. Van den Hof, K. R. Ramaswamy, A. G. Dankers, and G. Bottegal, “Local module identification in dynamic networks with correlated noise: the full input case,” in 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 5494–5499, IEEE, 2019.
  * [7] K. R. Ramaswamy and P. M. Vandenhof, “A local direct method for module identification in dynamic networks with correlated noise,” IEEE Transactions on Automatic Control, 2020.
  * [8] H. J. van Waarde, P. Tesi, and M. K. Camlibel, “Topology identification of heterogeneous networks of linear systems,” in 58th Conference on Decision and Control (CDC), pp. 5513–5518, IEEE, 2019.
  * [9] H. J. van Waarde, P. Tesi, and M. K. Camlibel, “Topology reconstruction of dynamical networks via constrained lyapunov equations,” IEEE Transactions on Automatic Control, vol. 64, pp. 4300–4306, 2019.
  * [10] A. S. Bazanella, M. Gevers, J. M. Hendrickx, and A. Parraga, “Identifiability of dynamical networks: which nodes need be measured?,” in IEEE 56th Annual Conference on Decision and Control (CDC), 2017.
  * [11] H. H. Weerts, P. M. Van den Hof, and A. G. Dankers, “Identifiability of linear dynamic networks,” Automatica, vol. 89, pp. 247–258, 2018.
  * [12] J. M. Hendrickx, M. Gevers, and A. S. Bazanella, “Identifiability of dynamical networks with partial node measurements,” IEEE Transactions on Automatic Control, vol. 64, no. 6, pp. 2240–2253, 2018.
  * [13] S. Shi, X. Cheng, and P. M. Van den Hof, “Excitation allocation for generic identifiability of a single module in dynamic networks: A graphic approach,” in Proc. 21st IFAC World Congress, 2020.
  * [14] X. Cheng, S. Shi, and P. M. Van den Hof, “Allocation of excitation signals for generic identifiability of dynamic networks,” in Proceedings of the IEEE Conference on Decision and Control, 2019.
  * [15] S. Shi, X. Cheng, and P. M. Van den Hof, “Generic identifiability of subnetworks in a linear dynamic network: the full measurement case,” arXiv preprint arXiv:2008.01495, 2020.
  * [16] H. H. Weerts, J. Linder, M. Enqvist, and P. M. Van den Hof, “Abstractions of linear dynamic networks for input selection in local module identification,” Automatica, vol. 117, p. 108975, 2020.
  * [17] H. J. van Waarde, P. Tesi, and M. K. Camlibel, “Topological conditions for identifiability of dynamical networks with partial node measurements,” IFAC-PapersOnLine, vol. 51, no. 23, pp. 319–324, 2018.
  * [18] H. J. Van Waarde, P. Tesi, and M. K. Camlibel, “Necessary and sufficient topological conditions for identifiability of dynamical networks,” IEEE Transactions on Automatic Control, 2019.
  * [19] A. S. Bazanella, M. Gevers, and J. M. Hendrickx, “Network identification with partial excitation and measurement,” in 58th Conference on Decision and Control (CDC), pp. 5500–5506, IEEE, 2019.
  * [20] S. Shi, X. Cheng, and P. M. Van den Hof, “Single module identifiability in linear dynamic networks with partial excitation and measurement,” arXiv preprint arXiv:2012.11414, 2020.
  * [21] A. Legat and J. M. Hendrickx, “Local network identifiability with partial excitation and measurement,” in 2020 59th IEEE Conference on Decision and Control (CDC), pp. 4342–4347, IEEE, 2020.
  * [22] A. Legat and J. M. Hendrickx, Identifiability test.  https://github.com/alegat/identifiable.
  * [23] E. J. Davison, “Connectability and structural controllability of composite systems,” Automatica, vol. 13, no. 2, pp. 109–123, 1977.
  * [24] J.-M. Dion, C. Commault, and J. Van der Woude, “Generic properties and control of linear structured systems: a survey,” Automatica, 2003.
  * [25] J. Van der Woude, “A graph-theoretic characterization for the rank of the transfer matrix of a structured system,” Mathematics of Control, Signals and Systems, vol. 4, no. 1, pp. 33–40, 1991.
  * [26] J. P. D’Angelo, An introduction to complex analysis and geometry, vol. 12.  American Mathematical Soc., 2010.

##  Appendix A APPENDIX

###  A-A Proofs of Section III.

Proof of Proposition III.1: We first prove that (7) holds if and only if
$\hat{K}(G,G^{\prime})$ has full rank: vectorizing (7) gives

$\displaystyle\left(B^{\top}T(G^{\prime})^{\top}\otimes
C\,T(G)\right)\operatorname{vec}(\Delta)=0\Rightarrow\operatorname{vec}(\Delta)=0,$

where $\operatorname{vec}$ stacks the columns of its
$|\mathcal{B}|\times|\mathcal{C}|$ matrix argument into a
$|\mathcal{B}|\cdot|\mathcal{C}|$ vector. Then, observing that
$\operatorname{vec}\left(\Delta\right)=I_{G^{\Delta}}\,\delta$ yields
$\hat{K}\delta=0\Rightarrow\delta=0\quad\forall\,\delta\in\mathbb{C}^{|E^{\Delta}|}$,
which gives that $\hat{K}(G,G^{\prime})$ has full column rank $|E^{\Delta}|$.

Then, we prove that the full-rankness of $\hat{K}(G,G^{\prime})$ is a generic
property. The following argument is inspired from the proof of Lemma B.1 in
[21]: the set on which the rank of $\hat{K}(G,G^{\prime})$ is not full can be
expressed as the intersection of zero sets of determinants of submatrices
[25], which are analytic in their entries $(G,G^{\prime})$. Since nonconstant
analytic functions vanish only on a lower-dimensional set [26], the set on
which $\hat{K}(G,G^{\prime})$ does not reach its maximal rank has lower
dimension. This extends the notion of generic rank introduced in [23].

Proof of Proposition III.2: Consider a generically locally identifiable
network. Take $G$, at which the network is locally identifiable. Equation (6)
holds, and it is decoupled-identifiable at $(G,G^{\prime}=G)$ by definition.
Proposition III.1 asserts that decoupled identifiability is a generic
property, thus the fact it holds for one variable $(G,G^{\prime}=G)$ implies
it holds for almost all $(G,G^{\prime})$.

Proof of Proposition III.3: From definitions 1 and 3, the decoupled network
$\hat{G}(G,G^{\prime})$ is generically identifiable if, for all
$\tilde{G}^{\Delta}$ with same zero entries as $G^{\Delta}$, there holds

$\displaystyle\hat{C}\,(I-\hat{\tilde{G}})^{-1}\,\hat{B}=\hat{C}\,(I-\hat{G})^{-1}\,\hat{B}\Rightarrow\tilde{G}^{\Delta}=G^{\Delta}$
(11)

for almost all $G^{\Delta}$ with same zeros, where
$\hat{\tilde{G}}=\begin{bmatrix}G&\tilde{G}^{\Delta}\\\
0&G^{\prime}\end{bmatrix}.$ Developing (11) yields

$\displaystyle
C\,T(G)\,\tilde{G}^{\Delta}\,T(G^{\prime})\,B=C\,T(G)\,G^{\Delta}\,T(G^{\prime})\,B\Rightarrow\tilde{G}^{\Delta}=G^{\Delta},$

and bringing out common terms gives

$\displaystyle
C\,T(G)\,\underbrace{(\tilde{G}^{\Delta}-G^{\Delta})}_{\triangleq\Delta}\,T(G^{\prime})\,B=0\Rightarrow\underbrace{\tilde{G}^{\Delta}-G^{\Delta}}_{=\Delta}=0,$

which matches Definition 4, and requiring this to hold for almost all
$(G,G^{\prime})$ completes the proof.

###  A-B Proofs of Section IV

Proof of Proposition IV.1: Consider a generically decoupled-identifiable
network. Corollary IV.1 yields that $\det\hat{K}\neq 0$ generically. Besides,
the Leibniz formula for the determinant gives (9):

$\displaystyle\det\hat{K}=\sum_{\sigma\in
S}\operatorname{sgn}(\sigma)\underbrace{\prod_{\alpha\in
E^{\Delta}}T^{\prime}_{\alpha,\sigma_{B}(\alpha)}T_{\sigma_{C}(\alpha),\alpha}}_{\triangleq\tau(\sigma)}.$

Since $\det\hat{K}\neq 0$ generically, there is at least one assignation
$\sigma$ with a nonzero term $\tau(\sigma)$. For this assignation, each factor
of the product must be generically nonzero, i.e.
$T^{\prime}_{\alpha,\sigma_{B}(\alpha)}\neq 0$ and
$T_{\sigma_{C}(\alpha),\alpha}\neq 0$ for all unknown edge $\alpha$. Lemma
IV.1 for $\mathcal{J}=\\{j\\},\mathcal{I}=\\{i\\}$ ensures that $T_{i,j}$ is
generically nonzero if and only if there is a path from $j$ to $i$. Hence, we
can construct a path going from $\sigma_{B}(\alpha)$ to $\sigma_{C}(\alpha)$,
including $\alpha$.

If more than one assignation $\sigma$ have a nonzero term $\tau(\sigma)$, the
condition is not sufficient since terms could cancel each other. But if there
is only one such assignation, then the condition is also sufficient.

The proof of Lemma IV.2 relies on Lemma A-B.1, which handles the signature
decomposition of the assignation.

###### Lemma A-B.1

Let $\sigma$ be a bijective assignation comprised of compatible $\sigma_{B}$
and $\sigma_{C}$. Its signature decomposes into

$\displaystyle\operatorname{sgn}(\sigma)=\operatorname{sgn}(\sigma_{B})\
\operatorname{sgn}(\sigma_{C}),$

where
$\operatorname{sgn}(\sigma_{B})\triangleq\prod_{c}\operatorname{sgn}(\sigma_{B,c}),\
\operatorname{sgn}(\sigma_{C})\triangleq\prod_{b}\operatorname{sgn}(\sigma_{C,b}),$
$\sigma_{C,b}$ is the sub-assignation of $\sigma_{C}$ for the edges $\alpha$
assigned to the same excitation $b$ in $\sigma_{B}$ (see Fig. 5 and equation
(13)), and $\sigma_{B,c}$ is defined analogously.

Proof: An assignation $\sigma$ is comprised of a sequence of transpositions
that swap two edges: denote them $\alpha_{1,1}$ and $\alpha_{2,2}$, initially
assigned to $(b_{1},c_{1})$ and $(b_{2},c_{2})$ respectively. Since $\sigma$
is bijective, either (i) $b_{1}=b_{2},c_{1}\neq c_{2}$, (ii) $b_{1}\neq
b_{2},c_{1}=c_{2}$, or (iii) $b_{1}\neq b_{2},c_{1}\neq c_{2}$. Besides, every
transposition of case (iii) can be decomposed into transpositions of type (i)
and (ii), as shown in Table I.

$(b_{1},c_{1})$ | $\boxed{\alpha_{1,1}}$ | $\xrightarrow[]{\text{(i)}}$ | $\alpha_{1,2}$ | $\xrightarrow[]{\text{(ii)}}$ | $\boxed{\alpha_{1,2}}$ | $\xrightarrow[]{\text{(i)}}$ | $\alpha_{2,2}$
---|---|---|---|---|---|---|---
| $\updownarrow$ |  |  |  | $\updownarrow$ |  |
$(b_{1},c_{2})$ | $\boxed{\alpha_{1,2}}$ |  | $\boxed{\alpha_{1,1}}$ |  | $\boxed{\alpha_{2,2}}$ |  | $\alpha_{1,2}$
$(b_{2},c_{1})$ | $\alpha_{2,1}$ |  | $\alpha_{2,1}\ \updownarrow$ |  | $\alpha_{2,1}$ |  | $\alpha_{2,1}$
$(b_{2},c_{2})$ | $\alpha_{2,2}$ |  | $\boxed{\alpha_{2,2}}$ |  | $\alpha_{1,1}$ |  | $\alpha_{1,1}$
TABLE I: Decomposition of a class-(iii) transposition into class-(i) and -(ii)
transpositions. The boxes and arrows indicate the swap to the next column.

Hence, any assignation $\sigma$ can be decomposed into a sequence of class-(i)
and -(ii) transpositions only. Take one such sequence, and denote its number
of transpositions $\pi(\sigma)$. Within this sequence, consider the
transpositions of class (i) with same fixed $b$. Combining these
transpositions gives the sub-assignation $\sigma_{C,b}$, defined in equation
(13) and illustrated in Fig. 5. Hence, the number of such transpositions gives
the parity of $\sigma_{C,b}$, and we denote it by $\pi(\sigma_{C,b})$. The
same holds for all $b$, and for the class-(ii) transpositions, for each $c$.
Enumerating all $\pi(\sigma)$ transpositions of the considered sequence allows
to rewrite it as a sum:

$\displaystyle\pi(\sigma)=\sum_{c}\pi(\sigma_{B,c})+\sum_{b}\pi(\sigma_{C,b}).$

Since $\operatorname{sgn}(\sigma)\triangleq(-1)^{\pi(\sigma)}$, this yields

$\displaystyle\qquad\quad\
\operatorname{sgn}(\sigma)=\prod_{c}(-1)^{\pi(\sigma_{B,c})}\
\prod_{b}(-1)^{\pi(\sigma_{C,b})}.\qquad\quad\ \rule{5.69054pt}{5.69054pt}$

Proof of Lemma IV.2: This proof is derived for $\sigma_{B}$, but an analogous
argument can be made for $\sigma_{C}$. Consider a generically decoupled-
identifiable network. Corollary IV.1 yields that $\det\hat{K}\neq 0$
generically. Besides, the Leibniz formula for the determinant gives (9):

$\displaystyle\det\hat{K}$ $\displaystyle=\sum_{\sigma\in
S}\operatorname{sgn}(\sigma)\prod_{\alpha\in
E^{\Delta}}T^{\prime}_{\alpha,\sigma_{B}(\alpha)}T_{\sigma_{C}(\alpha),\alpha}$
$\displaystyle\,\downarrow\quad\text{grouping $\sigma$ with same $\sigma_{B}$,
Lemma {\ref{lemma:signature}} for }\operatorname{sgn}$
$\displaystyle=\sum_{\sigma_{B}}\operatorname{sgn}(\sigma_{B})\sum_{\sigma_{C}\in
S_{C}^{\sigma_{B}}}\operatorname{sgn}(\sigma_{C})\prod_{\alpha}T^{\prime}_{\alpha,\sigma_{B}(\alpha)}T_{\sigma_{C}(\alpha)\alpha},$

where $S_{C}^{\sigma_{B}}$ is the set of all $\sigma_{C}$ compatible with
$\sigma_{B}$, i.e. such that $(\sigma_{B},\sigma_{C})$ is bijective: if two
edges are assigned to the same $b$ by $\sigma_{B}$, their $c$ assigned by
$\sigma_{C}$ must be different.

Bringing out common $T^{\prime}$ terms yields $\det\hat{K}=$

$\displaystyle\sum_{\sigma_{B}}\operatorname{sgn}(\sigma_{B})\prod_{\alpha}T^{\prime}_{\alpha,\sigma_{B}(\alpha)}\underbrace{\sum_{\sigma_{C}\in
S_{C}^{\sigma_{B}}}\operatorname{sgn}(\sigma_{C})\prod_{\alpha}T_{\sigma_{C}(\alpha),\alpha}}_{\triangleq\gamma_{\sigma_{B}}}$

We develop further $\gamma_{\sigma_{B}}$, the factor accounting for the
assignations $\sigma_{C}$ compatible with $\sigma_{B}$. We decompose the
product on all unknown edges $\alpha$ into subgroups of $\alpha$ assigned to
the same excitation $b$ in $\sigma_{B}$ (subgroups denoted as
$\sigma_{B}^{-1}(b)$):

$\displaystyle\gamma_{\sigma_{B}}=\sum_{\sigma_{C}\in
S_{C}^{\sigma_{B}}}\operatorname{sgn}(\sigma_{C})\prod_{b}\prod_{\alpha\in\sigma_{B}^{-1}(b)}T_{\sigma_{C}(\alpha),\alpha}$
(12)

The assignation $\sigma_{C}$ can be decomposed into $n_{B}$ sub-assignations
for the edges $\alpha$ assigned to the same excitation $b$ in $\sigma_{B}$
(i.e. the $\alpha$ belonging to $\sigma_{B}^{-1}(b)$, see Fig. 5):

Figure 5: The assignation $\sigma_{C}$ can be decomposed into $n_{B}$ sub-
assignations for the edges $\alpha$ assigned to the same excitation $b$ in
$\sigma_{B}$.

$\displaystyle\sigma_{C,b}:\sigma_{B}^{-1}(b)\rightarrow\mathcal{C}:\alpha\rightarrow\sigma_{C,b}(\alpha)$
(13)

Here is the key step of the proof: since those sub-assignations $\sigma_{C,b}$
do not depend on each other, (12) can be rewritten by regrouping in the same
factor the edges $\alpha$ assigned to the same excitation $b$ in the
assignation $\sigma_{B}$:

$\displaystyle\gamma_{\sigma_{B}}=\prod_{b}\sum_{\sigma_{C,b}\in
S_{C,b}^{\sigma_{B}}}\operatorname{sgn}(\sigma_{C,b})\prod_{\alpha\in\sigma_{B}^{-1}(b)}T_{\sigma_{C,b}(\alpha),\alpha},$
(14)

where $S^{\sigma_{B}}_{C,b}$ denotes the set of sub-assignations
$\sigma_{C,b}$ compatible with $\sigma_{B}$.

Furthermore, each sub-assignation $\sigma_{C,b}$ compatible with $\sigma_{B}$
is bijective: since each pair $(b,c)$ must be assigned to only one edge
$\alpha$ for compatibility, each $\alpha$ of $\sigma_{B}^{-1}(b)$ must be
assigned to a different measurement, i.e. $\sigma_{C,b}$ must be injective.
Besides, $\sigma_{B}^{-1}(b)$ is composed of $n_{C}$ edges, so $\sigma_{C,b}$
is bijective.

All compatible $\sigma_{C,b}$ must be taken in order to cover all assignations
$\sigma$. Thus $S^{\sigma_{B}}_{C,b}$ contains all possible bijections
$\sigma_{C,b}$. Hence, each factor of (14) is a determinant:

$\displaystyle\gamma_{\sigma_{B}}=\prod_{b}\det
T_{\mathcal{C},\sigma_{B}^{-1}(b)}.$

Putting it all together, we have

$\displaystyle\det\hat{K}=\sum_{\sigma_{B}}\operatorname{sgn}(\sigma_{B})\underbrace{\prod_{\alpha}T^{\prime}_{\alpha,\sigma_{B}(\alpha)}\prod_{b}\det
T_{\mathcal{C},\sigma_{B}^{-1}(b)}}_{\triangleq\tau_{B}(\sigma_{B})}$ (15)

Since $\det\hat{K}\neq 0$ generically, there is at least one assignation
$\sigma_{B}$ with a nonzero term $\tau_{B}(\sigma_{B})$. For this assignation,
each factor of the product must be nonzero, i.e. (i)
$T^{\prime}_{\alpha,\sigma_{B}(\alpha)}\neq 0$ for all unknown edge $\alpha$
and (ii) $\det T_{\mathcal{C},\sigma_{B}^{-1}(b)}\neq 0$ for all excitation
$b$.

Lemma IV.1 for single nodes $\mathcal{J}=\\{j\\},\mathcal{I}=\\{i\\}$ ensures
that $T_{i,j}$ is generically nonzero if and only if there is a path from $j$
to $i$, so (i) gives condition b.

Besides, (ii) is equivalent to
$\operatorname{rank}\,T_{\mathcal{C},\sigma_{B}^{-1}(b)}=n_{C}$ for all
excitation $b$, and Lemma IV.1 yields condition c.

If more than one assignation $\sigma$ have a nonzero term
$\tau_{B}(\sigma_{B})$, the condition is not sufficient since terms could
cancel each other. But if there is only one such assignation, the condition is
also sufficient.

Proof of Theorem IV.1: Consider a generically decoupled-identifiable network.
Corollary IV.1 yields that $\det\hat{K}\neq 0$ generically. Combine (15) with
its analogous version for the measurements:

$\displaystyle\left(\det\hat{K}\right)^{2}$
$\displaystyle=\left(\sum_{\sigma_{B}}\operatorname{sgn}(\sigma_{B})\
\tau_{B}(\sigma_{B})\right)\left(\sum_{\sigma_{C}}\operatorname{sgn}(\sigma_{C})\
\tau_{C}(\sigma_{C})\right)$
$\displaystyle=\sum_{\sigma_{B}}\sum_{\sigma_{C}}\operatorname{sgn}(\sigma_{B},\sigma_{C})\
\tau_{B}(\sigma_{B})\ \tau_{C}(\sigma_{C}),$

where $\tau_{B}(\sigma_{B})$ stands for the contribution of assignation
$\sigma_{B}$ in $\det\hat{K}$ as defined in (15), and $\tau_{C}(\sigma_{C})$
is defined analogously.

Since $\det\hat{K}$ is generically nonzero, so is $(\det\hat{K})^{2}$, and
there is at least one pair of assignations $(\sigma_{B},\sigma_{C})$ such that
$\tau_{B}(\sigma_{B})\neq 0$ and $\tau_{C}(\sigma_{C})\neq 0$. Hence, the
connectivity and vertex-disjoint path conditions follow from Lemma IV.2.

If more than one pair $(\sigma_{B},\sigma_{C})$ have a nonzero term
$\tau_{B}(\sigma_{B})\,\tau_{C}(\sigma_{C})$, the condition is not sufficient
since terms could cancel each other. But if there is only one such pair, the
condition is also sufficient.