# Feasible Region for DERs in Unbalanced Distribution Networks Considering
Uncertain Network Impedances
††thanks: The authors would like to acknowledge support from the CSIRO
Strategic Project on Network Optimisation & Decarbonisation under Project
Number: OD-107890.

Bin Liu Member, IEEE Energy Systems Program, Energy Centre
Commonwealth Scientific and Industrial Research Organisation
Mayfield West, NSW 2304, Australia
<EMAIL_ADDRESS>Jin Ma Member, IEEE School of Electrical and Information
Engineering
The University of Sydney
Camperdown, NSW 2006, Australia
<EMAIL_ADDRESS>
###### Abstract

The rapid development of distributed energy resources (DERs) has brought many
challenges to the operation of distribution networks in recent years, where
operating envelopes (OEs), as a key enabler to facilitate DER integration via
virtual power plants, have attracted increasing attention. Geometrically, OEs
are inherently linked to the feasible region (FR) concept, which depicts an
area that quantifies the permissible operational ranges for all DERs. This
paper studies how to calculate FR for DERs based on a linear unbalanced three-
phase power flow model while considering uncertainties from measured/estimated
network impedances, leading to the robust feasible region (RFR). With an
initialised FR, the proposed solution algorithm aims to seek the RFR by
solving a non-convex bilinear problem, formulated based on Motzkin
Transposition Theorem, and updating the RFR alternatively until reaching
convergence. Two cases, including a 2-bus illustrative network and a
representative Australian network, have been tested to demonstrate the
effectiveness of the proposed approach.

###### Index Terms:

Distributed energy resources (DERs), distribution network, network
uncertainty, operational envelopes, robust feasible region, unbalanced optimal
power flow

##  I Introduction

Renewable energy, particularly distributed energy resources (DERs), has been
rapidly growing in power systems worldwide in the past years. In Australia, an
estimated 3.04 million homes and businesses have installed rooftop PV systems
by the end of 2021, with more than 2.88 GW in 349,000 systems installed in the
year 2021 alone [1]. High penetration of DERs is a crucial enabler to
decarbonise the energy sector through more engagement of residential
customers. However, due to relatively low visibility and controllability, high
DER penetration could lead to a series of operational issues, such as
over/under-voltage and minimum demand issues [2, 3, 4]. Addressing such issues
to integrate DERs needs close coordination among transmission network
operators (TSOs), distribution network operators (DSOs), emerging DER
aggregators in operating virtual power plants (VPPs) and market operators,
potential future architectures of which were investigated in [5].

Operating envelope (OE), which specifies the operational range for each DER
device that is permissible within the network operational limits at a given
point in time and location [6, 7], is identified as a key enabler in potential
future power system architectures and has gained increasing interest from both
industry and academia [6, 7, 8, 9, 10]. OEs can be used by DSOs to manage the
export and import limits for DERs based on available real-time network hosting
capacity and facilitate DER participation in the electricity market to provide
network services [11]. Substantial advances have been made in developing
approaches to calculating OEs in recent years, including the unbalanced three-
phase power flow (UTPF)-based approach, where OEs are calculated iteratively
by solving a series of PF problems with varying export/import limits [12, 8,
9, 13], the unbalanced three-phase optimal power flow (UTOPF)-based approach,
where OEs of all DERs are calculated simultaneously subject network
constraints [14, 7, 15], and machine learning (ML)-based approach, which is
useful for networks with low visibility [16, 6].

Geometrically, OEs are inherently linked to the feasible region (FR), which
depicts an area for controllable variables, especially for the PF and OPF-
based approaches [17, 18, 19, 20]. For example, calculating OEs based on OPF
is equivalent to seeking an _operational point_ within the FR under a given
objective function. Moreover, the FR is mathematically a polyhedron if a
linear OPF model is employed. However, uncertainties naturally exist in power
grids and may come from load forecast errors, inaccurate estimations or
measurements for network impedances. Such uncertainties can potentially
undermine the accuracy of FR and, as a result, lead to unreliable and non-
robust OEs. To address this issue, this paper proposes an approach to
calculate the _robust_ FR (RFR) for DERs in an unbalanced distribution network
while considering uncertainties from network impedances based on a linear UTPF
model. Uncertainties from network impedances are formulated as box sets on
either self and mutual impedances directly or on the positive/negative and
zero sequence impedances of all types of line codes. With formulated
uncertainties and an initialised FR, seeking the RFR is achieved by
alternatively solving a _max-min_ problem, which is formulated based on
Motzkin Transposition Theorem (MTT), and updating the FR based on reported
solutions, until reaching convergence. Both a 2-bus illustrative network and a
representative Australian network aa tested to demonstrate the effectiveness
of the proposed approach.

It is noteworthy that the proposed approach only applies to cases based on
linear UTPF models, and we admit that errors inherently exist for linear UTPF
formulations and are unavoidable, which can further contribute to the
inaccuracy of both FR and RFR. Therefore, improving the accuracy of linear
UTPF will also help improve the effectiveness of the RFR, and calculating an
RFR based on non-convex UTOPF models is beyond the scope of this paper.

The remainder of the paper is organised as follows. Section II presents the
formulation of FR both with and without considering uncertainties from network
impedances. Section III discusses the solution algorithm to seeking the RFR,
followed by case studies in Section IV, and the paper is concluded in Section
V.

##  II Feasible Region for DERs

###  II-A UTPF with operational constraints

UTPF is an essential tool for the distribution network analysis to capture its
operational states [21, 22]. Although a typical low-voltage distribution
network may contain four wires, including the neutral conductor, this paper
assumes that the neutral conductor is well grounded and the line parameters
for the equivalent three-phase networks can be calculated via Carson’s
Equations [23].

Based on its current and voltage formulation, the typical UTPF model with
operational constraints can be formulated as (1).

$\displaystyle V^{\phi}_{i_{\text{ref}}}=V^{\phi}_{0}\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\phi$ (1a) $\displaystyle
V^{\phi}_{i}-V^{\phi}_{j}=\sum\nolimits_{\psi}{z_{ij}^{\phi\psi}I^{\phi}_{ij}}\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\phi,\forall ij$ (1b)
$\displaystyle\sum_{n:n\rightarrow i}{I^{\phi}_{ni}}-\sum_{m:i\rightarrow
m}{I^{\phi}_{im}}=\sum_{m}{\frac{\mu_{\phi,i,m}(P^{\text{}}_{m}-\text{j}Q^{\text{}}_{m})}{(V^{\phi}_{i})^{*}}}\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\phi,\forall i$ (1c) $\displaystyle
V^{\text{min}}_{i}\leq|V^{\phi}_{i}|\leq
V^{\text{max}}_{i}\leavevmode\nobreak\ \forall\phi,\forall i$ (1d)

where $i_{\text{ref}}$ is the index of the reference bus and $V^{\phi}_{0}$ is
its fixed voltage at phase $\phi$ (known parameter). $V^{\phi}_{i}$ is the
voltage of phase $\phi$ at node $i$. $I^{\phi}_{ni}$ is the current in phase
$\phi$ of line $ni$: flowing from bus $n$ to bus $i$. $P_{m}$ is the active
power demand of customer $m$ while $Q_{m}$ is the reactive power demand. For
simplicity, all $P_{m}$ and $Q_{m}$ are treated as variables in the
formulation. However, they will be fixed as the forecasted values if they are
uncontrollable. $\mu_{\phi,i,m}\in\\{0,1\\}$ is a parameter indicating the
phase connection of customer $m$ with its value being $1$ if it is connected
to phase $\phi$ of bus $i$ and being $0$ otherwise. $V^{\text{min}}_{i}$ and
$V^{\text{max}}_{i}$ are the lower and upper limits of $|V_{i}|$,
respectively.

In the formulation, (1a) specifies the voltage at the reference bus. (1b) is
the voltage drop equation in each line, and Kirchhoff’s current law is
formulated as (1c). Only constraints on voltage magnitudes are considered in
this paper, and they are expressed as (1d).

As discussed previously, a linear UTPF (LUTPF) model will be used, where the
physical characteristics of a real distribution network, including the reality
in most cases that differences of voltage angles in each phase are
sufficiently small [24, 25] and nodal voltages throughout the network are
close to $1.0\leavevmode\nobreak\ p.u.$, will be taken. Specifically, (1c) can
be linearised by fixing $V^{\phi}_{i}$ in the denominator on the right-hand
side, while (1d) will be approximated by a set of linear constraints. Details
of the linearisation approach can be found in [26] and are omitted here for
simplicity. For the convenience of later discussions, LUTPF can be expressed
in the following compact form,

$\displaystyle Ap+Bq+Cl=b$ (2a) $\displaystyle Dv+El=d$ (2b) $\displaystyle
Fv\leq f$ (2c)

where $p$ and $q$, which are independent variables, represent the active and
reactive powers; $l$ represents the vector consisting of state variables of
line currents; $v$ is a vector consisting of state variables related to nodal
voltages; $A,B,C,b,D,E,d,F$ and $f$ are constant parameters with appropriate
dimensions;

In the formulation, (2a) links back to (1c) after linearisation and represents
the relations between line currents and residential demands; (2b) represents
the linearised power flow equations that link the power demands and currents
running in all lines, i.e. (1a)-(1b), and (2c) represents all the operational
constraints after the linearisation.

###  II-B FR without network uncertainties

Noting that only $p$ and $q$ are independent variables, (2) defines the FR for
$(p,q)$. This paper focuses on deriving the feasibility region with $p$ as the
independent variable, as most residential homes consume active power loads.
Correspondingly, the FR depicts an area in which, if all realised values of
$p$ fall within, the operational security of the network can be guaranteed
[18, 19, 20]. After removing state variables $v$ and $l$, the FR for $p$ can
be expressed as (3).

$\displaystyle\mathcal{F}(q)=\\{p|Hp\leq
h\\}=\left\\{p\left|\begin{matrix}FD^{-1}EC^{-1}Ap\leq f\\\
-FD^{-1}(d-EC^{-1}b+EC^{-1}Bq)\\\ \end{matrix}\right.\right\\}$ (3)

where $H=FD^{-1}EC^{-1}A$ and $h=f-FD^{-1}(d-EC^{-1}b\\\ +EC^{-1}Bq)$.

Obviously, the $\mathcal{F}(q)$ is a function of $q$, and also a polyhedron if
$A,B,C,b,D,E,d,F$ and $f$ are fixed111For the convenience of later
discussions, parameters will not be specifically indicated that it is a
function of $q$ in following paragraphs.. However, uncertainties inherently
exist in real power systems and can arise from errors in forecasting load
profiles for passive customers, incorrectly recorded phase connectivity for
residential customers and network topology, and inaccuracy in measuring or
estimating network impedances. This paper will focus on uncertainty modelling
for network impedances and their impacts on the FR.

###  II-C FR considering uncertainties from network impedances

####  II-C1 Uncertainty modelling

In this paper, uncertainties on network impedances will be modelled as box
sets either on both self-impedances and mutual impedances directly, as
presented in (4), or on both the positive/negative and zero-sequence
impedances, as indicated by (5).

$\displaystyle z_{ij}^{\phi\psi,\text{min}}\leq z_{ij}^{\phi\psi}\leq
z_{ij}^{\phi\psi,\text{max}}\leavevmode\nobreak\ \forall
ij,\forall\phi,\forall\psi$ (4)

where $z_{ij}^{\phi\psi,\text{min}}$ and $z_{ij}^{\phi\psi,\text{max}}$ are
the low and upper bounds for the $z_{ij}^{\phi\psi}$, respectively.

$\displaystyle z_{ij}^{+,\text{min}}\leq z_{ij}^{+}=z_{ij}^{-}\leq
z_{ij}^{+,\text{max}}\leavevmode\nobreak\ \forall ij$ (5a) $\displaystyle
z_{ij}^{0,\text{min}}\leq z_{ij}^{0}\leq
z_{ij}^{0,\text{max}}\leavevmode\nobreak\ \forall ij$ (5b) $\displaystyle
z_{ij}^{\phi\psi}=(2z_{ij}^{+}+z_{ij}^{0})/3\leavevmode\nobreak\ \forall
ij,\forall\phi,\forall\psi\leavevmode\nobreak\ (\psi=\phi)$ (5c)
$\displaystyle
z_{ij}^{\phi\psi}=(-z_{ij}^{+}+z_{ij}^{0})/3\leavevmode\nobreak\ \forall
ij,\forall\phi,\forall\psi\leavevmode\nobreak\ (\psi\neq\phi)$ (5d)

where $z_{ij}^{+},z_{ij}^{-}$ and $z_{ij}^{0}$ are the positive, negative and
zero-sequence impedances for line $ij$, respectively, with their lower or
upper limits indicated by
$z_{ij}^{+,\text{min}},z_{ij}^{+,\text{max}},z_{ij}^{0,\text{min}}$ and
$z_{ij}^{0,\text{max}}$.

