# High efficiently numerical simulation of the TDGL equation with reticular
free energy in hydrogel

Jun Han  Hui Zhang  Zhengru Zhang School of Mathematical Sciences, Beijing
Normal University, Beijing, 100875, P.R. China .

###### Abstract

In this paper, we focus on the numerical simulation of phase separation about
macromolecule microsphere composite (MMC) hydrogel. The model equation is
based on Time-Dependent Ginzburg-Landau (TDGL) equation with reticular free
energy. We have put forward two $L^{2}$ stable schemes to simulate simplified
TDGL equation. In numerical experiments, we observe that simulating the whole
process of phase separation requires a considerably long time. We also notice
that the total free energy changes significantly in initial time and varies
slightly in the following time. Based on these properties, we introduce an
adaptive strategy based on one of stable scheme mentioned. It is found that
the introduction of the time adaptivity cannot only resolve the dynamical
changes of the solution accurately but also can significantly save CPU time
for the long time simulation.

Keywords: TDGL equation, $L^{2}$-norm stale, semi-implicit finite difference
schemes, adaptive time-stepping strategy

##  1 Introduction

Hydrogels, which are polymer with 3-D crosslinked hydrophilic network
structure, have increasingly extensive applications in industrial and
biomedical fields [8,9,22]. Therefore, hydrogels have recently received
tremendous attention in scientific communities due to their enhanced
properties [1,15]. Traditional hydrogels usually have poor mechanism as they
consist of water in a fragile polymer network and the cross-linking points in
the system are disordered leading to different chain length, which causes
shorter chains bearing more stress crack and the whole hydrogel breaks down.
However, in recent year, there has been a lot of novel methods to improve the
structure of hydrogels and thus significantly enhancing their mechanical
properties [8,11]. These methods primarily introduce special corss-linking
center in gel structure such as sliding crosslinking and multiple crosslinking
center, which can effectively disperse stress and thus the polymer is not easy
to fracture [9]. According to the structure of these novel hydrogels, they can
be classified as topological (TP) gels, nanocomposite (NC) hydrogels, double-
network (DN) gels, MMC hydrogels and so on [8,9]. TP gels have eight
crosslinkers that can slide along polymer chains [3] and can swells to about
500 times its original weight and can be stretched to nearly 20 times its
original length [14]. NC hydrogels, which utilize inorganic material to form
crosslinking junctions and adsorb strongly on surfaces of platelets with
enough chains attached to provide bridges, constitute a network structure
jielding unusual mechanical properties including very good toughness [7]. DN
gels generally consist of two independently cross-linked networks, one
consisting of a rigid polyelectrolyte and the other a flexible uncharged
polymer, and exhibit extremely high mechanical strength [4,6]. MMC hydrogels,
which were synthesized by Ting Huang in 2007 [8], have a unigue well-defined
microstructure and very high mechanical strength. MMC hydrogels are
environmentally sensitive and mainly possible to use in drug delivery and
other biomedical applications [20]. Much work in chemical structure and
dynamics simulation have been done about TP gels, NC hydrogels, and DN gels
and much progress has been made [5,21,10,11]. But there are relatively less
results about MMC hydrogels.

Meanwhile, many computational methods have been developed and applied in
modeling and simulation of phase separation of hydrogels. They range from
molecular scale (e.g., molecular dynamics, Monte Carlo), microscale (e.g.,
Brownian dynamics, dissipative particle dynamics, lattice Boltzmann, TDGL
method, dynamic density functional theory method) to mesoscale and macroscale
(e.g., micromechanics, equivalent-continuum and self-similar approaches,
finite element method) [21]. These methods are foundation to simulate phase
separation of hydrogels. However, observing the whole process of phase
separation is computationally expensive and costs a large amount of time when
these methods are used. Therefore, seeking for algorithms to decrease
computation cost is an urgent task. The adaptive time-stepping method is
introduced here. Adaptive time stepping has been well studied for solving
imitial value problems in ODEs. The review in [17,18] summarizes several step
control methods for local time adaptivity based on linear feedback theory. A
locally varying time step method is developed for solving hyperbolic
conservative PDEs where large time step is adopted on region with smooth
solution while small time step is taken in domain with nearly singular
solution within same time level [19]. In [16], the adaptive time-stepping
technique has been developed based on energy variation which is an important
physical quantity in the molecular beam epitaxy growth model. Following [16],
we improve the choice of adaptive time step according to properties of TDGL
equation.

In this paper, considering that MMC hydrogels have microstructure with network
crosslinking points, we exploit TDGL method to simulate phase evolution of MMC
hydrogels based on reticular free energy deduced in [22]. TDGL method is an
accurate and robust method to simulate structural evolution of phase
separation in polymer blends and copolymers and many nice results have been
gotten in numerical experiments [1,13,21]. There are two main objectives in
this paper. First, we put forward two $L^{2}$ stable second-order implicit
schemes to simulate the equation, which satisfy the property of total free
energy decreasing. One scheme needs to solve a high-dimensional matrix in each
step, which can be solved by conjugate gradient method. The other is a
nonlinear scheme where we need Newton iteration in each step, which will cost
more computation. Second, as the equation is unconditionally energy stale and
the schemes are $L^{2}$ norm stable, we introduce an adaptive strategy which
is vert effective for such evolutionary equations where total free energy
decays significantly in initial time and slightly later, called ”rough-smooth”
pattern. During the time with rapid free energy decay, small adaptive time
steps are chosen adaptively, while large time steps are used when free energy
is varied smoothly.

This paper is presented as follows. In section 2, TDGL equation is introduced
in detail with reticular free energy. Then, we employ variational principle to
derive and simplify the equation in a further form, which will be used in
following theoretical analysis and numerical simulation. In section 3,
analysis of the equation is given and two semi-implicit schemes are proposed.
In section 4, some numerical experiments are conducted to test one of the
numerical schemes. In section 5, based on observations, we put forward a
suitable adaptive time-stepping strategy and some concluding remarks are drawn
in final section.

##  2 TDGL equation

TDGL equation is a microscale method for simulating the structural evolution
of phase-separation in polymer blends and block copolymers. For an
incompressible polymer blend, the model equation can be written as follows
[1,13,21,22].

$\frac{\partial\phi(\textbf{x},t)}{\partial
t}=\nabla\cdot(M(\phi)\nabla\frac{\delta
U[\phi(\textbf{x},t)]}{\delta\phi(\textbf{x},t)})+\zeta(\textbf{x},t),$ (2.1)

where $\phi(\textbf{x},t)$ is order parameter, by definition $\phi\in(0,1)$,
the mobility $M$ may depend on $\phi$. $\textbf{x}\in{\Omega}$,
$\Omega=(0,2\pi)\times(0,2\pi)$. The total Flory-Huggins-de Gennes free energy
$U$ can be defined as:

$U[\phi(\textbf{x},t)]=k_{B}T\int_{\Omega}[\frac{F(\phi(\textbf{x},t))}{k_{B}T}+\kappa(\phi(\textbf{x},t))|\nabla\phi(\textbf{x},t)|^{2}]d\textbf{x},$
(2.2)

where $\kappa(\phi)=\frac{\sigma^{2}}{36\phi(1-\phi)}$, $k_{B}$ is Boltzmann
constant, $\sigma$ is the Huhn length of the polymer,$T$ is temperature
constant [2]. It is pointed out that for a binary polymer blend,
$F(\phi(\textbf{x},t))$ usually denotes the Flory-Huggins free energy.
However, as mentioned before, MMC hydrogels have network structure. Hence, in
[8], it is proven that the reticular free energy should be defined as follows.

$F=k_{B}T[\frac{\phi}{\tau}ln(\frac{\alpha\phi}{\tau})+\frac{\phi}{N}ln(\frac{\beta\phi}{\tau})+(1-\rho\phi)ln(1-\rho\phi)+\chi\rho\phi(1-\rho\phi)],$
(2.3)

where parameters $\tau$, $\alpha$, $\beta$, $N$, and $\rho$ are constant
variables [22]. $\chi$ is the enthalpic interaction parameter between two
polymer components. As derived in [22], it can be guaranteed that
$(1-\rho\phi)>0$.

The thermal noise in (2.1) is a random term with zero mean which should
satisfy following conditions:

$<\zeta(\textbf{x},t),\zeta(\textbf{x}^{\prime},t)>=\epsilon
M_{0}\nabla^{2}\delta(\textbf{x}-\textbf{x}^{\prime})\delta(t-t^{\prime}),$

where $\epsilon$ is the small magnitude of the fluctuation [12].

###  2.1 Analysis of energy stability

When the property of energy stability is considered, the noise term in (2.1)
should not be considered. Therefore, multiplying (2.1) (deleting noise term)
by $\frac{\delta U}{\delta\phi}$ and integrating over $\Omega$, we get

$\frac{dU}{dt}=\frac{1}{|\Omega|}\int_{\Omega}\frac{\delta
U(\phi)}{\delta\phi}\nabla\cdot(M(\phi)\nabla\frac{\delta
U(\phi)}{\delta\phi})d\textbf{x}.$

Here, $|\Omega|$ is the volume of $\Omega$. The above equality holds because
of the fact that $U$ is independent of variable x according to definition of
$U$ in (2.2). Then we can easily derive

$\frac{dU}{dt}=\frac{1}{|\Omega|}(\int_{\Omega}\nabla\cdot(\frac{\delta
U(\phi)}{\delta\phi}M(\phi)\nabla\frac{\delta
U(\phi)}{\delta\phi})d\textbf{x}-\int_{\Omega}M(\phi)\nabla\frac{\delta
U(\phi)}{\delta\phi}\cdot\nabla\frac{\delta U(\phi)}{\delta\phi}d\textbf{x}).$

By divergence theorem and considering periodic boundary condition for $\phi$,
the first term on the right side of above equality disappears. Hence, we get

$\frac{dU}{dt}=-\frac{1}{|\Omega|}\int_{\Omega}M(\phi)|\nabla\frac{\delta
U(\phi)}{\delta\phi}|^{2}d\textbf{x}\leq 0.$ (2.4)

Therefore, we know that the total free energy of whole system is decreasing
over time, i.e, the system is energy stable.

###  2.2 Deduction of the model equation

For any $\psi\in C^{1}(\Omega)$ and considering periodic boundary condition
for $\phi$ and $\psi$, we define

$H(\epsilon)=U[\phi+\epsilon\psi]=k_{B}T\int_{\Omega}[\frac{F(\phi(\textbf{x})+\epsilon\psi(\textbf{x}))}{k_{B}T}+\kappa(\phi(\textbf{x})+\epsilon\psi(\textbf{x}))|\nabla(\phi(\textbf{x})+\epsilon\psi(\textbf{x}))|^{2}]d\textbf{x}.$

Obviously, $H(\epsilon)\in C^{1}(\mathbb{R})$. Hence, $H^{\prime}(\epsilon)$
can be easily derived and when $\epsilon=0$, we get

$H^{\prime}(0)=k_{B}T\int_{\Omega}[\frac{F^{\prime}(\phi)\psi}{k_{B}T}+\kappa^{\prime}(\phi)\psi|\nabla\phi|^{2}+2\kappa(\phi)\nabla\phi\cdot\nabla\psi]d\textbf{x}.$
(2.5)

According to the Green’s formula and periodic boundary condition for $\phi$
and $\psi$, we derive

$\int_{\Omega}\kappa(\phi)\nabla\phi\cdot\nabla\psi
d\textbf{x}=\int_{\Omega}\nabla\cdot(\kappa(\phi)\psi\nabla\phi)d\textbf{x}-\int_{\Omega}\nabla\cdot(\kappa(\phi)\nabla\phi)\psi
d\textbf{x}$
$=-\int_{\Omega}[\kappa^{\prime}(\phi)\nabla\phi\cdot\nabla\phi+\kappa(\phi)\Delta\phi]\psi
d\textbf{x}.$

By the variational principle in (2.5), we get

