# Leveraging Neural Machine Translation for Word Alignment

firstname=Vilém uds surname=Zouhar uds institute=uds uds corresponding=yes uds
address=Institute of Formal and Applied Linguistics
Faculty of Mathematics and Physics,
Charles University
Malostranské náměstí 25
118 00 Praha 1, Czech Republic uds<EMAIL_ADDRESS>uds uds
firstname=Daria uds surname=Pylypenko uds institute=uds uds
<EMAIL_ADDRESS>uds uds

###### Abstract

The most common tools for word-alignment rely on a large amount of parallel
sentences, which are then usually processed according to one of the IBM model
algorithms. The training data is, however, the same as for machine translation
(MT) systems, especially for neural MT (NMT), which itself is able to produce
word-alignments using the trained attention heads. This is convenient because
word-alignment is theoretically a viable byproduct of any attention-based NMT,
which is also able to provide decoder scores for a translated sentence pair.

We summarize different approaches on how word-alignment can be extracted from
alignment scores and then explore ways in which scores can be extracted from
NMT, focusing on inferring the word-alignment scores based on output sentence
and token probabilities. We compare this to the extraction of alignment scores
from attention. We conclude with aggregating all of the sources of alignment
scores into a simple feed-forward network which achieves the best results when
combined alignment extractors are used.

Saarland University, Department of Language Science and Technology

##  1 Introduction

Although word alignment found its use mainly in phrase-based machine
translation (for generating phrase tables), it is still useful for many other
tasks and applications: boosting neural MT performance (Alkhouli et al.,
2016), exploring cross-linguistic phenomena (Schrader, 2006), computing
quality estimation (Specia et al., 2013), presenting quality estimation
(Zouhar and Novák, 2020) or simply highlighting matching words and phrases in
interactive MT (publicly available MT services).

The aim of this paper is to improve the word alignment quality and demonstrate
the capabilities of alignment based on NMT confidence. Closely related to this
is the section devoted to aggregating multiple NMT-based alignment models
together, which outperforms the individual models. This is of practical use
(better alignment) as well as of theoretical interest (word alignment
information encoded in NMT scores).

We first briefly present the task of word alignment, the metric and the used
tools and datasets. In Section 2 we introduce the soft word alignment models
based on MT scores and also several hard word alignment methods (extractors).
The models are evaluated together with other solutions (fast_align and
Attention) in Section 3. We then evaluate the models enhanced with new
features and combined together using a simple feed-forward neural network
(Section 4). In both cases, we explore the models’ behaviour on Czech-English
and German-English datasets.

All of the code is available open-source.111github.com/zouharvi/LeverageAlign

###  1.1 Word Alignment

Word alignment (also bitext alignment) is a task of matching two groups of
words together that are each other’s semantic translation. This is problematic
for non-content words which are specific for the given language but generally
one is able to construct a mapping as in the example in Figure 1. Word
alignment usually follows after sentence alignment. Even though it is called
word alignment, it usually operates on all tokens, including punctuation
marks.

Figure 1:  Illustration of English (top) to German (bottom) word alignment.
The token >>choose<< is aligned to two tokens >>wählen<< and >>Sie<< while the
token >>Option<< is left unaligned. The article >>die<< is mistakenly aligned
to two unrelated articles >>the<<.

Word alignment output can be formalized as a set containing tuples of source-
target words. For an aligner output $A$, a sure alignment $S$ and a possible
alignment $P$ ($S\subseteq P$),222Sure alignments can be treated as gold
alignments with very high confidence, while pairs marked with possible
alignments are still sensible to connect, but with the decision being much
less clear. The AER is designed not to penalize models by including more
possible alignments in the gold annotations. precision can be computed as
$\frac{|A\cap P|}{|A|}$ and recall as $\frac{|A\cap S|}{|S|}$. The most common
metric, Alignment Error Rate (AER), is defined as $1-\frac{|A\cap S|+|A\cap
P|}{|S|+|A|}$ (lower is better). Even though the test set is annotated with
two types of alignments, the aligner is expected to produce only one type.
These evaluation measures are described by Mihalcea and Pedersen (2003) and
Och and Ney (2003).

Traditionally word alignment models can be split into soft and hard alignment
parts. In soft alignment, the model produces a score for every source-target
pair. When producing hard alignment (extractors), the model makes decisions as
to which alignments to include in the output. For source sentence $S$ and
target sentence $T$, the output of soft alignment is a
$\mathbb{R}^{|S|\times|T|}$ matrix while hard alignment is a set $A\subseteq
S\times T$.

#### Symmetrization.

Assuming that we have access to bi-directional word alignment (in the context
of this paper to two MT systems of the opposite directions) we can compute
both the alignment from source to target ($X$) and target to source
($X^{\prime}$). Having access to both $X$ and $X^{\prime}$ makes it possible
to create a new alignment $Y$ with either higher precision through
intersection or higher recall through union (Koehn, 2009).

$\displaystyle X^{T}:=\\{(b,a):(a,b)\in X\\}$ $\displaystyle Y_{prec}=X\cap
X^{\prime T}\qquad Y_{rec}=X\cup X^{\prime T}$

We can make use of the fact that the models output soft alignment scores and
create new alignment scores in the following way using a simple linear
regression model. This allows us to fine-tune the relevance of each of the
directions as well as their interaction. However, it does not have the same
effects as the union or the intersection because it affects the soft alignment
and not hard alignment in contrast to the previous case.

$\displaystyle p^{sym}(s,t)=\beta_{0}\cdot p(s,t)+\beta_{1}\cdot
p^{r}(t,s)+\beta_{2}\cdot p(s,t)\cdot p^{r}(t,s)$

More complex symmetrization techniques have been proposed and implemented by
Och and Ney (2000); Junczys-Dowmunt and Szał (2011).

###  1.2 Relevant Work

Och and Ney (2003) introduce the word alignment task and systematically
compare the IBM word alignment models. The work of Li et al. (2019) is closely
related to this article as it examines the issue of word alignment from NMT
and proposes two ways of extracting it: prediction difference and explicit
model. They also show that without guided alignment in training, NMT systems
perform worse than fast_align baseline. Using attention for word alignment is
thoroughly discussed by Bahdanau et al. (2014) and Zenkel et al. (2019). Word
alignment based on static and contextualized word embeddings is explored by
the recent work of Sabet et al. (2020). Word alignment based on cross-lingual
(more than 2 languages) methods is presented by Wu et al. (2021). The work of
Chen et al. (2020b) focuses on inducing word alignments from glass-box NMT as
a replacement for using Transformer attention layers directly. Chen et al.
(2020a) document Mask Align, an unsupervised neural word aligner based on a
single masked token prediction.

Chen et al. (2016) propose guided attention, a mechanism that uses word
alignment to bias the attention during training. This improves the MT
performance on especially rare and unknown tokens. The usage of word alignment
in this work is, however, opposite to the goals of this paper. While for Chen
et al. (2016) the word alignment improved their MT system, here the MT system
improves the word alignment.

###  1.3 Tools

The experiments in this paper require an MT system capable of providing output
probabilities (decoder scores) and optionally also attention-based word
alignment. For comparison, we also use an IBM-model-based word aligner. This
tool is also used as an additional feature to the final aggregation model.

#### MarianNMT

(Junczys-Dowmunt et al., 2018a, b) is a popular (both in academia and in
deployment scenarios), actively developed and maintained system for fast
machine translation. It already contains options for producing word alignment,
output probabilities for words and sentences and also attention scores.

#### fast_align

(Dyer et al., 2013) is an unsupervised word aligner based on IBM Alignment
Model 2. It does not provide state of the art pre-neural performance but is
easy to build with modern toolchains and has low resource requirements (both
memory- and computational-wise).

###  1.4 Data

For training purposes, we make use mostly of the parallel corpora of
Czech–English word alignments by Mareček (2016), based on manually annotated
data. We also include a large Czech-English corpus by Kocmi et al. (2020) and
a large German–English corpus by Rozis and Skadiņš (2017), which are not word
aligned. From this corpus, 1M sentences were sampled randomly. A small
manually aligned German–English corpus by Biçici (2011) is included for
testing. An overview of the corpus sizes is displayed in Table 1.

CS/DE-English | Type | Domain | CS/DE Tok. | EN Tok. | Sent.
---|---|---|---|---|---
Czech Small | aligned | news, legal | 53k | 60k | 2.5k
Czech Big | unaligned | multi | 2618M | 3013M | 188M
German Small | aligned | legal | 1k | 1k | 0.1k
German Big | unaligned | tech, news, legal | 23M | 25M | 1M
Table 1: Used word aligned corpora with their sizes, domains and origin.

###  1.5 MT Models

We make use of the MT models made available333github.com/browsermt/students by
Germann et al. (2020) and Bogoychev et al. (2020). For both Czech-English and
English-German, CPU-optimized student models are used. They are transformer-
based (Vaswani et al., 2017) and were created by using knowledge distillation.
With WMT19 and WMT20 SacreBLEU (Post, 2018), the models achieve the following
BLEU scores: Czech-English (27.7), English-Czech (36.3) and English–German
(42.7).444BLEU+case.mixed+lang.cs-en+numrefs.1
+smooth.exp+test.wmt20+tok.13a+version.1.4.13
BLEU+case.mixed+lang.en-cs+numrefs.1
+smooth.exp+test.wmt20+tok.13a+version.1.4.13
BLEU+case.mixed+lang.en-
de+numrefs.1+smooth.exp+test.wmt19+tok.13a+version.1.4.8 Since the
English–German MT is only available in one direction, word alignment is
reported in this direction as well. Exceptions, such as word alignment using
an MT for the opposite direction, are explicitly mentioned.

##  2 Individual Models

In this section, we describe and evaluate the individual word alignment
models. All of the newly introduced models make use of the fact that NMT
systems can be viewed as language models and can produce translation
probabilities.