For the convenience of later discussion, $z$ is used to represent network
impedances, and its uncertainty is compactly formulated as $z^{\text{min}}\leq
z\leq z^{\text{max}}$. Noting that $z_{ij}^{\phi\psi}$ are coefficients only
for the _current_ terms in (1b)-(1c), in the compact form (2), such
uncertainties will be broadcasted to $E$ and also $H$ and $h$. Specifically,
taking $E=f(z)$ to represent the linear mapping from uncertain network
impedances $z$ to the matrix $E$, the uncertainty set for the pair $(H,h)$ can
be formulated as

$\displaystyle\mathcal{H}=\left\\{(H,h)\left|\begin{matrix}H=FD^{-1}f(z)C^{-1}A\\\
h=f-FD^{-1}(d-f(z)C^{-1}b+f(z)C^{-1}Bq)\\\ z^{\text{min}}\leq z\leq
z^{\text{max}}\\\ \end{matrix}\right.\right\\}$ (6)

####  II-C2 Mathematical formulation

With uncertain $E,H$ and $h$ in (3), to secure the network operation, it is
critical to find a more _conservative_ yet more reliable and robust FR for all
DERs’ generation, which will be a subset of $\mathcal{F}$. The more reliable
and robust FR, which is named _robust feasible region (RFR)_ in this paper,
can be mathematically defined as

$\displaystyle\mathcal{R}=\\{p|Gp\leq
g\\}=\cap_{(H,h)\in\mathcal{H}}\mathcal{F}=\\{p|Hp\leq h\leavevmode\nobreak\
\forall(H,h)\in\mathcal{H}\\}$ (7)

which is equivalent to

$\displaystyle\\{p|Gp\leq
g\\}\cap\\{p|H_{i,:}p>h_{i}\\}=\emptyset\leavevmode\nobreak\
\forall(H,h)\in\mathcal{H},\forall i$ (8)

where $H_{i,:}$ represents the $i-th$ row of $H$ and $h_{i}$ is the $i-th$
element in $h$.

Obviously, the $\mathcal{R}$ is also a polyhedron. However, noting that
uncertainties from $H$ and $h$ are depicted by a closed set, there are
indefinite constraints in (8) and a suitable solution algorithm needs to be
developed to get the exact formulation of the RFR.

##  III Solution algorithm

###  III-A Reformulation of (8) to an optimisation problem

The reformulation is based on the MTT [27], where we have the following
proposition for our case.

######  Proposition III.1.

The formulation (8) is equivalent to

$\displaystyle G^{T}x_{:,i}=H^{T}_{i,:},\leavevmode\nobreak\
\leavevmode\nobreak\ x_{:,i}\geq 0,\leavevmode\nobreak\ \leavevmode\nobreak\
g^{T}x_{:,i}\leq h_{i}\leavevmode\nobreak\ \forall i$ (9)

where $x$ is a matrix variable and and $x_{:,i}$ is its $i-th$ column.

###### Proof.

Given the matrix $A$, vectors $b,d$ and a number $c$, MTT indicates that only
one of the following two statements holds.

  1. 1.

The system $Ax\leq b,d^{T}x\leq c$ has a solution $x$;

  2. 2.

The system

$\displaystyle y\geq 0,\leavevmode\nobreak\ \leavevmode\nobreak\ z\geq
0,\leavevmode\nobreak\ \leavevmode\nobreak\ A^{T}y+z\cdot
d=0,\leavevmode\nobreak\ \leavevmode\nobreak\ b^{T}y+c\cdot z<0$ (10)

has a solution $(y,z)$.

With MTT, $\mathcal{R}=\\{p|Gp\leq g\\}$ being a subset of
$\mathcal{F}=\\{p|Hp\leq h\\}$ is equivalent to
$\mathcal{R}\cap\\{p|H_{i,:}p>h_{i}\\}=\emptyset\leavevmode\nobreak\ (\forall
i)$. Further, $\mathcal{R}$ and $H_{i,:}p>h_{i}$ can be reformulated as (11a)
and (11b), respectively.

$\displaystyle\left[G,\leavevmode\nobreak\ -G,\leavevmode\nobreak\
I_{s}\right]\left[\begin{matrix}p^{+}\\\ p^{-}\\\
s\end{matrix}\right]+1\cdot(-g)=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\left[\begin{matrix}p^{+}\\\ p^{-}\\\ s\end{matrix}\right]\geq 0$ (11a)
$\displaystyle\left[-H_{i,:},\leavevmode\nobreak\ H_{i,:},\leavevmode\nobreak\
0_{s}^{T}\right]\left[\begin{matrix}p^{+}\\\ p^{-}\\\
s\end{matrix}\right]+h<0$ (11b)

where $y=y^{+}-y^{-}$ and $s$ is an auxiliary vector variable.

After applying MTT, $\mathcal{R}\cap\\{p|H_{i,:}p>h_{i}\\}=\emptyset$ is
equivalent to

$\displaystyle G^{T}y=H_{i,:}^{T},\leavevmode\nobreak\ y\geq
0,\leavevmode\nobreak\ g^{T}y\leq h_{i}$ (12)

with $y$ being a vector variable and, accordingly, (8) can be equivalently
reformulated as

$\displaystyle G^{T}x_{:,i}=H^{T}_{i,:},\leavevmode\nobreak\
\leavevmode\nobreak\ x_{:,i}\geq 0,\leavevmode\nobreak\ \leavevmode\nobreak\
g^{T}x_{:,i}\leq h_{i}\leavevmode\nobreak\ \forall i$ (13)

where $x$ is now a matrix variable. ∎

Based on the proposition, if (8) holds for fixed $H$ and $h$, the optimal
objective value of (14), which is given below, must be 0.

$\displaystyle\min_{x,\Delta h}{1^{T}\Delta h}$ (14a) $\displaystyle
s.t.\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
G^{T}x=H^{T},\leavevmode\nobreak\ g^{T}x\leq h^{T}+\Delta h^{T}$ (14b)
$\displaystyle x\geq 0,\leavevmode\nobreak\ \Delta h\geq 0$ (14c)

Similarly, if (8) holds for all possible $H$ and $h$, the optimal objective
value of (15) will be 0.

$\displaystyle\max_{(H,h)\in\mathcal{H}}\min_{x,\Delta h}{1^{T}\Delta h}$
(15a) $\displaystyle s.t.\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ G^{T}x_{:,i}=H^{T}_{i}\leavevmode\nobreak\ \forall
i\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ (\alpha_{i})$
(15b) $\displaystyle x_{:,i}\geq 0,\leavevmode\nobreak\ \Delta h_{i}\geq
0\leavevmode\nobreak\ \forall i$ (15c) $\displaystyle g^{T}x_{:,i}\leq
h_{i}+\Delta h_{i}\leavevmode\nobreak\ \forall i\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ (\beta_{i})$ (15d)

where $\alpha_{i}$, as a vector, and $\beta$, as a scaler, are Lagrange
multipliers.

Based on strong duality theory, (15) is equivalent to

$\displaystyle\text{BiLP-Max:\leavevmode\nobreak\ \leavevmode\nobreak\
}\max_{(H,h)\in\mathcal{H},\alpha,\beta}{\sum\nolimits_{i}{(-\alpha_{i}^{T}H^{T}_{i,:}-\beta_{i}h_{i}})}$
(16a) $\displaystyle s.t.\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ 0\leq\beta_{i}\leq 1\leavevmode\nobreak\ \forall i$ (16b)
$\displaystyle\alpha_{i}^{T}G^{T}+\beta_{i}g^{T}\geq 0\leavevmode\nobreak\
\forall i$ (16c)

Obviously, (16) is a non-convex bilinear programming problem, where several
solution approaches are available. One of them is reformulating (16) into a
mixed-integer linear programming (MILP) problem noting that the optimal
solutions of $(H,h)$ must fall on the extreme points of $\mathcal{H}$.
However, such an approach could become computationally intractable if there is
a large number of uncertain variables [28, 18]. Another approach is based on
iteration-based linear programming (ITLP), where a series of linear
programming (LP) problems will be solved by alternatively fixing $(H,h)$ or
$(\alpha,\beta)$ in (16) until reaching convergence. Although the ITLP-based
approach cannot guarantee a global solution, it has been demonstrated with
good performance with sufficient and properly selected starting points on
$(H,h)$ [19, 28]. Noting that the number of uncertain parameters can be large
in our case, the ITLP-based approach will be employed to solve (15).

###  III-B Updating $\mathcal{R}$ and removing redundant constraints

After solving (16), a set of linear inequalities will be generated if the
objective value is positive and added to update $\mathcal{R}$. With reported
solutions, say $(H^{*},h^{*})$, from solving (16), the RFR will be updated as

$\displaystyle\text{U-RFR:\leavevmode\nobreak\ \leavevmode\nobreak\
}\mathcal{R}=\mathcal{R}\cap\\{p|H^{*}p\leq f^{*}\\}$ (17)

Further, to reduce the computational burden in solving (16), before conducting
the next iteration, redundant constraints in the updated $\mathcal{R}$ will be
removed first [29]. In summary, a schematic diagram showing the process of
seeking the RFR for DERs is given in Fig.1.

Figure 1: The schematic diagram of calculating the robust feasible region for
DERs.

Several remarks on calculating the RFR for DERs are given below.

  1. 1.

Noting that in a real network, the initial value of $\mathcal{R}$, denoted as
$\bar{\mathcal{R}}$, can be initialised after fixing $z$ as its measured or
estimated value.

  2. 2.

Although controllable reactive power $q$ can be regarded as a variable, it is
fixed as a constant when calculating the RFR in this paper.

  3. 3.

Properly selecting starting points for $(H,h)$ is critical to ensure the
quality of the solutions. In this paper, all projected points from
$(\widetilde{H},\widetilde{h})$ to its uncertainty set boundaries will be used
as starting points, which, as demonstrated in [19], can lead to a high-quality
solution for (16).

  4. 4.

After solving (15), several solutions, which link back to the various starting
point for $(H,h)$, can be reported. In such cases, $\mathcal{R}$ will be
updated by considering all solutions in this paper.

  5. 5.

All LP problems in this paper are formulated in Julia [30] with data parsed by
PowerModelsDistribution.jl [31], and solved by Gurobi [32] on a laptop with
Intel Core i7-1185G7 @ 3.00GHz and 16 GB RAM. It is noteworthy that parallel
computing techniques can potentially accelerate computational efficiency,
noting that (16) with different starting points can be solved simultaneously.
However, investigating how parallel computing techniques could improve
computational performance is beyond the scope of this paper and falls within
our future research interest.

##  IV Case Study

###  IV-A Case setup

In this section, two distribution networks, one of which is a 2-bus
illustrative network and the other one is a representative Australian network,
will be studied. For the illustrative distribution network, an ideal balanced
voltage source with the voltage magnitude being $1.0\leavevmode\nobreak\ p.u.$
is connected to bus 1. A three-phase line connects bus 1 and bus 2, and its
recorded impedance matrix is

$\bar{z}_{12}=\left[\begin{smallmatrix}0.3465+\text{j}1.0179&0.1560+\text{j}0.5017&0.1580+\text{j}0.4236\\\
0.1560+\text{j}0.5017&0.3375+\text{j}1.0478&0.1535+\text{j}0.3849\\\
0.1580+\text{j}0.4236&0.1535+\text{j}0.3849&0.3414+\text{j}1.0348\end{smallmatrix}\right]\Omega$

and, given a number $\delta\in[0,1]$ representing the uncertainty level, we
also assume the true network impedances are within the range
$[1-\delta,1+\delta]$ of its recorded value, i.e.

$\displaystyle(1-\delta)\bar{z}_{12}\leq z_{12}\leq(1+\delta)\bar{z}_{12}$
(18)

and $\delta$ is set as 20% for the 2-bus network.

Three single-phase customers, numbered 1 to 3, are connected to phase $b$, $a$
and $c$ of bus 2, respectively, where all reactive powers are fixed at 0 kVar,
and the active power of customer 2 is fixed at -2.0 kW. Moreover, both
customers’ default export/import limits are set as 20 kW. For this 2-bus
network, only voltage magnitude constraints are considered for bus 2, where
the lower and upper bounds for each phase are set as $0.95\leavevmode\nobreak\
p.u.$ and $1.05\leavevmode\nobreak\ p.u.$, respectively.

The topology of the representative Australian distribution network is given in
Fig. 2, and network parameters can be found in [33] for network J in OpenDSS
format [34], where the transformer has been changed to $Y_{n}/Y_{n}$
connection with $R$%=5 and $X$%=7. The voltage at the reference bus, grid in
Fig.2, is set as $[1.0,\leavevmode\nobreak\
1.0e^{-\text{j}\frac{2\pi}{3}},\leavevmode\nobreak\
1.0e^{\text{j}\frac{2\pi}{3}}]^{T}$ while the upper and lower limits on
voltage magnitudes of all the other buses are set as 1.05 $p.u.$ and 0.95
$p.u.$, respectively. Similar to the 2-bus network, the reactive powers of all
customers are set to 0 kVar, and the active powers of 10 customers (customers:
“1”, “10”, “11”, “12”, “13”, “14”, “15”, “16”, “17”, “18”), where default
export and import limits are set as 5 kW and 6 kW, respectively, form the
independent variables in calculating the FR. In contrast, demands of all the
other active customers are fixed at the values provided in [33]. In the
Australian network, network impedances are generated from 5 line codes, where
uncertainties are assumed on both their positive/negative and zero-sequence
impedances ($\delta$=10%).