$\frac{\delta
U[\phi]}{\delta\phi}=H^{\prime}(0)=F^{\prime}(\phi)-k_{B}T(\kappa^{\prime}(\phi)|\nabla\phi|^{2}+2\kappa(\phi)\Delta\phi).$
(2.6)

Combining (2.1) with (2.6), following [1,13,21,22] where the mobility $M$ is
equal to constant $M_{0}$, we get the model equation

$\frac{\partial\phi}{\partial t}=M_{0}\nabla\cdot(\nabla
F^{\prime}(\phi))-M_{0}k_{B}T\Delta(\kappa^{\prime}(\phi)|\nabla\phi|^{2})-2M_{0}k_{B}T\Delta(\kappa(\phi)\Delta\phi)+\zeta(\textbf{x},t).$
(2.7)

Straightforward computations in (2.3) jield $\nabla
F^{\prime}(\phi)=k_{B}T(\frac{1}{\tau\phi}+\frac{1}{N\phi}+\frac{\rho^{2}}{1-\rho\phi}-2\chi\rho^{2})\nabla\phi$.
Here, we denotes
$G(\phi)=k_{B}T(\frac{1}{\tau\phi}+\frac{1}{N\phi}+\frac{\rho^{2}}{1-\rho\phi}-2\chi\rho^{2})$.

In practical problems in materials, many chemists and physicists have taken
$\kappa(\phi)$ to be a constant parameter. Therefore, in practical simulation,
it is reasonable to assume that $\kappa(\phi)$ is constant, denoted by $K$.
Then, the model equation would be simplified as:

$\frac{\partial\phi}{\partial
t}=M_{0}G^{\prime}(\phi)|\nabla\phi|^{2}+M_{0}G(\phi)\Delta\phi-2M_{0}k_{B}TK\Delta^{2}\phi+\zeta(\textbf{x},t).$
(2.8)

##  3 $L^{2}$ stable semi-implicit finite difference schemes

As illustrated in [23], A system located at a metastable state initially will
fluctuate around this metastable state under thermal noise perturbation. When
the noise happens to be large enough to carry the system out of the basin of
attraction of this metastable state, the phase separation will occur. The role
of noise therm is to carry the system out of the basin of the metastable
state. Therefore, in theoretical analysis, we do not need to consider noise
term and the model equation can be rewritten:

$\frac{\partial\phi}{\partial
t}=M_{0}G^{\prime}(\phi)|\nabla\phi|^{2}+M_{0}G(\phi)\Delta\phi-2M_{0}k_{B}TK\Delta^{2}\phi.$
(3.9)

Before we construct finite difference schemes to solve equation (2.8), we
first prove that the analytic solution of (2.8) is $L^{2}$ stable if it
existed, which can be rewritten in following theorem.

THEOREM 3.1. If $\phi(\textbf{x},t)$ is a solution of Eq. (3.9), then
$\phi(\textbf{x},t)$ is $L^{2}$ stable with respect to the initial condition:

$\|\phi(\cdot,t)\|^{2}\leq\|\phi(\cdot,0)\|^{2},$ (3.10)

where $\|\cdot\|$ is the standard $L^{2}$-norm in $\Omega$.

Proof: Multiplying (3.9) by $\phi$ and integrating over $\Omega$, we have

$\frac{1}{2}\frac{d}{dt}\int_{\Omega}|\phi|^{2}d\textbf{x}=M_{0}\int_{\Omega}G^{\prime}(\phi)\phi|\nabla\phi|^{2}d\textbf{x}+M_{0}\int_{\Omega}G(\phi)\phi\Delta\phi
d\textbf{x}-2M_{0}k_{B}TK\int_{\Omega}\phi\Delta^{2}\phi d\textbf{x}.$

By the Green’s formula for the second term and third term on the right side of
the above equality and considering the periodic boundary condition for $\phi$,
we get

$\frac{1}{2}\frac{d}{dt}\int_{\Omega}|\phi|^{2}d\textbf{x}=M_{0}\int_{\Omega}G^{\prime}(\phi)\phi|\nabla\phi|^{2}d\textbf{x}-M_{0}\int_{\Omega}\nabla(G(\phi)\phi)\cdot\nabla\phi
d\textbf{x}-2M_{0}k_{B}TK\int_{\Omega}|\Delta\phi|^{2}d\textbf{x}.$

After simple calculations, we derive

$\frac{1}{2}\frac{d}{dt}\int_{\Omega}|\phi|^{2}d\textbf{x}=-M_{0}\int_{\Omega}G(\phi)|\nabla\phi|^{2}d\textbf{x}-2M_{0}k_{B}TK\int_{\Omega}|\Delta\phi|^{2}d\textbf{x}.$
(3.11)

According to the value of parameter $\rho,\chi,N$, and $\tau$, it is obvious
that $G(\phi)\geq 0$, Therefore, the following inequality holds:

$\frac{d}{dt}\int_{\Omega}|\phi|^{2}d\textbf{x}\leq 0.$

Hence, inequality (3.10) is proven.

###  3.1 A second-order semi-implicit linear scheme

In the following, we construct finite difference schemes to solve the Eq.
(3.9). First, we give a semi-implicit linear scheme with second-order as
follows.

$\frac{\phi_{i,j}^{n+1}-\phi_{i,j}^{n}}{\Delta
t}=M_{0}G^{\prime}(\phi_{i,j}^{n})\frac{|\nabla_{h}\phi_{i,j}^{n}|^{2}+\nabla_{h}\phi_{i,j}^{n}\cdot\nabla_{h}\phi_{i,j}^{n+1}}{2}+M_{0}G(\phi_{i,j}^{n})\frac{\Delta_{h}\phi_{i,j}^{n}+\Delta_{h}\phi_{i,j}^{n+1}}{2}$
$-2M_{0}k_{B}TK\frac{\Delta_{h}^{2}\phi_{i,j}^{n+1}+\Delta_{h}^{2}\phi_{i,j}^{n}}{2}.$
(3.12)

$\nabla_{h}$,$\Delta_{h}$ are central-scheme operators of finite difference
method, and h is the mesh size of the spatial discretizition. Obviously, in
each step, (3.11) can be solved by computing a linear system of equations.
Thus, we call it as a linear scheme.

THEOREM 3.2. let ${\phi_{i,j}^{n}}$ be the numerical solution of (3.9) at
$t_{n}$. Then we have

$\parallel\phi^{n}\parallel_{h}^{2}\leq\parallel\phi^{0}\parallel_{h}^{2},$
(3.13)

where $\parallel\cdot\parallel_{h}$ is the discrete $L^{2}$-norm.

Proof. Multiplying (3.12) with the rewritten form by
$\frac{\phi_{i,j}^{n+1}+\phi_{i,j}^{n}}{2}$, we get

$\frac{(\phi_{ij}^{n+1})^{2}-(\phi_{ij}^{n})^{2}}{2\Delta
t}=M_{0}G^{\prime}(\phi_{i,j}^{n})\nabla_{h}\phi_{i,j}^{n}\frac{\nabla_{h}\phi_{i,j}^{n}+\nabla_{h}\phi_{i,j}^{n+1}}{2}\frac{\phi_{i,j}^{n+1}+\phi_{i,j}^{n}}{2}$
$+M_{0}G(\phi_{i,j}^{n})\frac{\phi_{i,j}^{n+1}+\phi_{i,j}^{n}}{2}\frac{\Delta_{h}\phi_{i,j}^{n}+\Delta_{h}\phi_{i,j}^{n+1}}{2}-2M_{0}k_{B}TK\frac{\Delta_{h}^{2}\phi_{i,j}^{n+1}+\Delta_{h}^{2}\phi_{i,j}^{n}}{2}\cdot\frac{\phi_{i,j}^{n+1}+\phi_{i,j}^{n}}{2}.$

Summing for $i=1,2,3,...,N_{x}$ and $j=1,2,...,N_{y}$. By the Green’s formula
and considering periodic boundary condition for $\phi$, the above equality can
be rewritten:

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})=(M_{0}G^{\prime}(\phi^{n})\nabla_{h}\phi^{n}\frac{\nabla_{h}\phi^{n}+\nabla_{h}\phi^{n+1}}{2},\frac{\phi^{n+1}+\phi^{n}}{2})_{h}+$
$(M_{0}G(\phi^{n})\frac{\phi^{n+1}+\phi^{n}}{2},\frac{\Delta_{h}\phi^{n}+\Delta_{h}\phi^{n+1}}{2})_{h}-2M_{0}k_{B}TK\parallel\frac{\Delta_{h}(\phi^{n+1}+\phi^{n})}{2}\parallel_{h}^{2}.$

where $(\cdot,\cdot)_{h}$ denotes the discrete $L^{2}$ inner product. By
Green’s formula for second term on right side above and considering periodic
boundary condition for $\phi$, the above equality becomes

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})=(M_{0}G^{\prime}(\phi^{n})\nabla_{h}\phi^{n}\frac{\nabla_{h}\phi^{n}+\nabla_{h}\phi^{n+1}}{2},\frac{\phi^{n+1}+\phi^{n}}{2})_{h}$
$-(M_{0}\nabla_{h}(G(\phi^{n})\frac{\phi^{n+1}+\phi^{n}}{2}),\frac{\nabla_{h}\phi^{n}+\nabla_{h}\phi^{n+1}}{2})_{h}-2M_{0}k_{B}TK\parallel\frac{\Delta_{h}(\phi^{n+1}+\phi^{n})}{2}\parallel_{h}^{2}.$

The second term on right side above can be decomposed to two terms, one of
them can be varnished by the first term on right side. Thus, we can get

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})=-(M_{0}G(\phi^{n}),\frac{|\nabla_{h}\phi^{n+1}+\nabla_{h}\phi^{n}|^{2}}{2})_{h}$
$-2M_{0}k_{B}TK\parallel\frac{\Delta_{h}(\phi^{n+1}+\phi^{n})}{2}\parallel_{h}^{2}.$

Consequently, the right side above is less than 0, we get

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})\leq
0.$

Hence, inequality (3.13) holds.

###  3.2 A semi-implicit nonlinear scheme with second-order

Similar to the linear scheme, we put forward a nonlinear scheme where the
equation should be solved by the Newton’s iteration.

$\frac{\phi_{i,j}^{n+1}-\phi_{i,j}^{n}}{\Delta
t}=M_{0}G^{\prime}(\frac{\phi_{i,j}^{n+1}+\phi_{i,j}^{n}}{2})|\nabla_{h}\frac{\phi_{i,j}^{n}+\phi_{i,j}^{n+1}}{2}|^{2}+M_{0}G(\frac{\phi_{i,j}^{n+1}+\phi_{i,j}^{n}}{2})\frac{\Delta_{h}\phi_{i,j}^{n}+\Delta_{h}\phi_{i,j}^{n+1}}{2}$
$-2M_{0}k_{B}TK\frac{\Delta_{h}^{2}\phi_{i,j}^{n+1}+\Delta_{h}^{2}\phi_{i,j}^{n}}{2}.$
(3.14)

THEOREM 3.3. let ${\phi_{i,j}^{n}}$ be the numerical solution of (3.14) at
$t_{n}$. Then we have

$\parallel\phi^{n}\parallel_{h}^{2}\leq\parallel\phi^{0}\parallel_{h}^{2}.$
(3.15)

Proof. Similar to the proof of Th.3.2, multiply (3.14) by
$\frac{\phi_{i,j}^{n+1}+\phi_{i,j}^{n}}{2}$ and summing for
$i=1,2,3,...,N_{x}$ and $j=1,2,...,N_{y}$. By the Green’s formula and periodic
boundary condition for $\phi$ is considered, we get

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})=(M_{0}G^{\prime}(\frac{\phi^{n}+\phi^{n+1}}{2})|\frac{\nabla_{h}\phi^{n}+\nabla_{h}\phi^{n+1}}{2}|^{2},\frac{\phi^{n+1}+\phi^{n}}{2})_{h}+$
$(M_{0}G(\frac{\phi^{n}+\phi^{n+1}}{2})\frac{\phi^{n+1}+\phi^{n}}{2},\frac{\Delta_{h}\phi^{n}+\Delta_{h}\phi^{n+1}}{2})_{h}-2M_{0}k_{B}TK\parallel\frac{\Delta_{h}(\phi^{n+1}+\phi^{n})}{2}\parallel_{h}^{2}.$