###  2.1 Baseline Models

The first model is fast_align. The second is attention-based soft word
alignment extracted from MarianNMT (Attention), which was trained with guided
alignment during the distillation. For the rest of this subsection, we will
focus on models generating soft alignment scores (an unbounded real number
corresponding to the quality of a possible alignment between two tokens) and
not the alignments themselves.

#### One Token Translation ($M_{1}$).

The simplest approach to get alignment scores is to compute decoder
translation probability using the MT (function $m$) between every source and
target token $s_{i}$ and $t_{j}$ of the source and target sentences $S$ and
$T$. Single tokens are passed to the models as if they were a sentence pair.
The scores are not normalized which is not an issue in this case, since the
models working with these alignment scores (in Section 2.2) compare output
from sequences of the same length.

$\displaystyle\forall s_{i}\in S,t_{j}\in
T:p(s_{i},t_{j})=m(\\{s_{i}\\},\\{t_{j}\\})$

The produced values are in a log space $(-\inf,0]$. This approach requires
$|S|\cdot|T|$ of one-token translation scorings (decoder probability of the
target reference) for producing word alignments of a single sentence pair. On
a CPU,5558 threads 2.3GHz Ryzen 7 3700u, no RAM to disk swapping the models
average to $2.7$s per one sentence pair alignment.

#### Source Token Dropout ($M_{2}$).

A more refined approach was chosen by Zintgraf et al. (2017) in which the
alignment score is computed as the difference in target token probability when
the source token is unknown. The exact approach is too computationally
demanding (requires translation scorings with large amounts of replacement
words), and therefore we use a much simpler, yet conceptually similar method
by either omitting the token or replacing it with <unk>.666Even though
subword-based MT models do not need <unk>, SentencePiece reserves the token
<unk> for an unknown symbol. Assume $m_{j}(S,T)$ produces the log probability
of the $j$-th target token. The sentence $S^{a/b}_{i}$ with an obscured token
$s_{i}$ can be defined in two ways which leads to two versions of this model:
$M_{2}^{a}$ and $M_{2}^{b}$. Output is then possibly unbounded $(-\inf,\inf)$.

$\displaystyle\forall s_{i}\in S,t_{j}\in
T:p(s_{i},t_{j})=m_{j}(S,T)-m_{j}(S^{a/b}_{i},T)$

$\displaystyle\text{Word deletion}\,(M_{2}^{a}):$ $\displaystyle
S^{a}_{i}=s_{0},s_{1},\ldots,s_{i-1},s_{i+1},\ldots,s_{|S|}$
$\displaystyle\text{Word substitution}\,(M_{2}^{b}):$ $\displaystyle
S^{b}_{i}=s_{0},s_{1},\ldots,s_{i-1},\texttt{<unk>},s_{i+1},\ldots,s_{|S|}$

This requires $|S|$ translation scorings of source and target lengths $|S|$
and $|T|$, which is comparable to $M_{1}$. The models average to $1.5$s per
one sentence pair alignment.777 The running time is lower because in this case
it is $|S|$ scorings of length $|T|$, while in $M_{1}$ it is $|S|\times|T|$
scorings of length $1$.

#### Source and Target Dropout ($M_{3}$).

A very similar method would be to also dropout the target token and examine
how the sentence probability changes. Applying the two different ways of
dropout leads to four versions of this approach. Note that in this case we
compute the sentence probability (because the target word is hidden) and also
do not subtract from the base sentence probability, but rather use the new
sentence probability as it is. This probability should be highest if the
corresponding tokens are both obscured. The probability is in log space
$(-\inf,0]$.

$\displaystyle\forall s_{i}\in S,t_{j}\in
T:p(s_{i},t_{j})=m(S^{a/b}_{i},T^{a/b}_{j})$

$\displaystyle T^{a}_{j}$
$\displaystyle=t_{0},t_{1},\ldots,t_{j-1},t_{j+1},\ldots,t_{|T|}$
$\displaystyle T^{b}_{j}$
$\displaystyle=t_{0},t_{1},\ldots,t_{j-1},\texttt{<unk>},t_{j+1},\ldots,t_{|T|}$

Word deletion, deletion ($M_{3}^{aa}$) $\displaystyle S^{a}_{i},T^{a}_{j}$
Word deletion, substitution ($M_{3}^{ab}$) $\displaystyle S^{a}_{i},T^{b}_{j}$
Word substitution, deletion ($M_{3}^{ba}$) $\displaystyle S^{b}_{i},T^{a}_{j}$
Word substitution, substitution ($M_{3}^{bb}$) $\displaystyle
S^{b}_{i},T^{b}_{j}$

This approach requires $|S|\cdot|T|$ translation scorings of source and target
lengths of $|S^{a/b}|$ and $|T^{a/b}|$ for sentence $S$ translated to $T$
which is roughly $|T|$ times more than in $M_{1}$ and $M_{2}$. This makes it
it the most computationally demanding approach. On average it takes $46.1$s to
produce one sentence pair alignment on a CPU.

###  2.2 Direct Alignment from Baseline Models

All of the models (except for fast_align) are not producing the alignments
themselves, but soft alignment scores $p$ for each pair of tokens $(s,t)$ in
source $S$ $\times$ target $T$ sentence. The hard alignment itself can then,
for example, be computed in the following ways. The parameter $\alpha$ can be
estimated from the development set. The function $p$ is in general any soft
alignment function (e.g. attention scores or the alignment scores from IBM
model 1 EM algorithm).

  1. 1.

For every source token $s$ take the target tokens $t$ with the maximum score.

$\displaystyle A_{1}=\bigcup_{s\in S}\\{(s,t):p(s,t)=max_{r}\\{p(s,r)\\}\\}$

  2. 2.

For every source token $s$ take all target tokens $t$ with a high enough score
(above threshold). This method is used to control the density of alignments in
the work of Liang et al. (2006) and provides a parameter to tradeoff precision
and recall.

$\displaystyle A_{2}^{\alpha}=\bigcup_{s\in S}\\{(s,t):p(s,t)\geq\alpha\\}$

  3. 3.

For every source token $s$ take any target token which has a score of at least
$\alpha$ times the score of the best candidate. Special handling for negative
cases in the form of a division is required to make the formula work for the
whole $\mathbb{R}$. The motivation for this is $M_{2}$, which provides
possibly unbounded alignment scores. Assume $\alpha\in(0,1]$.

$\displaystyle A_{3}^{\alpha}=\bigcup_{s\in S}\\{(s,t):p(s,t)\geq
min\big{[}\max_{r}\ p(s,r)\cdot\alpha,\max_{r}\ p(s,r)\,/\,\alpha\big{]}\\}$

$A_{1}$ can then be expressed as $A_{3}^{1}$. Lower $alpha$ values lead to
lower precision and higher recall because the algorithm includes more, less
probable, alignments. A variation on this approach would be to subtract
$\alpha$ instead of multiplying it. The reason for choosing multiplication is
that it dynamically adapts to a wider range of intervals and bounds the
parameters between $0$ and $1$. This is not the case for substraction and
because of this, it would be harder to choose the right parameter.

  4. 4.

Similar approach is for $A_{3}$, but with the focus on the target side. For
every target token $t$ take any source token which has a score of at least
$\alpha$ times the score of the best candidate.

$\displaystyle A_{4}^{\alpha}=\bigcup_{t\in T}\\{(s,t):p(s,t)\geq
min\big{[}\max_{r}\ p(r,t)\cdot\alpha,\max_{r}\ p(r,t)/\alpha\big{]}\\}$

Similar reversal for $A_{2}$ would not make sense, because it already takes
all alignment above a threshold without any consideration for the direction.

  5. 5.

Similarly to $M_{3}$ and $M_{4}$ it is possible to create an extractor in
which instead of having a single dropout on the target side, there are a
multiple of them. This way, the score would not be between the source token
and the target token, but between the source token and a subset of all target
tokens. Formally, this would replace the (complete) weighted graph structure
with a (complete) hypergraph. Instead of just having a weight for
Choose–Wählen, there would also be a weight for {Choose}–{Wählen, Sie},
{Choose}-{Wählen, im, Popupmenü} etc. This would, however, lead to exponential
complexity in terms of target sentence length. The number of words
participating in an edge would then have to be limited to the number of
alignments to a single token that we can empirically expect of the given
language pair. Figure 1 suggests that for English-German this could be 3. Upon
computing the scores for all the edges in this hypergraph, a follow-up task
would be to find the maximum-weight matching.

#### Coverage.

The suggested greedy way of computing alignments from alignment scores is far
from perfect. In the scenario depicted in Figure 2, all but the last source
token (German) have been aligned with the target, each with different
alignment scores. Although the model may lack any lexical knowledge of the
word Übersetzung, it should consider the prior of a word being aligned to at
least one target token.

Figure 2:  Partial alignment from German (top) to Sesotho (bottom). The model
has no lexical knowledge about the alignment of >>Übersetzung<<, though
>>phetolelo<< is a good candidate because no other word aligned to it. Line
strength corresponds to the soft alignments produced by the model.

In this specific case, $A_{3}^{0.9}$ would probably include all alignments to
the word Übersetzung, since there is no single strong candidate (assume that
lines not visible depict soft alignments close to $0$). Similarly,
$A_{4}^{0.9}$ would also include most alignments of the word Übersetzung,
including the word phetolelo, since the alignment score with maschinelle is
weak and also close to $0$. Intersecting these two extractors $A_{3}^{0.9}\cap
A_{4}^{0.9}$ would yield the correct alignment Übersetzung–phetolelo. Other
tokens would not be aligned to either of these two words because they have
strong alignment scores with different tokens.

