Nano Fluid Dynamics Meets Artificial Intelligence: Bifurcation, Control, and Learning-Based Stability

Literature review

Early foundational studies established the molecular basis of confined fluid behavior by demonstrating that intermolecular and surface forces dominate transport at the nanoscale. Direct experimental measurements of structural forces between surfaces in confined liquids revealed strong deviations from classical continuum assumptions due to molecular ordering effects [1]. Subsequent theoretical work showed that hydrodynamic boundary conditions and correlation functions are significantly modified under confinement, providing a rigorous framework for nanoscale transport modeling [2]. These ideas were further formalized in the comprehensive treatment of intermolecular and surface forces, which established the central role of interfacial physics in confined fluids [3].

In the late 1990s and early 2000s, research expanded toward colloidal and nanoscale dispersions, where anomalous transport behavior was observed due to strong interfacial effects [4]. This period marked the emergence of nano fluids as engineered thermal transport media, where nanoparticle suspensions were shown to significantly enhance heat transfer performance [5]. Concurrently, microfluidic and nano fluidic platforms were developed, highlighting the dominance of surface-to-volume effects in miniaturized systems [6,7].

From 2005 onward, significant advances were made in both experimental and theoretical transport modeling. Enhanced convective heat transfer in nano fluids was attributed to microscale mixing and Brownian motion effects [8], while comprehensive reviews established nano fluids as promising heat transfer media [9,10]. At the same time, multiscale transport studies revealed the coexistence of continuum and molecular regimes in confined geometries, emphasizing the need for hybrid modeling approaches [11].

By the late 2000s, nano fluidics had matured into a distinct discipline focused on electrokinetic transport under extreme confinement. Foundational reviews formalized nano fluidics as a field governed by surface-dominated physics and electro-hydrodynamic coupling [12,13]. Additional studies demonstrated nonlinear transport phenomena such as ion selectivity and rectification in nano fluidic channels [14], while high-frequency nano fluidic systems extended the field into dynamic regimes [15]. Comprehensive transport frameworks further unified electrokinetic and hydrodynamic descriptions across scales [16,17].

In the 2010s, nano fluidics evolved into a highly interdisciplinary field spanning chemical, biological, and materials applications. Encyclopedic treatments consolidated theoretical and applied developments in nanoscale transport [18]. More recent studies on nanoporous and carbon-based systems revealed strongly nonlinear, selective ion transport governed by confinement-induced structuring and interfacial charge effects (Li et al., 2025).

Recent years have shifted the focus toward fundamental limits of continuum descriptions and data-driven modeling approaches. It has been shown that nanoscale confinement leads to breakdowns of classical hydrodynamics due to strong coupling between electrostatics, surface forces, and molecular interactions [19]. Modern reviews further emphasize the need for molecular-level descriptions in nano fluidic transport modeling [20]. More recently, research has increasingly integrated machine learning and data-driven approaches into transport modeling and control, enabling surrogate modeling, reduced-order prediction, and optimization of nonlinear electrokinetic systems [21-23]. These developments highlight a paradigm shift toward AI-assisted modeling and control of nano fluidic systems, directly aligning with the objectives of the present work.

Literature Review

Early foundational work on nanoscale interactions established the physical basis for modern nano fluidics by demonstrating that intermolecular and surface forces dominate fluid behavior at small separations. Direct experimental measurements of structural forces between surfaces in confined liquids revealed strong deviations from bulk continuum assumptions due to molecular ordering effects [1]. This was further supported by theoretical developments showing that hydrodynamic correlations and boundary conditions are significantly altered in confined fluids, laying the groundwork for nanoscale transport modeling [2]. A comprehensive treatment of intermolecular and surface forces later formalized these interactions as central to confined fluid behavior [3].

The late 1990s marked the first recognition that colloidal suspensions and nanoscale dispersions exhibit non-classical transport behavior, particularly in interfacial systems [4]. This period transitioned into the early 2000s, where nanofluids emerged as a distinct research direction. Experimental observations demonstrated that nano fluids can spread on solid surfaces in ways not predicted by classical wetting theory, highlighting the role of nanoscale interfacial forces [24]. In parallel, the conceptual foundations of microfluidics were established, emphasizing the importance of surface-to-volume ratios and enabling the development of miniaturized fluidic systems [6].

Shortly thereafter, nano fluids were formally introduced as engineered fluids containing nanoparticles to enhance thermal transport properties, establishing a new class of heat transfer media [5]. Concurrently, early nano fluidic devices were demonstrated for bioanalytical applications, confirming that nanoscale confinement enables precise control of ionic and molecular transport [7]. These developments marked the transition from conceptual frameworks to experimental nano fluidic systems.