By the Green’s formula for the second term on right side above and considering
the periodic boundary condition for $\phi$, the above equality becomes

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})=(M_{0}G^{\prime}(\frac{\phi^{n}+\phi^{n+1}}{2})|\frac{\nabla_{h}\phi^{n}+\nabla_{h}\phi^{n+1}}{2}|^{2},\frac{\phi^{n+1}+\phi^{n}}{2})_{h}$
$-(M_{0}\nabla_{h}(G(\frac{\phi^{n}+\phi^{n+1}}{2})\frac{\phi^{n+1}+\phi^{n}}{2}),\frac{\nabla_{h}\phi^{n}+\nabla_{h}\phi^{n+1}}{2})_{h}-2M_{0}k_{B}TK\parallel\frac{\Delta_{h}(\phi^{n+1}+\phi^{n})}{2}\parallel_{h}^{2}.$

Similarly, the second term on the right side above can be decomposed to two
terms, one of which is varnished by first term right side. Therefore, we
derive

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})=-(M_{0}G(\frac{\phi^{n+1}+\phi^{n}}{2}),\frac{|\nabla_{h}\phi^{n+1}+\nabla_{h}\phi^{n}|^{2}}{2})_{h}$
$-2M_{0}k_{B}TK\parallel\frac{\Delta_{h}(\phi^{n+1}+\phi^{n})}{2}\parallel_{h}^{2}.$

Therefore, for $\phi\in[0.5,1)$, the right side above is less than 0, we get

$\frac{1}{2\Delta
t}(\parallel\phi^{n+1}\parallel_{h}^{2}-\parallel\phi^{n}\parallel_{h}^{2})\leq
0.$

Hence, inequality (3.15) holds and the proof is complete.

##  4 Numerical simulation

In this section, we will conduct some numerical experiments. As the nonlinear
scheme (14) which requires Newton’s iteration in each step is computationally
a very demanding task, we exploit the second-order semi-implicit linear scheme
to solve equation (8) with difference scheme given in (9). The noise term is
just introduced in initial time for carrying the system out of metastable
state. After we have done simulation for some time, the noise term will not be
used, i.e., we employ the semi-implicit scheme (12) to solve equation (9) in
order to simulate phase separation of the system. The mesh used is square
grids of $N_{x}=64$ and $N_{y}=64$. We set $M_{0}=0.2$, $\chi=0.4$,
$\tau=10^{7}$, $N=800$ and $\rho=1.0$ as in [22]. The tolerance of the
conjugate gradient method is set to be $10^{-6}$. We will conduct experiments
under constant time step $\Delta t=0.001$. At first, we set $T=1$, and the
initial values $\phi$ are random numbers around 0.65. we get the phase
evolution in Figure 4.1.

Figure 4.1: phase evolution figures at t=0.1, 10, 50, 100, 500, 1000.

From Figure 4.1, we can see phase evolution and obviously observe the phase
separation. At first, hydrogel is distributed out of order in the system. As
the time evolves, hydrogel gradually has well-defined structure, which will
greatly improve its mechanical strength.

In the following, we change initial value $\phi$. $\phi$ are chosen as random
numbers around 0.35. We get the phase evolution in Figure 4.2.

Figure 4.2: phase evolution figures at t=0.1, 10, 100, 500, 1000, 2000.

In Figure 4.2, it is showed that the phase separation of hydrogel is observed.
At first, the initial value $\phi$ is considerably disordered, indicated in
first phase picture in Figure 4.2. In the end, it can still evolve to a steady
state with regular structure.

Meanwhile, we will change the temperature constant parameter to do some
numerical experiment. The initial values $\phi$ are the same as initial values
given in Figure 4.1. We just increase temperature constant to $T=20$ and T=50
respectively. The values of other parameters are same as values given above.
We have the steady phase state in Figure 4.3.

Figure 4.3: t=1000, phase figures when $T=20,50$.

Comparing Figure 4.1 to figure 4.3, we can find that even though we have the
same initial values for $\phi$, different temperature constants can derive
different steady phase figures. All of these experiments where clear phase
separation has been observed demonstrate that the scheme (LABEL:12)
constructed is reasonable. It can simulate the process of phase evolution over
long time simulation.

But during each step of numerical simulation, we should solve a linear system
of equations in very high dimensions, for example, in our simulation the
matrix is 4096 orders. Although the matrix is a sparse diagonally-dominant
matrix, which can be solved by the conjugate gradient method in order to
reduce the computation, it still costs a large amount of time when we will
conduct a long time simulation. Especially, as illustrated in [22], in
industrial applications, we have to enlarge the computational area. If we
enlarge the area to $(0,4\pi)\times(0,4\pi)$, the dimensional of matrix is 4
times than that above when we have to employ the same mesh size to ensure
accuracy. Therefore, the computation is unacceptable for long time simulation.
This stipulates us to seek method to decrease computational cost. During the
simulation, we find that the total free energy has a significant variation in
initial time and has a micro change in late time, which means the solution
also has such property. Hence, it allows us to introduce adaptive time step as
the scheme (3.12) keeps the property of energy decreasing and is $L^{2}$
stable [3].

##  5 An adaptive time stepping method

In section 2.1, we have demonstrated that the whole system is unconditionally
energy stale. And The semi-implicit linear scheme (3.12) and nonlinear scheme
are $L^{2}$ stable. These properties of stability allow us to introduce large
time step during numerical simulations [16, 17, 18]. For the sake of accuracy,
a very large time step is unacceptable except in the time intervals where
solutions change considerably little. In the following simulations, we will
conduct numerical simulations under different constant time steps, such as
$\Delta t=0.001$, $\Delta t=0.005$, $\Delta t=0.01$ and $\Delta t=0.05$. In
the following simulations, the initial values are the same values as given in
figure 1. In order to test the numerical accuracy, we take the numerical
solution obtained using $\Delta t=10^{-4}$ as ”exact” solution. We define
$L^{2}$ relative error (RE) as

$RE=(\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}(\phi^{\Delta
t}_{t}(i,j)-\phi^{\Delta t_{0}}_{t}(i,j))^{2}h_{x}h_{y})^{\frac{1}{2}},$
(5.16)

where $\phi^{\Delta t_{0}}_{t}(i,j)$ is the ”exact” solution $\phi$ we assume
at grid $(i,j)$ under $\Delta t=0.001$, and similarly, $\phi^{\Delta
t}_{t}(i,j)$ is the solution at $\Delta t$($\Delta t\geq 0.001$). Here
$h_{x}=\frac{2\pi}{N_{x}}$, $h_{y}=\frac{2\pi}{N_{y}}.$ We set constant
temperature T=1. The following figures show total free energy evolution over
time:

Figure 5.4: Total energy evolution over time at T=1 under different constant
time step.

From Figure 5.4, we can see total free energy under different constant time
steps are all decreasing over time. This indicates that numerical solutions
under different constant time steps are consistent with the property of energy
stability analyzed in section 2.1. Although there are errors at first, as time
evolves, solutions under different constant time steps are the same
ultimately. In initial time, solutions under larger constant time steps have
higher errors than solutions under small constant time steps. This
demonstrates that the systems are sensitive to time steps chosen at initial
time. However, we can also see from figures that the total free energy
decreases significantly at initial time and changes minor after some time. As
we have proven that the whole system is energy stable and the scheme 3.12 is
$L^{2}$ stable, this allows us to choose large time step without resulting in
much error.

At the same time we will conduct some numerical experiments under different
temperature to test the robustness of observations above. In the following
numerical simulations, we choose $T=50$.

Figure 5.5: Total energy evolution over time at T=1 under different constant
time step.

From Figure 5.5, we can see that when temperature increases, large constant
time step will deviate the energy curve from the energy curve obtained by
”exact” solution, which illustrate that the solution obtained by larger
constant time step is much higher than that of solution with small temperature
under the same initial value condition . This indicates that numerical
solutions are more sensitive to constant time steps. Therefore, we should be
more concerned about time steps when we numerically solve the equation.

In the following, we will introduce an adaptive time step to numerically solve
equations. The basic idea of choosing adaptive time step is that when the
total free energy decreases significantly, which mean solutions have a large
variation, we employ small time step and when total free energy has a micro
change, we will choose large time steps for it will not bring about large
error under this circumstance. In this way, small errors can be guaranteed
during the whole process.

Following [3], according to the choice of adaptive time step and definition of
total free energy in this equation, we improve the choice of adaptive time
step at time t as:

$\Delta t=\max(\Delta t_{min},\lambda(T,t)\frac{\Delta
t_{max}}{\sqrt{1+\mu|U^{\prime}(t)|^{2}}}).$

Here $\mu$ is a constant parameter, $U^{\prime}(t)$ is the derivative of
$U(t)$ and $\lambda(T,t)$ is introduced to ensure accuracy and efficiency. For
the same change of total free energy, we will choose small time step at first
for accuracy but it is safe to use a larger time step in late time for the
system evolves to a steady state. The total free energy is proportional to
temperature T. For the same change of $\phi$, the derivative of
$\mu|U^{\prime}(t)|^{2}$ at $T=50$ is $2500$ times than the derivative of it
at T=1. In this case, for the same change of $\phi$, the choice of adaptive
time step is considerably small. So we have to multiply $\lambda(T,t)$ to
ensure the time step is not too small.

In the following simulation, we choose $\Delta t_{min}=0.001,\Delta
t_{max}=0.1,\mu=1000$. and the definition of $\lambda(T,t)$ is as follows.

$\lambda(T,t)=\begin{cases}1,\quad if\quad 0\leq t\leq 100.\\\ 1.5,\quad
if\quad 100<t\leq 200.\\\ 2,\quad if\quad 200<t\leq 300.\\\ 3,\quad if\quad
300<t\leq 400.\\\ 4,\quad if\quad 400<t\leq 500.\\\ 5,\quad if\quad
t>500\end{cases}.$

First, we list the CPU time consumed at some selected time levels $t=1$,
$t=10$, $t=100$, $t=1000$ respectively under different time steps.

CPU time costs under different time steps at $t=1$, $t=10$, $t=100$, $t=500$.

Running time | t=1 | t=10 | t=100 | t=500
---|---|---|---|---
$\Delta t=0.05$ | 47.8 | 388.6 | 3416.9 | 14341.4
$\Delta t=0.01$ | 149.6 | 905.6 | 4190.2 | 5987.2
$\Delta t=0.005$ | 151.8 | 948.5 | 2335.1 | 3920.5
$\Delta t=0.001$ | 184.5 | 419.7 | 2125.2 | 10485.2
The adaptive time step | 110.1 | 202.7 | 668.8 | 2809.4

From the table above, we can see that the CPU time consumed under adaptive
time steps is much less than that of time under other constant time steps.
Meanwhile, we should consider the accuracy of solutions evolved. At first, we
see the

Figure 5.6: Total free energy evolution over time, blue:$\Delta t=0.001$, red:
the adaptive time step.

From Figure 5.6, we can see that the difference of total free energy evolution
is very small under constant time step $\Delta t=0.001$ and the adaptive time
step. The main reasons are as follows. At the initial time with significant
change of total free energy, the choice of time step in the adaptive method is
almost the same with that of time step in constant time step simulation
$\Delta t=0.001$. However, as the system evolves for some time, the free
energy changes little indicating that the change of $\phi$ is very small.
Under this circumstance, we can choose a large time step which might be $100$
times than $0.001$ without bringing about much error.

Also, we will give $L^{2}$-errors to test the accuracy of different time steps
at some selected time levels. In the following, we give comparison of
$L^{2}$-errors of different time steps both constant time steps and the
adaptive time step at time levels,$t=1$, $t=10$, $t=100$, $t=1000$
respectively.