This prior may not always be desirable. For this, intersecting with
$A_{2}^{\alpha}$ provides a limiting threshold. In an application where the
target token is erroneous, this prevents the alignment model from aligning the
two corresponding tokens. Inducing alignment based on graph properties is
examined by Matusov et al. (2004), though without the presence of NMT.

##  3 Evaluation of Individual Models

#### Baseline Models.

Figure 3 shows the results on Czech$\leftrightarrow$English data averaged from
both directions. Different models have different spans of their scores, and
therefore it is much harder to select the single best $\alpha$. The most basic
model, $M_{1}$, achieves the best performance (AER = $0.46$). The figure
serves as an illustration of the $A_{2}^{\alpha}$ landscape.

Figure 3: Precision, Recall and AER of individual models on
CS$\leftrightarrow$EN data extracted using $A_{2}$ (directions averaged)

The results on Czech$\leftrightarrow$English data averaged from both
directions with $A_{3}$ can be seen in Figure 4. The case of $\alpha=0$
corresponds to aligning everything with everything, while $\alpha=1$ means
aligning only the token with the highest score to the single source one (i.e.
$A_{1}$). The different model families behave similarly with respect to
Precision, Recall and AER. $M_{1}$ achieves again the best result (AER =
$0.34$), but with a smoother distinction between models.

Out of the model $M_{3}$ family, $M_{3}^{bb}$ outperformed the rest
significantly. In $A_{2}$ (Figure 3), the other models, $M_{3}^{aa}$,
$M_{3}^{ab}$ and $M_{3}^{ba}$, perform worse than $M_{2}^{a}$ and $M_{2}^{b}$.
This is reversed in case of using the $A_{3}$ extractor, as shown in Figure 4
and Figure 5. For the $M_{3}$ model family, models with mixed obscuring
functions ($M_{3}^{ab}$ and $M_{3}^{ba}$) perform worse than with the same
obscuring function on both the source and the target side ($M_{3}^{aa}$ and
$M_{3}^{bb}$).

Figure 4: Precision, Recall and AER of individual models on
CS$\leftrightarrow$EN extracted using $A_{3}$ (directions averaged)

The English$\rightarrow$German dataset proved to be more difficult. The AER,
that are shown in Figure 5, are higher than for Czech$\leftrightarrow$English.
The model $M_{1}$ again achieves the best results with AER = $0.43$. The model
ordering is preserved from Figure 4.

Figure 5: Precision, Recall and AER of individual models on EN$\rightarrow$DE
extracted using $A_{3}$

Figure 6 documents that different model types produce different number of
alignments per one token. It also shows that the performance rapidly decreases
with sentence length. The high AER in Figure 4 can be explained by the dataset
containing mostly longer sentences (21 tokens on average). The model $M_{1}$
is still better than $M_{3}^{bb}$ even on longer sentences despite the fact it
does not model the context.

Figure 6: AER for $\alpha=1$ (left) and average number of aligned tokens
(right) of individual baseline models on CS$\leftrightarrow$EN extracted using
$A_{3}$ (directions averaged)

The best results were achieved with $A_{4}^{1}$ using $M_{1}$: AER = $0.30$
for German$\rightarrow$English and AER = $0.31$ for
Czech$\leftrightarrow$English. The plots (not shown) are very similar to those
of $A_{3}$. Hence $M_{3}^{bb}$ follows up with AER = $0.38$ and AER = $0.36$
for German and Czech respectively using $A_{4}^{1}$.

Data | Precision | Recall | AER
---|---|---|---
Czech$\leftrightarrow$English Small | $0.54$ | $0.66$ | $0.41$
Czech$\leftrightarrow$English Big | $0.63$ | $0.64$ | $0.38$
German$\rightarrow$English Small | $0.49$ | $0.55$ | $0.48$
German$\rightarrow$English Small+Big | $0.63$ | $0.72$ | $0.34$
Table 2: Precision, Recall and AER of fast_align. Models were evaluated on the
respective annotated datasets part.

#### fast_align.

For comparison, the results of fast_align can be seen in Table 2. For both
language pairs, we use two models, trained on the Small and Big corpora. The
motivation for the latter is that the performance of fast_align on 5k sentence
pairs is unfairly low in comparison to the other methods because the used MT
system has had access to a much larger amount of data. This is shown by the
performance difference between these two models.

Data | Subword Aggregation | Precision | Recall | AER
---|---|---|---|---
Czech$\leftrightarrow$English Small | maximum | $0.64$ | $0.81$ | $0.29$
Czech$\leftrightarrow$English Small | average | $0.64$ | $0.81$ | $0.29$
German$\rightarrow$English Small | maximum | $0.69$ | $0.81$ | $0.26$
German$\rightarrow$English Small | average | $0.68$ | $0.80$ | $0.27$
Table 3: Precision, Recall and AER of attention-based word alignment extracted
using $A_{3}^{1}$

#### Attention Scores.

Extracting alignment from MT model attention using $A_{3}^{1}$ results in the
highest performance (Table 3). Since the attention scores are between subword
units from SentencePiece (Kudo and Richardson, 2018), we chose two methods of
aggregation to a single score between two tokens (two lists of subwords): (1)
taking the maximum probability between two subwords and (2) taking the average
probability. They, however, produce almost identical results with respect to
the word alignment quality. Scores are listed with $A_{3}^{1}$, but
$A_{2}^{0.25}$ achieved very close results.

Model | Method | Precision | Recall | AER
---|---|---|---|---
$M_{1}$ | reverse | $0.56$ | $0.82$ | $0.35$
$M_{1}$ | add | $0.59$ | $0.86$ | $0.31$
$M_{1}$ | intersect | $0.73$ | $0.77$ | $0.26$
Attention (avg) | reverse | $0.64$ | $0.81$ | $0.29$
Attention (avg) | multiply | $0.66$ | $0.83$ | $0.28$
Attention (avg) | intersect | $0.77$ | $0.70$ | $0.28$
Table 4: Average Precision, Recall and AER on Czech$\leftrightarrow$English
extracted using $A_{4}^{1}$ with symmetrization methods applied for $M_{1}$
and Attention (avg)

#### Symmetrization.

Results of symmetrization methods (akin to those described in Section 1.1) for
$M_{1}$ and Attention scores (attention scores aggregated by averaging) are
shown in Table 4. Each method is accompanied by an example formula; $p^{x}$
stands for either $M_{1}$ or Attention (avg) (in principle any function which
produces soft alignments). Similarly, $A_{4}^{1}$ could be replaced by other
extractors, even though this one worked the best. For reverse and add,
$A_{4}^{1}$ is applied on the final result, but for simplicity left out of the
formulas.

Method reverse consists of using TGT$\rightarrow$SRC translation direction to
get alignment scores but then transposing the soft alignment matrix so that
the scores are SRC$\rightarrow$TGT.

$\displaystyle p_{CS\rightarrow EN}^{\text{reverse}}(s,t)=p^{x}_{EN\rightarrow
CS}(t,s)$

Method add simply combines the original and reversed scores before alignments
are extracted. The scores of $M_{1}$ are in log space; therefore, addition is
used instead of multiplication. For attentions, multiplication is used, since
they are bounded by $[0,1]$.

$\displaystyle p_{CS\rightarrow EN}^{\text{add}}(s,t)=p^{x}_{CS\rightarrow
EN}(s,t)+p^{x}_{EN\rightarrow CS}(t,s)$ $\displaystyle p_{CS\rightarrow
EN}^{\text{mutliply}}(s,t)=p^{x}_{CS\rightarrow EN}(s,t)\cdot
p^{x}_{EN\rightarrow CS}(t,s)\hskip 4.83691pt$

Method intersect first extracts the alignments for the two directions and then
intersects the results (with one direction transposed). This method produces
the best results overall (AER = $0.26$), also surpassing $M_{1}$’s forward
direction and attention-based alignments.

$\displaystyle A_{4}^{1}(p_{CS\rightarrow
EN}^{\text{intersect}}(s,t))=A_{4}^{1}(p^{x}_{CS\rightarrow EN}(s,t))\cap
A_{4}^{1}(p^{x}_{EN\rightarrow CS}(t,s))$

In contrast to $M_{1}$, none of the other models, including attention-based,
improved rapidly. This is partly explained by the fact that in other models,
the precision-recall balance is shifted from recall to precision, while in
$M_{1}$ it became more balanced after intersection. The reversal also allowed
us to get significant results (AER = $0.27$) for the
English$\rightarrow$German direction using Attention (avg), for which we did
not have an MT system.

###  3.1 Extractor Limitations

Computing word alignments by taking the most probable target token
($A_{3}^{1},A_{4}^{1}$) has theoretical limitations to the AER because it
makes a faulty assumption that every token is aligned to at least one other
token. The Czech$\rightarrow$English dataset has $12\%$ of unaligned tokens
and an average of $1.16$ aligned target tokens per source tokens (excluding
non-aligned tokens).

Assuming access to a word alignment oracle ($0$ if not aligned, $1$ if
aligned), in case the token is not aligned to any other, all of the scores are
$0$. The extractor $A_{3}^{1}=A_{1}$ will then take all tokens with values
equal to the maximum, effectively aligning the in reality unaligned token to
every possible one. This extractor is then bound to have maximum recall, but
relatively poor accuracy.

The measured performance shows that the $A_{2}^{\alpha}$ is not the best
extraction method. However, it is objectively not prone to this issue because
it does not make any assumptions about the number of aligned tokens, and the
minimum possible AER is 0 ($A_{2}^{1}$ with an oracle). In the next section,
we will therefore make use of $A_{2}^{\alpha_{1}}\cap A_{3}^{\alpha_{2}}\cap
A_{4}^{\alpha_{3}}$, which provides better performance than individual
extractors.

##  4 Ensembling of Individual Models