Figure 2: Network topology of the real Australian network (Upstream grid
(reference bus) and its connected bus are marked in yellow).

###  IV-B The 2-bus illustrative distribution network

As only one line exists in the 2-bus illustrative network, there are 12
independent parameters with uncertainties, leading to a total number of 24
initial points by projecting $\bar{z}$ to the boundaries of its uncertainty
set for solving (16). The proposed approach takes 42.26 seconds in total after
three iterations to report the RFR, where the number of reported solutions and
the computational time for each iteration are presented in Fig.3. As shown in
the figure, the number of unique solutions reported from solving the BiLP-Max
problem becomes smaller while the objective value decreases along with
increasing iterations. Regarding the computational time, the solver takes most
of the time in solving BiLP-Max while removing redundant constraints is much
more efficient, averagely taking 0.2 seconds in each iteration.

Figure 3: Computational and solution information for the 2-bus illustrative
network.

The initial FR, the finally reported RFR, and FRs calculated with randomly
generated scenarios for $z$ are presented in Fig.4, where the calculated RFR
is a subset of both the randomly generated FRs and the initialised FR,
demonstrating its robustness.

Figure 4: FRs calculated for the 2-bus illustrative network.

Moreover, the percentage of random FRs containing the RFR and initialised FR
are also calculated, where the number is 100% for RFR and 0% for the
initialised FR. As a higher percentage implies higher reliability and
robustness of FR against uncertainties, such results further demonstrate the
effectiveness of the proposed approach.

###  IV-C The real Australian network

The Australian network is used to test the computational efficiency of the
proposed approach in a real network. Simulation results, including the number
of reported solutions, objective value and computational time for each
iteration, are presented in Fig.5.

Figure 5: Computational and solution information for the Australian network.

Similar to the 2-bus network, 100 FRs are also randomly generated to test the
effectiveness of the RFR, where the percentage of random FRs containing the
RFR and the initialised FR are 62% and 0%, respectively, where the number for
RFR being less than 100% is due to the local optimality in solving the BiLP-
Max problem by the ITLP approach.

Regarding the computational process, the proposed approach reports the RFR
within three iterations, where each iteration consumes 130-180 seconds.
Although the computational efficiency is acceptable if the RFR is used for
planning purposes, we admit that it may not be suitable for intra-day or real-
time operation purposes if the dispatch interval is less than 10 minutes.
However, the practicality of the proposed approach can be significantly
enhanced by the following measures.

  * •

Further investigating the uncertainty modelling of network impedances to
reduce the number of uncertain parameters and the conservativeness of
uncertainty set.

  * •

Applying parallel computing techniques (PCTs). Nothing that solving the BiLP-
Max problem with various starting points of $(H,h)$ in each iteration can be
paralleled, PCTs could considerably enhance the computational performance of
the proposed approach.

##  V Conclusion

The paper studies how to calculate feasible region for DERs in an unbalanced
distribution network while considering uncertainties from network impedances.
The proposed approach, with an initialised FR, tries to seek the optimal
solution by solving a non-convex bilinear optimisation problem and updating
the RFR successively until reaching convergence, which has been demonstrated
effective in two distribution networks. To make the proposed approach more
practical, several research directions are open in this area and fall within
our future research interest, including a) Improving computational efficiency
in networks with more types of uncertainties and more DERs, particularly in
properly modelling uncertainties and efficiently solving the BiLP-Max problem;
2) Investigating the implications of RFR on calculating operating envelopes.

## References

  * [1] “Solar Report,” Australian Energy Council, Melbourne, Australia, Tech. Rep., 2022. [Online]. Available: https://www.energycouncil.com.au/media/5wkkaxts/australian-energy-council-solar-report_-jan-2022.pdf
  * [2] M. M. Haque and P. Wolfs, “A review of high PV penetrations in LV distribution networks: Present status, impacts and mitigation measures,” _Renewable Sustainable Energy Rev._ , vol. 62, pp. 1195–1208, 2016.
  * [3] A. Kharrazi, V. Sreeram, and Y. Mishra, “Assessment techniques of the impact of grid-tied rooftop photovoltaic generation on the power quality of low voltage distribution network - A review,” _Renewable Sustainable Energy Rev._ , vol. 120, pp. 1–16, 2020.
  * [4] AEMO, “Negative electricity demand in South Australia ,” 2021. [Online]. Available: https://aemo.com.au/newsroom/news-updates/negative-electricity-demand-in-south-australia
  * [5] “Open Energy Networks Project - Energy Networks Australia Position Paper,” ENA and AEMO, Melbourne, Australia, Tech. Rep., 2020. [Online]. Available: https://tinyurl.com/4bd9tr5n
  * [6] “Dynamic Operating Envelopes working group outcomes report,” Distributed Energy Integration Program (DEIP), Melbourne, Australia, Tech. Rep. March, 2022\. [Online]. Available: https://arena.gov.au/assets/2022/03/dynamic-operating-envelope-working-group-outcomes-report.pdf
  * [7] K. Petrou, A. T. Procopiou, L. Gutierrez-Lagos, M. Z. Liu, L. F. Ochoa, T. Langstaff, and J. Theunissen, “Ensuring Distribution Network Integrity Using Dynamic Operating Limits for Prosumers,” _IEEE Trans. Smart Grid_ , vol. 12, no. 5, pp. 3877–3888, 2021.
  * [8] “Project EDGE,” https://aemo.com.au/en/initiatives/major-programs/nem-distributed-energy-resources-der-program/der-demonstrations/project-edge, accessed: 2022-04-12.
  * [9] “Project Symphony,” https://aemo.com.au/initiatives/major-programs/wa-der-program/project-symphony, accessed: 2022-04-12.
  * [10] “DER Implementation Plan – reform activities over three-year horizon,” Energy Security Board (ESB), Melbourne, Australia, Tech. Rep., 2021. [Online]. Available: https://tinyurl.com/yck2vxm6
  * [11] M. Z. Liu, L. N. Ochoa, S. Riaz, P. Mancarella, T. Ting, J. San, and J. Theunissen, “Grid and market services from the edge: using operating envelopes to unlock network-aware bottom-Up flexibility,” _IEEE Power and Energy Magazine_ , vol. 19, no. 4, pp. 52–62, 2021.
  * [12] L. Blackhall, “On the calculation and use of dynamic operating envelopes,” Zepben, Canberra, Australia, Tech. Rep., 2020. [Online]. Available: https://arena.gov.au/assets/2020/09/on-the-calculation-and-use-of-dynamic-operating-envelopes.pdf
  * [13] B. Liu and J. H. Braslavsky, “Sensitivity and Robustness Issues of Operating Envelopes in Unbalanced Distribution Networks,” _IEEE Access_ , vol. 10, no. September, pp. 92 789–92 798, 2022.
  * [14] K. Petrou, M. Z. Liu, A. T. Procopiou, and L. N. Ochoa, “Managing residential prosumers using operarting envelopes: An Australian case study,” in _CIRED Berlin Workshop_ , no. September, Berlin, Germany, 2020, pp. 1–6.
  * [15] M. Z. Liu, L. F. Ochoa, P. K. Wong, and J. Theunissen, “Using OPF based operating envelopes to facilitate residential DER services,” _IEEE Trans. Smart Grid_ , pp. 1–1, 2022.
  * [16] L. N. Ochoa, V. Bassi, D. Jaglal, T. Alpcan, and C. Leckie, “The Future of DER Hosting Capacity and Operating Envelopes,” University of Melbourne, Melbourne, Australia, Tech. Rep., 2022. [Online]. Available: https://www.researchgate.net/profile/Luisnando-Ochoa-2/publication/360513775_The_Future_of_DER_Hosting_Capacity_and_Operating_Envelopes/links/627b4387107cae29199bb105/The-Future-of-DER-Hosting-Capacity-and-Operating-Envelopes.pdf
  * [17] B. Liu and J. Braslavsky, “Technical Challenges in the Calculation and Implementation of Operating Envelopes for DER Grid Integration in Australia,” 2022. [Online]. Available: https://smartgrid.ieee.org/bulletins/august-2022/technical-challenges-in-the-calculation-and-implementation-of-operating-envelopes-for-der-grid-integration-in-australia
  * [18] W. Wei, F. Liu, and S. Mei, “Real-time dispatchability of bulk power systems with volatile renewable generations,” _IEEE Transs. Sustain. Energy_ , vol. 6, no. 3, pp. 738–747, 2015.
  * [19] ——, “Dispatchable Region of the Variable Wind Generation,” _IEEE Trans. Power Syst._ , vol. 30, no. 5, pp. 2755–2765, 2015.
  * [20] S. Riaz and P. Mancarella, “Modelling and Characterisation of Flexibility from Distributed Energy Resources,” _IEEE Transs. Power Syst._ , vol. 37, no. 1, pp. 38–50, 2022.
  * [21] T. Jen-Hao, “A direct approach for distribution system load flow solutions,” _IEEE Trans. Power Del._ , vol. 18, no. 3, pp. 882–887, 2003.
  * [22] H. Abdi, S. D. Beigvand, and M. L. Scala, “A review of optimal power flow studies applied to smart grids and microgrids,” _Renewable and Sustainable Energy Reviews_ , vol. 71, pp. 742–766, 2017.
  * [23] W. H. Kersting, _Distribution System Modeling and Analysis_ , 4th ed. Boca Raton: CRC Press, 2017.
  * [24] H. Ahmadi, J. R. Marti, and A. von Meier, “A linear power flow formulation for three-phase distribution systems,” _IEEE Trans. Power Syst._ , vol. 31, no. 6, pp. 5012–5021, 2016.
  * [25] L. Gan and S. H. Low, “Convex relaxations and linear approximation for opf in multiphase radial networks,” in _Proc. Power Systems Computation Conf._ , Dublin, Ireland, 2018.
  * [26] B. Liu, F. Geth, N. Mahdavi, and J. Zhong, “Load balancing in low-voltage distribution networks via optimizing residential phase connections,” in _2021 IEEE PES Innovative Smart Grid Technologies - Asia (ISGT Asia)_ , Brisbane, Australia, 2021, pp. 1–5.
  * [27] A. Ben-Israel, “Motzkin’s Transposition Theorem, And The Related Theorems Of Farkas, Gordan And Stiemke,” _Encyclopaedia of Mathematics_ , pp. 1–4, 2001\.
  * [28] B. Liu, F. Liu, W. Wei, K. Meng, Z. Y. Dong, and W. Zhang, “Probabilistic evaluation of a power system’s capability to accommodate uncertain wind power generation,” _IET Renewable Power Generation_ , vol. 13, no. 10, pp. 1780–1788, 2019.
  * [29] K. Fukuda, “Polyhedral Computation,” Department of Mathematics, and Institute of Theoretical Computer Science, ETH Zurich, Zurich, Switzerland, Tech. Rep., 2014. [Online]. Available: https://people.inf.ethz.ch/fukudak/lect/pclect/notes2014/PolyComp2014.pdf
  * [30] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” _SIAM review_ , vol. 59, no. 1, pp. 65–98, 2017.
  * [31] D. M. Fobes, S. Claeys, F. Geth, and C. Coffrin, “PowerModelsDistribution.jl: An open-source framework for exploring distribution power flow formulations,” _Electric Power Systems Research_ , vol. 189, p. 106664, 2020\.
  * [32] Gurobi Optimization, “Gurobi optimizer reference manual,” 2018. [Online]. Available: http://www.gurobi.com
  * [33] CSIRO, “National low-voltage feeder taxonomy study,” 2021, [accessed 11-April-2022]. [Online]. Available: https://near.csiro.au/assets/f325fb3c-2dcd-410c-97a8-e55dc68b8064
  * [34] “OpenDSS,” https://www.epri.com/pages/sa/opendss#:~:text=OpenDSS%20is%20an%20electric%20power,grid%20integration%20and%20grid%20modernization, accessed: 2022-04-29.


# REAL-TIME DENOISING AND DEREVERBERATION WTIH TINY RECURRENT U-NET

###### Abstract

Modern deep learning-based models have seen outstanding performance
improvement with speech enhancement tasks. The number of parameters of state-
of-the-art models, however, is often too large to be deployed on devices for
real-world applications. To this end, we propose Tiny Recurrent U-Net (TRU-
Net), a lightweight online inference model that matches the performance of
current state-of-the-art models. The size of the quantized version of TRU-Net
is 362 kilobytes, which is small enough to be deployed on edge devices. In
addition, we combine the small-sized model with a new masking method called
phase-aware $\beta$-sigmoid mask, which enables simultaneous denoising and
dereverberation. Results of both objective and subjective evaluations have
shown that our model can achieve competitive performance with the current
state-of-the-art models on benchmark datasets using fewer parameters by orders
of magnitude.

Index Terms—  real-time speech enhancement, lightweight network, denoising,
dereverberation