$L^{2}$ error under different constant time steps and the adaptive time step
at $t=1$, $t=10$, $t=100$, $t=200$.

$L^{2}$-error | t=1 | t=10 | t=100 | t=200
---|---|---|---|---
$\Delta t=0.05$ | 0.031 | 0.0286 | 0.0256 | 0.0376
$\Delta t=0.01$ | 0.0267 | 0.0152 | 0.0086 | 0.0083
$\Delta t=0.005$ | 0.0198 | 0.0058 | 0.0037 | 0.0056
The adaptive time step | 0.0095 | 0.0042 | 0.000987 | 0.000325

From the $L^{2}$ errors given above, we can see that larger constant time step
have higher errors which is consistent with analysis given above. The
$L^{2}$-error with adaptive time strategy is considerably small during the
whole process of simulation, which indicates that the adaptive method is not
only efficient but also accurate.

##  6 Concluding remarks

In this paper, we have demonstrated the simplified version of TDGL equation is
$L^{2}$ stable. At the same time, we put forward two second-order semi-
implicit schemes satisfying energy decreasing property, which are $L^{2}$
stale and consistent with the analytic property of $L^{2}$ stability. We also
have conducted some numerical experiments to test the robustness of the
schemes. Some numerical results are shown above using one of schemes
constructed. From the numerical results, we can clearly see the phase
evolution of MMC hydrogels and the process of phase separation has been
observed during the experiments, which indicates effectiveness of the
numerical scheme (3.12).

Meanwhile, during the simulation, we find that the total free energy of the
system changes significantly at initial time and varies little after some time
simulation, which mean that numerical solutions have a palpable variation at
first and have a mirror change in late time. As the system is energy stale and
the schemes are $L^{2}$ stable with energy decreasing property, which
guarantees that large time step can be introduced at late time without bring
about much error. Hence we implement the adaptive time-stepping method to
simulate the equation. The numerical results show the efficiency and accuracy
of the adaptive time strategy. These results indicates that simplified TDGL
equation is very suitable for simulating the phase evolution of MMC hydrogels
after introducing the adaptive time strategy.

Although we have gotten some reasonable results to show that scheme (3.12) is
a robust scheme to simulate the phase separation, much further work should be
done to test the accuracy of the scheme. As spectral method is high resolution
method to numerically solve equations with periodic boundary condition, we can
exploit the spectral method to compute (2.8) and make a comparison between
solutions obtained with scheme (3.12) and solutions with spectral method under
the same initial value condition. In our further effort, we will employ the
spectral method to solve the equation and investigate more about TDGL equation
in its initial form.

## References

  * [1] A. Chakrabarti, R. Toral, J. D. Gunton, et al; Dynamics of phase separation in a binary polymer blend of critical composition, J. Chem. Phys 1990. 92, 6899.
  * [2] M. Fermeglia, S. Pricl; Multiscale modeling for polymer systems of industrial interest, Progress in organic coatings, 2007, 58: 187-199
  * [3] H.W. Gibson, M.C. Bheda, P.T. Engen; Rotaxanes, catenanes, polyrotaxanes,polycatenanes and related materials. Prog Polym Sci 1994;19:843 C945.
  * [4] J.P. Gong, Y. Katsuyama, T. Kurokawa, et al; Double network hydrogels with extremely high mechanical strength. Adv Mater 2003;15:1155 C8.
  * [5] W.W. Graessley, D.S. Pearson; Stress Cstrain behavior in polymer networks containing nonlocalized junctions. J Chem Phys 1977;66:3363 C70.
  * [6] M. Huang, H. Furukawa, Y. Tanaka, et al; Importance of entanglements between first and second components in high-strength double network gels. Macromolecules 2007;40:6658 C64.
  * [7] K. Harag, T. Takehisa; Composite hydrogels: a unique organic Cinorganic network structure with extraordinary mechanical, optical, and swelling/deswelling properties. Advance Matter.2002;14:1120-4.
  * [8] T. Huang, H.G. Xu, K.X. Jiao, et al; A novel hydrogel with high mechanical strength: A macromolecular microsphere composite hydrogel, Advanced Materials, 2007, 19, 1622-1626.
  * [9] J.A. Johnson, N.J. Turro, J.T. Koberstein, et al; Some hydrogel having novel molecular structure, Progress in Polymer Science, 2010, 35:332-337.
  * [10] S.S. Jang, M. Yashaar, S. Kalani, et al; Mechanical and transport properties of the poly(ethylene oxide) Cpoly(acrylic acid) double network hydrogel from molecular dynamic simulations, J Phys Chem B 2007;111:1729 C37.
  * [11] S.S. Jang, M. Yashaar, S. Kalani, et al; Corrections to mechanical and transport properties of the poly(ethylene oxide) Cpoly(acrylic acid) double network hydrogel from molecular dynamic simulations. J Phys Chem B 2007;111,14440/1.
  * [12] K. Luo; The morphology and dynamics of polymerization of polymerization-induced phase separation, European Polymer Journal, 2006, 42: 1499-1505.
  * [13] D. Nwabunma, H.W. Chiu, T. Kyu; Theoretical investigation on dynamics of photopolymerization-induced phase separation and morphology development in nematic liquid crystal/polymer mixtures, J. Chem. Phys (2000). 113, 6429.
  * [14] Y. Okumura, K. Ito, The Polyrotaxane Gel: A Topological Gel by Figure-of-Eight Cross-links, Advanced matter. 2001, 13, 485.
  * [15] A. Okada, A. Usuki; Twenty years of polymer-clay nanocomposites. Macromol Mater Eng 2006;291:1449 C76.
  * [16] Z.H. Qiao, Z.R Zhang, T. Tang; An Adaptive time-stepping strategy for the Molecular Beam Epitaxy Models, SIAM J.SCI.COMPUT Vol. 33, No. 3, pp. 1395 C1414.
  * [17] G. Soderlind; Automatic control and adaptive time-stepping,Numer. Algorighms, 31 (2002),pp. 281 C310.
  * [18] G. Soderlind, L. Wang; Adaptive Time-Stepping and Computational Stability, Journal of Computational Methods,in Sciences and Engineering,vol. 2, no. 3, 2, pp. 1 C3.
  * [19] Z. J. Tan, Z. R. Zhang, Y. Q. Huang et al; Moving mesh methods with locally varying time steps, J. Comput. Phys., 200 (2004), pp. 347 C367.
  * [20] X.C. Xiao, L.Y. Chu, W.M. Chen, et al; Monodispersed thermosponsive hydrogel micropsheres with a volume phase transition driven by hydrogen boning[J]. Ploymer 2005, 46(9):3199-3209.
  * [21] Q.H. Zeng, A.B. Yu, G.Q. LU; Multiscale modeling and simulation of polymer nanocomposite, Progress in Polymer Science, 2008, 33: 199-269.
  * [22] D. Zhai, H. Zhang; Investigation on the application of the TDGL equation in macromolecular microsphere composite hydrogel, Soft Matter RSC Publishing, 2012. 10. 1039.
  * [23] W. Zhang, T.J. Li, P.W. Zhang; Numercial Study for the Nucleation of one-dimensional Stochastic Cahn-Hillard Dynamics, to appear in Commun. Math. Sci., (2012)


# Using Spectral Submanifolds for Nonlinear Periodic Control

