# Correlation of plastic events with local structure in jammed packings across
spatial dimensions

Sean A. Ridout  Jason W. Rocks  Andrea J. Liu  Department of Physics and
Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

###### Abstract

In jammed packings, it is usually thought that local structure only plays a
significant role in specific regimes. The standard deviation of the relative
excess coordination, $\sigma_{Z}/Z_{\mathrm{c}}$, decays like $1/\sqrt{d}$, so
that local structure should play no role in high spatial dimensions.
Furthermore, in any fixed dimension $d\geq 2$, there are diverging length
scales as the pressure vanishes approaching the unjamming transition, again
suggesting that local structure should not be sufficient to describe response.
Here we challenge the assumption that local structure does not matter in these
cases. In simulations of jammed packings under athermal, quasistatic shear, we
use machine learning to identify a local structural variable, softness, that
correlates with rearrangements in dimensions $d=2$ to $d=5$. We find that
softness - and even just the coordination number $Z$ \- are quite predictive
of rearrangements over a wide range of pressures, all the way down to
unjamming, in all $d$ studied. This result provides direct evidence that local
structure can play a role in higher spatial dimensions.

Supercooled liquids, glasses and jammed packings can relax via rearrangements,
in which particles experience large, sudden displacements, due to thermal
fluctuations or mechanical loading. In such systems in $d=2$ and $d=3$, where
$d$ is the spatial dimension, the details of local structure play an important
role, an idea supported by numerous studies that correlate local structure
with rearrangements [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. In
particular, machine learning has been used to identify a highly predictive
combination of local structural descriptors, named “softness,” which provides
insight into the underlying physics [15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 10]. For example, in thermal systems, each value of this continuous
structural variable is associated with a characteristic free energy barrier to
rearrangements that decreases with increasing softness [16, 18, 20, 21, 22,
23, 24].

On the other hand, a successful analytical approach to these systems has been
found in the form of a mean-field theory. These mean-field calculations, valid
in the limit $d=\infty$, seem not to include details of local structure, and
activated processes - which we associate with softness - are absent because
their free energy barriers become infinite. In spite of these simplifications,
these mean-field models predict many glass transition phenomena [26] and
Gardner transition phenomena [27] with qualitative accuracy, even in $d=2,3$,
and quantitatively predict aspects of jamming criticality in $d=2,3$ [26].

The successes of these two frameworks seem to be in tension with one another;
if the details of the local structure are empirically important in low
dimensions, how can mean-field theories which ignore them be so successful?
One pathway to reconciling these facts is to look for a smooth transition to
mean-field-like behaviour as the dimension is raised. It is useful to address
this question in the context of athermal, quasistatic shear to avoid the
existence of multiple relevant time scales governing rearrangements that are
inevitable at nonzero temperatures $T$. At $T=0$, there is reason to believe
that localized rearrangements correlated with local structural variables give
way with increasing $d$ to delocalized rearrangements that are poorly
correlated with local structure. At $T=0$, rattlers and bucklers, two types of
local structures associated with instabilities (rearrangements) as well as
deviation from mean field predictions in jammed packings, grow exponentially
more rare with increasing $d$ [28, 29]. Furthermore, as $d$ rises, a larger
fraction of the low-frequency vibrational modes become extended, and the
density of vibrational states smoothly approaches the form expected from mean
field theory, although in any finite dimension there still exist some
quasilocalized modes at sufficiently low frequency [30, 31]. These facts seem
to imply that local structure does matter in low dimensions, but becomes less
and less relevant as $d$ increases, smoothly converging to some “mean-field-
like” dynamics analogously to the density of states and distribution of
forces.

In light of these ideas, this paper quantifies directly the importance of
local structure as a function of spatial dimension. Here we examine the
prediction accuracy of machine-learned softness in jammed soft sphere packings
under athermal, quasistatic shear, in dimensions $d=2$ to $d=5$. In contrast
to the predominant belief that local structure only matters in low dimensions,
we find that as $d$ increases, local structure as quantified by softness
becomes more predictive of rearrangements. Moreover, we find that softness
increasingly coincides with simply the local number of interacting neighbors,
or coordination number, $Z$. This result suggests that the distribution of $Z$
values plays an important role and that $Z$ could be a valuable indicator of
mobility even in the mean field limit.

We prepare bi-disperse jammed packings of Hertzian particles at a range of
pressures in each spatial dimension. We study $N=4096$ particles in $2d$,
$N=8192$ in $3d$ and $4d$, and $N=16384$ in $5d$. Initially, particles are
assigned random positions ($T=\infty$), before quenching to a energy minimum,
producing a force-balanced, zero-temperature (jammed) state. The system is
sheared quasistatically by applying small strain steps
($\delta\gamma=10^{-4}$) and re-minimizing the energy $E$ with respect to
$\left\\{x_{i}\right\\}$ at each step. When the pressure $P$ is large, it is
relatively insensitive to shear
($\frac{1}{P}\partialderivative{P}{\gamma}\approx 0$), allowing us to minimize
with a fixed box volume $V$. However, when $P$ is small and the system is
finite (and thus slightly anisotropic), the pressure can fluctuate
substantially with shear, so we fix $P$ by minimizing the enthalpy, $H=E+PV$,
with respect to both $\left\\{x_{i}\right\\}$ and volume.When a plastic event
is detected in the form of a drop in the shear stress, the resulting state is
discarded, and the stress drop is approached with a shear step
$\delta\gamma/2$ [32, 33]. We repeat this procedure until the stress drop is
approached with a shear step of $\delta\gamma=10^{-12}$, at which point the
main numerical limitation to realizing true quasistatic shear is the force
tolerance in the energy minimization algorithm.

Next, we quantify local structure, by finding a combination of local
structural descriptors captures each particle’s propensity to rearrange. In
previous studies, populations of rearranging and non-rearranging particles
have been identified from a collection of configurations and then used to
train classifier to sort particles into these two groups [15]. Here we instead
utilize a more conceptually straightforward approach in which linear
regression is used to identify structures that correlate with local
fluctuations in the displacement field [25]. In addition to providing more
concise sets of structural descriptors, this method has the added advantage of
requiring far less data to construct a high quality training set as it can use
all particles in a data set, rather than a small subset of examples.

In this approach, we describe the local structure around each particle by
counting the number of contacts and gaps at each distance in a triangulation
of the packing. Once the Delaunay triangulation of the packing is constructed,
each particle $j$ in the triangulation is assigned a discrete distance
$d_{ij}$ relative to particle $i$ such that $d_{ii}=0$, $d_{ij}=1$ for
neighbours, $d_{ij}=2$ for particles which share a neighbor, etc. Each edge
$\left(j,k\right)$ in the triangulation is either a contact between particles
or a gap and may be assigned a distance from particle $i$, defined as
$d_{i,\left(j,k\right)}=d_{ij}+d_{jk}$. The number of gaps and contacts at
each distance up to $d_{i,\left(j,k\right)}=8$, together with the particle
radius, form a set of 17 structural descriptors that we use to encode the
local structural environment around each particle.

One complicating factor is that quasistatically sheared jammed packings
exhibit crackling noise in which one localized rearrangement triggers a
second, which triggers a third, and so forth, producing an avalanche at each
single stress drop or instability [34]. To focus on localized rearrangements
rather than avalanches, we follow what has become standard practice [35, 3,
10], designing our structural variable to predict the first rearrangement in
each avalanche. We capture this initial rearrangement by computing the lowest
eigenvector of the dynamical matrix immediately before the stress drops. This
eigenvector corresponds to the normal mode whose frequency vanishes at the
onset of the rearrangement, signaling the instability. Using this eigenvector,
we compute the quantity $D^{2}_{\min}$ for each particle, measuring its non-
affine displacement relative to its neighbours [36, 35].

Figure 1: The $S$ percentile $C_{S}$ (A) and $-Z$ percentile $C_{Z}$ (B) of
the particle with the highest $\Delta^{2}_{\min}$ in a plastic event, as a
function of dimension $d$ and pressure $P$ ($d=2$ black triangles, $d=3$ blue
squares, $d=4$ green pentagons, $d=5$ red hexagons). $C_{Z}$ clearly peaks at
a finite pressure, and decays at large pressure. $C_{S}$ appears to have a
weak pressure dependence which is difficult to resolve in our data. Except at
high pressure, $Z$ performs almost as well as $S$. Error bars are uncertainty
in mean value. Inset to (A) $R^{2}$ for fitting of $S$ to the local
displacement fluctuations $\Delta^{2}_{\min}$. Error bars show full range over
all pressures; all data are shown in the Supporting Information [37]. Although
the rearranging particle is predicted equally well in higher $d$, the
correlation to the rescaled mobility field used to fit $S$ is actually much
worse in higher $d$,. Inset to (B): Peak value over all pressures $C^{*}_{Z}$,
showing little dimension dependence.

Because the system is solid, a rearranging particle exerts a long-ranged
strain field on other particles in the system, causing a power-law decay with
distance of $D^{2}_{\min}$ averaged over directions. A particle near the
rearrangement core will move more than one far from it, regardless of whether
their respective structures predispose them to mobility or to rigidity. To
account for this effect, we rescale $D^{2}_{\min}$ for each particle by
dividing by the average $D^{2}_{\min}$ of its neighbors [25]; we denote this
quantity as $\Delta^{2}_{\min}$ and expect that it should better correlate
with local structure than $D^{2}_{\min}$.

Next, we perform linear regression to identify a linear combination of our
structural descriptors, the “softness” $S$, which best correlates with
$\Delta^{2}_{\min}$. Although the correlation coefficient $R^{2}$ is not high,
and becomes weaker in higher dimensions, we find that the prediction accuracy
on rearranging particles is just as high as previous classification-based
approaches [25]. This quantification of prediction accuracy is superior to
comparing classification accuracies on a training set or regression accuracies
because these numbers change if the definition of the training set is changed,
even if the resulting softness $S$ does not change. An additional issue is
that, as described above, our training variable $\Delta^{2}_{\min}$ does not
even directly correspond to rearranging and non-rearranging particles so it
does not make sense to calculate the prediction accuracy for rearranging
particles. An alternative approach might be to set a threshold on
$D^{2}_{\min}$ so that particles above this threshold are “rearranging”, and
then to ask which fraction of these particles have high softness [24].
However, this means of quantifying the prediction accuracy of softness is
problematic in our case because a different threshold $D^{2}_{\min}$ value
would have to be selected in different dimensions and pressures, making
meaningful comparisons across pressures and dimensionalities difficult.

Once training is complete, each particle in the system is assigned a value of
$S$. To evaluate the predictive power of $S$, or any structural indicator of
rearrangements, a standard approach is to check whether the particle most
associated with the center of each rearrangement also has a high value of the
softness $S$ [38, 10]. We identify this particle as that with the maximum
value of $D^{2}_{\min}$ in the critical mode. We them measure the average
percentile of $S$ for rearranging particles, which we call $C_{S}$, as our
measure of correlation between structure and dynamics. A value of $C=1$
corresponds to perfect prediction accuracy while $C=0.5$ corresponds to random
guessing.

In Fig. 1(A), we report $C_{S}$ for packings in all dimensions as a function
of pressure $P$. We compute error bars as an estimate of the uncertainty in
our measures due to sampling error calculated via cross-validation [25]. We
find that within uncertainty the predictive power of $S$ does not strongly
depend on either pressure nor dimension. At lower pressures in $d=5$, $C_{S}$
may be slightly less predictive, but it is otherwise difficult to see any
significant trends.

Although Fig. 1 shows that the local structure can predict the location of the
maximum of $D^{2}_{\min}$ equally well in all dimensions studied, we note that
softness is actually less predictive of the local fluctuations in
$D^{2}_{\min}$ which we use to train $S$ at higher $d$. The $R^{2}$ of the
linear regression, shown in the inset of the first panel, decays dramatically
$2$ and $5$. Thus, in higher dimensions, the local fluctuations in
$D^{2}_{\min}$ are actually very poorly correlated with the local structure,
but what little correlation remains is enough to identify a structural
variable which correlates strongly with the rearrangement itself.

Furthermore, we find that most of the predictive power of $C_{S}$ is
determined by the number of contacts at distance $d_{i,(j,k)}=1$, i.e. the
coordination number $Z$ of each particle. As a comparison, we compute the
predictive power of the coordination number $Z$ by computing the average
percentile of $-Z$ , which we denote as $C_{Z}$. In Fig. 1(B), we report
$C_{Z}$ for the same set of particles used to measure $C_{S}$. Except at high
pressures, we see that $C_{Z}$ is comparable to $C_{S}$; in $4d$ and $5d$ they
are equal within the uncertainty. Thus, the coordination in fact contains most
of the local structural information we are capable of finding using our method
– in high dimension it contains the only available information.

In order to more easily compare the properties of our structure and dynamics
across dimensions, we pick a single pressure to represent each dimension.
Comparing systems at the same pressure is not necessarily meaningful; pressure
has different units in different spatial dimensions. Alternative control
parameters include the difference in packing fraction from the jamming
transition $\Delta\phi$, the excess coordination number $\Delta Z$ needed for
rigidity, and the “prestress” $e\sim\Delta Z/d$ which has recently been
suggested as an appropriate control parameter when comparing between
dimensions [31]. However, we find that none of these choices results in a
clean collapse of the curves for $C_{Z}$ in Fig. 1(B). Therefore, whenever we
need to compare a property of $S$ or $Z$ across spatial dimension, we choose
the pressure $P^{*}(d)$ in each dimension which maximizes $C_{Z}$. In the
inset of Fig. 1(B), we report the predictive power of $Z$ at $P^{*}(d)$, which
we denote $C^{*}_{Z}$. We find that $C^{*}_{Z}$ appears to be approximately
constant between $d=3,4,5$ within sampling error.

Surprisingly, neither $C_{S}$ nor $C_{Z}$ shows a strong decrease at low
pressure. It is known that various length scales diverge at the jamming
transition, including the length scale corresponding to spatial correlations
of $Z$, as well as the size of a localized “core” where the response to a
perturbation deviates from linear elasticity [39, 40, 41]. Thus, one might
expect an approach based on simple local structural descriptors to decrease in
performance as the pressure approaches zero, or at least for the number of
necessary descriptors to increase strongly. Moreover, we find that $C_{S}$
shows no noticeable finite-size scaling at low pressures, further indicating
that the slight decrease in accuracy at low pressures is unlikely due to a
diverging correlation length [42] (see Supporting Information [37]).

Figure 2: Correlation lengths extracted from fits to initial decay of small-
small-particle correlation functions out to a distance of about $3$ particle
diameters, calculated at the pressure $P^{*}(d)$ . All correlation lengths
decrease slightly with increasing dimension $d$. Red squares: $D^{2}_{\min}$,
blue triangles: $\Delta^{2}_{\min}$, black circles: $S$.

In contrast, it is not suprising that structural predictors other than $Z$ do
not contribute to $S$ in higher $d$. In infinite dimensions, it suffices for
an exact treatment to truncate the virial expansion at second order, i.e. to
only include the effects of nearest interacting neighbors [43]. In high $d$,
each particle’s nearest neighbours are increasingly likely to also be
neighbours of one another, so structural descriptors other than particle type,
number of contacts, and the number of nearest-neighbour gaps should contain no
information in the limit $d\to\infty$.

The most interesting aspect of our findings is that rearrangements continue to
localize around low-$Z$ particles as the $d$ is increased, as evidenced by
$C_{\mathrm{Z}}$. Because the low-frequency vibrational at zero strain modes
become more extended in higher $d$, one might expect the initial rearrangement
would become less localized in higher dimensions, leading to a decrease of
$C_{\mathrm{Z}}$ [30, 31]. To determine the extent of localziation of the
initial rearrangement, we calculate the spatial correlations of $D^{2}_{\min}$
and $\Delta^{2}_{\min}$, along with the spatial correlations of softness $S$.
At long distances, the $D^{2}_{\min}$ of each rearrangement displays power-law
decay, consistent with continuum elasticity, but at short distances there is
an exponential decay in $D^{2}_{\min}$ correlations on the scale of the
particle diameter [37]. Using this decay, we fit an exponential function to
calculate an approximate correlation length. We do the same for the
correlations in $\Delta^{2}_{\min}$, for which small spatial gradients are
scaled out, and for the softness $S$.

Fig. 2 shows these correlation lengths evaluated at $P^{*}(d)$ as a function
of $d$. We observe that the correlation lengths all decrease with dimension,
indicating that rearrangements become more localized in space. Note, however,
that there are more nearby particles with increasing $d$ so the rearrangements
may be more extended in terms of the number of particles involved (e.g.,
higher participation ratio) even though they are more spatially localized. In
the Supporting Information, we show that the structural correlation lengths,
as well as that for $\Delta^{2}_{\min}$, depend only weakly on pressure in
each dimension, increasing slightly at higher pressures. On the other hand,
the $D^{2}_{\min}$ correlation length, grows as pressure is decreased as one
might expect, although it does not seem to diverge at unjamming [37]. In the
Supporting Information, we also provide the inverse participation ratio of the
critical modes for direct comparison with previous reports at zero strain [30,
31]. In our critical modes we do not see the same apparent delocalization in
IPR which was reported at zero strain, although we do not hold N fixed.

Figure 3: Probability of rearranging for small particles as a function of
$(Z-\langle Z\rangle)/\sigma_{Z}$ in varying spatial dimension ($d=2$ black
triangles, $d=3$ blue squares, $d=4$ green pentagons, $d=5$ red hexagons). As
dimension is varied, or pressure is varied (at pressures low enough for $Z$ to
still be strongly predictive), the dependence of the rearrangement probability
on this rescaled $Z$ is qualitatively similar, indicating that particles with
relatively low $Z$ are always more likely to rearrange, even in higher $d$
where they are not bucklers. Inset shows $\langle Z\rangle$ (black x) and
$\sigma_{Z}$ (blue circles) for small particles, indicating that at $P^{*}$
(where $\sigma_{Z}/\langle Z\rangle$ is smaller than at $P=0$) the asymptotic
$1/\sqrt{d}$ scaling of $\sigma_{Z}/\langle Z\rangle$ has yet to occur.

In Fig. 2, we also observe that the spatial correlations of
$\Delta^{2}_{\min}$ are the same as those of $S$, indicating that structural
correlations set the size of rearrangements. In earlier work [44], it was
concluded that the spatial correlations of $D^{2}_{\min}$ are the same as
those of $S$; here we find that spatial correlations of $D^{2}_{\min}$ are
noticeably longer in range in $d=2,3$. However, the systems studied in Ref.
[44] did not exhibit elastic correlations over as long a range as we find here
as a result of friction or thermal fluctuations. Therefore, it is likely that
the correlation lengths of $D^{2}_{\min}$ and $\Delta^{2}_{\min}$ in those
systems would be far more similar.

It has previously been noted that “bucklers”, particles with the minimum
number of contacts for local stability, are associated with deviations from
mean-field behaviour and become less and less common in higher-dimensional
packings [29]. Naively, one might imagine that this would result in a decrease
in $C_{Z}^{*}$ with dimension. To investigate this hypothesis, in Fig. 3 we
report the probability of the rearrangement being located at a given particle
as a function of $Z$ in each dimension at pressure $P^{*}(d)$. We collapse the
collapse the curves by multiplying probability by the total number of
particles $N$ and standardizing $Z$ by subtracting off its mean over all
particles $\expectationvalue*{Z}$ and dividing by its standard deviation
$\sigma_{Z}$. We find that particles that rearrange tend to have small $Z$
relative to the average, often by two or more standard deviations. This result
holds even as the $d$ increases and rattlers and bucklers become exceedingly
rare, demonstrating that it is a particle’s value of $Z$ relative to the rest
of the configuration that determines its propensity to rearrange, rather than
solely its status as a buckler or a rattler.

In the inset to Fig. 3, we show the scaling of $\langle Z\rangle$ and
$\sigma_{Z}$ evaluated at $P^{*}(d)$ for each $d$ studied. We see that
$\langle Z\rangle(P^{*}(d))$ scales linearly with $d$ as expected since
$\langle Z\rangle(P=0)=2d$. However, we find $\sigma_{Z}(P^{*}(d))\sim d$,
while we would expect $\sigma_{Z}(P=0)\sim\sqrt{d}$ [29]. We find that
$\sigma_{Z}/\langle Z\rangle$ decreases as a function of $P$, and thus we
would expect $\sigma_{Z}{\left(P^{*}\right)}\sim\sqrt{d}$ at large enough $d$
(since it is bounded above by $P=0$). Thus, although the collapse in Fig. 3
seems to indicate that lower $Z$ particles should dominate rearrangements in
all $d$, there is the caveat that $\sigma_{Z}$ has not yet reached its
asymptotic scaling in $d=5$.

In summary, our results suggest a refinement of the picture relating local
structure to rearrangements in high dimensions. Local structure remains
important in determining particle rearrangements. In fact, the spatial size of
rearrangements, as measured by the decay of $D^{2}_{\min}$, is somewhat
shorter in range. Nevertheless, the number of particles within that range
increases exponentially with $d$, so that rearrangements may become extended
in terms of the number of particles that participate. In high $d$, the contact
number $Z$ for each particle becomes a good structural predictor of
rearrangements, implying that rearrangements are controlled only by nearest
neighbors that interact directly with a particle. Even though the contact
number distribution narrows relative to the mean contact number $\langle
Z\rangle\sim d$, the standard deviation diverges, scaling as
$\sigma_{Z}\sim\sqrt{d}$ [29]. Thus, different particles can still have
different $Z$ in high dimensions as $d\rightarrow\infty$, and low $Z$
particles may still have a higher propensity to rearrange.

Defining $x=2Z/\langle Z\rangle\sim Z/d$, there are three simple scenarios
that could prevail in high $d$: (1)
$P_{\mathrm{R}}{\left(Z\right)}=f{\left(x\right)}$. (2) Just as the
distribution of contacts has a large-deviation form $P{\left(Z\right)}\sim
e^{-dB\left(x\right)}$, one might expect $P_{\mathrm{R}}{\left(Z\right)}\sim
e^{dA{\left(x\right)}}$ [29]. (3) The scenario of Fig. 3, where
$P_{\mathrm{R}}{\left(Z\right)}=f{\left((Z-\langle
Z)\rangle/\sigma_{Z}\right)}$. In the first scenario, the typical $x^{*}$ for
rearranging particles goes to $1$ in and $C_{Z}\to\frac{1}{2}$ (no
correlation) as $d\rightarrow\infty$. In the second, $1<x^{*}<2$: rearranging
particles in large dimension are not bucklers, but are still extremely
atypical, and $C_{\mathrm{Z}}\to 1$. In the third scenario, $x^{*}\to 2$ as
$d\to\infty$, but $2-x^{*}$ goes to zero at the same rate as
$\sigma_{\mathrm{Z}}$, so rearranging particles are still atypical and
$C_{\mathrm{Z}}$ is constant. We find the large deviation form of scenario (2)
does not collapse our data as well as in Fig. 3, but this could be because
$d=5$ is still low.

Either scenario (2) or (3) would connect the softness picture of dynamical
behavior associated with the glass transition in $d=3$ with mean field theory.

## Acknowledgements

We thank Ge Zhang, Sam Schoenholz, Robert Ivancic, François Landes, Eric
Corwin and Francesco Zamponi for helpful discussions. We thank Carl Goodrich
for providing useful code.

This research was supported by the US Department of Energy, Office of Basic
Energy Sciences, Division of Materials Sciences and Engineering under Award
DE-FG02-05ER46199 (J.W.R), the Natural Sciences and Engineering Research
Council of Canada under a Postgraduate Scholarship – Doctoral award (S.A.R),
and the Simons Foundation for the collaboration Cracking the Glass Problem via
award 454945 (J.W.R, S.A.R, and A.J.L) and Investigator Award 327939 (A.J.L.).

## References

  * [1] Asaph Widmer-Cooper, Heidi Perry, Peter Harrowell, and David R. Reichman.  Irreversible reorganization in a supercooled liquid originates from localized soft modes.  Nature Physics, 4(9):711–715, 2008.
  * [2] Asaph Widmer-Cooper, Heidi Perry, Peter Harrowell, and David R. Reichman.  Localized soft modes and the supercooled liquid’s irreversible passage through its configuration space.  The Journal of Chemical Physics, 131(19):194508, 2009.
  * [3] M. L. Manning and A. J. Liu.  Vibrational modes identify soft spots in a sheared disordered packing.  Physical Review Letters, 107(10):108302, aug 2011.
  * [4] S. S. Schoenholz, A. J. Liu, R. A. Riggleman, and J. Rottler.  Understanding plastic deformation in thermal glasses from single-soft-spot dynamics.  Phys. Rev. X, 4:031014, Jul 2014.
  * [5] Michel Tsamados, Anne Tanguy, Chay Goldenberg, and Jean-Louis Barrat.  Local elasticity map and plasticity in a model lennard-jones glass.  Phys. Rev. E, 80:026112, Aug 2009.
  * [6] Luka Gartner and Edan Lerner.  Nonlinear plastic modes in disordered solids.  Phys. Rev. E, 93:011001, Jan 2016.
  * [7] Jacques Zylberg, Edan Lerner, Yohai Bar-Sinai, and Eran Bouchbinder.  Local thermal energy as a structural indicator in glasses.  Proceedings of the National Academy of Sciences, 114(28):7289–7294, 2017.
  * [8] Alex Malins, Jens Eggers, C. Patrick Royall, Stephen R. Williams, and Hajime Tanaka.  Identification of long-lived clusters and their link to slow dynamics in a model glass former.  Journal of Chemical Physics, 138(12), March 2013.
  * [9] Hua Tong and Hajime Tanaka.  Revealing hidden structural order controlling both fast and slow glassy dynamics in supercooled liquids.  Phys. Rev. X, 8:011041, Mar 2018.
  * [10] D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J. L. Barrat, L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens, S. Sastry, D. Vandembroucq, E. Lerner, and M. L. Manning.  Predicting plasticity in disordered solids from structural indicators, 2020.
  * [11] Susana Marín-Aguilar, Henricus H. Wensink, Giuseppe Foffi, and Frank Smallenburg.  Tetrahedrality dictates dynamics in hard sphere mixtures.  Phys. Rev. Lett., 124:208005, May 2020.
  * [12] V. Bapst, T. Keck, A. Grabska-Barwińska, C. Donner, E. D. Cubuk, S. S. Schoenholz, A. Obika, A. W.R. Nelson, T. Back, D. Hassabis, and P. Kohli.  Unveiling the predictive power of static structure in glassy systems.  Nature Physics, 16(4):448–454, 2020.
  * [13] Emanuele Boattini, Susana Marín-Aguilar, Saheli Mitra, Giuseppe Foffi, Frank Smallenburg, and Laura Filion.  Autonomously revealing hidden local structures in supercooled liquids, 2020.
  * [14] Leevi Viitanen, Jonatan R. Mac Intyre, Juha Koivisto, Antti Puisto, and Mikko Alava.  Machine learning and predicting the time dependent dynamics of local yielding in dry foams, 2020.
  * [15] E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu.  Identifying structural flow defects in disordered solids using machine-learning methods.  Physical Review Letters, 114(10):108001, mar 2015.
  * [16] S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu.  A structural approach to relaxation in glassy liquids.  Nature Physics, 12(5):469–471, feb 2016.
  * [17] Daniel M. Sussman, Samuel S. Schoenholz, Ekin D. Cubuk, and Andrea J. Liu.  Disconnecting structure and dynamics in glassy thin films.  Proceedings of the National Academy of Sciences of the United States of America, 114(40):10601–10605, 2017.
  * [18] Tristan A. Sharp, Spencer L. Thomas, Ekin D. Cubuk, Samuel S. Schoenholz, David J. Srolovitz, and Andrea J. Liu.  Machine learning determination of atomic dynamics at grain boundaries.  Proceedings of the National Academy of Sciences of the United States of America, 115(43):10943–10947, 2018.
  * [19] Matt Harrington and Douglas J. Durian.  Anisotropic particles strengthen granular pillars under compression.  Phys. Rev. E, 97:012904, Jan 2018.
  * [20] Xiaoguang Ma, Zoey S. Davidson, Tim Still, Robert J. S. Ivancic, S. S. Schoenholz, A. J. Liu, and A. G. Yodh.  Heterogeneous activation, local structure, and softness in supercooled colloidal liquids.  Phys. Rev. Lett., 122:028001, Jan 2019.
  * [21] François P. Landes, Giulio Biroli, Olivier Dauchot, Andrea J. Liu, and David R. Reichman.  Attractive versus truncated repulsive supercooled liquids: The dynamics is encoded in the pair correlation function.  Physical Review E, 101(1):1–5, 2020.
  * [22] Ekin D. Cubuk, Andrea J. Liu, Efthimios Kaxiras, and Samuel S. Schoenholz.  Unifying framework for strong and fragile liquids via machine learning: a study of liquid silica, 2020.
  * [23] Indrajit Tah, Tristan A. Sharp, Andrea J. Liu, and Daniel M. Sussman.  Quantifying the link between local structure and cellular rearrangements using information in models of biological tissues, 2020.
  * [24] Ge Zhang, Sean Ridout, and Andrea J. Liu.  Interplay of rearrangements, strain, and local structure during avalanche propagation, 2020.
  * [25] Jason W. Rocks, Sean A. Ridout, and Andrea J Liu.  Learning-based approach to plasticity in athermal sheared amorphous packings: Improving softness, 2020.
  * [26] Patrick Charbonneau, Jorge Kurchan, Giorgio Parisi, Pierfrancesco Urbani, and Francesco Zamponi.  Glass and Jamming Transitions: From Exact Results to Finite-Dimensional Descriptions.  Annual Review of Condensed Matter Physics, 8(1):265–288, 2017.
  * [27] Ludovic Berthier, Giulio Biroli, Patrick Charbonneau, Eric I. Corwin, Silvio Franz, and Francesco Zamponi.  Gardner physics in amorphous solids and beyond.  Journal of Chemical Physics, 151(1), 2019.
  * [28] Patrick Charbonneau, Eric I. Corwin, Giorgio Parisi, and Francesco Zamponi.  Universal microstructure and mechanical stability of jammed packings.  Physical Review Letters, 109(20):1–5, 2012.
  * [29] Patrick Charbonneau, Eric I. Corwin, Giorgio Parisi, and Francesco Zamponi.  Jamming criticality revealed by removing localized buckling excitations.  Physical Review Letters, 114(12):1–5, 2015.
  * [30] Patrick Charbonneau, Eric I. Corwin, Giorgio Parisi, Alexis Poncet, and Francesco Zamponi.  Universal non-debye scaling in the density of states of amorphous solids.  Phys. Rev. Lett., 117:045503, Jul 2016.
  * [31] Masanari Shimada, Hideyuki Mizuno, Ludovic Berthier, and Atsushi Ikeda.  Low-frequency vibrations of jammed packings in large spatial dimensions, 2019.
  * [32] Merlijn S. van Deen, Johannes Simon, Zorana Zeravcic, Simon Dagois-Bohy, Brian P. Tighe, and Martin van Hecke.  Contact changes near jamming.  Phys. Rev. E, 90:020202, Aug 2014.
  * [33] Peter Morse, Sven Wijtmans, Merlijn van Deen, Martin van Hecke, and M. Lisa Manning.  Differences in plasticity between hard and soft spheres.  Phys. Rev. Research, 2:023179, May 2020.
  * [34] Shubha Tewari, Dylan Schiemann, Douglas J. Durian, Charles M. Knobler, Stephen A. Langer, and Andrea J. Liu.  Statistics of shear-induced rearrangements in a two-dimensional model foam.  Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, 60(4):4385–4396, 1999.
  * [35] Craig E. Maloney and Anaël Lemaître.  Amorphous systems in athermal, quasistatic shear.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 74(1):1–22, 2006.
  * [36] M L Falk and J S Langer.  Dynamics of viscoplastic deformation in amorphous solids M.  Physical Review E, 57(6):14, 1998.
  * [37] See Supplemental Material at [URL will be inserted by publisher] for data showing lack of finite size effects, regression $R^{2}$ as a function of pressure, and correlation lengths and $P{\left(R|Z\right)}$ as a function of pressure..
  * [38] Sylvain Patinet, Damien Vandembroucq, and Michael L. Falk.  Connecting Local Yield Stresses with Plastic Activity in Amorphous Solids.  Physical Review Letters, 117(4):045501, jul 2016.
  * [39] Andrea J. Liu and Sidney R. Nagel.  The jamming transition and the marginally jammed solid.  Annual Review of Condensed Matter Physics, 1:347–369, 2010.
  * [40] Edan Lerner, Eric Degiuli, Gustavo Düring, and Matthieu Wyart.  Breakdown of continuum elasticity in amorphous solids.  Soft Matter, 10(28):5085–5092, 2014.
  * [41] Daniel Hexner, Andrea J. Liu, and Sidney R. Nagel.  Two Diverging Length Scales in the Structure of Jammed Packings.  Physical Review Letters, 121(11):115501, 2018.
  * [42] Carl P. Goodrich, Simon Dagois-Bohy, Brian P. Tighe, Martin van Hecke, Andrea J. Liu, and Sidney R. Nagel.  Jamming in finite systems: Stability, anisotropy, fluctuations, and scaling.  Phys. Rev. E, 90:022138, Aug 2014.
  * [43] Giorgio Parisi, Pierfrancesco Urbani, and Francesco Zamponi.  Theory of Simple Glasses: Exact Solutions in Infinite Dimensions.  Cambridge University Press, New York, February 2020.
  * [44] E. D. Cubuk, R. J. S. Ivancic, S. S. Schoenholz, D. J. Strickland, A. Basu, Z. S. Davidson, J. Fontaine, J. L. Hor, Y.-R. Huang, Y. Jiang, N. C. Keim, K. D. Koshigan, J. A. Lefever, T. Liu, X.-G. Ma, D. J. Magagnosc, E. Morrow, C. P. Ortiz, J. M. Rieser, A. Shavit, T. Still, Y. Xu, Y. Zhang, K. N. Nordstrom, P. E. Arratia, R. W. Carpick, D. J. Durian, Z. Fakhraai, D. J. Jerolmack, Daeyeon Lee, Ju Li, R. Riggleman, K. T. Turner, A. G. Yodh, D. S. Gianola, and Andrea J. Liu.  Structure-property relationships from universal signatures of plasticity in disordered solids.  Science, 358(6366):1033–1037, 2017.

Supporing Information:
Correlation of plastic events with local structure in jammed packings across
spatial dimensions

Figure S1: $C$ for various $N$ in $2d$. This data shows no evidence of the
finite-size scaling seen in other properties associated with a diverging
correlation length at $P=0$, i.e., the low pressure at which accuracy
decreases does not become higher at smaller $N$.  Figure S2: $R^{2}$ for
regression against $\Delta^{2}_{\mathrm{min}}$ as defined in main text. In
addition to the decay with dimension shown in the main text, here the very
weak pressure dependence is made clear. Figure S3: Correlation lengths
extracted from fits to initial decay of small-small-particle correlation
functions in $2d$, showing weak dependence on pressure. The $D^{2}_{\min}$
correlation length depends more on pressure, but does not seem to diverge at
unjamming. Blue circles: $D^{2}_{\min}$, orange squares:
$\Delta^{2}_{\mathrm{min}}$, green triangles: $S$. Figure S4: Probability of
rearranging for small particles as a function of $(Z-\langle
Z\rangle)/\sigma_{Z}$ in several pressures in $3d$. ($P=1e-3$ black, $P=1e-2$
gray, $P=5e-2$ white). Leaving out pressures above $P^{*}$, the dependence of
the rearrangement probability on this rescaled $Z$ is qualitatively similar at
different $P$. Inset shows $\langle Z\rangle$ (black x’s) and $\sigma_{Z}$
(blue circles) for small particles. Since $\langle Z\rangle$ and $\sigma_{Z}$
are anticorrelated, rather than correlated as when $d$ is varied in the main
text, this provides further evidence that the rescaling used is the correct
one.  Figure S5:  Inverse participation ratio (IPR)
$\sum_{i=1}^{N}u_{i}^{4}/\sum_{i=1}^{N}u_{i}^{2}$ of critical modes (with
displacement $\mathbf{u}$). Black triangles: $2d$, $N=4096$. Blue squares:
$3d$, $N=8192$. Green pentagons: $4d$, $N=8192$. Red hexagons: $5d$,
$N=16384$. Because it is computationally intensive to study larger systems in
high $d$, we did not look at the dependence of the participation ratio on
system size, which is needed in order to determine whether the modes are
spatially localized or not


# State estimation of a carbon capture process through POD model reduction and
neural network approximation

Siyu Liu1,2, Xunyuan Yin3, Jinfeng Liu2 1Siyu Liu is with School of Internet
of Things Engineering, Jiangnan University, Wuxi 214122, China.
<EMAIL_ADDRESS>Liu and Siyu Liu are with Department of Chemical &
Materials Engineering, University of Alberta, Edmonton, AB, Canada, T6G 1H9
<EMAIL_ADDRESS>Yin is with School of Chemistry, Chemical
Engineering and Biotechnology, Nanyang Technological University, 62 Nanyang
Drive, Singapore, 637459<EMAIL_ADDRESS>
###### Abstract

This paper presents an efficient approach for state estimation of post-
combustion CO2 capture plants (PCCPs) by using reduced-order neural network
models. The method involves extracting lower-dimensional feature vectors from
high-dimensional operational data of the PCCP and constructing a reduced-order
process model using proper orthogonal decomposition (POD). Multi-layer
perceptron (MLP) neural networks capture the dominant dynamics of the process
and train the network parameters with low-dimensional data obtained from open-
loop simulations. The proposed POD-MLP model can be used as the basis for
estimating the states of PCCPs at a significantly decreased computational
cost. For state estimation, a reduced-order extended Kalman filtering (EKF)
scheme based on the POD-MLP model is developed. Our simulations demonstrate
that the proposed POD-MLP modeling approach reduces computational complexity
compared to the POD-only model for nonlinear systems. Additionally, the POD-
MLP-EKF algorithm can accurately reconstruct the full state information of
PCCPs while significantly improving computational efficiency compared to the
EKF based on the original PCCP model.

##  I INTRODUCTION

In recent years, post-combustion CO2 capture plants (PCCPs) have gained
significant attention due to their potential in reducing greenhouse gas
emissions and mitigate global warming. PCCPs are commonly used in power plants
and carbon-intensive industrial processes to separate CO2 from the flue gas
[1]. The operational safety, carbon capture efficiency and economic cost of
PCCPs are highly dependent on the performance of the advanced control systems
used for regulating the process operations [2].

Real-time information of the key quality variables of the PCCP is essential
for the advanced control system to make the most appropriate decisions for
safe and efficient process operation. However, measuring all quality variables
online through deploying hardware sensors is unrealistic. Therefore, it is
crucial to exploit state estimation to reconstruct the full state information
for PCCPs. Unfortunately, results on state estimation of PCCPs have been
limited. In [3], we made an initial attempt on estimating the states of the
absorber of a PCCP by developing a distributed moving horizon estimation
scheme. However, the other key quality variables associated with the
desorption unit and other physical units are not addressed within this
framework. Wang et al., presented a robust soft sensor using a neural network
and moving horizon estimator to monitor key operating parameters in the carbon
capture process [4]. In the context of nonlinear state estimation, there have
been some algorithms that have the potential to be leveraged for state
estimation of PCCPs, e.g., extended Kalman filtering [5], moving horizon
estimation methods [6], and particle filter [7]. However, most of nonlinear
estimation algorithms require accurate first-principles dynamic models of the
underlying nonlinear processes. Due to the large scales and complex structures
of PCCPs, it can be challenging to conduct first-principles modeling.
Additionally, even if a first-principles model is obtained, this type of model
that will involve partial differential equations that describe the dynamical
behaviors of the absorption and desorption columns will be computationally
expensive to simulate, especially considering the cases when optimization-
based estimation approaches are used.

Data-driven modeling using neural networks has been widely used as an
alternative to first-principles modeling for various nonlinear processes. For
example, Jeon et al. utilized the neural networks (NNs) to describe the
dynamics of a chemical reactor, which was further used to optimize its control
performance [8]. Cheng et al. employed NNs to model the ship motion and
improved its navigation accuracy [9]. Zhao et al., developed reduced-order
recurrent neural networks that capture the dominant dynamics of nonlinear
systems using an autoencoder [10]. However, these methods may face certain
limitations when considering possible applications to PCCPs. Specifically, the
complexity of the NN model, including the number of neurons and layers, and
the computational cost of training, may increase exponentially with the
dimensionality of the input and output data, commonly referred to as the
“curse of dimensionality”.

To address these limitations associated with state estimation of PCCPs, we
propose a solution that combines data-driven modeling using NNs with model
reduction techniques. Proper orthogonal decomposition (POD) has been widely
adopted in various engineering fields for reducing the dimensionality of high-
dimensional data sets while preserving dominant patterns or features of the
data. In the context of control systems, POD can be used to satisfactorily
approximate the dynamics of a large-scale nonlinear process in a lower-
dimensional state space [11, 12]. POD holds the promise to lower the
dimensionality of the originally high-dimensional state-space model for PCCP.
Specifically, our objective is to leverage POD to create a reduced-order
machine learning-based model with lower structural complexity. This model can
then facilitate reductions in the computational cost for model training and
implementation of state estimation based on the reduced-order model.

Motivated by the observations above, we propose a neural network-based state
estimation approach for PCCPs using POD reduced-order models. Specifically, we
normalize the data before POD and use the POD approach to obtain a reduced-
order model that accurately approximates the dynamics of the PCCP. Then, we
train a multi-layer perceptron (MLP) neural network to capture the dominant
dynamics of the POD reduced-order model from the low-dimensional data. The
resulting reduced-order POD-based MLP (POD-MLP) model is used as the basis of
state estimation. We develop a reduced-order extended Kalman filtering (EKF)
algorithm based on the POD-MLP model to estimate the states of the original
process. Our approach can effectively reduce the complexity of the NN model
and the computational workload required for training while preserving high
modeling accuracy. This approach provides a promising solution to state
estimation of PCCPs and can potentially be applied to other complex systems
with of high dimension for efficient state estimation through bypassing the
use of the original higher-dimensional dynamic process model. Our proposed
framework offers an advantage in that it can be employed when only data is
available. Specifically, in the absence of a physical model, we can still
leverage proposed method to build a reduced-order neural network and perform
efficient state estimation.

##  II MODEL DESCRIPTION

Figure 1 illustrates the post-combustion capture plant considered in this
paper. This diagram shows the four key components of the PCCP, which include
the absorption and desorption columns, lean-rich heat exchanger (LRHE), and
the reboiler. The flue gas, which contains a high concentration of CO2, is
introduced at the bottom of the absorber from the power plant and is then
mixed with a lean solvent having low CO2 levels. The 5M Monoethanolamine (MEA)
is used as the solvent in this study. The treated flue gas with a reduced
amount of CO2 leaves the absorption column, while the rich solvent with a high
concentration of CO2 is heated via the heat exchanger by exchanging heat with
the lean solvent coming from the reboiler. The rich solvent is then fed into
the top of the desorption column, where it is heated through contact with the
hot vapor from the reboiler. In the desorption column, the CO2 is stripped
from the rich solvent, which is then recycled back to the absorber. The
discharged CO2 gas, with a high concentration of CO2 (90-99$\%$), is obtained
from the desorption column.

Figure 1: A schematic of the CO2 capture plant.

The details of the PCCP model are briefly described in the following [13, 14]:

$\displaystyle\frac{\partial C_{L}(i)}{\partial
t}=\frac{F_{L}}{S_{c}}\frac{\partial C_{L}(i)}{\partial l}+(N(i)a^{I}),$ (1a)
$\displaystyle\frac{\partial C_{G}(i)}{\partial
t}=-\frac{F_{G}}{S_{c}}\frac{\partial C_{G}(i)}{\partial l}-(N(i)a^{I}),$ (1b)
$\displaystyle\frac{\partial T_{L}}{\partial
t}=\frac{F_{L}}{S_{c}}\frac{\partial T_{L}}{\partial
l}+\frac{(Q_{L}a^{I})}{\sum_{i=1}^{n}C_{L}(i)C_{p,i}},$ (1c)
$\displaystyle\frac{\partial T_{G}}{\partial
t}=-\frac{F_{G}}{S_{c}}\frac{\partial T_{G}}{\partial
l}+\frac{(Q_{G}a^{I})}{\sum_{i=1}^{n}C_{G}(i)C_{p,i}},$ (1d)
$\displaystyle\frac{dT_{tu}}{dt}=\frac{\dot{V}_{tu}}{V_{tu}}(T_{tu,in}-T_{tu,out})+\frac{\dot{Q}}{C_{p_{tu}}\rho_{tu}V_{tu}},$
(1e)
$\displaystyle\frac{dT_{sh}}{dt}=\frac{\dot{V}_{sh}}{V_{sh}}(T_{sh,in}-T_{sh,out})+\frac{\dot{Q}}{C_{p_{sh}}\rho_{sh}V_{sh}},$
(1f) $\displaystyle\rho
C_{p}V\frac{dT_{reb}}{dt}=f_{in}H_{in}-f_{V}H_{V,out}-f_{L}H_{L,out}+Q_{reb}.$
(1g)

The dynamic models for the absorption column and desorption column are
described by the partial differential equations (1a)–(1d), where
$i=CO_{2},MEA,H_{2}O,N_{2}$, and the subscripts $L$ and $G$ denote liquid and
gas phases, respectively. The dependent variables vary with time $t$ and axial
position $l$ of the column. It is assumed that each stage in the two columns
is well mixed. Their dynamic models are similar except for a few details like
the direction of reactions, temperature, and reaction rate constants. The
energy balance equations in (1e)–(1f) represent the dynamics of the lean-rich
heat exchanger, where $\dot{V}(\rm m^{3}/s)$ and $\dot{Q}(\rm kJ/s)$ represent
the volumetric flow and heat transfer rate, respectively, the subscripts $tu$,
$sh$, $in$ and $out$ denote the tube-side, shell-side, inlets and outlets of
the heat exchanger, respectively. It is assumed that the mass inside the heat
exchanger remains constant. Equation (1g) is the energy balance equation of
the reboiler, where $T_{reb}(\rm K)$ represents the temperature of the
reboiler, the subscripts $in$, $out$, $V$, and $L$ denote inlet, outlet,
vapour and liquid, respectively, and $Q_{reb}(\rm KJ/s)$ is the heat input.
The definitions of other variables and parameters of the PCCP can be found in
Table I. Physical property calculations of gas and liquid phases are necessary
for the model development, and they are estimated from seven nonlinear
algebraic correlations. Details of these calculations are not included in this
work but can be found in [13, 14].

TABLE I: Process variables of each unit of the PCCP. Notation | Definition | Unit
---|---|---
$C_{i}$ | Molar concentrations of component $i$ | $\rm mol/m^{3}$
$S_{c}$ | Cross-sectional area of the column | $\rm m^{2}$
$F$ | Volumetric flow | $\rm m^{3}/s$
$N_{i}$ | Mass transfer rate | $\rm kmol/m^{2}s$
$T$ | Temperature | $\rm K$
$l$ | Height of the column | $\rm m$
$C_{p}$ | Heat capacity | $\rm KJ/kmol$
$Q$ | Heat transfer rate | $\rm KJ/m^{2}s$
$a^{I}$ | Interfacial area | $\rm m^{2}/m^{3}$
$V$ | Volume | $\rm m^{3}$
$\rho$ | Average molar density | $\rm kmol/m^{3}$
$H$ | Enthalpy | $\rm KJ$
$f$ | Flow rate | $\rm mol/s$

The PCCP model consists of partial differential equations for the two columns,
and ordinary differential equations for the heat exchanger and reboiler, as
well as some algebraic equations for parameter calculations. As the variables
in the columns exhibit temporal and spatial distributions, the partial
differential equations are discretized using the method outlined in [13] to
convert them into ordinary differential equations, with the column length
divided into five stages. Therefore, the model presented in (1) is expressed
as a system of differential-algebraic equations (DAEs):

$\displaystyle\bm{x}(k+1)=\bm{F}(\bm{x}(k),\bm{a}(k),\bm{u}(k))+\bm{w}(k),$
(2a) $\displaystyle\bm{G}(\bm{x}(k),\bm{a}(k),\bm{u}(k))=\bm{0},$ (2b)
$\displaystyle\bm{y}(k)=\bm{H}(\bm{x}(k),\bm{u}(k))+\bm{v}(k),$ (2c)

where $\bm{x}(k)\in{\mathbb{R}}^{103}$ is the state vector, and the
definitions of state variables are defined in Table II,
$\bm{a}(k)\in{\mathbb{R}}^{7}$ is the algebraic state vector,
$\bm{u}(k)=[F_{L},Q_{reb},F_{G}]\in{\mathbb{R}}^{3}$ denotes the input vector:
solvent flow rate in L/s, reboiler heat in KJ/s, and flue gas flow rate m3/s,
and $\bm{w}(k)$ is the process noise.

TABLE II: Definition of the PCCP state variables at the $j$th discrete point ($j=1,2,\dots,5$). States | Def. | States | Def. | States | Def.
---|---|---|---|---|---
$x_{1-5}$ |  $C_{L}^{j}$(N2) | $x_{6-10}$ |  $C_{L}^{j}$(CO2) | $x_{11-15}$ |  $C_{L}^{j}$(MEA)
$x_{16-20}$ |  $C_{L}^{j}$(H2O) | $x_{21-25}$ | $T_{L}^{j}$ | $x_{26-30}$ |  $C_{G}^{j}$(N2)
$x_{31-35}$ |  $C_{G}^{j}$(CO2) | $x_{36-40}$ |  $C_{G}^{j}$(MEA) | $x_{41-45}$ |  $C_{G}^{j}$(H2O)
$x_{46-50}$ | $T_{G}^{j}$ | $x_{51-55}$ |  $C_{L}^{j}$(N2) | $x_{56-60}$ |  $C_{L}^{j}$(CO2)
$x_{61-65}$ |  $C_{L}^{j}$(MEA) | $x_{66-70}$ |  $C_{L}^{j}$(H2O) | $x_{71-75}$ | $T_{L}^{j}$
$x_{76-80}$ |  $C_{G}^{j}$(N2) | $x_{81-85}$ |  $C_{G}^{j}$(CO2) | $x_{86-90}$ |  $C_{G}^{j}$(MEA)
$x_{91-95}$ |  $C_{G}^{j}$(H2O) | $x_{96-100}$ | $T_{G}^{j}$ | $x_{101}$ | $T_{h1}$
$x_{102}$ | $T_{h2}$ | $x_{103}$ | $T_{reb}$ |  |

##  III POD AND ITS APPLICATION TO PCCP

In this section, we employ the proper orthogonal decomposition (POD) to derive
a reduced-order model that approximates the dynamics of the PCCP that
originally has 103 state variables. By projecting the high-dimensional state
variables onto a lower-dimensional subspace, that captures the dominant modes
of variability, we obtain a reduced-order model that accurately captures the
essential dynamics of the PCCP while significantly reducing its complexity.
This is accomplished by computing the singular value decomposition (SVD) of
the data matrix, which yields a set of orthonormal basis vectors that
represent the most significant patterns of variability in the data.

For general nonlinear systems described by (2), we obtain a state trajectory
by capturing and sampling the system’s response to a typical input trajectory
at fixed time intervals $\delta$. Then, we sample the resulting state
trajectory to construct a matrix of process states from time $0$ to $N$,
denoted as:

$\displaystyle\bm{\chi}=[\bm{x}(0)\ \bm{x}(1)\ \dots\
\bm{x}(N)]\in{\mathbb{R}}^{n\times(N+1)},$ (3)

where the snapshot matrix $\bm{\chi}$ is composed of the actual state at each
sampling interval, the number of state variables is denoted as $n$, and the
number of sampling intervals is represented by $N$. To ensure a sufficient
number of samples, we require $N$ to be much larger than $n$.

For PCCP, the magnitudes of different states vary greatly. To ensure that the
POD reduction method is not biased towards states with larger magnitudes in
PCCP, each state variable $x_{i}$ (where $i=1,2,\dots,103$) in the data matrix
$\bm{\chi}$ is normalized using (4) prior to performing SVD decomposition:

$\displaystyle x_{i,norm}=\frac{x_{i}-x_{i,min}}{x_{i,max}-x_{i,min}}.$ (4)

This normalization transforms each state variable $x_{i}$ to $x_{i,norm}$,
where $x_{i,min}$ and $x_{i,max}$ are the minimum and maximum values of
$x_{i}$ in the original dataset, respectively. The resulting normalized matrix
$\bm{\chi}_{norm}$ constructed by $x_{i,norm}$ has all states with magnitudes
between 0 and 1. SVD is then performed on the normalized matrix
$\bm{\chi}_{norm}$ as shown in the following,

$\displaystyle\bm{\chi}_{norm}=\bm{U}\bm{\varSigma}\bm{V}^{\tiny\mbox{\rm
T}},$ (5)

where $\bm{U}\in{\mathbb{R}}^{n\times n}$ and
$\bm{V}\in{\mathbb{R}}^{(N+1)\times(N+1)}$ are orthogonal matrices, the
rectangular matrix $\bm{\varSigma}\in{\mathbb{R}}^{n\times(N+1)}$ has non-
negative real values on its main diagonal. The diagonal entries $\sigma_{i}$
represent the singular values of the $\bm{\chi}$ matrix, where
$i\in{1,2,\dots,n}$. These values are sorted in descending order on the main
diagonal of $\bm{\varSigma}$.

To construct a reduced-order model, we select a positive integer $r$ that is
smaller than the number of states $103$, and truncate $\bm{\varSigma}$ at the
$r$th column and row to from the reduced-order matrix
$\bm{\varSigma}_{r}\in{\mathbb{R}}^{r\times r}$ using the first $r$ singular
values $\sigma_{i}$. Accordingly, we select the first $r$ columns of $\bm{U}$
and the first $r$ rows of $\bm{V}^{\tiny\mbox{\rm T}}$ to form the matrices
$\bm{U}_{r}$ and $\bm{V}_{r}^{\tiny\mbox{\rm T}}$, respectively. Using these
matrices, we obtain a reduced-order approximation of normalized process data,
given by

$\displaystyle\bm{\chi}_{norm}\approx\bm{U}_{r}\bm{\varSigma}_{r}\bm{V}_{r}^{\tiny\mbox{\rm
T}}.$ (6)

We define $\bm{\xi}\in{\mathbb{R}}^{r}$ as the state vector of the reduced-
order model, and set $\bm{\xi}(k):=\bm{U}_{r}^{\tiny\mbox{\rm T}}x(k)$. Using
the truncated SVD matrices, the original nonlinear model in (2) can be
expressed as a reduced-order model in state-space form:

$\displaystyle\bm{\xi}(k+1)=\bm{U}_{r}^{\tiny\mbox{\rm
T}}\bm{F}(\bm{U}_{r}\bm{\xi}(k),\bm{a}(k),\bm{u}(k))+\bm{U}_{r}^{\tiny\mbox{\rm
T}}\bm{w}(k),$ (7a)
$\displaystyle\bm{G}(\bm{U}_{r}\bm{\xi}(k),\bm{a}(k),\bm{u}(k))=\bf 0,$ (7b)
$\displaystyle\bm{y}(k)=\bm{H}(\bm{U}_{r}\bm{\xi}(k),\bm{u}(k))+\bm{v}(k).$
(7c)

The evolution of $\bm{\xi}(k)$ in the reduced-order model can be used to
approximate the actual state trajectory of the original nonlinear process
through the mapping $\bm{x}(k)\approx\bm{U}_{r}\bm{\xi}(k)$.

###### Remark 1

In this section, we improve the accuracy of the model approximation by
normalizing the data prior to applying the POD method to obtain a reduced-
order model. The effectiveness of the normalization for POD will also be
illustrated through simulations in Section VI.

##  IV APPROXIMATING REDUCED-ORDER MODEL WITH MLP NETWORKS

The POD technique is commonly applied in linear systems to decrease
computational costs by reducing the dimensionality of the problem. However, it
does not yield similar advantages for nonlinear systems due to the challenge
in explicitly expressing $\bm{U}_{r}^{\tiny\mbox{\rm
T}}\bm{F}(\bm{U}_{r}\bm{\xi},\bm{a},\bm{u})$ in terms of the reduced basis
$\bm{U}_{r}$. Consequently, evaluating the reduced-order model may require
more time than evaluating the original nonlinear function $\bm{F}$. Different
approaches have been adopted to address this issue. For instance, a linear
parameter varying model was utilized to approximate the reduced-order model,
however, the benefits were found to be insignificant [15].

To address this issue, we present a method to speed up the evolution of
reduced-order models for nonlinear systems such as (2). The method employs a
multi-layer perceptron (MLP) neural network to fit the reduced-order model and
hence decrease the computation time. The MLP neural network model of
$\bm{\xi}$ is given by:

$\displaystyle\bm{\xi}(k+1)=\bm{f}_{mlp}(\bm{\xi}(k),\bm{u}(k))+\bm{w}_{r}(k),$
(8)

where the vector function $\bm{f}_{mlp}\in{\mathbb{R}}^{r}$ approximates the
dynamic behavior of $\bm{\xi}$ in (7a), and $\bm{w}_{r}\in{\mathbb{R}}^{r}$
denotes the process noise and model error. Consequently, we do not have to
evaluate the vector function
$\bm{F}:{\mathbb{R}}^{103}\rightarrow{\mathbb{R}}^{103}$ of the full-order
model, leading to significant time savings.

MLP models are a type of artificial neural network widely used to approximate
complex and nonlinear problems. An MLP typically consist of an input layer,
one or more hidden layers, and an output layer. Each layer in an MLP contains
multiple fully connected neurons, which are connected to the next layer by
weights. In supervised learning, the weights are adjusted to approximate each
target value. The number of neurons in the input and output layers is
determined by the input and output variables, respectively. The computation
time required for MLP output is relatively short, as only a few matrix
multiplications, vector additions, and function evaluations are necessary.
Given these characteristics, MLP is a suitable choice for approximating the
vector function $\bm{f}_{mlp}$ in equation (2) after POD model reduction. The
basic model formulation of MLP is indicated below:

$\displaystyle\bm{z}^{(l)}=\sigma_{h}(\bm{w}^{(l)}\bm{z}^{(l-1)}+\bm{b}^{(l)}),\quad\bm{y}=\bm{w}^{o}\bm{z}+\bm{b}^{o}$
(9)

where $\bm{z}^{(l)}$ denotes the output vector of the $l$-th hidden layer,
obtained by applying an activation function $\sigma_{h}$ to the weighted sum
of the input vector $\bm{z}^{(l-1)}$ and the bias vector $\bm{b}^{(l)}$. The
weight matrix of the $l$-th hidden layer is denoted by $\bm{w}^{(l)}$. The
input vector $\bm{z}^{(0)}$ is the MLP input, and the output vector $\bm{y}$
is obtained by applying the weight matrix $\bm{w}^{o}$ to the output vector of
the last hidden layer and adding the bias vector $\bm{b}^{o}$. It is worth
noting that the MLP output layer is typically linear, while the activation
function of the hidden layer can be chosen from various options based on the
specific problem being solved.

The reduced-order PCCP model is enhanced by utilizing an MLP network, where
the input is $\bm{\xi}_{u}:=[\bm{u}^{\tiny\mbox{\rm
T}},\bm{\xi}^{\tiny\mbox{\rm T}}]^{\tiny\mbox{\rm T}}\in{\mathbb{R}}^{3+r}$,
and the output is $\hat{\bm{\xi}}\in{\mathbb{R}}^{r}$. The input layer has
$3+r$ neurons, and the output layer has $r$ neurons. The MLP model is trained
to minimize the mean-squared-error (MSE) between the predicted output
$\hat{\bm{\xi}}$ and the actual output $\bm{\xi}$:

$\displaystyle L=MSE(\bm{\xi},\hat{\bm{\xi}}).$ (10)

The model structure of the proposed POD-based MLP (POD-MLP) model is shown in
Figure 2.

Figure 2: The diagram of the POD-MLP model.

###### Remark 2

As the system dimension increases, the number of training samples required to
construct an MLP model that accurately approximates the system grows
exponentially. Therefore, direct training a large-scale nonlinear model can be
computationally expensive. Hence, performing POD dimension reduction for high-
dimensional PCCP models and then training the MLP model for the reduced-order
model is a highly practical approach.

###### Remark 3

The POD-MLP model not only reduces the order of the system, but also
eliminates the need for solving the DAEs in the original PCCP model. It
achieves this by extracting the dynamic information of the algebraic state
variable $\bm{a}$ from the data, and excluding $\bm{a}$ from the model.

##  V EKF USING the POD-MLP MODEL

In this section, we develop an extended Kalman filtering (EKF) based on the
POD-MLP model (the POD-MLP-EKF algorithm for short) to estimate the actual
process states. The POD-MLP model is summarized in the following:

$\displaystyle\bm{\xi}(k+1)=\bm{f}_{mlp}(\bm{\xi}(k),\bm{u}(k))+\bm{w}_{r}(k),$
(11a)
$\displaystyle\bm{y}(k)=\bm{H}(\bm{U}_{r}\bm{\xi}(k),\bm{u}(k))+\bm{v}(k).$
(11b)

Assuming that $\bm{w}_{r}(k)$ and $\bm{v}(k)$ are two mutually uncorrelated
Gaussian noise sequences with zero-mean, we further assume that they have
covariance matrices $\bm{Q}_{r}$ and $\bm{R}_{r}$.

Based on the above preparation, the POD-MLP-EKF algorithm is designed in the
following two steps:
Step 1: Prediction step:

$\displaystyle\hat{\bm{\xi}}(k+1|k)=\bm{f}_{mlp}(\hat{\bm{\xi}}(k|k),\bm{u}(k)),$
(12a) $\displaystyle\bm{P}(k+1|k)=\bm{A}(k)\bm{P}(k|k)\bm{A}^{\tiny\mbox{\rm
T}}(k)+\bm{Q}_{r},$ (12b)

where $\hat{\bm{\xi}}(k+1|k)$ is the prediction of the system state at time
$k+1$ based on the current state estimate $\hat{\bm{\xi}}(k|k)$ and the input
$\bm{u}(k)$, and $\bm{P}(k+1|k)$ contains a priori error covariance
information, incorporating the prediction error and the uncertainty associated
with the system dynamics through the covariance matrix $\bm{Q}_{r}$. The
matrix
$\bm{A}(k)=\frac{\partial\bm{f}_{mlp}(\bm{\xi},\bm{u})}{\partial\bm{\xi}}|_{\bm{\xi}=\hat{\bm{\xi}}(k|k)}$
is the Jacobian matrix of the MLP model with respect to the state vector
$\bm{\xi}$ evaluated at the predicted state estimate $\hat{\bm{\xi}}(k|k)$.
Step 2: Update step:

$\displaystyle\bm{K}(k+1)=\frac{\bm{P}(k+1|k)\bm{C}^{\tiny\mbox{\rm
T}}(k)}{\bm{C}(k)\bm{P}(k+1|k)\bm{C}^{\tiny\mbox{\rm T}}(k)+\bm{R}_{r}},$
(13a)
$\displaystyle\hat{\bm{\xi}}(k+1|k+1)=\hat{\bm{\xi}}(k+1|k)+\bm{K}(k+1)[y(k+1)$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\
-\bm{H}(\bm{U}_{r}\hat{\bm{\xi}}(k+1|k),\bm{u}(k))],$ (13b)
$\displaystyle\bm{P}(k+1|k+1)=[\bm{I}-\bm{K}(k+1)\bm{C}(k)]\bm{P}(k+1|k),$
(13c)

where the correction gain $\bm{K}(k+1)$ is computed based on the a priori
estimation error covariance $\bm{P}(k+1|k)$, the measurement error covariance
$\bm{R}_{r}$, and the observation matrix $\bm{C}(k)$, which maps the predicted
state $\hat{\bm{\xi}}(k+1|k)$ to the measurement space. The state estimate is
then updated to $\hat{\bm{\xi}}(k+1|k+1)$ using the correction gain and the
measurement innovation
$y(k+1)-\bm{H}(\bm{U}_{r}\hat{\bm{\xi}}(k+1|k),\bm{u}(k))$, which represents
the difference between the actual measurement and the predicted measurement
based on the predicted state. Finally, the a posteriori estimation error
covariance matrix $\bm{P}(k+1|k+1)$ is computed based on the updated state
estimate and the correction gain, which reflects the reduced uncertainty in
the estimated state after incorporating the measurement information.

Then, we can obtain the state estimates of the actual PCCP states, denoted by
$\hat{\bm{x}}$, by utilizing the reduced-order state estimate
$\hat{\bm{\xi}}(k+1|k+1)$ and the linear mapping $\bm{U}_{r}$, as follows:

$\displaystyle\hat{\bm{x}}(k+1)=\bm{U}_{r}\hat{\bm{\xi}}(k+1|k+1).$ (14)

##  VI SIMULATION RESULTS

To perform the model order reduction, the input vector
$\bm{u}(k)=[F_{L},Q_{reb},F_{G}]$ is used to excite the PCCP, which are
constrained as 0.48 L/s $\leq F_{L}\leq$ 0.66 L/s, 0.14 KJ/s $\leq
Q_{reb}\leq$ 0.20 KJ/s and 0.8 m3/s $\leq F_{G}\leq$ 1.2 m3/s. The dynamic
model of the PCCP is discretized at a sample interval of $\Delta=30$s. Pseudo-
random multi-level signals (PRMSs) are commonly used as excitation signals for
identifying nonlinear systems. For example, a ten-level PRMS used in the PCCP
is shown in Figure 3 (Samples 1 to 3000). The switching times are randomly
chosen from a uniform distribution ranging between 900s and 3000s (30 to 100
sampling times). The system states are sampled over a duration of 100 hours to
construct the snapshot matrix $\bm{\chi}$ for POD reduction. With 12,000
sampling intervals, the requirement $N\gg n$ is satisfied.

Figure 3: The three input signals for generating system state trajectories in
model reduction.

Next, we apply the SVD to the matrix $\bm{\chi}$ which has dimensions 103 by
12000. This results in a unitary matrix $U$ that can be used for coordinate
transformation. To evaluate the reduced-order model accuracy, we use the root-
mean-square error (RMSE) to denote the model errors, which is defined as
follows:

$\displaystyle{\rm RMSE}$
$\displaystyle:=\sqrt{\frac{\sum^{N}_{j=0}\sum^{103}_{i=1}(x_{i,norm}(j)-\hat{x}_{i,norm}(j))^{2}}{N}},$

where $\hat{x}_{i,norm}:=\bm{U}_{r}\xi_{i}$ is the $i$th approximated state
obtained from a reduced order model. To validate the accuracy of a reduced-
order model, we should use input trajectories that are different from the ones
used in POD model reduction. Additional 600 sampling data are used for
validation. Based on the actual state trajectories and reduced-order model
state trajectories, the RMSE is calculated for each model. The values of
log(RMSE) at different orders $r=20,\dots,90$ under the POD with normalization
and the POD without normalization are shown in Figure 4. It shows that the
degree of model mismatch increases with the decrease in the model order for
both cases. Moreover, the RMSE value of POD with normalization is smaller than
that without normalization. This is because for the POD method without
normalization of the data matrix $\bm{\chi}$, the approximation of the state
with small numerical values is not accurate when the model is reduced, and the
model approximation error is added to each state in the form of absolute
value. While the method that normalizes the data matrix $\bm{\chi}$ reflects
the model approximation error in the form of relative values for each state.
Therefore, RMSE values of reduced-order models at different orders obtained
from POD with normalization are smaller than those obtained from POD. This can
also be further demonstrated from Fig. 5, which shows the trajectories of some
of states based on the original model and the reduced model with order 30
using POD with normalization and POD methods respectively. From the figure, it
can be seen that for states with large numerical values ($x_{17}$, $x_{101}$),
both methods provide very good approximation effects, but for states with
small numerical values or small fluctuations ($x_{3}$, $x_{11}$, $x_{26}$,
$x_{31}$), the approximation of POD with normalization is significantly better
than that of POD.

Figure 4: The values of $log$(RMSE) at different orders. Figure 5: The
trajectories of some of state based on the original model and the reduced
model using POD with normalization and POD without normalization.

Next, we will evaluate the open-loop prediction performance of the POD-MLP
model. Using the PRMS input signal, we generate 100,000 samples for each state
by utilizing the first-principle model with $\bm{U}_{r}$. These data are split
into training (70$\%$), validation (20$\%$), and testing (10$\%$) datasets.
The developed POD-MLP model consists of 3 hidden layers with 128 neurons, the
input of the POD-MLP is the normalized 33-dimensional vector $\bm{\xi}_{u}$
and the output is the normalized state vector $\bm{\xi}$. The activation
function for the hidden layers has been selected to be Tanh, and a linear
activation function is used in the output layer. Figure 6 demonstrates the
testing performance of the POD-MLP model in multi-step open-loop prediction
for the actual state trajectory of the PCCP. It can be seen that the fit is
accurate.

Figure 6: Multi-step prediction for the PCCP. Figure 7: The actual states and POD-MLP-EKF state estimates for the PCCP. TABLE III: Running time. | EKF | POD-EKF | POD-MLP-EKF
---|---|---|---
Model prediction | 8s | 8s | 3s
Discretization | 115s | 117s | 3s
Other steps for EKF | 2s | 2s | 0.4s
Total time (s) | 125s | 127s | 6.4s

In the following analysis, we focus on the state estimation of the PCCP using
the POD-MLP-EKF algorithm. The PCCP is subject to process disturbances, and
the output measurements are corrupted by random noise. Specifically, each
process disturbance sequence associated with the $i$th state $x_{i}$ is
generated following normal distribution with zero mean and a standard
deviation 0.01$x_{i,s}$, where $x_{i,s}$ is the steady-state value of $x_{i}$.
Random noise is added to each measurement $y_{i}$ as Gaussian white noise with
zero mean and a standard deviation 0.01$y_{i,s}$, where $y_{i,s}$ is the value
of $y_{i}$ at steady-state. As a result, the covariance matrices of process
noise and measurement noise are $\bm{Q}={\rm diag}((0.01\bm{x}_{s})^{2})$ and
$\bm{R}={\rm diag}((0.01\bm{y}_{s})^{2})$, respectively. The tuning parameters
in the POD-MLP-EKF are $\bm{Q}_{r}=\bm{U}_{r}^{\tiny\mbox{\rm
T}}\bm{Q}\bm{U}_{r}$ and $\bm{R}_{r}=\bm{U}_{r}^{\tiny\mbox{\rm
T}}\bm{R}\bm{U}_{r}$. The initial guess for the normalized states is set to
$0.5\bf 1_{103}$. Figure 7 shows some of the state estimates and the actual
states. The proposed estimation scheme provides accurate state estimates.

We compare the computational efficiency of the proposed POD-MLP-EKF algorithm
based on model reduction and neural network, a centralized EKF design directly
based on the PCCP model, and an EKF for POD model, denoted as the POD-EKF.
Specifically, we evaluate the average computation time for 600 steps required
by the three algorithms, which is 125 s, 127 s, and 6.4 s, respectively. The
results demonstrate that the combination of POD reduction and neural network
significantly reduces the computing time.

## ACKNOWLEDGMENT

The first author, S.Y. Liu, is a visiting Ph.D. student in the Department of
Chemical and Materials Engineering at the University of Alberta from March
2021 to February 2023. She acknowledges the financial support from the China
Scholarship Council (CSC) during this period.

## References

  * [1] N. MacDowell, N. Florin, A. Buchard, J. Hallett, A. Galindo, G. Jackson, C. S. Adjiman, C. K. Williams, N. Shah, and P. Fennell, “An overview of CO2 capture technologies,” Energy & Environmental Science, vol. 3, no. 11, pp. 1645–1669, 2010.
  * [2] N. A. Manaf, A. Qadir, and A. Abbas, “Efficient energy management of CO2 capture plant using control-based optimization approach under plant and market uncertainties,” Journal of Process Control, vol. 74, pp. 2–12, 2019\.
  * [3] X. Yin, B. Decardi-Nelson, and J. Liu, “Distributed monitoring of the absorption column of a post-combustion CO2 capture plant,” International Journal of Adaptive Control and Signal Processing, vol. 34, no. 6, pp. 757–776, 2020.
  * [4] Q. Wang, C. Zheng, X. Wu, and M. Wang, “Robust monitoring of solvent based carbon capture process using deep learning network based moving horizon estimation,” Fuel, vol. 321, p. 124071, 2022.
  * [5] S. Liu, X. Zhang, L. Xu, and F. Ding, “Expectation–maximization algorithm for bilinear systems by using the Rauch–Tung–Striebel smoother,” Automatica, vol. 142, p. 110365, 2022.
  * [6] X. Yin and J. Liu, “Distributed moving horizon state estimation of two-time-scale nonlinear systems,” Automatica, vol. 79, pp. 152–161, 2017\.
  * [7] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” IEEE Transactions on signal processing, vol. 50, no. 2, pp. 174–188, 2002\.
  * [8] P. R. Jeon, M. S. Hong, and R. D. Braatz, “Compact neural network modeling of nonlinear dynamical systems via the standard nonlinear operator form,” Computers & Chemical Engineering, vol. 159, p. 107674, 2022.
  * [9] X. Cheng, G. Li, R. Skulstad, S. Chen, H. P. Hildre, and H. Zhang, “A neural-network-based sensitivity analysis approach for data-driven modeling of ship motion,” IEEE Journal of Oceanic Engineering, vol. 45, no. 2, pp. 451–461, 2020.
  * [10] T. Zhao, Y. Zheng, J. Gong, and Z. Wu, “Machine learning-based reduced-order modeling and predictive control of nonlinear processes,” Chemical Engineering Research and Design, vol. 179, pp. 435–451, 2022.
  * [11] Y. Liang, H. Lee, S. Lim, W. Lin, K. Lee, and C. Wu, “Proper orthogonal decomposition and its applications—part i: Theory,” Journal of Sound and vibration, vol. 252, no. 3, pp. 527–544, 2002.
  * [12] G. Kerschen, J.-c. Golinval, A. F. Vakakis, and L. A. Bergman, “The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview,” Nonlinear dynamics, vol. 41, pp. 147–169, 2005.
  * [13] B. Decardi-Nelson, S. Liu, and J. Liu, “Improving flexibility and energy efficiency of post-combustion CO2 capture plants using economic model predictive control,” Processes, vol. 6, no. 9, p. 135, 2018.
  * [14] N. Harun, T. Nittaya, P. L. Douglas, E. Croiset, and L. A. Ricardez-Sandoval, “Dynamic simulation of MEA absorption process for CO2 capture from power plants,” International Journal of Greenhouse Gas Control, vol. 10, pp. 295–309, 2012.
  * [15] A. Zhang and J. Liu, “Economic mpc of wastewater treatment plants based on model reduction,” Processes, vol. 7, no. 10, p. 682, 2019.