##  1 Introduction

In this paper, we focus on developing a deep learning-based speech enhancement
model for real-world applications that meets the following criteria: 1\. a
small and fast model that can reduce single-frame real-time-factor (RTF) as
much as possible while keeping competitive performance against the state-of-
the-art deep learning networks, 2. a model that can perform both the denoising
and derverberation simultaneously.

To address the first issue, we aim to improve a popular neural architecture,
U-Net [1], which has proven its superior performance on speech enhancement
tasks [2, 3, 4]. The previous approaches that use U-Net on source separation
applications apply convolution kernel not only on the frequency-axis but also
on the time-axis. This non-causal nature of U-Net increases computational
complexity because additional computations are required on past and future
frames to infer the current frame. Therefore, it is not suitable for online
inference scenarios where the current frame needs to be processed in real-
time. In addition, the time-axis kernel makes the network computationally
inefficient because there exists redundant computation between adjacent frames
in both the encoding and decoding path of U-Net. To tackle this problem, we
propose a new neural architecture, Tiny Recurrent U-Net (TRU-Net), which is
suitable for online speech enhancement. The architecture is designed to enable
efficient decoupling of the frequency-axis and time-axis computations, which
makes the network fast enough to process a single frame in real-time. The
number of parameters of the proposed network is only 0.38 million (M), which
is small enough to deploy the model not only on a laptop but also on a mobile
device and even on an embedded device combined with a quantization technique
[5]. The details of TRU-Net is described more in section 2.

Next, to suppress the noise and reverberation simultaneously, we propose a
phase-aware $\beta$-sigmoid mask (PHM). The proposed PHM is inspired by [6],
in which the authors propose to estimate phase by reusing an estimated
magnitude mask value from a trigonometric perspective. The major difference
between PHM and the approach in [6] is that PHM is designed to respect the
triangular relationship between the mixture, the target source, and the
remaining part, hence the sum of the estimated target source and the remaining
part is always equal to the mixture. We extend this property into a
quadrilateral by producing two different PHMs simultaneously, which allows us
to effectively deal with both denoising and dereverberation. We will discuss
PHM in further details in section 3.

##  2 Tiny Recurrent U-Net

Fig. 1: The network architecture of TRU-Net

###  2.1 PCEN feature as an input

A spectrogram is perhaps the most popular input feature for many speech
enhancement models. Per-channel energy normalization (PCEN) [7] combines both
dynamic range compression and automatic gain control together, which reduce
the variance of foreground loudness and supress background noise when applied
to a spectrogram [8]. PCEN is also suitable for online inference scenarios as
it includes a temporal integration step, which is essentially a first-order
infinite impulse response filter that depends solely on a previous input
frame. In this work, we employ the trainable version of PCEN.

###  2.2 Network architecture

TRU-Net is based on U-Net architecture, except that the convolution kernel
does not span the time-axis. Therefore, it can be considered a frequency-axis
U-Net with 1D Convolutional Neural Networks (CNNs) and recurrent neural
networks in the bottle-neck layer. The encoder is composed of 1D Convolutional
Neural Network (1D-CNN) blocks and a Frequency-axis Gated Recurrent Unit
(FGRU) block. Each 1D-CNN block is a sequence of pointwise convolution and
depthwise convolution similar to [9], except the first layer, which uses the
standard convolution operation without a preceding pointwise convolution. To
spare the network size, we use six 1D-CNN blocks, which downsample the
frequency-axis size from 256 to 16 using strided convolutions. This results in
a small receptive field (1,750Hz) which may be detrimental to the network
performance. To increase the receptive field, we use a bi-directional GRU
layer [10] along the frequency-axis instead of stacking more 1D-CNN blocks.
That is, the sequence of 16 vectors from 1D-CNN blocks is passed into the bi-
directional GRU to increase the receptive field and share the information
along the frequency-axis. We call this frequency-axis bi-directional GRU layer
an FGRU layer. A pointwise convolution, batch normalization (BN), and
rectified linear unit (ReLU) are used after the FGRU layer, composing an FGRU
block. We used 64 hidden dimensions for each forward and backward FGRU cell.

The decoder is composed of a Time-axis Gated Recurrent Unit (TGRU) block and
1D Transposed Convolutional Neural Network (1D-TrCNN) blocks. The output of
the encoder is passed into a uni-directional GRU layer to aggregate the
information along the time-axis. We call this GRU layer a TGRU layer. While
one can apply different GRU cells to each frequency-axis index of the encoder
output, we shared the same cell on each frequency-axis index to save the
number of parameters. A pointwise convolution, BN, and ReLU follow the TGRU
layer, composing a TGRU block. We used 128 hidden dimensions for the TGRU
cell. Finally, 1D-TrCNN blocks are used to upsample the output from the TGRU
block to the original spectrogram size. The 1D-TrCNN block takes two inputs -
1. a previous layer output, 2. a skipped tensor from the encoder at the same
hierarchy - and upsamples them as follows. First, the two inputs are
concatenated and projected to a smaller channel size (192 $\rightarrow$ 64)
using a pointwise convolution. Then, 1D transposed convolution is used to
upsample the compressed information. This procedure saves both the number of
parameters and computation compared to the usual U-Net implementation where
the two inputs are concatenated and upsampled immediately using the transposed
convolution operation. Note that we did not use depthwise convolution for
1D-TrCNN block as we empirically observed that it drops the performance
significantly when used in the decoding stage.

Every convolution operation used in the encoder and decoder is followed by BN
and ReLU. We denote the convolution configurations as follows, $l$-th:
($\kappa$, $s$, $c$) , where $l$, $\kappa$, $s$, $c$ denotes layer index,
kernel size, strides, and output channels, respectively. The detailed
configurations of the encoder and decoder are as follows, EncoderConfig =
{1-th: (5,2,64), 2-th: (3,1,128), 3-th: (5,2,128), 4-th: (3,1,128), 5-th:
(5,2,128), 6-th: (3,2,128)}, DecoderConfig = {1-th: (3,2,64), 2-th: (5,2,64),
3-th: (3,1,64), 4-th: (5,2,64), 5-th: (3,1,64), 6-th: (5,2,10)}. Note that the
pointwise convolution operations share the same output channel configuration
with the exception that $\kappa$ and $s$ are both 1. The overview of TRU-Net
and the number of parameters used for 1D-CNN blocks, FGRU block, TGRU block,
and 1D-TrCNN blocks are shown in Fig. 1.

##  3 Single-stage Denoising and Dereverberation

A noisy-reverberant mixture signal $\bm{x}$ is commonly modeled as the sum of
additive noise $\bm{y}^{(n)}$ and reverberant source $\tilde{\bm{y}}$, where
$\tilde{\bm{y}}$ is a result of convolution between room impulse response
(RIR) $\bm{h}$ and dry source $\bm{y}$ as follows,

$\bm{x}=\tilde{\bm{y}}+\bm{y}^{(n)}=\bm{h}\circledast\bm{y}+\bm{y}^{(n)}$ (1)

More concretely, we can break down $\bm{h}$ into two parts. First, the direct
path part $\bm{h}^{(d)}$, which does not include the reflection path, and
second, the rest of the part $\bm{h}^{(r)}$ including all the reflection paths
as follows,

$\bm{x}=\bm{h}^{(d)}\circledast\bm{y}+\bm{h}^{(r)}\circledast\bm{y}+\bm{y}^{(n)}=\bm{y}^{(d)}+\bm{y}^{(r)}+\bm{y}^{(n)},$
(2)

where $\bm{y}^{(d)}$ and $\bm{y}^{(r)}$ denotes a direct path source and
reverberation, respectively. In this setting, our goal is to separate $\bm{x}$
into three elements $\bm{y}^{(d)}$, $\bm{y}^{(r)}$, and $\bm{y}^{(n)}$. Each
of the corresponding time-frequency $(t,f)$ representations computed by short-
time Fourier transform (STFT) is denoted as $X_{t,f}\in\mathbb{C}$,
$Y^{(d)}_{t,f}\in\mathbb{C}$, $Y^{(r)}_{t,f}\in\mathbb{C}$,
$Y^{(n)}_{t,f}\in\mathbb{C}$, and the estimated values will be denoted by the
hat operator $\hat{\,\cdot\,}$.

###  3.1 Phase-aware $\beta$-sigmoid mask

The proposed phase-aware $\beta$-sigmoid mask (PHM) is a complex-valued mask
which is capable of systemically restricting the sum of estimated complex
values to be exactly the value of mixture, $X_{t,f}=Y^{(k)}_{t,f}+Y^{(\lnot
k)}_{t,f}$. The PHM separates the mixture $X_{t,f}$ in STFT domain into two
parts as one-vs-rest approach, that is, the signal $Y^{(k)}_{t,f}$ and the sum
of the rest of the signals $Y^{(\lnot k)}_{t,f}=X_{t,f}-Y^{(k)}_{t,f}$, where
index $k$ could be one of the direct path source ($d$), reverberation ($r$),
and noise ($n$) in our setting, $k\in\\{d,r,n\\}$. The complex-valued mask
$M^{(k)}_{t,f}\in\mathbb{C}$ estimates the magnitude and phase value of the
source of interest $k$.

Computing PHM requires two steps. First, the network outputs the magnitude
part of two masks $\lvert M^{(k)}_{t,f}\rvert$ and $\lvert M^{(\lnot
k)}_{t,f}\rvert$ with sigmoid function $\sigma^{(k)}(\bm{z}_{t,f})$ multiplied
by coefficient $\beta_{t,f}$ as follows, $\lvert
M^{(k)}_{t,f}\rvert=\beta_{t,f}\cdot\sigma^{(k)}(\bm{z}_{t,f})=\beta_{t,f}\cdot(1+e^{-(z^{(k)}_{t,f}-z^{(\lnot
k)}_{t,f})})^{-1}$, where $z^{(k)}_{t,f}$ is the output located at $(t,f)$
from the last layer of neural-network function $\psi^{(k)}(\phi)$, and $\phi$
is a function composed of network layers before the last layer. $\lvert
M^{(k)}_{t,f}\rvert$ serves as a magnitude mask to estimate source $k$ and its
value ranges from 0 to $\beta_{t,f}$. The role of $\beta_{t,f}$ is to design a
mask that is close to optimal values with a flexible magnitude range so that
the values are not bounded between 0 and 1, unlike the commonly used sigmoid
mask. In addition, because the sum of the complex valued masks $M^{(k)}_{t,f}$
and $M^{(\lnot k)}_{t,f}$ must compose a triangle, it is reasonable to design
a mask that satisfies the triangle inequalities, that is, $\lvert
M^{(k)}_{t,f}\rvert+\lvert M^{(\lnot k)}_{t,f}\rvert$ $\geq 1$ and
$\mathinner{\\!\left\lvert\lvert M^{(k)}_{t,f}\rvert-\lvert M^{(\lnot
k)}_{t,f}\rvert\right\rvert}\leq 1$. To address the first inequality we
designed the network to output $\beta_{t,f}$ from the last layer with a
softplus activation function as follows,
$\beta_{t,f}=1+\texttt{softplus}((\psi_{\beta}(\phi))_{t,f})$, where
$\psi_{\beta}$ denotes an additional network layer to output $\beta_{t,f}$.
The second inequality can be satisfied by clipping the upper bound of the
$\beta_{t,f}$ by
$1\mathbin{/}\lvert\,\sigma^{(k)}(\bm{z}_{t,f})-\sigma^{(\lnot
k)}(\bm{z}_{t,f})\rvert$.

Once the magnitude masks are decided, we can construct a phase mask
$e^{j\theta_{t,f}^{(k)}}$. Given the magnitudes as three sides of a triangle,
we can compute the cosine of the absolute phase difference
$\Delta\theta_{t,f}^{(k)}$ between the mixture and source $k$ as follows,
$\cos(\Delta\theta_{t,f}^{(k)})=\nicefrac{{(1+{\lvert
M_{t,f}^{(k)}\rvert}^{2}-{\lvert M_{t,f}^{(\lnot k)}\rvert}^{2})}}{{(2\,\lvert
M_{t,f}^{(k)}\rvert)}}$. Then, the rotational direction
$\xi_{t,f}\in\\{1,-1\\}$ (clockwise or counterclockwise) for phase correction
is estimated for the phase mask as follows,
$e^{j\theta_{t,f}^{(k)}}=\cos(\Delta\theta_{t,f}^{(k)})+j\xi_{t,f}\sin(\Delta\theta_{t,f}^{(k)})$.
Two-class straight-through Gumbel-softmax estimator was used to estimate
$\xi_{t,f}$ [11]. $M^{(k)}_{t,f}$ is defined as follows, $M^{(k)}_{t,f}=\lvert
M^{(k)}_{t,f}\rvert\cdot e^{j\theta_{t,f}^{(k)}}$. Finally, $M^{(k)}_{t,f}$ is
multiplied with $X_{t,f}$ to estimate the source $k$ as follows,
$\hat{Y}^{(k)}_{t,f}=M^{(k)}_{t,f}\cdot X_{t,f}$.

###  3.2 Masking from the perspective of a quadrilateral