Florian Mahlknecht1, John Irvin Alora1, Shobhit Jain2,
Edward Schmerling1, Riccardo Bonalli3, George Haller2, Marco Pavone1 1
Department of Aeronautics and Astronautics, Stanford University, Stanford, CA,
94305, USA {florianm, jjalora, schmrlng<EMAIL_ADDRESS>Institute of
Mechanical Systems, ETH Zurich, 8092 Zurich, Switzerland {shjain,
<EMAIL_ADDRESS>of Signals and Systems, University of Paris-
Saclay, CNRS, CentraleSupélec<EMAIL_ADDRESS>
###### Abstract

Very high dimensional nonlinear systems arise in many engineering problems due
to semi-discretization of the governing partial differential equations, e.g.
through finite element methods. The complexity of these systems present
computational challenges for direct application to automatic control. While
model reduction has seen ubiquitous applications in control, the use of
nonlinear model reduction methods in this setting remains difficult. The
problem lies in preserving the structure of the nonlinear dynamics in the
reduced order model for high-fidelity control. In this work, we leverage
recent advances in Spectral Submanifold (SSM) theory to enable model reduction
under well-defined assumptions for the purpose of efficiently synthesizing
feedback controllers.

##  I INTRODUCTION

Automatic control of complex, infinite-dimensional systems (i.e., dynamically
evolving continua) such as soft robots as well as aircraft and underwater
vehicles with coupled fluid-structure interactions remains challenging. Model
reduction provides a principled approach to reduce model complexity while
capturing the essential physics required for controller synthesis. In optimal
control, we are interested in optimizing over a set of control inputs to track
a desired trajectory or stabilize around an operating point. In these
settings, working with the full-order model (FOM) is computationally
intractable. The need to utilize high-fidelity models to control these
challenging systems has resulted in significant research efforts towards
application of reduced-order models (ROM) for controller design.

In this work we explore recent developments in Spectral Submanifold (SSM)
theory [haller2016nonlinear] for nonlinear model reduction and control. SSMs
are the smoothest invariant manifolds that act as nonlinear continuations of
the eigenspaces from linearization of a system at a fixed point. This
nonlinear continuation is tangent to a corresponding spectral subspace of the
linearized system. The additional structure given by this continuation allows
us to capture highly-nonlinear behavior outside the vicinity of the linear
approximation. Under certain conditions, these nonlinearities can be
approximated arbitrarily well without ever increasing the size of the ROM.

Contributions: Motivated by the established theory on SSMs and their
successful application to model reduction [jain2021compute], we propose the
adaptation of SSM theory to automatic control. In particular, we aim to
synthesize optimal, nonlinear, periodic, feedback controllers on the reduced-
order SSM to exploit the computational speed-up while retaining the fidelity
of the high-order model.

Our contributions are threefold:

  1. (i)

We present, for the first time, this novel model reduction technique in the
context of control, enabling us to track (quasi-)periodic trajectories. In
this setting we restrict ourselves to a special case of trajectory tracking
and synthesize optimal control laws that guarantee the existence and
persistence of lower-dimensional SSMs on which the true system trajectory
lies.

  2. (ii)

We demonstrate the utility of SSMs in accurately capturing the nonlinear “slow
dynamics” of a system while neglecting its “fast dynamics”. This gives a
natural setting in which our control effort is focused on the dynamics that
persist.

  3. (iii)

We illustrate our approach on a pedagogical example, highlight its advantages
and disadvantages, and motivate its application to high-dimensional models.

Related work:  Most of the applications of model reduction for control exploit
projection-based methods. They involve a data-driven procedure to identify a
linear subspace from simulation rollouts of the FOM. This approach has been
successfully leveraged in literature for real-time control of infinite-
dimensional systems using Model Predictive Control (MPC). The work in
[LorenzettiMcClellanEtAl2021] considers the effect of proper orthogonal
decomposition (POD) on the closed-loop error dynamics and the authors develop
a constraint-tightening scheme to ensure satisfaction of safety constraints in
an MPC framework. While similar works [ghiglieri2014optimal,
alla2015asymptotic, altmuller2014model] adopt POD-based MPC schemes for the
control of certain PDE-classes, others explore different combinations of
subspace identification and optimal control schemes. In [alla2017error], the
authors investigate the sub-optimality of LQR due to the projection error
introduced by POD, while [antil2010domain] considers linear quadratic optimal
control using balanced truncation.

Direct application of projection-based methods for nonlinear systems is
difficult since evaluation of the nonlinear terms results in a more expensive
procedure than evaluating the full model directly, as the change of reference
frame involves high-dimensional matrix multiplications [farhat2015structure].
To overcome this limitation, much of the literature involves construction of
locally approximating linear ROMs for which standard linear control techniques
can be applied. In [huang2020balanced], the authors propose an iterative LQR
scheme combined with balanced truncation to control the 1D Burger’s Equation.
The authors in [TonkensLorenzettiEtAl2021] apply POD in a piecewise-affine
fashion by reducing linear approximations of the high-fidelity model. They
then evaluate the nonlinearities through interpolation of the linear
approximations and apply an MPC framework to control a soft robot.

While these approaches have been demonstrated to work well on various real-
world systems, their performance and theoretical guarantees are limited to
linear ROMs [LorenzettiMcClellanEtAl2021]. For highly nonlinear systems
operating in less constrained workspaces, linear ROMs can result in low-
fidelity surrogates that exhibit poor closed-loop performance and even
instability. This necessitates the need to capture the structure of the
nonlinearities in a more direct way and we propose a new direction for
addressing nonlinear model reduction for control.

Organization: In Section II we introduce our notation and definitions used in
this work. Section III defines the optimal control problem where we introduce
the tracking error in the periodic orbit of the FOM. In the methodology in
Section IV, we describe how we achieve model reduction using SSM in our
setting and how we leverage the reduced representation of the dynamics to
optimize the tracking error previously introduced in the full-order state-
space. LABEL:sec:validation showcases the application of the theory on an
illustrative example and provides the insights needed for tackling higher-
dimensional examples. We conclude this work in LABEL:sec:conclusions,
highlighting the most promising future avenues.

##  II PRELIMINARIES

This section provides the preliminaries that contextualize our approach. We
first describe the system dynamics model in Section II-B and then define
necessary notions in Section II-C to lay the groundwork for SSM theory.

###  II-A Notation

The set of integers and reals are denoted by $\mathbb{Z}$ and $\mathbb{R}$,
with their non-negative counterparts denoted by $\mathbb{Z}_{+}$ and
$\mathbb{R}_{+}$. The complex numbers are denoted by $\mathbb{C}$.
$S^{1}=\mathbb{R}^{1}/\left(2\pi\mathbb{Z}^{1}\right)$ represents the circle
on the real line. $C^{k}$ represents the space of $k$ continuously-
differentiable functions and $C^{a}$ represents the space of analytic
functions. $L^{2}(V,W)$ is the space of square integrable functions from a
complete vector space $V$ to $W$. $\mathcal{O}(\cdot)$ represents the standard
big-O notation. $\otimes$ is the tensor product, where $z^{\otimes 3}=z\otimes
z\otimes z$.

###  II-B System Model

####  II-B1 Full Order Model

Consider the following continuous-time, control-affine, nonlinear dynamics
with equilibrium point at the origin

$\displaystyle\begin{cases}\begin{aligned}
\dot{x}(t)&=Ax(t)+f_{0}(x(t))+\sum_{i=1}^{m}f_{i}(x(t))\,u_{i}(t),\\\
y(t)&=Hx(t),\end{aligned}\end{cases}$ (1)

where the state $x(t)\in\mathbb{R}^{N}$ is high-dimensional, i.e. $N$ is
large; $A\in\mathbb{R}^{N\times N}$ is the stability matrix; $f_{0}(x)$ are
the nonlinearities of the uncontrolled system;
$f_{i}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ for $i=1,...,m$ are nonlinear
functions that describe the state-dependence of the control effort via an
$m$-dimensional control input $u(t)\in\mathbb{R}^{m}$; the observed output of
the system is denoted as $y(t)\in\mathbb{R}^{o}$; and $H\in\mathbb{R}^{o\times
N}$ is the selection matrix of output variables, where $o\ll N$. In this work,
we assume that the performance and output variables are the same and that they
are perfectly observable. We introduce the following assumption on the form of
$A$.

###### Assumption 1

$A$ is negative definite, i.e. $A\prec 0$.

In other words, we assume that the origin $\bar{x}=0$ is a locally
asymptotically stable equilibrium point. Many physical systems and phenomena
of interest such as soft robots and fluid structure interactions satisfy this
assumption (possibly up to a shift in origin).

In addition, we introduce the following assumption on the form of $f_{i}$.

###### Assumption 2

The functions $f_{0},\dots,f_{m}\in C^{a}$. 111We make this assumption for
ease of exposition. In general, the right-hand side of System (1) can have
finite smoothness, infinite smoothness, or be analytic; correspondingly the
spectral submanifold defined in Section IV-A is as smooth as the right-hand
side.

We remark that this assumption is not particularly limiting since many
physical systems (e.g. soft robots) generically satisfy this assumption and we
are only interested in controlling smooth behavior.

###  II-C Spectral Subspace

Consider the uncontrolled part of System (1)

$\displaystyle\dot{x}(t)$ $\displaystyle=Ax(t)+f_{0}(x(t)),$ (2)

whose linearization around the origin is given by

$\displaystyle\dot{x}(t)$ $\displaystyle=Ax(t).$ (3)

For any eigenvalue $\lambda_{j}$ of $A$, there exists an eigenspace
$\mathbb{E}_{j}\subset\mathbb{R}^{N}$ spanned by the (generalized)
eigenvectors of $A$. These eigenspaces are invariant subspaces of the
linearized system (3).

###### Definition 3

A spectral subspace $\mathbb{E}_{j_{1},...,j_{n}}$ of System (1) is defined as
the direct sum of an arbitrary collection of eigenspaces of $A$ i.e.

$\displaystyle\mathbb{E}\vcentcolon=\mathbb{E}_{j_{1},...,j_{n}}=\mathbb{E}_{j_{1}}\oplus\mathbb{E}_{j_{2}}\oplus...\oplus\mathbb{E}_{j_{n}}.$

By linearity of System (3), any spectral subspace is an invariant subspace of
$A$. In projection-based methods, ROMs are constructed by projecting the
dynamics onto a nested hierarchy of the slowest $k$ spectral subspaces i.e.
$\mathbb{E}^{k}=\mathbb{E}_{1,...,k}$ where
$\mathbb{E}^{1}\subset\mathbb{E}^{2}\subset\mathbb{E}^{3}\subset\dots\subset\mathbb{E}^{k}$
and $k\ll N$.

However, such projections of the governing equations to spectral subspaces can
be guaranteed to work only for linear systems and do not capture the effects
of the nonlinear terms and control inputs of the FOM. To find a faithful
reduction of System (1), it is necessary to reason about how the additional
nonlinear terms and time dependent forcing influence the structure of the
spectral subspace. To this end, we propose using SSMs and their reduced
dynamics for reducing the following nonlinear control problem.

##  III PROBLEM STATEMENT

In this section we provide a formal problem definition of the full-order,
periodic optimal control problem in Section III-A.

###  III-A Periodic Optimal Control Problem

In this work we design periodic orbits, minimizing the mean distance to some
desired trajectory $z^{\star}(\omega\,t)\in\mathbb{R}^{o}$, where
$\omega\in\mathbb{R}_{+}$ is the frequency of the reference trajectory. Our
approach is to formulate the following optimal control problem

$\displaystyle\min_{u(\cdot)}\quad$
$\displaystyle\frac{1}{T}\int_{0}^{T}\left\|z^{\star}(\omega\,t)-y(t)\right\|_{2}dt$
subj. to System (1) (4) $\displaystyle\quad x(0)=x(T).$

In Equation (III-A), we minimize over a class of periodic feedback control
laws of the form

$\displaystyle u(t)=\kappa(y(t),\omega t),$ (5)

where $\kappa\in L^{2}(\mathbb{R}^{o}\times S^{1},\mathbb{R}^{m})$. For ease
of notation, throughout the rest of the paper we denote $\varphi=\omega t$.

Informally, we minimize the mean-squared trajectory error between our system’s
periodic orbit and the desired trajectory, after its fast dynamics have
sufficiently decayed. We emphasize that in this work, we are interested in
synthesizing control laws that neglect transients and control for a periodic
orbit.

We remark that while we consider the case of periodic control laws in this
paper, our approach generalizes to the quasi-periodic setting.

##  IV METHODOLOGY

###  IV-A Spectral Submanifold Preliminaries

An SSM serves as the unique nonlinear continuation of a nonresonant spectral
subspace $\mathbb{E}$ for the nonlinear system (2) and is defined as follows
[haller2016nonlinear].


# PCA-RF: An Efficient Parkinson’s Disease Prediction Model based on Random
Forest Classification

Ishu Gupta*
Cloud Computing Research Center
Department of Computer Science and Engineering
National Sun Yat-sen University
Kaohsiung, Taiwan
<EMAIL_ADDRESS>
&Vartika Sharma
Department of Computer Applications
National Institute of Technology
Kurukshetra, India
136119
<EMAIL_ADDRESS>&Sizman Kaur
Department of Computer Applications
National Institute of Technology
Kurukshetra, India
136119
<EMAIL_ADDRESS>& Ashutosh Kumar Singh
Department of Computer Applications
National Institute of Technology
Kurukshetra, India
136119
<EMAIL_ADDRESS>

###### Abstract

In this modern era of overpopulation disease prediction is a crucial step in
diagnosing various diseases at an early stage. With the advancement of various
machine learning algorithms, the prediction has become quite easy. However,
the complex and the selection of an optimal machine learning technique for the
given dataset greatly affects the accuracy of the model. A large amount of
datasets exists globally but there is no effective use of it due to its
unstructured format. Hence, a lot of different techniques are available to
extract something useful for the real world to implement. Therefore, accuracy
becomes a major metric in evaluating the model. In this paper, a disease
prediction approach is proposed that implements a random forest classifier on
Parkinson’s disease. We compared the accuracy of this model with the Principal
Component Analysis (PCA) applied Artificial Neural Network (ANN) model and
captured a visible difference. The model secured a significant accuracy of up
to 90%.

_Keywords_ Parkinson’s disease $\cdot$ Prediction $\cdot$ Learning $\cdot$
Feature selection $\cdot$ Random Forest Classification (RFC) $\cdot$ Principal
component analysis (PCA) $\cdot$ Pyspark $\cdot$ Accuracy $\cdot$ Training
$\cdot$ Testing $\cdot$ Artificial neural network (ANN)

##  1 Introduction

According to the fact sheets provided by the World Health Organization (WHO),
the major cause that leads to disability and death all over the world is
chronic diseases. Non-communicable diseases (NCD) kill around 40 million
people each year, which is equal to 71% of all deaths worldwide. Heart
diseases account for most of the NCD deaths, followed by cancers, respiratory
diseases, and diabetes. This led to economic output loss of 47 US trillion
dollars in the previous two decades. This loss represents globally around 75%
of gross domestic product (GDP) in 2010 [1]. A study showed that in India
around 25% of families with a member of cardiovascular disease and 50% of
families suffering from cancer experience severe disastrous expenses and 10%
and 25%, respectively, are affected due to poverty. Most of the estimations
proved that the NCDs in India account for an overall economic loss in the
range of 5-10% of GDP, which is a significant number and is thus slowing down
GDP causing a loss in development [2].

People in India have to spend more on treatment due to the lack of medical
facilities and limited access to health insurance [3, 4]. A lot of people do
not buy insurance at the early stage of their life and then repent at an older
age [5, 6]. That’s why it is very important to have an emergency fund saved
separately. The healthcare issue of these types of diseases is crucial in
other parts of the world as well [7, 8]. In America, over $10,000 per person
is spent out of the pocket annually on health care, more than any other
country in the world [9, 10]. The Partnership to Fight Chronic Disease
estimates that around 83 million people in the U.S. will have 3 or more
chronic health problems by 2030. Some behavioral factors including unhealthy
eating, lack of exercise, excessive smoking, and use of alcohol are the main
reasons for these chronic diseases [11, 12]. Therefore, it is essential to
conduct risk analysis for chronic diseases. With a great increase in the
world’s population, it would decline our life quality if we still depend on
traditional systems of healthcare service [13, 14]. Such traditional systems
of observing the real-time health condition of the person provided real quick
assistance but with the rise in the number of data sources, there has been a
flood of data in the healthcare sector [15, 16].

Due to such a big size of data records, it becomes very difficult to use the
existing traditional techniques for effective results [17, 18]. Since Big data
is a recent technology in the real world that can bring large benefits to
business organizations, it becomes really necessary that various types of
challenges and problems associated with adopting this technology are brought
into the light [19, 20]. Detection of chronic diseases on early-stage helps in
the early inception of preventive measures and on time-effective treatment at
an initial stage has always been a better help for patients [21]. Currently,
maintaining clinical data sets has become a crucial step in the medical field
[22, 23]. The patient data that holds varied features and diagnostics related
to a specific disease should be inserted in the dataset with the utmost care
to provide the best quality results and services [24]. The data entered
manually in medical data sets can contain a great amount of incomplete data
and redundant values, mining such healthcare data becomes time-consuming [25].
As it can change the output of mining, it is important to incorporate good
data preparation and data extraction prior to applying data mining algorithms.
Prediction of the disease becomes faster and more efficient if data is
accurate, consistent, and free from noises [26].

In this era of data explosion, a large amount of medical data is generated and
updated daily [27]. The healthcare record, as well as any electronic variant
of the traditional files, contains a proper identification of the patient [28,
29]. Medical data includes Electronic Health Records (EHR) which consists of
clinical reports of patients, test reports from diagnosis, a prescription from
a doctor, information from a pharmacist, information related to a person’s
health insurance, social media posts such as blogs, tweets [30, 31, 32]. While
the healthcare sector is moving towards using even more data sets into the
daily routine diagnosis of patients suffering from chronic diseases, big data
has truly become an effective tool in improving the patient’s treatment
experience. Basic machine learning techniques with the assistance of big data
technology have brought forward the concept of prediction in the medical
sector [33]. Machine Learning capabilities have changed healthcare in various
ways, improving diagnosis of treatment choices [34, 35, 36]. The predictive
analysis implies doctors focus more on services and patient care. Prediction
using classical models required data set of patients with disease and initial
disease risk models used supervised machine learning techniques for training
data [37, 38]. These models are helpful in clinical situations and are still
studied globally [39, 40]. The features affecting a particular disease were
selected through the doctor’s experience but they could not satisfy the
changes in the disease later [41, 42].

In this paper, we proposed a Principal Component Analysis based on Random
Forest (PCA-RF) model where we applied PCA on the medical data to extract the
relevant features from the data set and along with that random forest
algorithm is implemented as a classifier model. Principal Component Analysis
(PCA) is an unsupervised learning algorithm employed for dimensionality
reduction in machine learning. It is a statistical process that converts the
observations of correlated features into a set of linearly uncorrelated
features with the help of orthogonal transformation. Further, the accuracy of
the model along with specificity and sensitivity is calculated. Finally, a
comparison between the PCA-based Artificial Neural Network (ANN) and PCA-RF
model is performed.

##  2 Related Work

The medical data contents can be maintained within the electronic format
including images and may contain the patient’s identifiable details, like –
films, digital imprints, or written conclusive outlines or interpretable
findings. Medical data is generally represented in images or texts such as
ECG, prescriptions, and MRI. Various researchers have suggested different
technologies for feature selection from medical data for solving various
problems like classification, regression, and retrieval. The prediction has
various applications such as workload [43, 44, 45], security [46, 47, 48] and
much more. In [49], the projected system used delivers a deep learning
technique for efficient prediction of multiple distinctive diseases occurring
in vulnerable areas which hold the high frequency of diseases. “It
experimented with the modified estimate models over real-life medical data
gained. It used a latent factor prototype to complete the missing data. It
experimented on a regional chronic illness of cerebral infarction. It used
machine learning Decision Tree algorithm and Map Reduce algorithm for data
partitioning on structured and unstructured data taken from the
hospital.”Compared to various crucial estimate algorithms, the accuracy of
this proposed model reached 94.8%.

Various researchers have used different machine learning algorithms for
feature extraction, classification of other chronic diseases such as Alzheimer
disease [50, 51], heart disease [52, 53, 54], diabetes [55, 56]. In this
paper, researchers used Indian Diabetes data taken from the UCI Repository.
The system was implemented in Matlab. The Pima Indian Diabetes data set
comprises approximately 768 real-life instances. The dataset contains the
patient’s summary and history and the output was predicted either as positive
or negative. By analyzing the performance, it was perceived that among all the
algorithms that were used for training, the Levenberg-Marquardt Algorithm gave
the best results based on the epochs. In paper [57], a model was proposed to
prove the outcome of a deep learning technique and to diagnose the heart
disease the data set attained from the University of California, Irvine was
used. The deep learning model was analyzed on the basis of performance and
further was compared with the other four efficient machine learning algorithms
for predicting the status of disease from data containing a record of 566
patients from a two-record set taken from the UCI database. This learning
model accomplished an accuracy score of approximately 94% and an AUC score
obtaining 0.964 in comparison to other models. The performance of this model
and the non-linear machine learning algorithms was way better when compared to
linear machine learning prediction models on increasing the data set size.

In the paper, [58], a hybrid learning prediction model was implemented with
the help of missing value imputation (HPM-MI) that analyzes different
techniques using K-means clustering and then applied with the most optimal
imputation to a data set. This model was the first used amalgamation of
K-means clustering along with Multi-layer Perceptron. To validate class labels
for given data before applying classifier Kmeans clustering is used. This
proposed system had greatly improved the quality of data by the use of the
best imputation technique after the study of different eleven approaches. “The
model is evaluated on the basis of prediction and classification system and is
investigated on three criteria of medical data such as Pima Indians Diabetes,
Wisconsin Breast Cancer, and Hepatitis dataset from the UCI Repository.”With
respect to the performance matrix, the specificity, accuracy, sensitivity;
kappa statistics, and the area under ROC were evaluated for best possible
outcomes.

In paper [59], the investigation was done on the basis of the performance of
different classification models. The breast cancer dataset which had 683
instances and 10 features were used for testing, based on classification
accuracy. The breast cancer data found from the Wisconsin data set from the
UCI repository was analyzed with the intention of developing an efficient
prediction classifier for breast cancer disease with the help of various data
mining approaches. In this observation, there were three classification
methods that were implemented and comparison results showed that the
Sequential Minimal Optimization (SMO) had a greater accuracy i.e. 96.2% than
IBK and BF methods. There are a lot of other techniques used by researchers
for predicting Parkinson’s disease [60, 61, 62, 63, 64, 65]. In this paper, we
proposed a PCA-based decision tree model for the diagnosis of chronic disease
at an early stage. Then, in the end, we have compared the ANN and Random
Forest techniques for the same data set.

##  3 Proposed System

The proposed PCA-RF model uses a decision tree-based random forest algorithm
in the core. In addition to that, it also exploits PCA that combines a large
number of features into a comparatively small set of new principal components.
Further, we pass these reduced features into our classifier that predicts the
occurrence of disease. Finally, a comparative study is done using two
classifying techniques – PCA along with ANN and PCA along with Random Forest
based on various parameters. Fig. 1 represents the overall process of the
proposed technique.

Figure 1: The workflow of the proposed disease prediction scheme

###  3.1 Feature Extraction

More features can decrease the accuracy of the model since there is more data
that needs generalization. To reduce the model’s complexity and avoid over-
fitting of data, feature selection or feature extraction can be used. Feature
extraction deals with deriving information from the data set containing a
large number of features to construct a new feature sub-space. Here, the PCA
algorithm is used for feature extraction. PCA finds the direction of greatest
variance in higher dimension data set and projects it onto a new sub-space
with fewer dimensions than the original one. Mathematically, the covariance
matrix for the available data is computed from the mean-subtracted data matrix
in Eq. (1). Then, the Eigen vector-matrix V (must be unit Eigen-vectors) which
diagonalizes C is calculated in Eq. (2), where D represents the diagonal
matrix for Eigen-values of C. The Eigen-vector with the maximum Eigen-value is
the principal component of the data. The corresponding Eigen-vectors give the
components in order of their importance. The reduced set of principal
components is passed to the classifier for prediction.

$C=\frac{1}{n-1}B*B$ (1) $V^{-1}CV=D$ (2)

###  3.2 Network Training

Random Forest is a prototype consisting of an ample number of decision trees.
This model makes use of two major theories suggested the name “random ”.
First, a random specimen of trained data points is done on constructing trees.
The next random divisions of features are taken into consideration when nodes
split. During training, every tree in a random forest prototype understands
from an arbitrary sample of the data points. The conclusive predictions of
random forest are done by taking the average of the predictions of each
distinctive tree. In order for random forest to make decisions like humans, it
needs to learn things and for this learning, it needs to have a humongous
amount of information in its training set. After learning from the training
set, it is able to process the information and they can classify the given set
of data into a predefined class.

###  3.3 Model Representation

The brick of a random forest is a decision tree which is an intuitive model.
For example, a sequence of yes-no queries asked about the data gradually
leading to a predicted class is a decision tree. It is an interpretable model
because it constructs categorizations almost as humans do. This tree is
constructed by shaping the questions via splitting of nodes as in Fig. 2. The
answers from nodes lead to a high order reduction in Gini Impurity. Decision
trees form nodes comprising a higher proportion of data points from a unit
class by taking out values in the features that vividly divide the data into
different classes.

Figure 2: Predictive model’s network representation

###  3.4 Gini Impurity

The weighted average of Gini Impurity keeps on decreasing as we proceed
further to the lower levels of the tree. The Gini Impurity of a node is the
prospect that a sample preferred randomly in a node would be wrongly
categorized if it got categorized by the allocation of experiments in the
node. Eq. (3) defines the formula for finding a Gini Impurity of any node $n$.
At every node, the tree looks out for all the characteristics for finding the
value to split on that can result in the severe drop in Gini Impurity. Then it
reiterates the split process in a greedy, recursive manner until it reaches a
maximum depth or each node has only samples from one class.

$I_{G}(n)=1-\Sigma^{J}_{i=1}(p_{i})^{2}$ (3)

##  4 Performance Evaluation

The implementation of the proposed model is described in this section followed
by accuracy results.

###  4.1 Experimental Set-up and Benchmark Dataset

For implementing the model, the Pyspark platform is used with python language.
Pyspark is used for structured and semi-structured data. It can also read data
from various data sources having different file formats with the help of an
API. The dataset is attained from the deep learning repository which is from
the University of California(UCI)[66]. The patient’s ages were ranging from 33
to 87 who was suffering from Parkinson’s disease. The task is of
classification and the dataset is of multivariate characteristics. The number
of attributes is 754 and it has no missing values. The class values are 0 and
1, depicting the occurrence of disease or not. The database contains
information like age etc. Forecast Accuracy, sensitivity, specificity,
precision, and F1 Score are evaluated for the random forest classifier to
validate the performance. The evaluation of these metrics for the proposed
scheme is done using a confusion matrix.

###  4.2 Experimental Results

To evaluate the model, a 2×2 confusion matrix is used. The confusion matrix of
size 2 × 2 is the measure of correct and wrong predictions which are abridged
with count values. The model is evaluated using the following various
performance metrics that are computed from the confusion matrix where $TN$
reflects a true negative case, where $TP$ reflects truly-positive cases, $FP$
shows false-positive cases and $FN$ holds the false-negative cases.

  * •

Sensitivity: Sensitivity can be defined as the model’s ability to accurately
detect the patients who actually have the disease.

$Sensitivity=\frac{TP}{TP+FN}$ (4)

  * •

Specificity: Specificity can be defined as the model’s ability to accurately
identify the people who are healthy and who don’t have the condition.

$Specificity=\frac{TN}{TN+FP}$ (5)

  * •

Accuracy: Accuracy is delineated as to the frequency of correct result
predictions out of the total predictions from the dataset.

$Accuracy=\frac{\text{\\#Number of correct predictions}}{\text{\\#Total number
of predictions}}$ (6)

  * •

Precision: Precision tells us how precise the prototype is out of the
predicted positive and what number of those are really positive.

$Precision=\frac{TP}{TP+FP}$ (7)

  * •

F1 Score: This score is the harmonic mean between two vital stats which are
sensitivity and precision.

$F1Score=\frac{2TP}{2TP+FP+FN}$ (8)

  * •

ROC Curve: Receiver’s Operating characteristics curve is a performance metric
for evaluating a classification model’s performance. It is a probability curve
that is plotted as TPR on the vertical dimension Y-axis and FPR on the
horizontal dimension X-axis.

The performance in terms of accuracy, sensitivity, specificity, precision, and
F1 Score of PCA-RF model without PCA and with PCA is depicted in Table 1. It
can be seen that the PCA-RF model secures 89.9% and 76.7% accuracy, 70.2% and
55.6% sensitivity, 96.5% and 80.6% specificity, 70.2% and 35.1% precision,
77.7%, and 43% F1 Score without and with PCA respectively. It is observed that
the performance without PCA is better compared to with PCA for every parameter
accuracy, sensitivity, specificity, precision, and F1 Score since the
performance greatly relies on the amalgamation of the feature reduction
technique with the proposed classifier model.

Table 1: Performance metrıcs for PCA-RF model Performance Metrics | Without PCA | With PCA
---|---|---
Accuracy | 89.867 | 76.651
Sensitivity | 70.175 | 55.555
Specificity | 96.470 | 80.628
Precision | 70.175 | 35.087
F1 Score | 77.669 | 43.010

we have plotted the ROC curve for Parkinson’s disease that is depicted in
Fig.3. The figure shows the deflection of the model curve from the baseline
where the horizontal axis (x-axis) demonstrates a false positive rate and the
vertical axis (y-axis) depicts the true positive rate. It can be seen that the
true positive rate of the proposed model is high which validates its
performance.

Figure 3: ROC curve of the proposed model

###  4.3 Comparative Analysis

Random forest classifier, when compared with Artificial Neural Network, shows
a great deflection in the accuracy. Random forest works well with tabular data
and being an intuitive model, it works on decision trees. At every step, each
attribute plays a significant role in arriving at the desired result. Whereas
Artificial Neural Network is an advanced algorithm of machine learning where a
group of attributes is selected that participate in the classification process
to predict the result. ANN has a large number of layers, each layer’s output
being input to the next layer. The accuracy of both the models is demonstrated
in Table 2 when applied with and without a feature reduction technique. Table
2 shows the performance in terms of accuracy, sensitivity, specificity,
precision, F1 score, and the comparison among ANN and Random Forest with PCA
and without PCA respectively. The depicted results show a diverge deflection
when applied with and without PCA. When applied with PCA the random forest
classifier has low performance compared to ANN. Whereas, in the case of
Artificial Neural Network when PCA has applied the accuracy greatly increases.
While random forest classifier has high performance whereas ANN as
comparatively low performance in case of without PCA. It depends on the
importance of features and the dimension of the dataset. Dimensionality
reduction is required for removing the redundancy in data that reduces time
and storage. It becomes easy to visualize the data by reducing it to a low
dimension.

Table 2: Comparison among ANN and Random Forest Metrics | with PCA | without PCA
---|---|---
ANN | Random Forest | ANN | Random Forest
Accuracy | 97.354 | 76.651 | 79.470 | 89.534
Sensitivity | 95.454 | 55.555 | 64.516 | 70.175
Specificity | 97.931 | 80.628 | 83.333 | 96.470
Precision | 93.333 | 35.087 | 50.000 | 70.175
F1 Score | 94.382 | 43.010 | 56.338 | 77.669

##  5 Conclusion and Future Scope

The work elaborated in this paper is based on the evaluation of Random forest
classification on the particular dataset which is a high dimensional data
containing 754 attributes. In terms of tabular data, random forest works well
and so far the particular data we have opted for this classifier is the most
appropriate selection. The paper also deals with how accuracy deflects when we
pass the dataset to the proposed random forest model and when passed to the
artificial neural network along with PCA. There is a visible difference
between both models. Accuracy is high for the ANN model since a feature
reduction technique is used along with it. It depends on the classifier model
used along with the feature reduction technique that ensures good accuracy.
Therefore, the accuracy of the proposed model is approximately 90%. The paper
deals with the prediction of only one type of disease, so further diseases can
be predicted simultaneously for early detection. Furthermore, the accuracy of
the various models can be contrasted and the best classifier model can be
selected for that particular disease.

## References

  * [1] Kundu, M.K., Hazra, S., Pal D., Bhattacharya M.: A review on Noncommunicable Diseases (NCDs) burden, Indian Journal of Public Health, 62:4:302-4 (2018)
  * [2] P. Godha, S. Jadon, A. Patle, I. Gupta, B. Sharma, and A. K. Singh, “Flooding and Forwarding Based on Efficient Routing Protocol,” in _International Conference on Innovative Computing and Communications_ , vol. 1166. Singapore: Springer Singapore, 2021, pp. 215–223, advances in Intelligent Systems and Computing.
  * [3] Treating chronic diseases requires a lot more than hospitalization: Know the cost, The Indian Express, Nov 30 (2019)
  * [4] K. Kaur, I. Gupta, and A. K. Singh, “A Comparative Study of the Approach Provided for Preventing the Data Leakage,” _International Journal of Network Security & Its Applications_, vol. 9, no. 5, pp. 21–33, 2017.
  * [5] P. Tiwari, S. Mehta, N. Sakhuja, I. Gupta, and A. K. Singh, “Hybrid Method in Identifying the Fraud Detection in the Credit Card,” in _Evolutionary Computing and Mobile Sustainable Networks_ , vol. 53. Singapore: Springer Singapore, 2021, pp. 27–35, data Engineering and Communications Technologies.
  * [6] I. Gupta, S. Mittal, A. Tiwari, P. Agarwal, and A. K. Singh, “TIDF-DLPM: Term and Inverse Document Frequency based Data Leakage Prevention Model,” 2022.
  * [7] I. Gupta and A. K. Singh, “A Holistic View on Data Protection for Sharing, Communicating, and Computing Environments: Taxonomy and Future Directions,” 2022\.
  * [8] V. Sharma, S. Jalwa, A. R. Siddiqi, I. Gupta, and A. K. Singh, “A Lightweight Effective Randomized Caesar Cipher Algorithm for Security of Data,” in _Evolutionary Computing and Mobile Sustainable Networks_ , vol. 53. Singapore: Springer Singapore, 2021, pp. 411–419, data Engineering and Communications Technologies.
  * [9] Waters, H., Graf, M.: Chronic diseases are taxing our healthcare system and our economy, May 31 (2018)
  * [10] I. Gupta and A. K. Singh, “Dynamic Threshold based Information Leaker Identification Scheme,” _Information Processing Letters_ , vol. 147, pp. 69 – 73, 2019.
  * [11] A. K. Singh and I. Gupta, “Online Information Leaker Identification Scheme for Secure Data Sharing,” _Multimedia Tools and Applications_ , vol. 79, no. 41, pp. 31 165–31 182, November 2020.
  * [12] A. Acharya, H. Prasad, V. Kumar, I. Gupta, and A. K. Singh, “Host Platform Security and Mobile Agent Classification: A Systematic Study,” in _Computer Networks and Inventive Communication Technologies_ , vol. 58. Singapore: Springer Singapore, 2021, pp. 1001–1010, data Engineering and Communications Technologies.
  * [13] I. Gupta, T. K. Madan, S. Singh, and A. K. Singh, “HISA-SMFM: Historical and Sentiment Analysis Based Stock Market Forecasting Model,” 2022.
  * [14] A. Kesharwani, A. Nag, A. Tiwari, I. Gupta, B. Sharma, and A. K. Singh, “Real-Time Human Locator and Advance Home Security Appliances,” in _Evolutionary Computing and Mobile Sustainable Networks_ , vol. 53. Singapore: Springer Singapore, 2021, pp. 37–49, data Engineering and Communications Technologies.
  * [15] I. Gupta and A. K. Singh, “An Integrated Approach for Data Leaker Detection in Cloud Environment,” _Journal of Information Science and Engineering_ , vol. 36, pp. 993–1005, Sep. 2020.
  * [16] R. Verma, V. Gautam, C. P. Yadav, I. Gupta, and A. K. Singh, “A Survey on Data Leakage Detection and Prevention,” in _Proceedings of the International Conference on Innovative Computing & Communications (ICICC) 2020_. SSRN, Elsevier, May 2020, pp. 1–7.
  * [17] I. Gupta and A. K. Singh, “GUIM-SMD: Guilty User Identification Model using Summation Matrix-based Distribution,” _IET Information Security_ , vol. 14, pp. 773–782, November 2020.
  * [18] S. Jalwa, V. Sharma, A. R. Siddiqi, I. Gupta, and A. K. Singh, “Comprehensive and Comparative Analysis of Different Files Using CP-ABE,” in _Advances in Communication and Computational Technology_ , vol. 668. Singapore: Springer Singapore, 2021, pp. 189–198, electrical Engineering.
  * [19] I. Gupta, H. Mittal, D. Rikhari, and A. K. Singh, “MLRM: A Multiple Linear Regression based Model for Average Temperature Prediction of A Day,” 2022.
  * [20] K. Kaur, I. Gupta, and A. K. Singh, “E-Mail Protection System to Prevent Data Leakage,” _Vigyan Prakash_ , vol. 16, pp. 30–36, 2018.
  * [21] I. Gupta, P. K. Yadav, S. Pareek, S. Shakeel, and A. K. Singh, “Auxiliary Informatics System: an Advancement towards a Smart Home Environment,” 2022.
  * [22] I. Gupta and A. K. Singh, “A Framework for Malicious Agent Detection in Cloud Computing Environment,” _International Journal of Advanced Science and Technology (IJAST)_ , vol. 135, pp. 49–62, Feb 2020.
  * [23] K. Kaur, I. Gupta, and A. K. Singh, “Data Leakage Prevention: E-Mail Protection via Gateway,” _Journal of Physics: Conference Series_ , vol. 933, p. 012013, jan 2018.
  * [24] P. Godha, S. Jadon, A. Patle, I. Gupta, B. Sharma, and A. Kumar Singh, “Architecture, an Efficient Routing, Applications, and Challenges in Delay Tolerant Network,” in _2019 International Conference on Intelligent Computing and Control Systems (ICCS)_. IEEE, 2019, pp. 824–829.
  * [25] I. Gupta and A. K. Singh, “A Hybrid Technique for the Detection of Data Leakage in Cloud computing Environment,” in _Ist International Conference on Science in Hindi_ , August 2017, vigyan Prakash.
  * [26] K. N. Kaur, Divya, I. Gupta, and A. K. Singh, “Digital Image Watermarking Using (2, 2) Visual Cryptography with DWT-SVD Based Watermarking,” in _Computational Intelligence in Data Mining_ , vol. 711. Singapore: Springer Singapore, 2019, pp. 77–86, advances in Intelligent Systems and Computing.
  * [27] Jain, D., Singh, V.: Feature selection and classification systems for chronic disease prediction: A review, Egyptian Informatics Journal, 19 (3), 179-189 (2018)
  * [28] A. Acharya, H. Prasad, V. Kumar, I. Gupta, and A. K. Singh, “MACI: Malicious API Call Identifier Model to Secure the Host Platform,” in _Proceedings of the Seventh International Conference on Mathematics and Computing_. Singapore: Springer Singapore, 2022, pp. 309–320.
  * [29] K. Gupta and I. Gupta, “A Comprehensive Study on Architecture, Security issues and Challenges in Cloud Computing,” _International Journal of Scientific & Engineering Research_, vol. 7, no. 12, pp. 128–131, Dec. 2016.
  * [30] Kumar, M., N., Manjula R.: Role of Big data analytics in rural health care – a step towards svasth bharath, (IJCSIT) International Journal of Computer Science and Information Technologies, Vol. 5 (6) , 7172-7178 (2014)
  * [31] P. Agarwal, S. Mittal, A. Tiwari, I. Gupta, A. K. Singh, and B. Sharma, “Authenticating Cryptography over Network in Data,” in _2019 International Conference on Intelligent Computing and Control Systems (ICCS)_. IEEE, 2019, pp. 632–636.
  * [32] A. Nag, A. Kesharwani, B. Sharma, I. Gupta, A. Tiwari, and A. K. Singh, “Potential and Extention of Internet of Things,” in _Second International Conference on Computer Networks and Communication Technologies (ICCNCT)_ , vol. 44. Cham: Springer International Publishing, 2020, pp. 542–551.
  * [33] A. K. Singh, I. Gupta, R. Verma, V. Gautam, and C. P. Yadav, “A Survey on Data Leakage Detection and Prevention,” in _Proc. Int. Conf. Innov. Comput. Commun._ , 2020.
  * [34] I. Gupta, N. Singh, and A. Singh, “Layer-based Privacy and Security Architecture for Cloud Data Sharing,” _Journal of Communications Software and Systems (JCOMSS)_ , vol. 15, no. 2, 2019.
  * [35] I. Gupta and A. K. Singh, “SELI: Statistical Evaluation based Leaker Identification Stochastic Scheme for Secure Data Sharing,” _IET Communications_ , vol. 14, pp. 3607–3618, December 2020.
  * [36] Khushbu, P. Nishad, V. Kashyap, and I. Gupta, “A Classification and Distribution Model for Data Leakage Prevention and Detection,” _International Research Journal of Modernization in Engineering Technology and Science_ , vol. 3, no. 2, pp. 348–354, Feb. 2021.
  * [37] G. Batra, H. Singh, I. Gupta, and A. K. Singh, “Best Fit Sharing and Power Aware (BFSPA) Algorithm for VM Placement in Cloud Environment,” in _2017 3rd International Conference on Advances in Computing,Communication & Automation (ICACCA) (Fall)_. IEEE, 2017, pp. 1–4.
  * [38] U. Arora, S. Verma, I. Gupta, and A. K. Singh, “Implementing Privacy using Modified Tree and Map Technique,” in _2017 3rd International Conference on Advances in Computing,Communication & Automation (ICACCA) (Fall)_. IEEE, 2017, pp. 1–5.
  * [39] I. Gupta, “A Comparative Study of the Approach Provided for Preventing the Data Leakage,” _Other Topics Engineering Research eJournal_ , vol. 9, no. 5, September 2017.
  * [40] I. Gupta and K. Gupta, “Evaluation of Intrusion Detection Schemes in Wireless Sensor Network,” _IOSR Journal of Computer Engineering_ , vol. 18, no. 2, pp. 60–63, Mar-Apr. 2016.
  * [41] Khushbu, P. Nishad, V. Kashyap, I. Gupta, and A. K. Singh, “An Organized Study on Data Divulge Elimination and Discernment,” in _Computer Networks and Inventive Communication Technologies_. Singapore: Springer Singapore, 2021, pp. 569–578.
  * [42] I. Gupta and K. Gupta, “Review on Intrusion Detection System Architectures in WSN,” _International Journal of Scientific & Engineering Research_, vol. 7, no. 12, pp. 111–115, Dec. 2016.
  * [43] I. Gupta and A. K. Singh, “A Probabilistic Approach for Guilty Agent Detection using Bigraph after Distribution of Sample Data,” _Procedia Computer Science_ , vol. 125, pp. 662 – 668, 2018.
  * [44] I. Gupta, R. Gupta, A. K. Singh, and R. Buyya, “MLPAM: A Machine Learning and Probabilistic Analysis Based Model for Preserving Security and Privacy in Cloud Environment,” _IEEE Systems Journal_ , vol. 15, no. 3, pp. 4248–4259, 2021.
  * [45] I. Gupta and A. K. Singh, “A Probability based Model for Data Leakage Detection using Bigraph,” in _Proceedings of 7th International Conference on Communication and Network Security (ICCNS)_ , ser. ICCNS 2017. New York, NY, USA: Association for Computing Machinery (ACM), 2017, p. 1–5.
  * [46] D. Saxena, I. Gupta, J. Kumar, A. K. Singh, and X. Wen, “A Secure and Multiobjective Virtual Machine Placement Framework for Cloud Data Center,” _IEEE Systems Journal_ , pp. 1–12, 2021.
  * [47] I. Gupta and A. K. Singh, “A Confidentiality Preserving Data Leaker Detection Model for Secure Sharing of Cloud Data using Integrated Techniques,” in _2019 7th International Conference on Smart Computing Communications (ICSCC)_. Curtin University, Sarawak Malaysia: IEEE, 2019, pp. 1–5.
  * [48] K. Kaur, I. Gupta, and A. K. Singh, “A Comparative Evaluation of Data Leakage/Loss Prevention Systems (DLPS),” in _Proc. 4th International Conference Computer Science & Information Technology_, 2017, pp. 87–95.
  * [49] Vinitha, S., Sweetlin, S., Vinusha, H., Sajini, S.:Disease Prediction Using Machine Learning Over Big Data, Computer Science & Engineering: An International Journal (CSEIJ), Vol.8, No.1, February (2018)
  * [50] A. Alberdi et al., :Smart Home-Based Prediction of Multidomain Symptoms Related to Alzheimer’s Disease, IEEE Journal of Biomedical and Health Informatics, vol. 22, no. 6, pp. 1720-1731, Nov., doi: (2018)
  * [51] Kruthika, K., R., Rajeswari, Maheshappa, H., D.:Multistage classifier-based approach for Alzheimer’s disease prediction and retrieval, Informatics in Medicine Unlocked, Volume 14, pp. 34-42, ISSN 2352-9148(2019)
  * [52] Panda, D., Dash, S.R.:Predictive System: Comparison of Classification Techniques for Effective Prediction of Heart Disease, In: Satapathy S., Bhateja V., Mohanty J., Udgata S. (eds) Smart Intelligent Computing and Applications. Smart Innovation, Systems and Technologies, vol 159. Springer, Singapore (2020)
  * [53] Ching-seh, W., Mustafa, B., Bhagwat, V.:Heart Disease Prediction Using Data Mining Techniques, 7-11. (2019)
  * [54] Singh, P., Singh, S., Pandi-Jain, G., S.:Effective heart disease prediction system using data mining techniques, International Journal of Nanomedicine, 15 March Volume 2018:13(T-NANO 2014 Abstracts) Pages 121—124.(2018)
  * [55] Jerjawi, E., Samer, N., Naser, A., Samy S.:Diabetes Prediction Using Artificial Neural Network, International Journal of Advanced Science and Technology, 121:54-64 (2018)
  * [56] Marie-Sainte S., L., Aburahmah, Almohaini, Saba,.:Current Techniques for Diabetes Prediction: Review and Case Study. Applied Sciences. 9 (2019)
  * [57] Saji, S., A., Balachandran, K.:Performance analysis of training algorithms of multilayer perceptrons in diabetes prediction, 2015 International Conference on Advances in Computer Engineering and Applications, Ghaziabad, pp. 201-206 (2016)
  * [58] Purwar, A., Singh, S., K.:Hybrid prediction model with missing value imputation for medical data, Expert Systems with Applications,Volume 42, Issue 13, 5621-5631 (2015)
  * [59] Chaurasia, V., Pal, S.: Novel Approach for Breast Cancer Detection Using Data Mining Techniques, International Journal of Innovative Research in Computer and Communication Engineering (An ISO 3297: 2007 Certified Organization) Vol. 2, Issue 1, January (2014)
  * [60] Rastegar, A., D., Ho, N., Halliday, G., M.: Parkinson’s progression prediction using machine learning and serum cytokines. npj Parkinsons Dis. 5, 14 (2019)
  * [61] Grover, S., Bhartia, S., Akshama, Yadav, A., Seeja K., R.: Predicting Severity Of Parkinson’s Disease Using Deep Learning, Procedia Computer Science, Volume 132, Pages 1788-1794 (2018)
  * [62] Challa, K., N., R., Pagolu, V., S., Panda, G., Majhi, B.:An improved approach for prediction of Parkinson’s disease using machine learning techniques, 2016 International Conference on Signal Processing, Communication, Power and Embedded System (SCOPES), Paralakhemundi, pp. 1446-1451 (2016)
  * [63] Agarwal, A., Chandrayan, S., Sahu, S., S.:Prediction of Parkinson’s disease using speech signal with Extreme Learning Machine, 2016 International Conference on Electrical, Electronics, and Optimization Techniques (ICEEOT), Chennai, pp. 3776-3779 (2016)
  * [64] Gokul S., Sivachitra, M., Vijayachitra, S.:Parkinson’s disease prediction using machine learning approaches, 2013 Fifth International Conference on Advanced Computing (ICoAC), Chennai, 2013, pp. 246-252 (2013)
  * [65] Naghavi, N., Wade, E..:Prediction of Freezing of Gait in Parkinson’s Disease Using Statistical Inference and Lower–Limb Acceleration Data. IEEE Transactions on Neural Systems and Rehabilitation Engineering. 27. 947 - 955. (2019)
  * [66] Parkinson disease dataset from UCI repository : https://archive.ics.uci.edu/ml/datasets/Parkinson%27s+Disease+Classification


Master Thesis

_Author:_

_Supervisors:_

A thesis submitted in fulfilment of the requirements
for the
in the

January, 2016

###### Abstract

In the recent years, there has been a tremendous increase in the amount of
video content uploaded to social networking and video sharing websites like
Facebook and Youtube. As of result of this, the risk of children getting
exposed to adult and violent content on the web also increased. To address
this issue, an approach to automatically detect violent content in videos is
proposed in this work. Here, a novel attempt is made also to detect the
category of violence present in a video. A system which can automatically
detect violence from both Hollywood movies and videos from the web is
extremely useful not only in parental control but also for applications
related to movie ratings, video surveillance, genre classification and so on.

Here, both audio and visual features are used to detect violence. MFCC
features are used as audio cues. Blood, Motion, and SentiBank features are
used as visual cues. Binary SVM classifiers are trained on each of these
features to detect violence. Late fusion using a weighted sum of
classification scores is performed to get final classification scores for each
of the violence class target by the system. To determine optimal weights for
each of the violence classes an approach based on grid search is employed.
Publicly available datasets, mainly Violent Scene Detection (VSD), are used
for classifier training, weight calculation, and testing. The performance of
the system is evaluated on two classification tasks, Multi-Class
classification, and Binary Classification. The results obtained for Binary
Classification are better than the baseline results from MediaEval-2014.

I, , declare that this thesis titled, ‘’ and the work presented in it are my
own. I confirm that:

  * $\blacksquare$

This work was done mainly while in candidature for a masters degree at this
University.

  * $\blacksquare$

Where I have consulted the published work of others, this is always clearly
attributed.

  * $\blacksquare$

Where I have quoted from the work of others, the source is always given. With
the exception of such quotations, this thesis is entirely my own work.

  * $\blacksquare$

I have acknowledged all main sources of help.

Signature:

Date:

“Satisfaction lies in the effort, not in the attainment, full effort is full
victory.”

Mahatma Gandhi

Abstract

###### Acknowledgements.

  * •

First of all, I would like to express my gratitude to my supervisor Dr.
Christian Schulze for his support.

  * •

I would like to thank my professors Prof. Andreas Dengel for introducing me to
this topic.

  * •

Furthermore, I would like to thank my loved ones, who have supported me
through out my journey.

  * •

I will be grateful to you all, for all your care and support.

ll ANP & Adjective Noun Pair
AP Average Precision
BPM Blood Probability Map
EER Equal Error Rate
HSV Hue Saturation Value
MAP Mean Average Precision
MFCC Mel Frequency Cepstral Coefficients
MoSIFT Motion Scale Invariant Feature Transform
RBF Radial Basis Function
RGB Red Green Blue
ROC Receiver Operating Characteristic
SIFT Scale Invariant Feature Transform
STIP Space Time Interest Points
SVM Support Vector Machines
ViF Violent Flows
VSD Vioent Scene Dataset
VSO Visual Sentiment Ontology
XML EXtended Markup Language
ZCR Zero Cross Rate

This thesis is dedicated to my father, for his love, support, and
encouragement…