In the previous section, we saw that multiple methods with different
properties achieved good results, but were sensitive to the method used to
induce hard alignment. This section combines them together in a small feed-
forward neural network, which can be trained on a small amount of data.

###  4.1 Model

The ensemble neural network itself is a regressor:
$\mathbf{F}\rightarrow(0,1)$, where $\mathbf{F}$ is the set of feature vectors
for every pair of source and target tokens.888A completely different approach
would be to simply use (pretrained) word embeddings as an input to the
network. This is, however, not possible due to the low amount of gold
alignment data. By applying sigmoid to the output and establishing a threshold
value for the positive class, the network would become a classifier. This
behaviour can, however, be simulated using $A_{2}^{\alpha}$. We work with the
threshold explicitly and use the network for computing alignment score and not
for the alignment itself. For the hard alignment, we use $A_{2}^{0.001}\cap
A_{3}^{1}\cap A_{4}^{1}$, which we found to work the best with this ensemble
on the training data.

#### Additional Features.

Apart from $M_{1}$, $M_{2}^{b}$, $M_{3}^{aa}$, $M_{3}^{bb}$ and Attention with
averaging aggregation (Individual), we also include the output of fast_align
as one of the features. Moreover, four other manually crafted features
(Manual) are added. The motivation for the first two manual features is that
the position and token length help in determining the alignment in some cases.
The last two are specifically targeted at named entities, which have sparse
occurrences in the data, and also at non-word tokens, such as full stops,
delimiters and quotation marks. We list Pearson’s correlation coefficient with
true alignments on Czech$\leftrightarrow$English data (the two directions
averaged).

  * •

Difference in sentence positions:
$\rho=\text{-}\ 0.18$, $abs(\,i/|S|-j/|T|\,)$

  * •

Difference in token lengths:
$\rho=\text{-}\ 0.11$, $abs(\,|s_{i}|-|t_{j}|\,)$

  * •

Difference in subword unit counts:
$\rho=\text{-}\ 0.03$, $abs(\,|\text{subw}(s_{i})|-|\text{subw}(t_{j})|\,)$

  * •

Normalized token case-insensitive Levenshtein distance:
$\rho=\text{-}\ 0.30$, $lev(s_{i},t_{j})/max(|s_{i}|,|t_{j}|)$

  * •

Number of subword units which are present in both tokens:999Normalized version
of this feature had slightly lower correlation coefficient: $0.30$.
$\rho=\,\,\,0.32$, $|\,\text{subw}(s_{i})\cap\text{subw}(t_{j})\,|$

  * •

Token string case-insensitive equality (equal to zero Levenshtein distance):
$\rho=\,\,\,0.28$, $I_{s_{i}\simeq t_{j}}$

#### Architecture.

For every model, the epoch with the lowest AER on the validation dataset is
used for the test dataset. This extractor was found to work best across all
ensemble models. The training was done with cross-entropy loss. The model was
composed of series of hidden linear layers, each with biases and Tanh as the
activation function with dropouts around the innermost layer:

$\displaystyle L_{|\text{Input}|}^{\text{Tanh}}\circ L_{32}^{\text{Tanh}}\circ
D_{0.2}\circ L_{16}^{\text{Tanh}}\circ D_{0.2}\circ L_{16}^{\text{Tanh}}\circ
L_{8}^{\text{Tanh}}\circ L_{1}^{\text{Softmax}}$

###  4.2 Data

The Czech$\leftrightarrow$English dataset contains 1.5M source-target pairs,
out of which $2.64\%$ is of a positive class (aligned). For
German$\leftrightarrow$English Small these quantities are 22k and $5.61\%$
respectively. This could be an issue for a simple classifier network and would
need e.g. oversampling of the positive or undersampling of the negative class.

For Czech$\leftrightarrow$English, we used $10\%$ and $10\%$ (250 sentences
each) for validation and test data and the rest for training. Samples were
split on sentence boundaries. The English$\rightarrow$German was used solely
for testing, due to its small size.

###  4.3 Evaluation

The averaged results of each ensemble on Czech$\leftrightarrow$English are in
Table 5. We also show the results of $M_{1}$, but without $A_{2}$. Due to the
range of $M_{1}$’s values, it is difficult to establish a cut-off threshold.
Attention uses $A_{3}^{1}$, since intersection with other extractors did not
improve the performance, as described in Section 3. The results demonstrate
that adding any feature improves the overall ensemble. All features combined
together improve on the best individual model by $-0.11$ AER.101010Performed
by Student’s t-test on 10 runs with $p<0.001$.

Model / Features | Precision | Recall | AER
---|---|---|---
$M_{1}$ ($A_{3}^{1}\cap A_{4}^{1}$) | $0.75$ | $0.78$ |  $0.25$ $\star$
Attention (max, $A_{3}^{1}$) | $0.64$ | $0.81$ | $0.29$
fast_align Small | $0.54$ | $0.66$ | $0.41$
fast_align Big | $0.63$ | $0.64$ | $0.38$
Manual features | $0.55$ | $0.46$ | $0.50$
Individual ($M_{1}$, $M_{2}^{b}$, $M_{3}^{aa}$, $M_{3}^{bb}$, attention) | $0.84$ | $0.73$ | $0.23$
Manual + Indiv. | $0.85$ | $0.79$ | $0.19$
Manual + Indiv. + fast_align | $0.86$ | $0.79$ | $0.18$
Manual + Indiv. + fast_align + Attention | $0.85$ | $0.84$ | $0.16$
Manual + Indiv. + fast_align + Attention + M1 rev. | $0.86$ | $0.86$ |  $0.14$ $\star$
Table 5: Average Precision, Recall and AER of $M_{1}$ (best individual) and
different ensemble models (using $A_{2}^{0.001}\cap A_{3}^{1}\cap A_{4}^{1}$)
on Czech$\leftrightarrow$English data (averaged)

#### Transfer.

The best models on Czech$\leftrightarrow$English (one for each direction) were
then used on the English$\rightarrow$German dataset, resulting on average in
AER = $0.18$. This is higher than for Czech but still significantly lower (by
a margin of $-0.08$)111111Performed by Student’s t-test with $p<0.001$. than
for the best individual model, Attention (max). This suggests that the
features generalize well and models can be trained even on other language
data. Furthermore, since the alignment datasets come from different origins,
there may be systematic biases, which lower the performance of the transfer.

##  5 Summary

This paper explored and compared different methods of inducing word alignment
from trained NMT models. Despite its simplicity, estimating scores with single
word translations (combined with reverse translations) appears to be the
fastest and the most robust solution, even compared to word alignment from
attention heads. Ensembling individual model scores with a simple feed-forward
network improves the final performance to AER = $0.14$ on
Czech$\leftrightarrow$English data.

#### Future work.

Section 2.1 presented but did not explore an idea of target dropout with
multiple tokens in order to better model the fact that words rarely map 1:1.
We then used neural MT for providing alignment scores but then used a
primitive extractor algorithm for obtaining hard alignment. More sophisticated
approaches which consider the soft alignment origin (NMT), could vastly
improve the performance.

Although it was possible to use any alignment extractor to get hard alignments
out of soft ones, we found that the choice of the mechanism and also the
parameters had a considerable influence on the performance. These alignment
extractors are, however, not bound to alignment from NMT and their ability to
be used with other soft alignment models and other symmetrization techniques
should be examined further.

Finally, we did not explore the possible effects of fine-tuning the
translation model on the available data or training it solely on this data.
Similarity based on word embeddings could be used as yet another soft-
alignment feature.

## Acknowledgements

We would like to thank Bernadeta Griciūtė and Nikola Kalábová for their
reviews and insightful comments.

This article has used services provided by the LINDAT/CLARIAH-CZ Research
Infrastructure (lindat.cz), supported by the Ministry of Education, Youth and
Sports of the Czech Republic (Project No. LM2018101). The work described also
used trained models made available by the H2020-ICT-2018-2-825303 (Bergamot)
grant.