Fig. 2: The illustration of PHM masking method on a quadrilateral

Because we desire to extract both the direct and reverberant source, two pairs
of PHMs are used for each of them. The first pair of masks, $M^{(d)}_{t,f}$
and $M^{(\lnot d)}_{t,f}$, separate the mixture into the direct source and the
rest of the components, respectively. The second pair of masks,
$M^{(n)}_{t,f}$ and $M^{(\lnot n)}_{t,f}$, separate the mixture into the noise
and the reverberant source. Since PHM guarantees the mixture and separated
components to construct a triangle in the complex STFT domain, the outcome of
the separation can be seen from the perspective of a quadrilateral, as shown
in Fig 2. In this setting, as the three sides and two side angles are already
determined by the two pairs of PHMs, the fourth side of the quadrilateral,
$M^{(r)}_{t,f}$, is uniquely decided.

###  3.3 Multi-scale objective

Recently, a multi-scale spectrogram (MSS) loss function has been successfully
used in a few audio synthesis studies [12, 13]. We incorporate this multi-
scale scheme not only in the spectral domain but also in the waveform domain
similar to [14].

Learning to maximize cosine similarity can be regarded as maximizing the
signal-to-distortion ratio (SDR) [2]. Cosine similarity loss $C$ between the
estimated signal $\hat{\bm{y}}^{(k)}\in\mathbb{R}^{N}$ and the ground truth
signal $\bm{y}^{(k)}\in\mathbb{R}^{N}$ is defined as follows,
$C(\bm{y}^{(k)},\hat{\bm{y}}^{(k)})=-\frac{<\bm{y}^{(k)},\hat{\bm{y}}^{(k)}>}{\mathinner{\\!\left\lVert\bm{y}^{(k)}\right\rVert}\mathinner{\\!\left\lVert\hat{\bm{y}}^{(k)}\right\rVert}}$,
where $N$ denotes the temporal dimensionality of a signal and $k$ denotes the
type of signal ($k\in\\{d,r,n\\}$). Consider a sliced signal
$\bm{y}^{(k)}_{[\frac{N}{M}(i-1)\mathrel{\mathop{\ordinarycolon}}\frac{N}{M}i]}$,
where $i$ denotes the segment index and $M$ denotes the number of segments. By
slicing the signal and normalizing it by its norm, each sliced segment is
considered a unit for computing $C$. Therefore, we hypothesize that it is
important to choose a proper segment length unit $\frac{N}{M}$ when computing
$C$. In our case, we used multiple settings of segment lengths
$g_{j}=\frac{N}{M_{j}}$ as follows,

$\mathcal{L}_{wav}^{(k)}=\sum_{j}\frac{1}{M_{j}}\sum_{i=1}^{M_{j}}C(\bm{y}^{(k)}_{[g_{j}(i-1)\mathrel{\mathop{\ordinarycolon}}g_{j}i]},\hat{\bm{y}}^{(k)}_{[g_{j}(i-1)\mathrel{\mathop{\ordinarycolon}}g_{j}i]}),$
(3)

where $M_{j}$ denotes the number of sliced segments. In our case, the set of
$g_{j}$'s was chosen as follows, $g_{j}\in\\{4064,2032,1016,508\\}$.

Next, the multi-scale loss on spectral domain is defined as follows,

$\mathcal{L}_{spec}^{(k)}=\sum_{i}\mathinner{\\!\left\lVert\,{\lvert\text{STFT}_{i}(\bm{y}^{(k)})\rvert}^{0.3}-{\lvert\text{STFT}_{i}(\hat{\bm{y}}^{(k)})}^{0.3}\rvert\right\rVert}^{2},$
(4)

where $i$ denotes the FFT size of $\text{STFT}_{i}$. The only difference to
the original MSS loss is that we replaced the log transformation into the
power-law compression, as it has been successfully used in previous speech
enhancement studies [15, 16]. We used the FFT sizes of STFT, (1024, 512, 256),
with 75% overlap. The final loss function is defined by adding all the
components as follows,
$\mathcal{L}_{\text{final}}=\sum_{k\in\\{d,r,n\\}}\mathcal{L}_{wav}^{(k)}+\mathcal{L}_{spec}^{(k)}$.

##  4 Experiments

|  |  | Synthetic without Reverb | Synthetic with Reverb
---|---|---|---|---
Methods | Size(M/MB) | RT | PESQ1 | PESQ2 | CBAK | COVL | CSIG | SI-SDR | STOI | PESQ1 | PESQ2 | CBAK | COVL | CSIG | SI-SDR | STOI
Noisy | - | - | 2.45 | 1.58 | 2.53 | 2.35 | 3.19 | 9.07 | 91.52 | 2.75 | 1.82 | 2.80 | 2.64 | 3.50 | 9.03 | 86.62
NSnet [17] | 1.27/4.84 | ✓ | 2.68 | 1.81 | 2.00 | 2.24 | 2.78 | 12.47 | 90.56 | 2.45 | 1.52 | 1.94 | 1.95 | 2.52 | 9.18 | 82.15
DTLN [18] | 0.99/3.78 | ✓ | 3.04 | - | - | - | - | 16.34 | 94.76 | 2.70 | - | - | - | - | 10.53 | 84.68
ConvTasNet [19] | 5.08/19.38 | ✗ | - | 2.73 | 3.64 | 3.41 | 4.07 | - | - | - | 2.71 | 3.67 | 3.47 | 4.21 | - | -
PoCoNet1 [3] | 50/190.73 | ✗ | - | 2.71 | 3.02 | 3.29 | 3.85 | - | - | - | 2.83 | 3.21 | 3.35 | 3.83 | - | -
PoCoNet2 [3] | 50/190.73 | ✗ | - | 2.75 | 3.04 | 3.42 | 4.08 | - | - | - | - | - | - | - | - | -
DCCRN-E [4] | 3.7/14.11 | ✓ | 3.27 | - | - | - | - | - | - | 3.08 | - | - | - | - | - | -
DCCRN-CL [4] | 3.7/14.11 | ✗ | 3.26 | - | - | - | - | - | - | 3.10 | - | - | - | - | - | -
TRU-Net (FP32) | 0.38/1.45 | ✓ | 3.36 | 2.86 | 3.66 | 3.55 | 4.21 | 17.55 | 96.32 | 3.35 | 2.74 | 3.62 | 3.48 | 4.17 | 14.87 | 91.29
TRU-Net (INT8) | 0.38/0.36 | ✓ | 3.35 | 2.84 | 3.62 | 3.53 | 4.18 | 17.23 | 96.12 | 3.31 | 2.70 | 3.56 | 3.45 | 4.16 | 14.47 | 91.01

Table 1: Objective evaluation results on DNS-challenge synthetic development
sets. PoCoNet2 denotes the model with partial dereverberation described in
[3], and PoCoNet1 is the model trained without it. We denote the network size
(Size) in two aspects, the number of parameters in million (M) and the actual
model size in megabyte (MB). The models with real-time (RT) capability are
marked with ✓, otherwise ✗.

###  4.1 Implementation details

Since our goal is to perform both denoising and dereverberation, we used
pyroomacoustics [20] to simulate an artificial reverberation with randomly
sampled absorption, room size, location of source and microphone distance. We
used 2 seconds of speech and noise segments, and mixed them with a uniformly
distributed source-to-noise ratio (SNR) ranging from -5 dB to 25 dB. Input
features were used as a channel-wise concatenation of log-magnitude
spectrogram, PCEN spectrogram, and real/imaginary part of demodulated phase.
We used AdamW optimizer [21] and the learning rate was halved when the
validation score did not improve for three consecutive epochs. The initial
learning rate was set to 0.0004. The window size and hop size were set to 512
(32 ms) and 128 (8 ms), respectively.

We also quantized the proposed model into INT8 format and compared the model
size with prior works. The purpose of our quantized model experiments is to
reduce the model size and computational cost for embedded environments. We
adopted the computation flow using quantized numbers suggested in [5] to
quantize the neural network. In addition, the uniform symmetric quantization
scheme [22], which uses uniform quantization and restricts zero-point to 0,
was applied for efficient hardware implementation. In the experiments, all the
layers in the neural network are processed using quantized weights,
activations, and inputs; only bias values are represented in full precision.
Other processing steps such as feature extraction and masking are computed in
full precision. For encoder and decoder layers, we observe the scale
statistics of intermediate tensors during training. Then, during inference, we
fix the scales of activations using the average of the observed minimum and
maximum values. Only GRU layers are dynamically quantized during the inference
time due to the large dynamic range of internal activations at each time step.

###  4.2 Ablation study

In order to confirm the effect of PCEN, multi-scale objective, and FGRU block,
we trained and validated the model using the CHiME2 training set and
development set, respectively. An ablation study was conducted on the CHiME2
test set. TRU-Net-A denotes the proposed method. TRU-Net-B denotes the model
trained without multi-scale objective. TRU-Net-C denotes the model trained
without the PCEN feature. TRU-Net-D denotes the model trained without FGRU
block. We used the original SDR [23] to compare our model with other models.
The results are shown in Table 2. It is clearly observable that all the
proposed methods are contributing to performance improvement. Note that FGRU
block contributes significantly on the performance. We also compared the
proposed model with other models using the CHiME2 test set. The proposed model
showed better performance than not only the recent lightweight model Tiny-LSTM
(TLSTM) and its pruned version (PTLSTM) [24], but also the large-sized model
[16].

|  | Input SNR
---|---|---
Methods | Size (M/MB) | -6 | -3 | 0 | 3 | 6 | 9 | Avg.
TLSTM (FP32) [24] | 0.97/3.70 | 10.01 | 11.54 | 13.08 | 14.23 | 15.85 | 17.46 | 13.70
PTLSTM (FP32) [24] | 0.52/1.97 | 10.07 | 11.59 | 13.10 | 14.31 | 15.89 | 17.50 | 13.74
PTLSTM (INT8) [24] | 0.61/0.58 | 9.82 | 11.37 | 12.91 | 14.20 | 15.74 | 17.44 | 13.58
PTLSTM (INT8) [24] | 0.33/0.31 | 9.33 | 10.91 | 12.46 | 13.79 | 15.46 | 17.16 | 13.18
Wilson et al. [16] | 65/247.96 | 12.17 | 13.44 | 14.70 | 15.83 | 17.30 | 18.78 | 15.37
TRU-Net-A (FP32) | 0.38/1.45 | 12.36 | 13.62 | 15.08 | 16.21 | 17.70 | 19.39 | 15.73
TRU-Net-B (FP32) | 0.38/1.45 | 12.21 | 13.39 | 14.91 | 16.09 | 17.53 | 19.24 | 15.56
TRU-Net-C (FP32) | 0.38/1.45 | 11.96 | 13.24 | 14.69 | 15.97 | 17.47 | 19.18 | 15.42
TRU-Net-D (FP32) | 0.31/1.18 | 11.83 | 13.14 | 14.63 | 15.85 | 17.28 | 18.97 | 15.28
TRU-Net-A (INT8) | 0.38/0.36 | 12.35 | 13.62 | 15.03 | 16.18 | 17.62 | 19.30 | 15.68
TRU-Net-B (INT8) | 0.38/0.36 | 12.23 | 13.40 | 14.91 | 16.08 | 17.51 | 19.21 | 15.56
TRU-Net-C (INT8) | 0.38/0.36 | 11.96 | 13.20 | 14.64 | 15.94 | 17.42 | 19.11 | 15.38
TRU-Net-D (INT8) | 0.31/0.30 | 11.79 | 13.13 | 14.56 | 15.78 | 17.19 | 18.85 | 15.22

Table 2: Objective evaluation results on the CHiME2 test set.

###  4.3 Denoising results

We further checked the denoising performance of our model by training the
model on the large scale DNS-challenge dataset [25] and internally collected
dataset. It was tested on two non-blind DNS development sets, 1) synthetic
clips without reverb (Synthetic without Reverb) and 2) synthetic clips with
reverb (Synthetic with Reverb). We compared our model with the recent models
[3, 4, 17, 18, 19] submitted to the previous 2020 Interspeech DNS-challenge. 6
evaluation metrics, PESQ, CBAK, COVL, CSIG, SI-SDR, and STOI [26, 27, 28, 29],
were used. Note that although it is recommended to use ITU-T P862.2 wide-band
version of PESQ (PESQ2), a few studies reported their score using ITU-T P862.1
(PESQ1). Therefore, we used both PESQ versions to compare our model with other
models. The results are shown in Table 1. We can see that TRU-Net shows the
best performance in the Synthetic without Reverb set while having the smallest
number of parameters. In the Synthetic with Reverb set, TRU-Net showed
competitive performance using orders of magnitude fewer parameters than other
models.

###  4.4 Dereverberation results

The performance of simultaneous denoising and dereverberation was tested on
min subset of WHAMR dataset, which contains 3,000 audio files. The WHAMR
dataset is composed of noisy-reverberant mixtures and the direct sources as
ground truth. TRU-Net models (FP32 and INT8) in Table 1 were used for the
test. We show the denoising and dereverberation performance of our model in
Table 3 along with two other models that were tested on the same WHAMR
dataset. Our model achieved the best results compared to the other baseline
models, which shows the parameter efficiency of TRU-Net on simultaneous
denoising and dereverberation task.