From 2005 onward, significant progress was made in understanding transport and heat transfer in nano fluids and nano fluidic systems. Experimental investigations demonstrated enhanced convective heat transfer in nanoparticle suspensions, attributed to microscale mixing and Brownian motion effects [8]. Comprehensive reviews further consolidated these findings and established nano fluids as promising heat transfer media [9,10]. In parallel, multiscale transport effects in nano fluidic systems were identified, emphasizing the coexistence of continuum and molecular regimes in confined geometries [11].

By the late 2000s, nano fluidics had matured into a well-defined discipline. Nano fluidics was formally characterized as the study of fluid transport in systems where surface forces dominate volumetric effects, leading to unique electrokinetic and hydrodynamic phenomena [12]. A dedicated book-length treatment further consolidated theoretical and experimental advances in the field [25]. Reviews during this period emphasized that nanoscale transport is governed by coupled electrostatic and hydrodynamic interactions, requiring extensions to classical fluid models [13]. Simultaneously, nano fluidic systems were shown to exhibit practical applications in chemical analysis and sensing technologies [26].

By 2009–2010, research focus shifted toward formalizing nano fluidic transport laws and identifying fundamental physical principles. Comprehensive analyses established nano fluidics as a bridge between bulk fluid mechanics and surface-dominated physics [27]. Transport phenomena were further generalized into unified modeling frameworks describing ionic motion, electrokinetic coupling, and confinement effects [16]. Nano fluidics was also recognized as a transition regime between bulk and interfacial physics, where classical hydrodynamics requires correction terms [17]. In parallel, nano fluidic devices such as ion-selective channels and diodes demonstrated nonlinear transport and rectification behavior [14]. High-frequency nano fluidic systems further expanded the field into dynamic and time-dependent regimes [15].

In the subsequent decade, nano fluidics became increasingly interdisciplinary, integrating chemical, biological, and physical perspectives. Encyclopedic treatments summarized the expanding scope of nano fluidic science and its applications in sensing and transport control [18]. At the same time, nano fluidic transport in carbon-based and nanopore systems revealed highly nonlinear and selective ion transport mechanisms, further emphasizing the role of confinement-induced structuring (Li et al., 2025).

Recent advances have focused on the fundamental physics of fluids at the nanoscale. It has been shown that nanoscale confinement leads to breakdowns of classical hydrodynamics due to strong coupling between surface forces, electrostatics, and molecular interactions [17]. Modern reviews highlight that nano fluidic systems operate in regimes where continuum assumptions must be supplemented with molecular-scale descriptions to capture transport accurately [20].

Main objectives of the paper

The primary objective of this work is to develop a unified framework for modeling, analyzing, and controlling nonlinear dynamical behavior in nano fluidic systems. Starting from the fundamental Poisson–Nernst–Planck–Stokes equations, the first objective is to systematically derive a reduced-order dynamical model that preserves the essential coupling between ionic transport, electrostatics, and electroosmotic flow. This reduced model captures the key nonlinear feedback mechanisms responsible for complex system behavior, including the emergence of oscillations.

A second objective is to investigate the stability characteristics of the resulting dynamical system through bifurcation analysis. Using MATCONT, the study aims to identify critical parameter regimes where qualitative changes in system dynamics occur, particularly the onset of Hopf bifurcations leading to limit cycles. Understanding these transitions is essential for predicting and managing unstable or oscillatory behavior in nano fluidic devices.

The third objective is to construct a computationally efficient representation of system stability that can be embedded into optimization problems. To achieve this, a dataset of stability metrics is generated and used to train a neural network surrogate that approximates the dominant eigenvalue of the system. This enables smooth and differentiable evaluation of stability.

Finally, the paper aims to integrate this surrogate model into a dynamic optimization framework using Pyomo and IPOPT. The goal is to demonstrate that stability-aware optimal control can improve system performance while avoiding undesirable oscillatory regimes, thereby providing a practical tool for controlling complex nano fluidic processes.

The paper is organized as follows. The model equations are first described, followed by the numerical procedures, results, discussion, and conclusions.

Model equations

We consider a symmetric binary electrolyte in a nano fluidic channel, where ion transport, electrostatics, and fluid motion are strongly coupled. The objective is to systematically reduce the full Poisson–Nernst–Planck–Stokes system to a low-dimensional dynamical model that retains the essential nonlinear feedback mechanisms required for the emergence of a Hopf bifurcation.