## References

  * Alkhouli et al. (2016) Alkhouli, Tamer, Gabriel Bretschner, Jan-Thorsten Peter, Mohammed Hethnawi, Andreas Guta, and Hermann Ney.  Alignment-based neural machine translation.  In _Proceedings of the First Conference on Machine Translation: Volume 1, Research Papers_ , pages 54–65, 2016.
  * Bahdanau et al. (2014) Bahdanau, Dzmitry, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate.  _arXiv preprint arXiv:1409.0473_ , 2014.
  * Biçici (2011) Biçici, Ergun.  _The Regression Model of Machine Translation_.  PhD thesis, Koç University, 2011.  Supervisor: Deniz Yuret.
  * Bogoychev et al. (2020) Bogoychev, Nikolay, Roman Grundkiewicz, Alham Fikri Aji, Maximiliana Behnke, Kenneth Heafield, Sidharth Kashyap, Emmanouil-Ioannis Farsarakis, and Mateusz Chudyk.  Edinburgh’s Submissions to the 2020 Machine Translation Efficiency Task.  In _Proceedings of the Fourth Workshop on Neural Generation and Translation_ , pages 218–224, Online, July 2020. Association for Computational Linguistics.  URL https://www.aclweb.org/anthology/2020.ngt-1.26.
  * Chen et al. (2020a) Chen, Chi, Maosong Sun, and Yang Liu.  Mask-Align: Self-Supervised Neural Word Alignment.  _arXiv preprint arXiv:2012.07162_ , 2020a.
  * Chen et al. (2016) Chen, Wenhu, Evgeny Matusov, Shahram Khadivi, and Jan-Thorsten Peter.  Guided alignment training for topic-aware neural machine translation.  _arXiv preprint arXiv:1607.01628_ , 2016.
  * Chen et al. (2020b) Chen, Yun, Yang Liu, Guanhua Chen, Xin Jiang, and Qun Liu.  Accurate Word Alignment Induction from Neural Machine Translation.  _arXiv preprint arXiv:2004.14837_ , 2020b.
  * Dyer et al. (2013) Dyer, Chris, Victor Chahuneau, and Noah A Smith.  A simple, fast, and effective reparameterization of ibm model 2.  In _Proceedings of the 2013 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 644–648, 2013.
  * Germann et al. (2020) Germann, Ulrich, Roman Grundkiewicz, Martin Popel, Radina Dobreva, Nikolay Bogoychev, and Kenneth Heafield.  Speed-optimized, Compact Student Models that Distill Knowledge from a Larger Teacher Model: the UEDIN-CUNI Submission to the WMT 2020 News Translation Task.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 190–195, Online, Nov. 2020. Association for Computational Linguistics.  URL https://www.aclweb.org/anthology/2020.wmt-1.17.
  * Junczys-Dowmunt and Szał (2011) Junczys-Dowmunt, Marcin and Arkadiusz Szał.  Symgiza++: symmetrized word alignment models for statistical machine translation.  In _International Joint Conferences on Security and Intelligent Information Systems_ , pages 379–390. Springer, 2011.
  * Junczys-Dowmunt et al. (2018a) Junczys-Dowmunt, Marcin, Roman Grundkiewicz, Tomasz Dwojak, Hieu Hoang, Kenneth Heafield, Tom Neckermann, Frank Seide, Ulrich Germann, Alham Fikri Aji, Nikolay Bogoychev, et al.  Marian: Fast neural machine translation in C++.  _arXiv preprint arXiv:1804.00344_ , 2018a.
  * Junczys-Dowmunt et al. (2018b) Junczys-Dowmunt, Marcin, Kenneth Heafield, Hieu Hoang, Roman Grundkiewicz, and Anthony Aue.  Marian: Cost-effective high-quality neural machine translation in C++.  _arXiv preprint arXiv:1805.12096_ , 2018b.
  * Kocmi et al. (2020) Kocmi, Tom, Martin Popel, and Ondrej Bojar.  Announcing czeng 2.0 parallel corpus with over 2 gigawords.  _arXiv preprint arXiv:2007.03006_ , 2020.
  * Koehn (2009) Koehn, Philipp.  _Statistical machine translation_.  Cambridge University Press, 2009.
  * Kudo and Richardson (2018) Kudo, Taku and John Richardson.  Sentencepiece: A simple and language independent subword tokenizer and detokenizer for neural text processing.  _arXiv preprint arXiv:1808.06226_ , 2018.
  * Li et al. (2019) Li, Xintong, Guanlin Li, Lemao Liu, Max Meng, and Shuming Shi.  On the word alignment from neural machine translation.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 1293–1303, 2019.
  * Liang et al. (2006) Liang, Percy, Ben Taskar, and Dan Klein.  Alignment by agreement.  In _Proceedings of the Human Language Technology Conference of the NAACL, Main Conference_ , pages 104–111, 2006.
  * Mareček (2016) Mareček, David.  Czech-English Manual Word Alignment, 2016.  URL http://hdl.handle.net/11234/1-1804.  LINDAT/CLARIAH-CZ digital library at the Institute of Formal and Applied Linguistics (ÚFAL), Faculty of Mathematics and Physics, Charles University.
  * Matusov et al. (2004) Matusov, Evgeny, Richard Zens, and Hermann Ney.  Symmetric word alignments for statistical machine translation.  01 2004.  10.3115/1220355.1220387.
  * Mihalcea and Pedersen (2003) Mihalcea, Rada and Ted Pedersen.  An evaluation exercise for word alignment.  In _Proceedings of the HLT-NAACL 2003 Workshop on Building and using parallel texts: data driven machine translation and beyond_ , pages 1–10, 2003.
  * Och and Ney (2000) Och, Franz Josef and Hermann Ney.  Improved statistical alignment models.  In _Proceedings of the 38th annual meeting of the association for computational linguistics_ , pages 440–447, 2000.
  * Och and Ney (2003) Och, Franz Josef and Hermann Ney.  A systematic comparison of various statistical alignment models.  _Computational linguistics_ , 29(1):19–51, 2003\.
  * Post (2018) Post, Matt.  A Call for Clarity in Reporting BLEU Scores.  In _Proceedings of the Third Conference on Machine Translation: Research Papers_ , pages 186–191, Belgium, Brussels, Oct. 2018. Association for Computational Linguistics.  URL https://www.aclweb.org/anthology/W18-6319.
  * Rozis and Skadiņš (2017) Rozis, Roberts and Raivis Skadiņš.  Tilde MODEL-multilingual open data for EU languages.  In _Proceedings of the 21st Nordic Conference on Computational Linguistics_ , pages 263–265, 2017.
  * Sabet et al. (2020) Sabet, Masoud Jalili, Philipp Dufter, and Hinrich Schütze.  Simalign: High quality word alignments without parallel training data using static and contextualized embeddings.  _arXiv preprint arXiv:2004.08728_ , 2020.
  * Schrader (2006) Schrader, Bettina.  How does morphological complexity translate? A cross-linguistic case study for word alignment.  In _Proceedings of Linguistic Evidence Conference_ , pages 189–191, 2006.
  * Specia et al. (2013) Specia, Lucia, Kashif Shah, José GC De Souza, and Trevor Cohn.  QuEst-A translation quality estimation framework.  In _Proceedings of the 51st Annual Meeting of the Association for Computational Linguistics: System Demonstrations_ , pages 79–84, 2013.
  * Vaswani et al. (2017) Vaswani, Ashish, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention Is All You Need, 2017.
  * Wu et al. (2021) Wu, Di, Liang Ding, Shuo Yang, and Dacheng Tao.  SLUA: A Super Lightweight Unsupervised Word Alignment Model via Cross-Lingual Contrastive Learning.  _arXiv preprint arXiv:2102.04009_ , 2021.
  * Zenkel et al. (2019) Zenkel, Thomas, Joern Wuebker, and John DeNero.  Adding interpretable attention to neural translation models improves word alignment.  _arXiv preprint arXiv:1901.11359_ , 2019.
  * Zintgraf et al. (2017) Zintgraf, Luisa M, Taco S Cohen, Tameem Adel, and Max Welling.  Visualizing deep neural network decisions: Prediction difference analysis.  _arXiv preprint arXiv:1702.04595_ , 2017.
  * Zouhar and Novák (2020) Zouhar, Vilém and Michal Novák.  Extending Ptakopět for Machine Translation User Interaction Experiments.  _The Prague Bulletin of Mathematical Linguistics_ , (115):129–142, 2020.


# Model fusion for efficient learning of nonlinear dynamical systems

Vatsal Kedia  Vivek S. Pinnamaraju  Dinesh Patil  Electrical Engineering
Department, Indian Institute of Technology Bombay, Mumbai, India (e-mail:
vatsalkedia@ ee.iitb.ac.in).  ABB Corporate Research, Bangalore, India
(e-mail: vivek.pinnamaraju@ in.abb.com)  ABB Corporate Research, Bangalore,
India
(e-mail: dinesh.patil@ in.abb.com)

###### Abstract

In the context of model-based control of industrial processes, it is a common
practice to develop a data-driven linear dynamical model around a specified
operating point. However, in applications involving wider operating
conditions, representation of the dynamics using a single linear dynamic model
is often inadequate, requiring either a nonlinear model or multiple linear
models to accommodate the nonlinear behaviour. While the development of the
former suffers from the requirements of extensive experiments spanning
multiple levels, significant compromise in the nominal product quality and
dealing with unmeasured disturbances over wider operating conditions, the
latter faces the challenge of model switch scheduling and inadequate
description of dynamics for the operating regions in-between. To overcome
these challenges, we propose an efficient approach to obtain a parsimonious
nonlinear dynamic model by developing multiple linear models from data at
multiple operating points, lifting the data features obtained from individual
model simulations to adequately accommodate the underlying nonlinear behaviour
and finally, sparse optimization techniques to obtain a parsimonious model.
The performance and effectiveness of the proposed algorithm is demonstrated
through simulation case studies.

###### keywords:

Nonlinear System Identification, Learning for Control, Advanced Process
Control

##  1 Introduction

Data-driven learning of linear dynamical systems is very popular owing to its
ease of interpretation, strong theoretical guarantees, mathematical
tractability and wide range of applications (Ljung (1998); Ljung (2010)). In
the context of process industries, linear dynamic models are easy to develop
as they need simple experiments involving use of step tests or pseudo random
binary signal (PRBS) excitation to generate informative data (Zhu (2001)). It
is well known that linear models provide accurate approximation around an
operating point but as we deviate away from the operating point the prediction
becomes worse since the linearization approximation no longer holds (Ren et
al. (2022)). Therefore linear models are not well suited for non-linear
systems due to less flexibility and hence are unable to capture the all
nonlinearities of the system.

While it is reasonable to develop linear models around a particular operating
point, a single linear model to capture the underlying dynamics is inadequate
in applications involving wider operating regions. In order to ensure good
control performance, it is necessary to describe the dynamics either using a
nonlinear model or multiple linear models. For directly identifying nonlinear
model we need input excitation covering the whole operating region of interest
to have informative data since step and/or PRBS tests are inadequate (Sriniwas
et al. (1995)). Hence we need to perform staircase or generalized multiple-
level noise (GMN) test (Zhu (2001)) which disturbs the normal plant operation
affecting the product quality and expensive due to prolonged test duration.
There are several popular nonlinear model structures which can accommodate the
nonlinear behavior such as NARX, NARMAX, Neural networks, etc. However,
learning the underlying non-linear dynamics using these models poses
significant practical challenges such as lack of parsimony, requirement of
extensive experiments, complexity in parameter estimation etc. (see Schoukens
and Ljung (2019) and references therein). Further, the neural networks
approach to model dynamical systems includes recurrent neural networks (RNNs),
LSTMs, etc. (see Pillonetto et al. (2023) and references therein for survey on
the topic) necessitate larger datasets compared to first-principles
(parsimonious) models. Moreover, these models often lack extrapolation
capabilities and can be challenging to interpret due to non-parsimonious model
structure. An additional limitation of these models lies in their mathematical
structure, which poses challenges when deriving methods for analysis of
closed-loop stability and performance (Nikolakopoulou et al. (2020)).