Method | Size (M/MB) | PESQ1 | SI-SDR | STOI
---|---|---|---|---
Noisy | - | 1.83 | -2.73 | 73.00
NSnet [17] | 1.27/4.84 | 1.91 | 0.34 | 73.02
DTLN [18] | 0.99/3.78 | 2.23 | 2.12 | 80.40
TRU-Net (FP32) | 0.38/1.45 | 2.51 | 3.51 | 81.22
TRU-Net (INT8) | 0.38/0.36 | 2.49 | 3.03 | 80.56

Table 3:  Objective evaluation of simultaneous denoising and dereverberation
results on the WHAMR dataset.

###  4.5 Listening test results

Using the proposed model (TRU-Net (FP32)) in Table 1, we participated in 2021
ICASSP DNS Challenge Track 1 [25]. For better perceptual quality, we mixed the
estimated direct source and reverberant source at 15 dB, and applied a zero-
delay dynamic range compression (DRC). The average computation time to process
a single frame (including FFT, iFFT, and DRC) took 1.97 ms and 1.3 ms on 2.7
GHz Intel i5-5257U and 2.6 GHz Intel i7-6700HQ CPUs, respectively. The
lookahead of TRU-Net is 0 ms. The listening test was conducted based on ITU-T
P.808. The results are shown in Table 4. The model was tested on various
speech sets including singing voice, tonal language, non-English (includes
tonal), English, and emotional speech. The results show that TRU-Net can
achieve better performance than the baseline model, NSnet2 [30].

Method | Size (M/MB) | Singing | Tonal | Non-English | English | Emotional | Overall
---|---|---|---|---|---|---|---
Noisy | - | 2.96 | 3.00 | 2.96 | 2.80 | 2.67 | 2.86
NSnet2 [30] | 2.8/10.68 | 3.10 | 3.25 | 3.28 | 3.30 | 2.88 | 3.21
TRU-Net | 0.38/1.45 | 3.08 | 3.38 | 3.43 | 3.41 | 2.88 | 3.32

Table 4:  MOS results on the DNS-challenge blind test set

##  5 Relation to prior works

Recently, there has been increasing interest in phase-aware speech enhancement
because of the sub-optimality of reusing the phase of the mixture signal.
While most of these works tried to estimate the clean phase by using a phase
mask or an additional network, the absolute phase difference between mixture
and source can be actually computed using the law of cosines [31]. Inspired by
this, [6] proposed to estimate a rotational direction of the absolute phase
difference for speech separation.

The FGRU and TGRU used in TRU-Net is similar to the work in [32]. They used
bidirectional long short-term memory (bi-LSTM) networks on the frequency-axis
and the time-axis combined with 2D-CNN-based U-Net. The difference is that bi-
LSTM was utilized to increase performance in [32], whereas we employ FGRU and
uni-directional TGRU to better handle the online inference scenario combined
with the proposed lightweight 1D-CNN-based (frequency-axis) U-Net.

##  6 Conclusions

In this work, we proposed TRU-Net, which is an efficient neural network
architecture specifically designed for online inference applications. Combined
with the proposed PHM, we successfully demonstrated a single-stage denoising
and dereverberation in real-time. We also showed that using PCEN and multi-
scale objectives improves the performance further. Experimental results
confirm that our model achieve comparable performance with state-of-the-art
models having a significantly larger number of parameters. For future work, we
plan to employ modern pruning techniques on an over-parameterized model to
develop a big-sparse model which may provide better performance than a small-
dense model with the same number of parameters.

## References

  * [1] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in Proc. MICCAI, 2015, pp. 234–241.

  * [2] Hyeong-Seok Choi, Jang-Hyun Kim, Jaesung Huh, Adrian Kim, Jung-Woo Ha, and Kyogu Lee,  “Phase-aware speech enhancement with deep complex u-net,”  arXiv preprint arXiv:1903.03107, 2019.
  * [3] Umut Isik, Ritwik Giri, Neerad Phansalkar, Jean-Marc Valin, Karim Helwani, and Arvindh Krishnaswamy,  “Poconet: Better speech enhancement with frequency-positional embeddings, semi-supervised conversational data, and biased loss,”  in Proc. INTERSPEECH, 2020.
  * [4] Yanxin Hu, Yun Liu, Shubo Lv, Mengtao Xing, Shimin Zhang, Yihui Fu, Jian Wu, Bihong Zhang, and Lei Xie,  “Dccrn: Deep complex convolution recurrent network for phase-aware speech enhancement,”  in Proc. INTERSPEECH, 2020.
  * [5] Benoit Jacob, Skirmantas Kligys, Bo Chen, Menglong Zhu, Matthew Tang, Andrew Howard, Hartwig Adam, and Dmitry Kalenichenko,  “Quantization and training of neural networks for efficient integer-arithmetic-only inference,”  in Proc. CVPR, 2018, pp. 2704–2713.
  * [6] Zhong-Qiu Wang, Ke Tan, and DeLiang Wang,  “Deep learning based phase reconstruction for speaker separation: A trigonometric perspective,”  in Proc. ICASSP, 2019, pp. 71–75.
  * [7] Yuxuan Wang, Pascal Getreuer, Thad Hughes, Richard F Lyon, and Rif A Saurous,  “Trainable frontend for robust and far-field keyword spotting,”  in Proc. ICASSP, 2017, pp. 5670–5674.
  * [8] Vincent Lostanlen, Justin Salamon, Mark Cartwright, Brian McFee, Andrew Farnsworth, Steve Kelling, and Juan Pablo Bello,  “Per-channel energy normalization: Why and how,”  IEEE Signal Processing Letters, vol. 26, no. 1, pp. 39–43, 2018\.
  * [9] Andrew G Howard, Menglong Zhu, Bo Chen, Dmitry Kalenichenko, Weijun Wang, Tobias Weyand, Marco Andreetto, and Hartwig Adam,  “Mobilenets: Efficient convolutional neural networks for mobile vision applications,”  arXiv preprint arXiv:1704.04861, 2017.
  * [10] Kyunghyun Cho, Bart van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshuas Bengio,  “Learning phrase representations using RNN encoder–decoder for statistical machine translation,”  in Proc. EMNLP, 2014, pp. 1724–1734.
  * [11] Eric Jang, Shixiang Gu, and Ben Poole,  “Categorical reparameterization with gumbel-softmax,”  in Proc. ICLR, 2017.
  * [12] Xin Wang, Shinji Takaki, and Junichi Yamagishi,  “Neural source-filter-based waveform model for statistical parametric speech synthesis,”  in Proc. ICASSP, 2019, pp. 5916–5920.
  * [13] Jesse Engel, Lamtharn (Hanoi) Hantrakul, Chenjie Gu, and Adam Roberts,  “Ddsp: Differentiable digital signal processing,”  in Proc. ICLR, 2020.
  * [14] Jian Yao and Ahmad Al-Dahle,  “Coarse-to-Fine Optimization for Speech Enhancement,”  in Proc. INTERSPEECH, 2019, pp. 2743–2747.
  * [15] Hakan Erdogan and Takuya Yoshioka,  “Investigations on data augmentation and loss functions for deep learning based speech-background separation.,”  in INTERSPEECH, 2018, pp. 3499–3503.
  * [16] Kevin Wilson, Michael Chinen, Jeremy Thorpe, Brian Patton, John Hershey, Rif A Saurous, Jan Skoglund, and Richard F Lyon,  “Exploring tradeoffs in models for low-latency speech enhancement,”  in IWAENC, 2018, pp. 366–370.
  * [17] Yangyang Xia, Sebastian Braun, Chandan KA Reddy, Harishchandra Dubey, Ross Cutler, and Ivan Tashev,  “Weighted speech distortion losses for neural-network-based real-time speech enhancement,”  in Proc. ICASSP, 2020, pp. 871–875.
  * [18] Nils L Westhausen and Bernd T Meyer,  “Dual-signal transformation lstm network for real-time noise suppression,”  in Proc. INTERSPEECH, 2020.
  * [19] Yuichiro Koyama, Tyler Vuong, Stefan Uhlich, and Bhiksha Raj,  “Exploring the best loss function for dnn-based low-latency speech enhancement with temporal convolutional networks,”  arXiv preprint arXiv:2005.11611, 2020.
  * [20] Robin Scheibler, Eric Bezzam, and Ivan Dokmanić,  “Pyroomacoustics: A python package for audio room simulation and array processing algorithms,”  in Proc. ICASSP, 2018, pp. 351–355.
  * [21] Sashank J. Reddi, Satyen Kale, and Sanjiv Kumar,  “On the convergence of adam and beyond,”  in Proc. ICLR, 2018.
  * [22] Raghuraman Krishnamoorthi,  “Quantizing deep convolutional networks for efficient inference: A whitepaper,”  arXiv preprint arXiv:1806.08342, 2018.
  * [23] Emmanuel Vincent, Rémi Gribonval, and Cédric Févotte,  “Performance measurement in blind audio source separation,”  IEEE transactions on audio, speech, and language processing, vol. 14, no. 4, pp. 1462–1469, 2006.
  * [24] Igor Fedorov, Marko Stamenovic, Carl Jensen, Li-Chia Yang, Ari Mandell, Yiming Gan, Matthew Mattina, and Paul N Whatmough,  “Tinylstms: Efficient neural speech enhancement for hearing aids,”  in Proc. INTERSPEECH, 2020.
  * [25] Chandan KA Reddy, Harishchandra Dubey, Vishak Gopal, Ross Cutler, Sebastian Braun, Hannes Gamper, Robert Aichner, and Sriram Srinivasan,  “Icassp 2021 deep noise suppression challenge,”  arXiv preprint arXiv:2009.06122, 2020.
  * [26] ITU-T Recommendation,  “Perceptual evaluation of speech quality (pesq): An objective method for end-to-end speech quality assessment of narrow-band telephone networks and speech codecs,”  Rec. ITU-T P. 862, 2001.
  * [27] Philipos C Loizou,  Speech enhancement: theory and practice,  CRC press, 2013.
  * [28] Jonathan Le Roux, Scott Wisdom, Hakan Erdogan, and John R Hershey,  “Sdr–half-baked or well done?,”  in Proc. ICASSP, 2019, pp. 626–630.
  * [29] Cees H Taal, Richard C Hendriks, Richard Heusdens, and Jesper Jensen,  “A short-time objective intelligibility measure for time-frequency weighted noisy speech,”  in Proc. ICASSP, 2010, pp. 4214–4217.
  * [30] Sebastian Braun and Ivan Tashev,  “Data augmentation and loss normalization for deep noise suppression,”  in International Conference on Speech and Computer, 2020, pp. 79–86.
  * [31] Pejman Mowlaee, Rahim Saeidi, and Rainer Martin,  “Phase estimation for signal reconstruction in single-channel source separation,”  in Thirteenth Annual Conference of the International Speech Communication Association, 2012.
  * [32] Tomasz Grzywalski and Szymon Drgas,  “Using recurrences in time and frequency within u-net architecture for speech enhancement,”  in Proc. ICASSP, 2019, pp. 6970–6974.


# TMM-Fast: A Transfer Matrix Computation Package for Multilayer Thin-Film
Optimization

Alexander Luce
Max Planck Institute for the Science of Light
Friedrich-Alexander-Universität
Erlangen-Nürnberg
OSRAM Opto Semiconductors GmbH
Regensburg
<EMAIL_ADDRESS>
&Ali Mahdavi
OSRAM Opto Semiconductors GmbH
Regensburg
Florian Marquardt
Max Planck Institute for the Science of Light
Friedrich-Alexander-Universität
Erlangen-Nürnberg
&Heribert Wankerl
Universität Regensburg
OSRAM Opto Semiconductors GmbH
Regensburg

###### Abstract

Achieving the desired optical response from a multilayer thin-film structure
over a broad range of wavelengths and angles of incidence can be challenging.
An advanced thin-film structure can consist of multiple materials with
different thicknesses and numerous layers. Design and optimization of complex
thin-film structures with multiple variables is a computationally heavy
problem that is still under active research. To enable fast and easy
experimentation with new optimization techniques, we propose the Python
package TMM-Fast which enables parallelized computation of reflection and
transmission of light at different angles of incidence and wavelengths through
the multilayer thin-film. By decreasing computational time, generating
datasets for machine learning becomes feasible and evolutionary optimization
can be used effectively. Additionally, the sub-package TMM-Torch allows to
directly compute analytical gradients for local optimization by using PyTorch
Autograd functionality. Finally, an OpenAi Gym environment is presented which
allows the user to train reinforcement learning agents on the problem of
finding multilayer thin-film configurations.

_K_ eywords Multilayer Thin-Film $\cdot$ Transfer Matrix Method $\cdot$
Parallelization $\cdot$ Optimization $\cdot$ Machine Learning

##  1 Introduction