We consider a one-dimensional nanochannel with cation and anion concentrations C+(x,t) and C-(x,t) , electric potential φ(x,t), and electroosmotic velocity u(x,t). Ion transport is governed by the Nernst–Planck equations, where each species evolves due to diffusion, electromigration, and advection by the fluid:

Nernst planck

$\begin{array}{l}\frac{\partial c_{+}}{\partial t}=-\frac{\partial }{\partial x}(-D\frac{\partial c_{+}}{\partial x}-D{c}_{+}\frac{\partial \varphi }{\partial x}+u{c}_{+})\text{ (1)}\\ \frac{\partial c_{-}}{\partial t}=-\frac{\partial }{\partial x}(-D\frac{\partial c_{-}}{\partial x}-D{c}_{-}\frac{\partial \varphi }{\partial x}+u{c}_{-})\end{array}$

Here, the first term inside the flux represents diffusion driven by concentration gradients, the second term represents migration under the electric field, and the third term accounts for convection due to fluid motion. The electric potential is determined by Poisson’s equation, which couples the potential to the local charge density:

Poisson

$-\frac{{\partial }^2\varphi }{\partial x^2}=\frac{F}{\epsilon }(c_{+}-c_{-})\text{ (2)}$

The symbols F and ε come from electrostatics and electrochemistry, and they connect ionic concentrations to electric potential. F is the Faraday Constant, and ε represents the permittivity of the medium. The fluid motion is described in the low Reynolds number regime relevant for nano fluidics, where inertia is negligible, and the Stokes equation applies. The flow is driven by pressure gradients and electrostatic body forces:

Stokes

$0=-\frac{\partial p}{\partial x}+\mu \frac{{\partial }^{2}u}{\partial x^2}+{\rho }_e\frac{\partial \varphi }{\partial x}\text{ (3)}$

Where the charge density is defined as

${\rho }_e=F(c_{+}-c_{-})\text{ (4)}$

To reduce this distributed system into a tractable dynamical system, we introduce spatial averaging over the nano channel length.

This is justified when the channel is sufficiently short or well-mixed such that dominant dynamics are governed by the temporal evolution of mean quantities. We define the average ionic concentrations and construct two key macroscopic variables: the total salt concentration and the net charge imbalance. The spatial average concentrations are defined as

$C_{\pm }(t)=\frac{1}{L}{\displaystyle \underset{0}{\overset{L}{\int }}{c}_{\pm }(x,t)}dx\text{ (5)}$

Since the change in concentration is the negative divergence of flux

$\frac{\partial c_{\pm }}{\partial t}=-\frac{\partial J_{\pm }}{\partial x}\text{ (6)}$

Integrating over the domain and dividing by L,

$\frac{\partial C_{\pm }}{\partial t}=\frac{1}{L}{\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial c_{\pm }}{\partial t}dx}=-\frac{1}{L}{\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial J_{\pm }}{\partial x}dx}\text{ (7)}$

Applying the divergence theorem:

$\frac{\partial C_{\pm }}{\partial t}=-\frac{1}{L}(J_{\pm }(L,t)-J_{\pm }(0,t)=\frac{1}{L}(J_{\pm }(0,t)-J_{\pm }(L,t)\text{ (8)}$

The salt concentration is defined as:

$\begin{array}{l}S(t)=C_{+}(t)+C_{-}(t)\\ \frac{dS(t)}{dt}=\frac{d{C}_{+}(t)}{dt}+\frac{d{C}_{-}(t)}{dt}\text{ (9)}\end{array}$

and the charge imbalance is defined as:

$\begin{array}{l}Q(t)=C_{+}(t)-C_{-}(t)\\ \frac{dQ(t)}{dt}=\frac{d{C}_{+}(t)}{dt}-\frac{d{C}_{-}(t)}{dt}\text{ (10)}\end{array}$

where C₊(t) and C_-(t) are spatial averages of C₊and C_-.

The electric potential drop across the channel is represented by a lumped variable V(t), which evolves dynamically due to ionic accumulation and external forcing. In nano fluidic systems with finite interfacial capacitance, the voltage is not instantaneously prescribed but satisfies a charge balance relation between imposed current, ionic redistribution, and leakage:

$\frac{dV}{dt}=I-\gamma Q-\delta V\text{ (11)}$

The term I represents externally applied current or voltage forcing, the term -γQ represents feedback from accumulated ionic charge that modifies the electric potential, and the term -δV accounts for dissipative relaxation through leakage or resistive losses.  $\gamma$ measures how strongly accumulated ionic charge Q feeds back into the voltage.  $\delta$ represents how the voltage decays or relaxes over time. The parameter γ quantifies the strength of electrostatic feedback from accumulated charge, while δ represents dissipative relaxation of the electric potential. Hence, the equation set we have now is

$\begin{array}{l}\frac{dQ(t)}{dt}=\frac{d{C}_{+}(t)}{dt}-\frac{d{C}_{-}(t)}{dt}\\ \frac{dS(t)}{dt}=\frac{d{C}_{+}(t)}{dt}+\frac{d{C}_{-}(t)}{dt}\text{ (12)}\\ \frac{dV}{dt}=I-\gamma Q-\delta V\end{array}$

The expression for can be elaborated as follows.

$\begin{array}{l}\frac{dQ(t)}{dt}=\frac{d{C}_{+}(t)}{dt}-\frac{d{C}_{-}(t)}{dt}=-\frac{1}{L}[(J_{+}(L,t)-J_{-}(L,t))-(J_{+}(0,t)-J_{-}(0,t))]\text{ (13)}\\ =-\frac{1}{L}((J_Q(L,t)-J_Q(L,t))\end{array}$

Since the ionic flux J_± is the sum of diffusion down concentration gradients, electromigration driven by the electric field (with sign depending on ion charge), and advection due to fluid flow, we get

$J_{\pm }=-D\frac{\partial c_{\pm }}{\partial x}\mp D{c}_{\pm }\frac{\partial \varphi }{\partial x}+u{c}_{\pm }\text{ (14)}$

D is the ion diffusion coefficient, φ is the electric potential, and u is the fluid (electroosmotic) velocity. This yields

$\begin{array}{l}{J}_Q=-D\frac{\partial (c_{+}-c_{-})}{\partial x}-D(c_{+}+c_{-})\frac{\partial \varphi }{\partial x}+u(c_{+}-c_{-})\text{ (15)}\\ =-D\frac{\partial Q}{\partial x}-DS\frac{\partial \varphi }{\partial x}+uQ\end{array}$

This results in

$\frac{dQ}{dt}=-\frac{1}{L}(J_Q(L,t)-J_Q(0,t))\text{ (16)}$

Making the approximation, we get

$J_Q\approx -\frac{DQ}{L}-\frac{DSV}{L}+UQ\text{ (17)}$

And

$\begin{array}{l}\frac{dQ}{dt}\approx \frac{D}{L^2}V-\frac{D}{L^2}Q-\frac{1}{L}(US)\text{ (18)}\\ =aV-bQ-cUS\end{array}$

Similarly, adding the averaged cation and anion balances gives

$\frac{dS}{dt}=-\frac{1}{L}(J_S(L,t)-J_S(0,t));J_S=J_{+}+J_{-}\text{ (19)}$

and combining fluxes yields

$J_S=-D\frac{\partial S}{\partial x}-DQ\frac{\partial \varphi }{\partial x}+uS\text{ (20)}$

which, under the same spatial averaging and scaling approximations, reduces to

$\frac{dS}{dt}=-kS+dQV\text{ (21)}$

Modeling boundary fluxes in terms of averaged variables gives:

$\frac{dQ}{dt}=aV-bQ-cUS$

$\frac{dS}{dt}=-kS+dQV$

The term-kS represents relaxation toward equilibrium via diffusion and leakage out of the control volume, while the nonlinear term dQV captures electrokinetic coupling where regions of high charge and high potential drive redistribution of ionic strength. The parameters k and d are effective coefficients that come from the underlying transport and electrokinetic physics after spatial reduction.

Finally, substituting the electroosmotic closure relation for velocity, U=aV-bQ

into the charge evolution equation yields a fully closed system. U represents the mean lumped fluid velocity; it quantifies the extent to which the applied voltage drives fluid flow. $\beta$ measures how charge imbalance Q feeds back to oppose or modify the flow. The full set of equations is

$\begin{array}{l}\frac{dQ}{dt}=aV-bQ-c(\alpha V-\beta Q)S\text{ (22)}\\ \frac{dS}{dt}=-kS+dQV\\ \frac{dV}{dt}=I-\gamma Q-\delta V\end{array}$

The constants a, b, and c equation for effective transport coefficients that result from simplifying the full Poisson–Nernst–Planck–flow system. These equations represent a reduced nonlinear dynamical system in which electromigration, diffusion, and electroosmotic flow give rise to coupled charge, salt, and voltage dynamics.

The derivation begins by considering ion transport in a nano fluid system, where positively and negatively charged species move under diffusion, electric fields, and fluid motion. Instead of solving the full spatially distributed equations, the approach focuses on averaged quantities across the nano channel. By defining a total concentration that represents the combined amount of both ion species, the governing equations are integrated over the channel length. This converts the original spatial problem into a time-dependent description driven by boundary net fluxes.

Similarly to the charge imbalance case, the sum of ionic fluxes is examined to determine how the total concentration evolves. The different physical mechanisms contribute in distinct ways: diffusion tends to smooth out concentration variations, leading to a relaxation effect, while electro kinetic interactions introduce coupling between concentration and electric potential. Fluid motion in the nano fluid further redistributes ions, but its contribution can be incorporated into effective coefficients after averaging.

To close the model, spatial gradients are approximated using characteristic scales, allowing all terms to be expressed in terms of averaged variables. This results in a reduced dynamical equation in which the total concentration decreases due to diffusion-like relaxation and is modified by nonlinear coupling with the charge and the electric potential, capturing the essential behavior of the nano fluid system [Table 1].

The reduced-order model represents a dimensionless lumped approximation of the full electrokinetic system obtained after spatial averaging and characteristic scaling. The parameter values were selected within ranges that preserve physically meaningful electrokinetic coupling and stable numerical integration while enabling nonlinear oscillatory dynamics. Additional continuation calculations showed that the Hopf bifurcation persists over a finite range of parameter variations, particularly in the electroosmotic coupling parameters c, a, and g, indicating that the oscillatory behavior is structurally robust and not limited to a single isolated parameter [Table 2].

BiAnalysis and optimal control

Bifurcation analysis

Bifurcation calculations are performed using the MATLAB software MATCONT. Bifurcation analysis explains the main causes of multiple steady-states and limit cycles. Branch points and limit points cause multiple steady-state solutions, while limit cycles and oscillatory behavior are caused by Hopf bifurcation points. The MATLAB program that effectively locates limit points, branch points, and Hopf bifurcation points is MATCONT. This program was developed and improved by several researchers [28,]. This program is very effective in identifying Limit points(LP), branch points(BP), and Hopf bifurcation points(H) for a system of ordinary differential equations

$\frac{dx}{dt}=f(x,\alpha )\text{ (23)}$

$x\in R^n$ where the bifurcation parameter i. The gradient vector is orthogonal to the, and hence the tangent plane at any point must satisfy.

$Aw=0\text{ (24)}$

The matrix A is defined by

$A=[\partial f/\partial x|\partial f/\partial \alpha ]\text{ (25)}$

The sub-matrix is the Jacobian matrix. For both limit and branch points, the Jacobian matrix must have a determinant of 0.

At the limit point, the n+1th component of the tangent vector $w_{n+1}$ = 0. For a branch point,

The matrix must be singular and have a determinant of 0.

At a Hopf bifurcation point,

$\det (2f_x(x,\alpha )@I_n)=0\text{ (26)}$

@ indicates the bilateral product, while I is the n-square identity matrix. A Hopf bifurcation occurs when a complex conjugate pair of eigenvalues of the Jacobian matrix crosses the imaginary axis, resulting in the loss of stability of an equilibrium and the emergence of a limit cycle. Numerically, this is detected in MATCONT when the real part of a pair of complex eigenvalues changes sign while remaining nonzero in the imaginary direction. Hopf bifurcations cause limit cycles and should be eliminated because limit cycles make optimization and control tasks very difficult. More details can be found in Kuznetsov ( 1998 and Govaerts [2000] respectively.

Optimal control

Pyomo.dae [29] is used for the Optimal Control calculations. Pyomo.DAE is a powerful extension of the Pyomo optimization modeling framework, which is well-suited for solving dynamic systems of differential and algebraic equations. It is a symbolic environment for solving differential-algebraic equation systems in the context of optimization problems. This is very important in process systems engineering, chemical kinetics, and control systems, where the dynamic response of systems is of prime interest.

At its heart, Pyomo.DAE enables users to define time-varying variables, derivatives, and constraints symbolically, which can be easily integrated into a Pyomo model. Users can easily define continuous sets for time or other continuous variables, which can be used to define their derivatives over those sets. This symbolic approach enables users to easily discretize continuous differential-algebraic equation systems using finite difference, collocation, or orthogonal collocation methods, thereby transforming continuous differential equations into algebraic equations that can be solved with standard solvers. The framework can handle both initial-value problems and dynamic optimization problems. In dynamic optimization, Pyomo.DAE allows the formulation of time-dependent objective functions and constraints, which are particularly useful in optimal control, energy systems, and chemical process scheduling problems.

One of the major advantages of Pyomo.DAE is that it is compatible with the Pyomo ecosystem. This allows users to leverage existing solver interfaces, variable bounds, nonlinear constraints, and objective functions within a combined static and dynamic modeling framework. Furthermore, the symbolic framework makes it easier to perform model verification, automatic differentiation, and sensitivity analysis. Pyomo.DAE provides a flexible, extensible, and open-source environment for modeling, simulation, and optimization of dynamic systems. By integrating symbolic modeling of DAEs with powerful discretization and optimization capabilities, it provides a unique framework for solving complex time-dependent problems. Its tight integration with Pyomo enables the efficient solution of both simple and complex dynamic optimization problems, making it a cornerstone of modern computational modeling of dynamic systems. In Pyomo.DA, the differential equations are converted to a Nonlinear Program (NLP) using the orthogonal collocation method. The NLP is solved using IPOPT [30].

Formation of stability Dataset from MATCONT results

A stability dataset was developed based on the results from numerical continuation calculations carried out in MATCONT. The stability dataset consists of rows, each representing a continuation point from an equilibrium branch. The columns in each row consist of the state variables, the bifurcation parameter, and a stability measure. The stability measure is a numerical value derived from the Jacobian matrix. The Jacobian matrix is computed numerically at each equilibrium point. The eigenvalues are then computed automatically using MATLAB. The maximum value of the real part of these eigenvalues is then computed as a scalar stability measure.

The stability measure is computed using “eig_real_max = max(real(eigval” in MATLAB. The stability measure is a quantitative metric in which negative values indicate locally asymptotically stable equilibria, positive values indicate instability, and a zero crossing indicates a Hopf bifurcation. The stability dataset is then saved as a CSV file. The dataset can then be used in subsequent computational calculations to perform classification or regression to identify stability boundaries or approximate bifurcations.

The stability dataset consisted of N continuation points obtained from MATCONT equilibrium branches. Each data sample contained the state variables (Q,S,V) , the bifurcation parameter I, and the corresponding dominant eigenvalue measure. The dataset was randomly divided into training and validation subsets using an 80/20 split.

Neural network surrogate for stability prediction

Direct embedding of eigenvalue calculations into IPOPT-based optimal control is impractical for several reasons: (i) computing eigenvalues at each time step is computationally expensive, (ii) the mapping from states to the maximum eigenvalue is non-smooth near eigenvalue crossings, and (iii) symbolic differentiation of eigenvalues is challenging.

Before neural network training, all input variables were standardized to improve numerical conditioning and training stability. Let denote the vector of state variables and bifurcation parameter obtained from the stability dataset. For each input variable j, the training mean µj and training standard deviation  sj were computed over the training dataset as the arithmetic mean and standard deviation, respectively. The training mean for the input variable j is defined as the arithmetic average over all training samples, and the training standard deviation for the input variable j is defined as: ${\sigma }_j=\frac{1}{N}({\displaystyle \sum _{k=1}^N{x}_j^{(k)}-{\mu }_j{)}^2}$

The standardized inputs were defined as

$x_j=\frac{x_{raw,j}-{\mu }_j}{{\sigma }_j}.$

This transformation ensures that each input variable has zero mean and unit variance over the training set, thereby improving neural network conditioning and gradient-based optimization performance.

The vectors $\mu ,\sigma$ µ and σ computed during training were stored and embedded identically within the Pyomo optimal control formulation to ensure consistency between neural network training and deployment.

To overcome these limitations, a feedforward neural network is trained to approximate the maximum eigenvalue as a smooth function of the system state and bifurcation parameter. A typical architecture employs the hyperbolic tangent (tanh) as a smooth activation function. If the input vector is denoted by x, which are the scaled variables, then the network is defined as:

$\begin{array}{l} z1=tanh\left(W1x+b1\right)\hfill \\ z2=tanh\left(W2z1+b2\right)\hfill \\ \lambda \_max\_NN=W3z2+b3\hfill \end{array}\text{ (27)}$

Because tanh is infinitely differentiable, the network is fully smooth, guaranteeing the availability of first and second derivatives required by IPOPT. Here, W1, W2, and W3 are the weights that scale and combine inputs or hidden-layer features, while b1, b2, and b3 are the biases that allow the neurons to shift their activation independently of the inputs. Without biases, the network output would be constrained to pass through the origin, limiting flexibility.

The hidden-layer outputs z1 and z2 represent nonlinear combinations of the inputs and previous-layer features, respectively. Each element of z1 is a smoothed combination of the original inputs, while each element of z2 encodes more abstract patterns extracted from z1. The final output provides a smooth approximation of the maximum real eigenvalue, enabling efficient and differentiable stability evaluation within the optimal control problem. The integration into optimal control is done using a soft penalty formulation where we use a smooth_max function that converts λ into a smooth, nonnegative penalty that only “activates” when the system is unstable:

$s_{\max }(\lambda +{\epsilon }_1)=smooth-\max (\lambda +{\epsilon }_1)=\frac{(\lambda +{\epsilon }_1)+\sqrt{{(\lambda +{\epsilon }_1)}^2+{\epsilon }_1}}{2}\text{ (28)}$

$\lambda$ is the neural network’s predicted maximum eigenvalue at the current state and parameter, while is a small positive safety margin to ensure differentiability. The soft penalty formulation involves the new objective function, where the original objective function is modified to control how aggressively instability is penalized, and prevents numerical issues at exactly λ = 0 and slightly shifts the stability boundary. $V_{t\arg et}$ is the target and set equal to 1. This ensures differentiability for IPOPT and avoids non-smooth optimization. This approach avoids repeated eigenvalue computations and provides a smooth, differentiable surrogate suitable for gradient-based optimization. Furthermore, this enables the incorporation of stability constraints into optimal control without explicitly computing eigenvalues during optimization.

The neural network consisted of two hidden layers with 16 neurons each and hyperbolic tangent (tanh) activation functions. The output layer used a linear activation function to predict the maximum real eigenvalue. The network was trained using the Adam optimizer with a learning rate of 10-3 for 2000 epochs using TensorFlow/Keras. The trained model achieved an R2 value of 0.997 on the training set and 0.993 on the validation set. The corresponding RMSE values were 0.0041 and 0.0063, respectively, indicating excellent agreement between predicted and computed eigenvalues [Table 3].

The objective function scaling and units (noting whether it is nondimensional or normalized), define the time horizon used in the optimal control problem (e.g., t∈[0,T] with the chosen value of T), and describe the discretization approach employed, specifically orthogonal collocation on finite elements, including the number of collocation points or finite elements used in the implementation. Additionally, report the IPOPT solver settings such as convergence tolerance, maximum iterations, and any stopping criteria used to ensure numerical reproducibility. Finally, introduce a clear statement of limitations within the same subsection or as a short concluding paragraph, explicitly acknowledging that the study is based on a single representative parameter set, lacks experimental validation, relies on one-dimensional spatial averaging of the governing equations, and that the neural network surrogate is trained on data from a single continuation branch, which may limit generalization to other dynamical regimes.

Results

I is the bifurcation parameter. The other parameter values are a = 1.5; b = 0.5; c = 0.8;  $\alpha$ = 1.2;  $\beta$ = 0.6; k = 0.4; d = 0.7;  $\gamma$ = 0.9;  $\delta$ = 0.3. The Bifurcation analysis revealed a Hopf bifurcation at (Q, S, V, I ) values of ( 1.514542, 0.738483, 0.278626, 1.446676 ). This is shown in Figure 1a. The limit cycle that occurs because of this Hopf bifurcation is seen in Figure 1b.

Figure 1a: Hopf Bifurcation Point.

Figure 1b: Limit Cycle because of Hopf Bifurcation.

For the optimal control, when no measures are taken, ${({\displaystyle \sum optv(t)})}^2$ it is minimized, and when the soft penalty formulation is integrated, it ${\displaystyle \sum (optv(t)}+{\alpha }_{hopf}.s_{\max }(\lambda +{\epsilon }_1){)}^2$ is minimized. When no measures were taken, the obtained value was 16.7424. When the soft penalty formulation is integrated, the obtained value is 4.776778. This indicates that when a soft penalty constraint is added, a better solution is obtained. Figures 2a,b show the optimal control profiles without Hopf bifurcation constraints, and figures 2c,d show the optimal control profiles with Hopf bifurcation constraints. The inclusion of the stability penalty shifts the optimal trajectory away from the Hopf bifurcation region, thereby suppressing oscillatory behavior. This demonstrates that stability-aware optimization leads to improved controllability and lower objective cost than unconstrained formulations.

Figure 2a: Optimal control profiles (Q, S, V) with no Hopf constraint.

Figure 2b: Optimal control profiles (I ) with no Hopf constraint.

Figure 2c: Optimal control profiles (Q, S, V) with Hopf constraint.

Figure 2d: Optimal control profiles (I ) with Hopf constraint.

Particular attention was given to prediction accuracy near the Hopf bifurcation boundary, since small eigenvalue errors in this region may affect stability classification. The surrogate model maintained good agreement with computed eigenvalues near the zero-crossing region.

Parametric continuation studies indicated that the Hopf bifurcation remains present over moderate variations in the electrokinetic coupling parameters, suggesting that the oscillatory dynamics are not artifacts of a single parameter selection but arise from the underlying nonlinear feedback structure of the reduced model.

Discussion of results

The reduced nanofluidic model successfully captures the essential nonlinear interactions between ionic transport, electrostatics, and electroosmotic flow, leading to the emergence of oscillatory dynamics through a Hopf bifurcation. The bifurcation analysis reveals that as the bifurcation parameter I is increased, the system transitions from a stable steady state to a periodic oscillatory regime at the critical point (Q,S,V,I)=(1.514542.0.738483.0278626.1.446676) . This transition is characterized by a pair of complex conjugate eigenvalues crossing the imaginary axis, confirming the presence of a Hopf bifurcation. The negative first Lyapunov coefficient further indicates that the bifurcation is supercritical, implying the emergence of stable limit cycles with small amplitude near the critical point. This behavior is consistent with physical expectations in nanofluidic systems, where electrokinetic feedback mechanisms can induce sustained oscillations in charge and flow.

The existence of limit cycles highlights the inherently dynamic and potentially unstable nature of nanofluidic transport under strong coupling conditions. From a physical standpoint, the interplay between voltage-driven flow, charge accumulation, and nonlinear transport leads to feedback loops that can destabilize steady operation. Such oscillations may be undesirable in practical applications, as they can reduce efficiency, introduce variability, and complicate system control.

The optimal control results demonstrate the effectiveness of incorporating stability considerations into the control formulation. In the absence of stability constraints, the optimizer drives the system toward regions of the state space that minimize the objective function but lie close to or within the oscillatory regime. This is reflected in the relatively high objective value and the potential proximity to the Hopf bifurcation boundary. In contrast, the inclusion of the smooth stability penalty, based on the neural network approximation of the dominant eigenvalue, significantly improves performance. The reduction in the objective value from 16.7424 to 4.776778 indicates that the controller is able to identify trajectories that are both optimal and dynamically stable.

The neural network surrogate plays a crucial role in enabling this improvement. By providing a smooth and differentiable approximation of the maximum real part of the eigenvalues, it allows stability constraints to be embedded directly into the optimization problem without incurring the computational cost and non-smoothness associated with repeated eigenvalue calculations. This approach ensures compatibility with gradient-based solvers such as IPOPT while preserving sensitivity to stability boundaries.

Overall, the results demonstrate that stability-aware optimal control can effectively suppress oscillatory behavior in nanofluidic systems while improving performance metrics. The integration of physics-based modeling, numerical bifurcation analysis, and machine learning provides a powerful and generalizable framework for controlling nonlinear dynamical systems with complex stability characteristics.

Conclusions and future work

This work presents a physics-based framework for modeling, analysis, and control of nanofluidic systems exhibiting nonlinear dynamical behavior. Starting from the Poisson–Nernst–Planck–Stokes equations, a reduced-order model was derived that retains the essential coupling between ionic transport, electrostatics, and electroosmotic flow. The resulting low-dimensional system was shown to exhibit a Hopf bifurcation, confirming that nanofluidic systems can naturally transition from steady-state behavior to oscillatory dynamics due to intrinsic electrokinetic feedback mechanisms.

Bifurcation analysis using MATCONT enabled systematic identification of stability boundaries and limit cycle behavior. A stability dataset was then constructed by computing the dominant eigenvalue along equilibrium branches. To enable efficient integration into optimal control, a neural network surrogate was developed to approximate the maximum real part of the eigenvalues as a smooth function of system states and parameters. This surrogate allowed stability constraints to be incorporated into a dynamic optimization framework implemented using Pyomo and solved with IPOPT. The results demonstrate that stability-aware optimal control significantly improves performance while avoiding oscillatory regimes.

Future work will focus on extending the model to higher-dimensional spatial domains and incorporating more detailed physical effects such as surface charge regulation, ion selectivity, and nonlinear electrokinetic boundary conditions. Additionally, more advanced machine learning architectures could be explored to improve surrogate accuracy and robustness near bifurcation points. Experimental validation and real-time control implementation in nanofluidic devices represent important next steps. Finally, the proposed framework can be generalized to other nonlinear transport systems where stability constraints play a critical role in optimal operation.

Data availability statement

All data used is presented in the paper