In this regard, multiple linear model framework to explain the nonlinear
behavior over wider operating region has gained attention (Johansen and Foss
(1999)). For more details on this approach, the readers are encouraged to
refer the review paper by Adeniran and El Ferik (2016) and references therein.
As the name suggests this approach involves developing different local linear
models at different operating points. The models used for control are switched
depending on the process operating condition. This approach faces several
challenges such as deciding on number of local linear models to be developed,
switching strategy for an intermediate operating region where none of the
individual models are adequate to explain the underlying dynamics etc.

Learning non-linear models poses several challenges as mentioned above. To
overcome these challenges, we propose an efficient model fusion algorithm. The
proposed approach utilizes simulated data from locally learned linear models
at different operating points, followed by lifting techniques to generate the
features necessary to explain the nonlinear behavior and finally, sparse
optimization techniques to obtain a parsimonious model. The key idea is while
each model simulation generates locally linear data, the fusion across
operating points essentially results in nonlinear data, which is appropriate
for nonlinear model identification. Since the non-linearity of the underlying
dynamic model is unknown, we postulate several non-linear features generated
from the simulated data using the ideas of lifting. While this results in
over-parameterized models, to obtain parsimonious models, we further use
sparse regularization on the model parameters to enable feature selection and
extract only the features that are most relevant in explaining the underlying
nonlinear behavior (Zou and Hastie (2005)). Furthermore, the models can be re-
identified using selected features by employing simple least squares to
enhance the parameter estimates and consequently the model predictions.

In Sun and Braatz (2020), the authors have also used ideas of lifting and
sparse optimization which has been applied on the actual nonlinear process
data. However this approach requires to deal with the effect of unmeasured
disturbances and measurement noise in the outputs. On the contrary our
approach differs in the following aspects: (i) Development of local linear
models across different operating points (sufficiently far) using traditional
input excitation signals. In the present work, these models are assumed to be
available. (ii) Simulation of noise-free data from multiple linear models
(iii) Model fusion of multiple linear models utilizing the transformed
simulation data and postulated nonlinear model, in combination with sparse
optimization techniques (iv) Model re-identification post feature selection to
enhance robustness. The major advantages of the proposed method in comparison
to existing works is as follows:

  * •

Simplicity of experiment design (Step/PRBS) for development of multiple linear
models.

  * •

Reliability on multi-model simulation data in comparison to actual data
enables effective rejection of measurement noise and unmeasured disturbance
while developing the nonlinear model.

  * •

Ease of nonlinear feature selection due to minimal noise in the features
unlike the actual data.

In the context of control, the closed-loop performance for a model-based
control strategy such as model predictive control (MPC) heavily depends on the
goodness of model predictions. Our proposed approach essentially develops a
reliable nonlinear model applicable to a wider operating region, which can be
easily deployed in a general nonlinear MPC framework. In a special case where
none of the selected features have any kind of non-linearity wrt. inputs, a
linear MPC can also be deployed.

To the best of the authors’ knowledge, model fusion to obtain a global
(applicable for relatively wider operating regions) nonlinear model using
multiple local linear models has not been studied and considered in the
literature. The applicability of the developed model is over a wider operating
region, including the intermediate transit regions for which a local linear
model is not available. As a result, significant experimental efforts in the
intermediate regions can be avoided for the sake of model development.

##  2 Preliminaries and Problem formulation

In this section, we review a variant of nonlinear auto-regressive with
exogenous input (NARX) model structure, namely, the polynomial NARX (p-NARX)
model. This model structure can accommodate mild nonlinear behaviour, which is
adequate in many of the process applications.

###  2.1 Polynomial NARX model (p-NARX)

The standard NARX model has the following form:

$y_{k}=f(y_{k-1},y_{k-2},\ldots,y_{k-l_{y}},u_{k-l_{d}},\ldots,u_{k-l_{u}-l_{d}})$
(1)

where $f(.)$ is some nonlinear function, while $l_{y}$, $l_{u}$ and $l_{d}$
represent the number of output lags (model order), input lags and the process
delay, respectively. Restricting $f(.)$ to multi-variable polynomial model
results in a linear-in-parameters form: $f(.)=\Psi(.)\alpha$, where $\Psi$
contains non-linear transformations of lagged inputs and outputs and $\alpha$
comprises of weights for all these features. This class of models are referred
to as p-NARX. For instance, if $\Psi(.)$ contains $\\{1,x,x^{2}\\}$ as the
polynomial basis functions, assuming $l_{u}=l_{y}=2$ and $l_{d}=0$, the
resulting p-NARX model is given by,

$\displaystyle y_{k}$
$\displaystyle=a_{1}y_{k-1}+a_{2}y_{k-2}+b_{1}u_{k-1}+b_{2}u_{k-2}+a_{11}y_{k-1}^{2}$
(2)
$\displaystyle+a_{22}y_{k-2}^{2}+b_{11}u_{k-1}^{2}+b_{22}u_{k-2}^{2}+a_{12}y_{k-1}y_{k-2}$
$\displaystyle+b_{12}u_{k-1}u_{k-2}+{ab}_{11}y_{k-1}u_{k-1}+{ab}_{22}y_{k-2}u_{k-2}$
$\displaystyle+{ab}_{12}y_{k-1}u_{k-2}+{ab}_{21}y_{k-2}u_{k-1}$

where $\alpha$ contains $a_{(.)}$ and $b_{(.)}$ coefficients.

###### Note 2.1

The non-linear transformation basis function, namely, $\Psi(.)$ is not
restricted to a polynomial basis. It can also accommodate other non-linear
basis functions such as transcendental (logarithmic, exponential) and
trigonometric functions. These can also be populated based on the underlying
process knowledge.

In this article, we restrict our focus to learning nonlinear polynomial NARX
model (p-NARX) as these are known to adequately accommodate process non-linear
behavior in wide applications. The basis function becomes
$\\{1,x,x^{2},\ldots,x^{q}\\}$, where $q\in\mathbb{Z}$. For instance, in
Sriniwas et al. (1995) the authors have modelled the distillation column as
p-NARX model considering upto $3^{rd}$-order interaction. Intuitively, for any
given nonlinear function $f(.)$, first-order approximation using Taylor’s
series around an operating point can explain only linear dynamics whereas as
we increase the order of approximation, the approximated function would be
able to explain non-linear behavior as well. It is to be noted that all
higher-order terms may not be necessary as we need only sufficient terms that
can explain the underlying nonlinear behaviour over the operating region of
interest. Despite truncation to a finite order, the number of terms for any
given order would be extremely high due to the different polynomial
combinations possible between the lagged inputs and outputs.

Since we do not have any information on actual nonlinear features, we need to
consider all potential features, which often results in over-fitting and
multicollinearity issues. To address these, we add additional regularization
terms on the parameters to the prediction error minization cost function while
learning the model. It is well known that $l_{1}$ regularizer (LASSO) promotes
sparsity-based solution, hence acts as feature selector (Tibshirani (1996))
while Ridge-regression ($l_{2}$ regularization) can take care of
multicollinearity in the data. Therefore, we use Elastic-net (Zou and Hastie
(2005)) as this can address the above-mentioned issues. Throughout this
article the following assumptions are made as described below.

###### Assumption 1

We assume the following:

  1. 1.

Mild non-linear behaviour over the operating region of interest.

  2. 2.

Local linear models are reliable locally around operating point.

The mild non-linear assumption is applicable to many real industrial processes
which may comprise of sub-units such as tanks, reactors, columns etc. (Xavier
et al. (2024)). In all process industries, the experiments are conducted
locally around some operating point using step or PRBS test. This leads to
reliable linear models locally around that operating point (see Assumption 1)
under assumption of sufficient number of input-output samples.

###  2.2 Problem formulation

Let us assume we have $l$-local linear models at distinct operating points
satisfying Assumption 1. Further let the local linear model corresponding to
$i^{th}$ operating point be denoted as $M_{i}\hskip 5.69046pt\forall
i\in\\{1,2,\ldots,l\\}$. Let $Y$ be vector containing outputs at various
operating points, similarly let $Z$ be the input regressor matrix (containing
linear and nonlinear features) and $\beta$ be the vector containing unknown
co-efficients.

$Y:=\begin{bmatrix}Y^{(1)}\\\ Y^{(2)}\\\ \vdots\\\ Y^{(l)}\end{bmatrix}\hskip
14.22636ptZ:=\begin{bmatrix}Z^{(1)}\\\ Z^{(2)}\\\ \vdots\\\
Z^{(l)}\end{bmatrix}$ (3)

Then the problem can be formulated as linear regression,

$\min_{\beta}\big{\lVert}Y-Z\beta\big{\rVert}_{2}^{2}+\lambda\bigg{(}\gamma\lVert\beta\rVert_{1}+\frac{1-\gamma}{2}\lVert\beta\rVert_{2}^{2}\bigg{)}$
(4)

where $\gamma\in(0,1)$ specifying the trade-off between LASSO and Ridge
regression.

###### Note 2.2

The local linear models are either already available (historical model) if not
then step or PRBS tests can be performed on the real process plant to obtain
these models.

###### Problem 2.3

Efficient learning of nonlinear dynamic process (solving (4)) given linear
models at different operating points.

The efficiency in above statement means that we do not need to perform any
additional experiment on the actual process plant (ie. without disturbing the
process) to collect the input-output data.