Although the existence of a globally optimal multilayer thin-film was proven
[1, 2, 3], the optimization of multilayer thin-films concerning reflectivity
and transmittivity over wavelength and incidence angle remains a challenging
task for the scientific community [4, 5, 6]. As such it gained attention in
many publications [7, 8, 9, 10, 11] including various methods based on
reinforcement learning [12, 13] and machine learning [14, 15]. Naturally, the
need for a standardized, encompassing environment arises to make future
research more trustful, comparable, and consistent. Here, we propose a
comprehensive Python package that provides functionality to researchers to
design and optimize multilayer thin-films. At its core, the proposed TMM-Fast
package (Transfer Matrix Method - Fast) is a parallelized and re-implemented
revision of the existing TMM code that was initially published by S. Byrnes
[16]. It implements Abèles transfer matrix method [17] in Python to calculate
transmission and reflection of an incident plane wave through a slab of
layered thin-film materials with different thicknesses. Given a broad,
discretized spectrum of light irradiating a thin-film under particular
incident angles, the transmission and reflection must be calculated for each
contributing wavelength at each angle of incident. The coefficients of
transmission and reflection for a specific wavelength and angle of incidence
depend on the thicknesses and dispersive and dissipative refractive index of
the layers. Those coefficients can deviate significantly at wavelengths which
differ only slightly for the same multilayer thin-film. An intuitive approach
to increase the computational speed of the response of a multilayer thin-film
is to parallelize the pair-wise independent computations regarding wavelengths
and angles. The parallelization is implemented via matrix operations based on
the open-source numeric algebra package Numpy [18] and the thread-management
package Dask. The parallelization showed to reduce the computational time by
$\sim 100\times$ with an additional factor on the order of available CPU cores
by using Dask. Additionally, by using PyTorch Autograd [19], analytical
gradients of the multilayer thin-film layer thicknesses and refractive indices
can be computed for local optimization. Since the implementation is done via
PyTorch methods, it can be integrated into advanced Neural Networks. Moreover,
the parallelized TMM-Fast code is used to implement an OpenAI Gym environment
111https://github.com/MLResearchAtOSRAM/gym-multilayerthinfilm that can be
used by physicists and reinforcement learning researchers. In the environment,
the optimization of multilayer thin-films is introduced as a sequence
generation process and thus considered as a parameterized Markov decision
process [20], shown in section 6. An overview of the scope of the package is
shown in Figure 1. All code is available on
GitHub222https://github.com/MLResearchAtOSRAM/tmm_fast under MIT license.

Figure 1: Schematic overview of the contents of the TMM-Fast package. The
package consists of three subdivisions. Part a) of Figure 1 contains the core
functionality of the TMM-package, computing the optical response of a
multilayer thin-films quickly via the transfer matrix method over a broad
range of wavelengths and incident angles. Another part of the core
functionality of solving multilayer thin-films is the possibility to generate
huge datasets with >1e5 data samples for machine learning models.
Part b) encompasses an OpenAI Gym Environment which allows easy comparisons of
different reinforcement learning agents. Here, the environment state is given
by the total layer number, layer thicknesses, and material choice. The agent
takes an action that specifies the material and thickness of the layer to
stack next. The environment implements the multilayer thin-film generation as
consecutive conduction of actions and assigns a reward to a proposed multi-
layer thin film based on how close the actual (solid orange line) fulfils a
desired (dashed orange line) characteristic - in this case, a one-dimensional
optical characteristic, e.g. reflectivity over wavelength for perpendicular
incidence. The experience accumulated by the taken actions of the agent is
used to adapt subsequent actions in order to increase the reward and thus
generate more and more sophisticated multilayer thin-films. An investigation
of a reinforcement learning agent which was trained with the proposed Gym
Environment is given by Wankerl et. al. [13].
Part c) encompasses an implementation of the transfer matrix method via
PyTorch functions which allow for backpropagation through the entire
computation. This enables easy differentiation of the parameters of interest
and allows gradient based and quasi-newton optimization methods to be used for
optimization of multilayer thin-films. In the example shown, a random thin-
film is optimized with respect to the given optimization target, weighted by
the mean squared error. Gradients via autograd are used with the BFGS
algorithm [21]. The iterative evolution of the layer thicknesses and the
correspondingly decreasing figure of merit are shown below. For the given
initial values the optimization converged in about 25 iteration steps.

##  2 Physical Background of the Transfer Matrix Method (TMM)

Light of a particular wavelength passing from one into another material
experiences a sudden change of the refractive index $n_{1}$ to $n_{2}$ that
results in reflection and transmission of the incoming wave. A light wave with
a wavevector k which is not parallel to the surface normal experiences
refraction. The ratio of angle of incidence $\theta_{1}$ and angle of
refraction $\theta_{2}$ is given by Snell’s law
$\frac{\sin{\theta_{1}}}{\sin{\theta_{2}}}=\frac{n_{2}}{n_{1}}$. Given the
Fresnel equations, reflection and transmission coefficients $r$ and $t$ can be
computed, where a distinction is made between s and p polarization by using a
subscript

$\displaystyle
t_{s}=\frac{2n_{1}\cos{\theta_{1}}}{n_{1}\cos{\theta_{1}}+\frac{\mu_{r1}}{\mu_{r2}}n_{2}\cos{\theta_{2}}}\hskip
14.22636ptr_{s}=\frac{n_{1}\cos{\theta_{1}}-\frac{\mu_{r1}}{\mu_{r2}}n_{2}\cos{\theta_{2}}}{n_{1}\cos{\theta_{1}}+\frac{\mu_{r1}}{\mu_{r2}}n_{2}\cos{\theta_{2}}}$
(1) $\displaystyle
t_{p}=\frac{2n_{1}\cos{\theta_{1}}}{\frac{\mu_{r1}}{\mu_{r2}}n_{2}\cos{\theta_{1}}+n_{1}\cos{\theta_{2}}}\hskip
14.22636ptr_{p}=\frac{\frac{\mu_{r1}}{\mu_{r2}}n_{2}\cos{\theta_{1}}-n_{1}\cos{\theta_{2}}}{\frac{\mu_{r1}}{\mu_{r2}}n_{2}\cos{\theta_{1}}+n_{1}\cos{\theta_{2}}}$
(2)

and $\mu_{r}$ being the respective magnetic permeability. Note that for
$\theta_{1}$ = 0∘, i.e. vertical incidence, $r_{s}=r_{p}$ and $t_{s}=t_{p}$.
Based on the reflection and transmission coefficients, the reflectivity and
transmittivity R and T can be computed as shown in Equation 3.

$\displaystyle R_{\lambda,\theta_{1}}=r_{i}^{2};\hskip
14.22636ptT_{\lambda,\theta_{1}}=\frac{n_{2}\cos{\theta_{2}}}{n_{1}\cos{\theta_{1}}}t_{i}^{2}$
(3)

where the subscript $i$ indicates the polarization. Note that, for
absorptionless materials $T+R=1$ holds.

Now, consider a multilayer thin-film with $L\in\mathds{N}$ layers and the
individual layers are denoted by $l\leq L$. The light enters from an injection
layer of semi-infinite thickness $l=0$ with a relative amplitude of 1 and
exits the multilayer thin-film in the outcoupling layer of semi-infinite
thickness with $l=L+1$. From the outcoupling layer, no light enters the thin-
film. The transmitted part of the light in layer $l$, that travels in
"forward" direction, ie. towards the layer with $l=l_{\text{current}\,+\,1}$
is denoted by $v_{l}$ and the reflected part that travels in "backward"
direction is given by $w_{l}$. By using the reflection and transmission
coefficients r and t, the response of any layer can be written as

$\displaystyle\begin{pmatrix}v_{l}\\\
w_{l}\end{pmatrix}=\begin{pmatrix}e^{-i\delta_{l}}&0\\\
0&e^{i\delta_{l}}\end{pmatrix}\begin{pmatrix}1&r_{l,l+1}\\\
r_{l,l+1}&1\end{pmatrix}\frac{1}{t_{l,l+1}}\begin{pmatrix}v_{l+1}\\\
w_{l+1}\end{pmatrix}=M_{l}\begin{pmatrix}v_{l+1}\\\ w_{l+1}\end{pmatrix}.$ (4)

where $\delta=d_{l}\,k_{z}$ is the accumulated phase of the light wave when
travelling through a layer with a specific thickness $d_{l}$ and with wave
vector $k_{l}$.

Eventually, the total characteristic matrix of the multilayer thin-film is
given by

$\displaystyle\tilde{M}=\prod_{i=0}^{L-1}M_{i}.$ (5)

Finally, to compute the reflection and transmission of the entire multilayer
thin-film, one needs to evaluate

$\displaystyle\begin{pmatrix}1\\\ r\end{pmatrix}=\tilde{M}\begin{pmatrix}t\\\
0\end{pmatrix}.$ (6)

The transmission and reflection coefficients are separated

$\displaystyle
r_{\lambda,\vartheta}=\frac{\tilde{M}_{10}}{\tilde{M}_{00}}\hskip
28.45274ptt_{\lambda,\vartheta}=\frac{1}{\tilde{M}_{00}}.$ (7)

and allow to compute the reflectivity and transmittivity via Equation 3. Note
that one can easily calculate the partial transmission and reflection
coefficients by only multiplying Equation 5 up to $L-l$. The initial angle of
incidence in the first layer is given by $\vartheta$.

##  3 Implementation in Numpy

The key contribution of TMM-Fast is the parallelized handling of the
characteristic matrix that reduces computational time. The matrix $M_{l}$
consists of three separate matrices, matrix $A$ which encompasses the
accumulated phase, and the two matrices holding the coefficients of reflection
and transmission, respectively. They are of shape
$\left[N_{\lambda},N_{\vartheta},L,2,2\right]$, where $N_{\lambda}$ and
$N_{\theta}$ are the number of wavelengths and incident angles, respectively.
To get the characteristic matrix $M_{l}$, Numpy’s einsum333Documentation:
https://numpy.org/doc/stable/reference/generated/numpy.einsum.html method
allows to specify multiplication and contractions of different dimensions
easily:

⬇

1M_l = np.zeros((num_lambda, num_angles, num_layers, 2, 2),

2 dtype=complex)

3F = r_list[:, 1:]

4M_l[:,:,1:-1,0,0] = np.einsum(’hji,ji->jhi’, 1/A, 1/t_list[:, 1:])

5M_l[:,:,1:-1,0,1] = np.einsum(’hji,ji->jhi’, 1/A, F/t_list[:, 1:])

6M_l[:,:,1:-1,1,0] = np.einsum(’hji,ji->jhi’, A, F/t_list[:, 1:])

7M_l[:,:,1:-1,1,1] = np.einsum(’hji,ji->jhi’, A, 1/t_list[:, 1:])

8

9Mtilde = np.empty((num_angles, num_lambda, 2, 2), dtype=complex)

10Mtilde[:, :] = make_2x2_array(1, 0, 0, 1, dtype=complex)

11for i in range(1, num_layers-1):

12 Mtilde = np.einsum(’ijkl,ijlm->ijkm’, Mtilde, M_l[:,:,i])

Finally $\tilde{M}$ is computed by multiplying out the thickness dimensions.

The entire function is called coh_tmm_fast or by coh_tmm_fast_disp. Both
methods differ since the former assumes dispersionless materials whereas the
latter accepts dispersive materials. An Example is give in Appendix A.1.

##  4 Core functionality speedup and dataset generation

To verify the speedup that the TMM-Fast package provides compared to the
native implementation of Byrnes, ten multilayer thin-films are generated with
21 layers are generated and evaluated in a spectral range from 400 nm to 700
nm at 100 equally spaced points. The angles of incidence range from 0° to 90°
based on 20 equally spaced supporting points. The original TMM method requires
a computation time of $17.5\pm 0.131$ s for the evaluation while our proposed
method computes the coefficients of the multilayer thin-film in $0.168\pm
0.002$ s which corresponds to an acceleration of $\sim$ 100$\times$.

Finally, the Python package Dask444https://dask.org/ manages parallel threads
of the CPU to distribute the computation of different independent computations
on all available CPU cores. The application is straightforward: by calling the
coh_tmm_fast function implicitly inside the delayed() function of Dask, a list
of all necessary computations is created. Then, the entire list is executed
implicitly to create a computational graph for Dask, which is required to
orchestrate the parallel threads efficiently. Lastly, the compute() method
triggers the actual computation and the result is returned. By running
coh_tmm_fast on all available threads, an additional speedup of the dataset
creation on the order of the number of available CPU cores is possible.
However, the benefit might decrease for very large computational clusters with
many cores since the management of the parallel threads creates computational
overhang. By calling multithread_coh_tmm from the TMM-Fast package, the
computation is easily started. A Dataset with one million samples of 9-layer
multilayer thin-films at 100 wavelengths and 10 angles of incidence can be
created in approximately 40 min on an eight-core machine. Computing more
layers, wavelengths points, and angles can increase the computational time
significantly. An example is shown in the Appendix subsection A.2.

##  5 TMM-Torch: Multilayer thin-film gradients via Autograd