##  3 Model Fusion: Learning nonlinear dynamics

In this section we propose a novel model fusion algorithm to estimates p-NARX
model. The p-NARX model is an polynomial model wherein we have linear as well
well as non-linear features. The major steps of the proposed algorithm is
discussed in the next subsections.

###  3.1 Linear features

The local linear models are simulated to obtain noise-free input-output data
corresponding to each operating point. Let the simulated input-output data for
$i^{th}$ operating point is denoted by
$\\{u_{j}^{(i)},y_{j}^{(i)}\\}_{j=1}^{N}$ where $N$ is total number of
samples, $u_{j}\in\mathbb{R}^{m}$ and $y_{j}\in\mathbb{R}^{p}$. Since the
input-output lags is unknown therefore we assume $n_{y}>l_{y}$ and
$n_{u}>l_{u}$. The total number of linear features is denoted by
$n_{l}:=n_{y}+n_{u}$. For simplicity of understanding we assume $m=p=1$ ie.
SISO case which can be easily extended for MIMO case.

###  3.2 Appending non-linear features and feature selection

In this article we have considered polynomial non-linearity (p-NARX model) so
we make use of linear features obtained from simulated linear models to
generate the polynomial non-linear features. Using basis function $\Psi$, the
linear features can be transformed into non-linear features. Suppose we are
interested upto quadratic non-linearity, then the number of possible non-
linear features will be $n_{nl}=\frac{n_{l}(n_{l}+1)}{2}$. The regression data
matrix (Z) has $n_{tot}=n_{l}+n_{nl}=\frac{n_{l}^{2}+3n_{l}}{2}$ features.
Therefore the regression equation for $i^{th}$ operating point is given by:

$\displaystyle Y^{(i)}$ $\displaystyle=Z^{(i)}\beta$
$\displaystyle=\begin{bmatrix}Z_{l}^{(i)}&Z_{nl}^{(i)}\end{bmatrix}\beta$ (5)

where $Y^{(i)}\in\mathbb{R}^{N}$ is output data matrix corresponding to
$i^{th}$ operating point, $Z^{(i)}\in\mathbb{R}^{N\times n_{tot}}$ is
regressor matrix, $Z^{(i)}_{l}\in\mathbb{R}^{N\times n_{l}}$ is linear
regressor matrix, $Z^{(i)}_{nl}\in\mathbb{R}^{N\times n_{l}}$ is non-linear
regressor matrix and $\beta\in\mathbb{R}^{n_{tot}}$ is co-efficient matrix.
Stacking together for all operating points we get:

$\underbrace{\begin{bmatrix}Y^{(1)}\\\ Y^{(2)}\\\ \vdots\\\
Y^{(l)}\end{bmatrix}}_{=:Y}=\underbrace{\begin{bmatrix}Z^{(1)}\\\ Z^{(2)}\\\
\vdots\\\ Z^{(l)}\end{bmatrix}}_{=:Z}\beta$ (6)

In compact form we can write above equation as $Y=Z\beta$ where
$Y\in\mathbb{R}^{lN}$ and $Z\in\mathbb{R}^{lN\times n_{tot}}$. Since the
features affecting the dynamics of the system is unknown so we have taken all
possible combination of the features which may lead to over-fitting and data
colinearity issues. To discover the important and independent features we use
Elastic net as regularization method. Then the problem can be formulated as
linear regression,

$\min_{\beta}\big{\lVert}Y-Z\beta\big{\rVert}_{2}^{2}+\lambda\bigg{(}\gamma\lVert\beta\rVert_{1}+\frac{1-\gamma}{2}\lVert\beta\rVert_{2}^{2}\bigg{)}$
(7)

where $\gamma\in(0,1)$ specifying the trade-off between LASSO and Ridge
regression. Let $n_{s}$ be the number of selected features by solving the
above problem and let $\hat{\beta_{s}}\in\mathbb{R}^{n_{s}}$ denotes the co-
efficient of selected features. A straightforward way would be to use solution
of above equation $\hat{\beta}_{s}$ as final estimate. But it turns-out that
due to the regularization term, the actual importance of the co-efficient gets
re-distributed amongst other co-efficients (features) . Therefore it is good
practice to perform model re-identification post feature selection.

###  3.3 Model re-identification and validation

The next step is to perform model re-identification using selected features.
Let $Z_{f}\in\mathbb{R}^{lN\times n_{s}}$ denotes the final regessor matrix
corresponding to the selected features and $\beta_{f}\in\mathbb{R}^{n_{s}}$
denotes the final re-identified co-efficient vector corresponding to the
selected features. The model re-identification can be done by solving
following equation:

$\min_{\beta_{f}}\big{\lVert}Y-Z_{f}\beta_{f}\big{\rVert}_{2}^{2}$ (8)

The above problem is solved using p-NARX method to obtain the final model.
Next we perform model validation on test dataset. We use prediction error as
the validation metric given by:

$MSE=\frac{1}{N_{v}}\sum_{k=1}^{N_{v}}\big{(}y_{k}-\hat{y}_{k}\big{)}^{2}$ (9)

where $N_{v}$ is number of samples in test data, $y_{(.)}$ is the actual
output and $\hat{y}_{(.)}$ is the predicted output. The steps involved in the
proposed algorithm is presented in Algorithm 1.

Algorithm 1 Model fusion (p-NARX model)

Input: Input-output data at various operating points ie.
$\\{u_{j}^{(i)},y_{j}^{(i)}\\}_{j=1}^{N}$ $\forall i\in\\{1,2,\ldots,l\\}$.
Output: Global nonlinear p-NARX model.

  1. 1.

Choose $n_{y}>l_{y}$ and $n_{u}>l_{u}$.

  2. 2.

Generate linear features using given data with $n_{y}$ and $n_{u}$ lags and
formulate $Z_{l}$ matrix.

  3. 3.

Using linear features generate all possible polynomial (nonlinear) features
and formulate $Z_{nl}$ matrix.

  4. 4.

Append linear and nonlinear features together and choose hyperparameters
$\lambda$ and $\gamma$ for feature selection.

  5. 5.

Solve equation (7) for feature selection.

  6. 6.

Perform model re-identification to obtain p-NARX model using selected
features.

  7. 7.

Validate model using test data (unseen by the model during training).

##  4 Case Studies

In this section we verify the effectiveness of the proposed algorithm by
identifying several nonlinear dynamical systems. For simplicity we have
considered only SISO systems which can be easily extended for MIMO case. In
all cases, the comparison was done with model predictions from available
linear model at the nearest operating points.

###  4.1 Toy example

Let us consider nonlinear polynomial model with quadratic non-linearity given
by:

$\displaystyle y(k)$ $\displaystyle=0.5y(k-1)+0.25y(k-2)-0.5y^{2}(k-2)$
$\displaystyle+0.1y(k-1)u(k-1)+u(k-1)+0.25u(k-2)$

The local linear models were given at two ($l=2$) steady-sate operating points
namely $OP1-(0.1,0.3146)$ and $OP2-(0.3,0.6735)$. Let the models at given
operating points are denoted by $M_{1}$ and $M_{2}$ respectively. From Fig. 1,
it is clear that around operating point 1 $(OP1)$ ie. $0.3146$ model response
is similar to first-order while at operating point 2 $(OP2)$ ie. $0.6735$ the
model behaviour is similar to second-order dynamics. So the model is
exhibiting different behaviour at two different operating points. Intuitively,
using one model to predict behaviour at in-between operating point may lead to
inferior prediction. Hence we use our proposed algorithm (Algorithm 1) to
estimate a reliable nonlinear model.

Figure 1: Toy example: Different behaviour at two operating points

Since the model order is not known we chose $n_{u}=3>l_{u}$, $n_{y}=3>l_{y}$.
Therefore, total features $n_{tot}=27$. Total number of training samples for
each operating point is $N=448$ and number of test data samples $N_{v}=155$.
The regressor matrix is of size $Z\in\mathbb{R}^{698\times 27}$. The hyper-
parameters chosen for this case are $\gamma=0.5$ (equal weightage to address
over-fitting and colinearity) while $\lambda$ was chosen using 3-fold cross
validation (see (7)). The optimization problem given by equation (7), selected
nine features ie. $n_{s}=9$. The model was re-identified as global p-NARX
model using the selected features and then validated using test data (see
Algorithm 1). In order to show the efficiency of the proposed method we chose
an intermediate steady-state operating point ie. $(0.2,0.5214)$. The model
predictions were plotted for (a) given linear models at $OP1$ and $OP2$ (b)
global p-NARX model - proposed method. From Fig. 2 we can deduce that neither
of the linear models ($M_{1}$ and $M_{2}$) have good prediction capability at
intermediate operating point $(0.2,0.5214)$ while the proposed model fusion
approach works very well. The MSE at intermediate operating point for proposed
algorithm was found to be more than 60 times lesser than that of linear model
predictions (see Table 1). To further investigate the effectiveness of the
proposed method we compare the model predictions at various operating points.
From Table 1, it can be deduced that the proposed model fusion approach works
better than individual linear model predictions.

Figure 2: Comparison of model predictions at new operating point (us = 0.2)
Table 1: Performance comparison for toy example at various operating points

| OP
---
(us)
| MSE
---
(Model fusion)
| MSE
---
($M_{1}$)
| MSE
---
($M_{2}$)
| MSE ratio
---
($M_{1}$/$MF$)
| MSE ratio
---
($M_{2}$/$MF$)
0.05 | 4.237e-06 | 5.41e-04 | 0.0026 | 127.66 | 619.95
0.1 | 2.537e-06 | 3.286e-05 | 0.0012 | 12.95 | 49.29
0.2 | 9.839e-06 | 0.0015 | 6.65e-04 | 154.84 | 67.6
0.3 | 2.969e-05 | 0.0069 | 5.779e-04 | 232.96 | 19.46
0.35 | 4.855e-05 | 0.0115 | 5.77e-04 | 237.1 | 11.88

###  4.2 Conical tank process

Conical tanks find wide applications in various process industries, including
hydro-metallurgical, food processing, and wastewater treatment industries.
Moreover understanding the nonlinear dynamics for conical tank can provide
valuable insights for comprehending other types of nonlinear processes in
process industries.

Dynamics of nonlinear conical tank system is given by (Prasad et al. (2020)):

$\frac{dh}{dt}=\frac{\alpha(q_{i}-C_{d}\sqrt{h})}{h^{2}}$

where $\alpha=\frac{4H^{2}}{\pi D^{2}}$ is constant, $h$ denotes water level
(height), $q_{i}$ is inlet flow-rate, $C_{d}$ is valve constant, $D$ is
maximum diameter and $H$ is maximum height.

Figure 3: Conical tank: Different behaviour at two operating points

In this case the local linear models were given at two ($l=2$) operating
points which is steady-state height ie. $h_{s}=5cm$ and $h_{s}=10cm$. Similar
to previous case of toy example, the model behaviour is different at different
operating point (refer Fig. 3). At operating point 1 (OP1) ie. $h_{s}=5cm$ the
dynamics is faster as compared to OP2 ie. $h_{s}=10cm$ (see Fig. 3). Similar
process was followed as in toy example. The parameters used in this case are
$D=30cm$, $H=62cm$, $C_{d}=1$, $n_{u}=n_{y}=3$, $n_{tot}=27$, $N=7109$,
$N_{v}=3048$, $\alpha=0.75$ and $\lambda$ was chosen based on 5-fold cross-
validation. The optimization problem selected 12 features ie. $n_{s}=12$.
Model re-identification to obtain p-NARX model was done using selected
features and validation was performed using test data. In order to show the
efficiency of the proposed method we chose an intermediate operating point ie.
$h_{s}=7.5cm$. Model predictions were plotted for (a) linear models at given
operating point and (b) global p-NARX model fusion (proposed method). From
Fig. 4 we can deduce that neither of the linear models ($M_{1}$ and $M_{2}$)
have good prediction capability at intermediate operating point $h_{s}=7.5cm$
while the proposed model fusion approach works very well. From Table 3, MSE at
the intermediate operating point for the proposed algorithm is atleast 10
times better as compared to the individual linear models. To further
investigate the effectiveness of the proposed method we compare the model
predictions at various operating points (see Table 2), it can be deduced that
model fusion approach works better than individual linear model predictions.
Moreover, the performance at operating points where local models were given
was found to be better than that of the proposed method. This is intuitive
because linear models are very good approximation at operating points where it
is developed, hence better performance.

Figure 4: Comparison of model predictions at new operating point (7.5 cm)
Table 2: Performance comparison for conical tank at various operating points

| OP
---
($h_{s}$)
| MSE
---
(Model fusion)
| MSE
---
($M_{1}$)
| MSE
---
($M_{2}$)
| MSE ratio
---
($M_{1}$/$MF$)
| MSE ratio
---
($M_{2}$/$MF$)
4 | 0.0336 | 0.0495 | 0.4279 | 1.4732 | 12.7351
5 | 0.0061 | 5.36e-06 | 0.2561 | 8.79e-04 | 41.9836
7.5 | 0.0042 | 0.1140 | 0.0443 | 27.1429 | 10.5476
8.5 | 0.0058 | 0.1728 | 0.0144 | 29.7931 | 2.4828
10 | 0.0082 | 0.2477 | 1.143e-04 | 30.2073 | 0.0139
11 | 0.0147 | 0.2881 | 0.0027 | 19.5989 | 0.1837

###  4.3 Hammerstein-Wiener models

Hammerstein-Wiener models is represented by a linear transfer function and
capture the nonlinearities using nonlinear functions of inputs and outputs of
the linear system. The configuration looks like series connection of static
nonlinear blocks with a dynamic linear block. It has several applications such
as modeling electromechanical system (Paduart et al. (2012)) and radio
frequency components (Taringou et al. (2010)), predictive control of chemical
processes like Bio-reactors (Luo and Song (2018)), etc. Let us consider the
following example:

$\displaystyle s(t)=\frac{A(z)}{B(z)}x(t)+v(t)$ $\displaystyle
A(z)=1-0.45z^{-1}-0.35z^{-2}$ $\displaystyle B(z)=0.5z^{-1}-0.25z^{-2}$
$\displaystyle x(t)=f(u(t))=0.5u(t)-0.18u^{2}(t)$ $\displaystyle
y(t)=g(s(t))=\frac{1}{0.3}\bigg{(}-1+\sqrt{1+0.6s(t)}\bigg{)}$

where $f(.)$ and $g(.)$ are static nonlinear function acting on input ($u(t)$)
and intermediate output ($s(t)$) respectively. In this case, local linear
models were given at two ($l=2$) operating points ie. $M_{1}(u_{s}=0)$ and
$M_{2}(u_{s}=1)$. Similar analysis was performed as in previous examples with
$n_{u}=n_{y}=3$. From Fig. 5 it can be clearly deduced that individual linear
model predictions are worse than the proposed model fusion approach.

Figure 5: Comparison of model predictions at new operating point (us = 0.5)

##  5 Conclusion

We proposed a practical and efficient model fusion approach from multiple
linear models to obtain a global non-linear reliable model over a wider
operating region. The efficacy of the proposed approach has been demonstrated
in three simulation examples. The predictions of the developed nonlinear
models are observed to be close to the real process responses even for the
intermediate operating regions for which there is no local linear model
available without the requirement of additional experimental efforts. The
results show promising steps towards developing models in unknown transition
regions which is necessary in applications involving frequent random grade
changes over wider operating regions and also in ensuring bump-less transfer
while switching from one grade to the another. Future work considers an
extension of these ideas for multi-input multi-output systems, accounting for
variable delays and model fusion in state space framework.

## References

  * Adeniran and El Ferik (2016) Adeniran, A.A. and El Ferik, S. (2016).  Modeling and identification of nonlinear systems: A review of the multimodel approach—part 1.  _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , 47(7), 1149–1159.
  * Johansen and Foss (1999) Johansen, T.A. and Foss, B.A. (1999).  Multiple model approaches to modelling and control.  _International journal of control_ , 72(7-8), 575–575.
  * Ljung (1998) Ljung, L. (1998).  System identification.  In _Signal analysis and prediction_ , 163–173. Springer.
  * Ljung (2010) Ljung, L. (2010).  Perspectives on system identification.  _Annual Reviews in Control_ , 34(1), 1–12.
  * Luo and Song (2018) Luo, X.S. and Song, Y.D. (2018).  Data-driven predictive control of hammerstein–wiener systems based on subspace identification.  _Information Sciences_ , 422, 447–461.
  * Nikolakopoulou et al. (2020) Nikolakopoulou, A., Hong, M.S., and Braatz, R.D. (2020).  Feedback control of dynamic artificial neural networks using linear matrix inequalities.  In _2020 59th IEEE Conference on Decision and Control (CDC)_ , 2210–2215. IEEE.
  * Paduart et al. (2012) Paduart, J., Lauwers, L., Pintelon, R., and Schoukens, J. (2012).  Identification of a wiener–hammerstein system using the polynomial nonlinear state space approach.  _Control Engineering Practice_ , 20(11), 1133–1139.
  * Pillonetto et al. (2023) Pillonetto, G., Aravkin, A., Gedon, D., Ljung, L., Ribeiro, A.H., and Schön, T.B. (2023).  Deep networks for system identification: a survey.  _arXiv preprint arXiv:2301.12832_.
  * Prasad et al. (2020) Prasad, G.M., Kedia, V., and Rao, A.S. (2020).  Multi-model predictive control (mmpc) for non-linear systems with time delay: an experimental investigation.  In _2020 First IEEE International Conference on Measurement, Instrumentation, Control and Automation (ICMICA)_ , 1–5. IEEE.
  * Ren et al. (2022) Ren, Y.M., Alhajeri, M.S., Luo, J., Chen, S., Abdullah, F., Wu, Z., and Christofides, P.D. (2022).  A tutorial review of neural network modeling approaches for model predictive control.  _Computers & Chemical Engineering_, 107956.
  * Schoukens and Ljung (2019) Schoukens, J. and Ljung, L. (2019).  Nonlinear system identification: A user-oriented road map.  _IEEE Control Systems Magazine_ , 39(6), 28–99.
  * Sriniwas et al. (1995) Sriniwas, G.R., Arkun, Y., Chien, I.L., and Ogunnaike, B.A. (1995).  Nonlinear identification and control of a high-purity distillation column: a case study.  _Journal of Process Control_ , 5(3), 149–162.
  * Sun and Braatz (2020) Sun, W. and Braatz, R.D. (2020).  Alven: Algebraic learning via elastic net for static and dynamic nonlinear model identification.  _Computers & Chemical Engineering_, 143, 107103.
  * Taringou et al. (2010) Taringou, F., Hammi, O., Srinivasan, B., Malhame, R., and Ghannouchi, F.M. (2010).  Behaviour modelling of wideband rf transmitters using hammerstein–wiener models.  _IET circuits, devices & systems_, 4(4), 282–290.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 58(1), 267–288.
  * Xavier et al. (2024) Xavier, J., Patnaik, S., and Panda, R.C. (2024).  Nonlinear system identification in coherence with nonlinearity measure for dynamic physical systems—case studies.  _Nonlinear Dynamics_ , 1–27.
  * Zhu (2001) Zhu, Y. (2001).  _Multivariable system identification for process control_.  Elsevier.
  * Zou and Hastie (2005) Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the elastic net.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 67(2), 301–320.