For optimization, computing gradients enable gradient-based optimization
algorithms such as gradient descent [21]. These gradient based optimization
methods generally converge to local minima with fewer iterations than other
non-gradient based algorithms such as the Nelder-Mead downhill simplex
algorithm [22]. The Python package PyTorch555https://PyTorch.org/ implements
matrix multiplication methods which allow parameters to be automatically
differentiated (Autograd) via the chain rule [19]. Autograd enables the user
to compute the gradients of the input parameters without the necessity of
deriving the gradient analytically and can be dynamically adapted to the
problem. By using the TMM-Torch subpackage of TMM-Fast, the gradients for a
multilayer thin-film can be readily computed. An example is shown in the
Appendix A.3. The application of the transfer matrix takes slightly longer by
using the PyTorch routines. Therefore, using the regular TMM-Fast algorithms
is still advised to reduce computational time if gradients are not necessary.

##  6 Environment for reinforcement learning

Reinforcement learning [23] is an area of machine learning concerned with how
intelligent agents ought to take actions in an environment in order to
maximize a notion of reward. The proposed code implements such an environment,
where agents can stack, characterize and optimize multilayer thin-films.
Therefore, the generation of multilayer thin-films is considered as a
parameterized Markov decision processes [20] and thereby implemented as a
sequence generation process: Beginning from $l=1$, an agent subsequently
executes parameterized actions $a_{l}=\left(d_{l},m_{l}\right)$ that specify
the thickness $d_{l}$ and material index $m_{l}$ of the $l$-th layer. These
actions determine which material of which thickness to stack next, thereby
consecutively forming a multilayer thin-film as illustrated in Figure 1. The
stacked $l$ layers and - optional - the optical characteristics of these
intermediate multilayer thin films are provided to the agent as the
environmental state $s$. Given this state, the agent takes the next action
$a_{l+1}$ until the pre-defined maximum number of layers is reached or the
agent decides to terminate stacking. The optical characteristic of the final
proposed multilayer thin-film of $L$ layers, e.g. regarding reflectivity
$R_{\lambda,\vartheta}\left(\mathbf{I},\mathbf{d}\right)$ over wavelength
$\lambda$ and angle $\vartheta$ of incidence, is computed via the proposed
transfer-matrix method (TMM-Fast). Here, $\mathbf{I}\in\mathbb{C}^{L\times
N_{\lambda}}$ refers to the matrix of (dispersive and dissipative) refractive
indices of the material of each layer. Each material is identified by the
material identifier index $m_{l}$. $\mathbf{d}\in\mathbb{R}^{L}$ denotes the
vector of layer thicknesses. When the agent starts to stack a thin-film,
$\mathbf{d}$ and $\mathbf{I}$ are initialized with zeros and get filled
according to the taken actions with the layer thicknesses and (dispersive)
refractive indices, respectively. The observed reflectivity is compared to a
user-defined, desired reflectivity $R^{\text{target}}_{\lambda,\vartheta}$, in
order to derive a notion of numeric reward, e.g. an inverted reconstruction
error

$-\sum_{\lambda,\vartheta}|R^{\text{target}}_{\lambda,\vartheta}-R_{\lambda,\vartheta}\left(\mathbf{I},\mathbf{d}\right)|.$

Based on this reward, the agent learns — for example, based on Q-learning [12,
24] — to distinguish between good and bad actions and thus derive an optimal
thin-film design.

Whereas the contained physical methods are well-studied and known for decades,
the contribution of the code available at
https://github.com/MLResearchAtOSRAM/gym-multilayerthinfilm lies in the
implementation of an OpenAI gym-related environment. Here, the intention is to
enable AI researchers without optical expertise to solve the corresponding
parameterized Markov decision processes based on common code in the future.

##  7 Conclusion

In this technical note, the comprehensive Python package TMM-Fast for
multilayer thin-film computation, optimization and reinforcement learning is
presented. At its core, the package is comprised of revised and speedup
transfer matrix method code from the original TMM package[16]. TMM-Fast
enables the user to compute multilayer thin-films with Numpy and PyTorch
methods and gives full control over the data to the user. Since the TMM-Fast
package evaluates multilayer thin-films especially fast, it can be used to
generate datasets for machine learning. The reduced computational time also
enables evolutionary optimization to be executed in a reasonable amount of
time. The TMM-Torch implementation allows the user to compute analytical
gradients via automatic differentiation. Quasi-Newton or gradient based
optimization algorithms such as gradient descent can be used and converge
faster to local minima by using an analytical gradient. Finally, an OpenAi Gym
environment is proposed which allows researchers to easily test and experiment
with new reinforcement agents on solving the multilayer thin-film problem. All
code proposed in this paper is open-source and available at
https://github.com/MLResearchAtOSRAM/ under the MIT-licence.

## References

  * [1] A. V. Tikhonravov and J. A. Dobrowolski.  Quasi-optimal synthesis for antireflection coatings: a new method.  Appl. Opt., 32(22):4265–4275, 1993.
  * [2] Alexander V. Tikhonravov.  Some theoretical aspects of thin-film optics and their applications.  Appl. Opt., 32(28):5417–5426, 1993.
  * [3] Mahdi Ebrahimi and Mohsen Ghasemi.  Design and optimization of thin film polarizer at the wavelength of 1540 nm using differential evolution algorithm.  Optical and Quantum Electronics, 50(192), 2018.
  * [4] Stephan W. Anzengruber, Esther Klann, Ronny Ramlau, and Diana Tonova.  Numerical methods for the design of gradient-index optical coatings.  Appl. Opt., 51(34):8277–8295, 2012.
  * [5] Hans Becker, Diana Tonova, Michael Sundermann, Henrik Ehlers, Stefan Günster, and Detlev Ristau.  Design and realization of advanced multi-index systems.  Appl. Opt., 53(4):A88–A95, 2014.
  * [6] Heather Mary Liddell and H. G. Jerrard.  Computer-Aided Techniques for the Design of Multilayer Filters.  A. Hilger, 1981.
  * [7] C. P. Chang, Y. H. Lee, and S. Y. Wu.  Optimization of a thin-film multilayer design by use of the generalized simulated-annealing method.  Optics Letters, 1990.
  * [8] Wojciech Paszkowicz.  Genetic algorithms, a nature-inspired tool: A survey of applications in materials science and related fields: Part ii.  Materials and Manufacturing Processes, 28, 2013.
  * [9] Chenying Yang, Liang Hong, Weidong Shen, Yueguang Zhang, Xu Liu, and Hongyu Zhen.  Design of reflective color filters with high angular tolerance by particle swarm optimization method.  Opt. Express, 21, 2013.
  * [10] X. Guo, H. Y. Zhou, S. Guo, X. X. Luan, W. K. Cui, Y. F. Ma, , and L. Shi.  Design of broadband omnidirectional antireflection coatings using ant colony algorithm.  Optics Express, 22, 2014.
  * [11] S. Martin, J. Rivory, and M. Schoenauer.  Synthesis of optical multilayer systems using genetic algorithms.  Applied Optics, 34, 1995.
  * [12] Anqing Jiang, Yoshie Osamu, and Liangyao Chen.  Multilayer optical thin film design with deep q learning.  Sci Rep, 10, 2020.
  * [13] Heribert Wankerl, Maike L. Stern, Ali Mahdavi, Christoph Eichler, and Elmar W. Lang.  Parameterized reinforcement learning for optical system optimization.  J. Phys. D: Appl. Phys, 54(30), 2021.
  * [14] Ravi S. Hedge.  Accelerating optics design optimizations with deep learning.  Optical Engineering, 58, 2019.
  * [15] John Roberts and Evan W. Wang.  Modeling and optimization of thin-film optical devices using a variational autoencoder.  Technical report, Stanford University, 2018.
  * [16] Steven J. Byrnes.  Multilayer optical calculations, 2019.
  * [17] F. Abelès.  La théorie générale des couches minces.  Journal De Physique Et Le Radium, 11:307–309, 1950.
  * [18] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [19] Adam Paszke, S. Gross, Soumith Chintala, Gregory Chanan, E. Yang, Zach DeVito, Zeming Lin, Alban Desmaison, L. Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  2017\.
  * [20] Warwick Masson, Pravesh Ranchod, and George Konidaris.  Reinforcement learning with parameterized actions.  In Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, page 1934–1940, 2016.
  * [21] Yann Lecun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  In Proceedings of the IEEE, pages 2278–2324, 1998.
  * [22] J. Nelder and R. Mead.  A simplex method for function minimization.  Comput. J., 7:308–313, 1965.
  * [23] Sutton and Barto.  Introduction to reinforcement learning.  Cambridge: MIT Press, 1998.
  * [24] Christopher John Cornish Hellaby Watkins.  Learning from Delayed Rewards.  PhD thesis, King’s College, Cambridge, UK, May 1989.

##  Appendix A Appendix

###  A.1 TMM-Fast Example

To demonstrate the functionality, a minimal example is given for a multilayer
thin-film with random thicknesses and refractive indices.

⬇

1import tmm_fast as tmmf

2L = 12 # number of layers

3d = np.random.uniform(20, 150, L)*1e-9 # thicknesses of the layers

4d[0] = d[-1] = np.inf # set first and last layer as injection layer

5n = np.random.uniform(1.2, 5, L) # random constant refractive index

6n[-1] = 1 # outcoupling into air

7

8wl = np.linspace(500, 900, 301)*1e-9

9theta = np.deg2rad(np.linspace(0, 90, 301))

10# here s and p polarization is computed and averaged to simulate incoherent
light

11result = (tmmf.coh_tmm_fast(’s’, n, d, theta, wl)[’R’]

12 + tmmf.coh_tmm_fast(’p’, n, d, theta, wl)[’R’])/2

The result of the computation can then be further evaluated. By using the
plot_stacks function from the TMM-Fast package, the multilayer thin-film can
be visualised. An example is shown in Figure 2.

Figure 2: A multilayer thin-film with random layer thicknesses and refractive
index under unpolarized illumination. In this example, the materials are
dispersion- and dissipationless. The injection region, which is below the
thin-film in this depiction, possesses a refractive index of $n=2$. The
outcoupling region above the thin-film possesses a refractive index of $n=1$.
The reflectivity is computed over a wavelength range of 500 nm to 900 nm and
from 0∘ to 90∘ on a 300$\times$300 grid.

###  A.2 Dataset Generation

A dataset containing tens of thousands of datasamples can be created by using
the multithread_coh_tmm function. The external package tqdm is used to display
a progress bar in order to keep track of the generation progress. In this
example, a dataset with 1e6 thin-films at 100 wavelengths should be created
for vertical incidence.

⬇

1import numpy as np

2import tqdm

3

4n_samples = 1e6

5n_lambda = 100

6n_layers = 12

7wavelength = np.linspace(1000, 1700, n_lambda)*1e-9

8

9theta = np.array([0]) # vertical incidence

10stack_layers = np.random.uniform(5, 180, (n_samples, n_layers))*1e-9

11stack_layers[:,0] = stack_layers[:,-1] = np.inf # injection and outcoupling
layer

12optical_index = np.array([2.5] +[2.0, 1.4]*5 +[1.])

13optical_index = np.tile(optical_index, (n_samples, 1))

14n = 10000 # the dataset is computed in steps of 1e5

15dataset = np.empty((n_samples, n_lambda))

16for i in tqdm.tqdm(np.array(range(n_samples))[::n]):

17 tmm_res = np.empty((n, n_lambda))

18 tmm_res = multithread_coh_tmm(’s’, optical_index[i:i+n],
stack_layers[i:i+n], theta, wavelength, TorR=’R’).squeeze()

19 dataset[i:i+n] = tmm_res

Now the dataset can be saved for example by using .hdf or Numpy’s .npz format.

###  A.3 TMM-Torch Example

Here, an example how to compute gradients via automatic differentiation with
the TMM-Torch package is shown.

⬇

1import torch

2import tmm_fast_torch as tmmt

3n_layers = 12 # number of layers

4stack_layers = np.random.uniform(20, 150, n_layers)*1e-9 # thicknesses of the
layers

5stack_layers[0] = stack_layers[-1] = np.inf # set first and last layer as
injection layer

6optical_index = np.random.uniform(1.2, 5, n_layers) # random constant
refractive index

7optical_index[-1] = 1 # outcoupling into air

8

9stack_layers = torch.tensor(stack_layers, requires_grad=True)

10

11wl = np.linspace(500, 900, 301)*1e-9

12theta = np.deg2rad(np.linspace(0, 90, 301))

13

14result = tmmt.coh_tmm_fast(’s’, optical_index, stack_layers, theta,
wavelength)[’R’]

15mse = torch.nn.MSELoss()

16error = mse(result, torch.zeros_like(result)

17error.backward()

18

19gradients = stack_layers.grad

Since the error is computed with respect to zero reflectivity for all
wavelengths at all incidences, the gradient points towards the steepest
descent for a broadband anti-reflection coating. By using the minimize
function from the Python package Scipy and a gradient based optimization
algorithm, for example
"L-BFGS-B"666https://docs.scipy.org/doc/scipy/reference/optimize.minimize-
lbfgsb.html, a local optimum is easily found.