Chaotic behavior and vibration suppression of a ship model via NDF controller

Abstract

It is essential to investigate the nonlinear phenomena associated with rolling ships to enhance marine operations, reliability, and stability under challenging circumstances. Ship design, operating protocols, and control strategies are improved by comprehending and anticipating these patterns of behavior. A dynamic analysis of the nonlinear rolling ship model under a parametric force is presented in this work. The nonlinear derivative feedback (NDF) control is used to enhance stability and lessen unwanted vibrations, particularly in the resonance zone. To obtain the approximate solutions (AS), the multiple-scales process (MSP) develops and solves the differential formula regarding the nonlinear rolling ship to the proper order. The correctness of the AS is verified by comparing it to the calculated numerical solutions (NS) using the fourth-order Runge–Kutta process (RK-4). Solvability requirements and resonance examples are explored to acquire the modulation formulas. MATLAB programming is used to display the frequency responses and time histories of the obtained solutions graphically. Furthermore, the temporal histories of the obtained solutions are analyzed graphically both with and without control. Stability analysis and steady-state results are also investigated using resonance curves. To determine the crucial factors that result in significant fluctuations in the structure’s behavior, the nonlinear dynamics of the model are further examined using bifurcation analysis, as represented by bifurcation diagrams. Periodic as well as chaotic oscillations are studied using the Poincaré diagram, which sheds light on the long-term equilibrium of the structure and the genesis of chaotic processes. These findings facilitate a better understanding and mitigation of dynamic disruptions in flexible structures across diverse operating environments. The dynamical analysis of the rolling ship model plays a vital role in enhancing ship stability, safety, and performance. It aids in the design of anti-roll systems, ensures cargo security, and informs safe operational guidelines under varying sea conditions. Additionally, it supports simulation-based crew training and helps predict critical transitions to chaotic motion, reducing the risk of capsizing and improving overall maritime safety.

Introduction

The rolling motion of a ship refers to its oscillatory rotation about the longitudinal axis due to external disturbances such as ocean waves. This motion is inherently nonlinear, influenced by factors like hull geometry, damping, and restoring moments resulting from buoyancy and gravity. Understanding the dynamics of ship roll motion is crucial for stability analysis, control strategies, and optimizing ship design to mitigate excessive rolling and ensure safe navigation.^1–3 Recent studies extend traditional ship stability analysis by considering dynamic instabilities in waves. Numerical simulations were employed to investigate transient roll responses to sudden heeling, with a focus on wave effects, damping, and nonlinear restoring forces, in order to identify potential instability. ⁴ Moreover, ⁵ derived kinematic models to unify maneuvering and seakeeping analyses using a robotics-based matrix approach, enabling the transformation of seakeeping equations into body-fixed coordinates. Furthermore, efficient computational models for real-time ship motion simulation were presented, accommodating various ship types, engines, and sea conditions. A simulation system was developed for educational and entertainment purposes, allowing users to interactively adjust ship and environmental parameters. The models, based on Newtonian mechanics, capture six degrees of freedom, ensuring high-fidelity motion simulations. ⁶ The authors in ⁷ studied the significant impact of ship rolling on the marine gear transmission system, highlighting increased vibrations, quasi-periodic and chaotic behavior, and a 23.81% rise in response amplitudes. By incorporating rolling effects into the marine gear transmission model, the research provides critical insights for optimizing marine transmission systems and navigation strategies to enhance stability and performance. On the other hand, the work of ⁸ explored the potential of harvesting wave-induced energy on ships using linear generators. Results show that while insufficient for propulsion, the generated power (1–2 kW per ton) could supplement auxiliary systems, offering a practical marine energy solution.

In the engineering curriculum, vibration control is crucial for minimizing negative impacts. A modified integral resonant control method was employed by ⁹ to enhance the bandwidth of a weakly damped resonant structure. Subsequently, a method for reducing the controller’s directive via a selected input by choice was created, which deviated from the conventional approach of this controller. Additionally, controller settings have been theoretically established to offer the maximum detection bandwidth. The periodic solution of the nonlinear oscillator, which describes the extended Rayleigh formula, was examined by. ¹⁰ Two distinct time delays−one for location and one for proportional velocity—were employed to regulate the vibration of the Van der Pol-Duffing-Rayleigh oscillator. The authors in ¹¹ assessed the immediate standard using the method of harmonic balance and MSP. Regarding the study of the unforced and undamped Duffing oscillator, it has been demonstrated that these techniques provide sufficient precision for analyzing backbone shapes when the amplitude fluctuation is modest. To control the bifurcation occurring in a fractional predator-prey mechanism, ¹² employed delayed feeding. In, ¹³ an investigation of the dynamically stable and unstable performance of a set of linked van der pol oscillators was performed. Additionally, it is analytically stated that the oscillator’s amplitude grows if the stability conditions are not met.

To increase the stability of the one-degree-of-freedom giant magnetostrictive actuator apparatus and manage its irregular vibratory properties, such as its principal frequency response, chaotic position, and regulating phase intensity, a feedback algorithm for displacement as well as velocity involving a time lag is devised in. ¹⁴ Giant magnetostrictive materials have gained widespread applications in various fields, including electro-hydraulic servo valves, aircraft, and gasoline infusion systems operated under electrically controlled conditions.^15–20 The researchers in ²¹ developed a regressive variable differential equation for a cable-stayed beam, regulated through time-delayed responses. The outcomes showed that displacement or acceleration feedback regulation outperformed velocity input control. In, ²² a longitudinal vibration structure for an electromechanical link transmission mechanism was created. The zone of disturbance shrank as the feedback gain rose, causing the motion to change from chaotic to periodic.

In, ²³ the MSP was utilized to mitigate unwanted residual harmonics in both conventional and unconventional Duffing systems, employing a two-scale feed-forward control technique. One basic method is the linear negative displacement feedback control system, in which the regulating force acts directly against the displacement.^24,25 This controller exhibits special stability properties in time-delayed situations, particularly when controlling various types of vibration.^26,27 In, ²⁸ the transversal time-delay loop control strategy was used to mitigate sub-harmonic and super-harmonic oscillations in a suspended cable. The Galerkin procedure was employed to derive the nonlinear differential equations of the planner motion, which were then transformed into delayed differential equations. Although time delay in the loop of feedback may decrease the efficacy of velocity feedback (VFB), investigations have demonstrated that it still performs well, regardless of various contextual circumstances. ²⁹ The intrinsic nonlinearities in cable structures are addressed by the adverse cubic VFB controller. ³⁰ Liu et al. ³¹ investigated the nonlinear oscillation characteristics of micro-nano-harmonic oscillators combined with delay time VFB, applying the multi-scale technique to determine the critical parameters for the chaotic motion of the structure.

The micro-oscillator subsystem’s chaotic motion was recorded by Bassinello et al., ³² who also examined the effects associated with each influencing factor and effectively managed the chaotic oscillation of the entire apparatus. Shang and Wen ³³ examined a micro-electromechanical system resonator with delayed feedback regarding position. It was shown that using delayed control via feedback efficiently suppresses the chaos in the micromechanical resonator. The stability and responsiveness of the structure to small voltage and vibration loads were investigated in ³⁴ concerning delayed VFB. They also demonstrated how a negative time delay feedback control gain can lead to an unstable reaction, while an excessive time lag can effectively enhance the system’s stability. Siewe ³⁵ examined the influence of three distinct active controllers on a chaotic system and found that a negative VFB can effectively control the vibration of the structure.

A dynamic scheme with a multi-excitation force procedure was proposed in ³⁶ for the shearer’s permanently magnet compact driven breaking apparatus. The intricacies of this structure were examined and addressed in this investigation. Additionally, the harmful impacts of resonance occurrences were examined by the deployment of a cubic velocity feedback (CVF) controller and a VFB controller. The relationship between the VFB and CVF regulators was clarified, demonstrating that roughly 99.95% of the control influence can be attributed to the VFB controller, while around 96.7% is related to the CVF controller. When both parametric and excitation forces were applied to the structure, the VFB controller demonstrated exceptional results. In, ³⁷ a mathematical representation of a revolving beam at different speeds was studied, the nonlinear system of differential equations was analyzed under the influence of the MSP, and the approximate solutions for the behavior of the system during the conditions of resonance were assessed, including the use of a proportional-derivative (PD) controller to enable delayed control of displacement and velocity. Reference 38 reported an investigation on a coupled pitch-roll ship model with CVF control in response to parametric excitation. The authors in^39–41 examined various dynamical models, obtained their approximate solutions using the MSP, and compared these with those obtained from the numerical algorithm. Along with an examination of the system’s chaotic dynamics using bifurcation diagrams, Poincaré maps, and phase portraits, control strategies were also incorporated into these models to lessen detrimental vibrations.

A semi-active nonlinear controller utilizing fuzzy logic and a chaotic fruit fly algorithm was developed to enhance the vibration isolation of marine engines. Integrated with a two-stage isolation system, the controller significantly improved vibration reduction, especially in low-frequency resonance regions, outperforming passive methods. ⁴² Moreover, in, ⁴³ an Adaptive Neuro-Fuzzy Control System-based semi-active control system using magnetorheological dampers was proposed to reduce ship engine vibrations, achieving over 20% improvement in vibration suppression compared to passive methods, with validated performance based on real-world data.

The dynamics of a ship model under a parametric force are examined in this paper. The NDF control is used to enhance consistency and minimize undesired vibrations, particularly in the resonance zone. To obtain an AS, the MSP addresses the mathematical representation of the ship model’s nonlinear ordinary differential equations. The correctness of the AS is confirmed by comparison with the NS, which is obtained using the RK-4. Solvability criteria and resonance scenarios are analyzed to derive the modulation formulas. A comparative graphical analysis of the temporal histories under controlled and uncontrolled settings is presented, along with the frequency responses and time histories of the obtained solutions. The system’s stability and steady-state properties are examined using resonance and the amplitude-frequency curves. The nonlinear dynamics of the model are further investigated using bifurcation analysis, illustrated through bifurcation diagrams, to identify the critical elements that lead to exploratory fluctuations in the structure’s performance. The Poincaré map is applied to study periodic and chaotic vibrations, providing insight into the system’s long-term stability and the emergence of chaotic dynamics. Engineers are better able to assess the impacts of both static and dynamic forces, such as weight, winds, or movement loads, on constructions as they investigate and construct ship constructions. We provide a feedback control approach based on the NDF approach, which may considerably lower the peak values of the ship model’s amplitude-frequency characteristic curve. In addition to providing an innovative technique and viewpoint for the active vibration reduction of rotating machinery, this scheme enhances oscillations, lessening performance.

Model’s description perturbation processes

In this section, we present the mathematical modeling of ship roll motion. Previous studies on the mentioned model^44,45 have shown that its dynamics can be described by a second-order differential equation. ⁴⁶

Here, we present the equation of motion that simulates the rolling angle $\phi \left(t\right)$ of the ship, incorporating the NDF controller to enhance stability and mitigate unwanted vibrations, as follows

$\ddot{\varphi }+{\omega }_0^2\varphi +\epsilon q\varphi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t+\epsilon \left(2{\mu }_1+{\gamma }_1{\varphi }^2+{\gamma }_2{\dot{\varphi }}^2\right)\dot{\varphi }+\epsilon \left({\delta }_1+{\delta }_2{\varphi }^2\right){\varphi }^3=\epsilon G_1\dot{y},\ddot{y}+{\omega }_1^{2}y+\epsilon {\mu }_2\dot{y}=-\epsilon G_2\dot{\varphi },$

(1)

where all utilized symbols and abbreviations are defined in Table 1.

Table 1.

Description of the utilized symbols and abbreviations.

The symbols and abbreviations	Description
$\phi \left(t\right)$	Ship rolling angle
${\mu }_1$ and ${\mu }_2$	Linear damping coefficients
${\gamma }_1$ and ${\gamma }_2$	Nonlinear damping coefficients
${\omega }_0$	Natural frequency of the model
$q$ and $\Omega$	Excitation amplitude and frequency
${\delta }_1$ and ${\delta }_2$	Coefficients of the nonlinear restoring force
$y$ and ${\omega }_1$	Controller amplitude and its internal frequency
$G_1$ and $G_2$	Controller feedback gains
$\epsilon$	Small parameter
NDF	Nonlinear derivative feedback
AS	Approximate solutions
MSP	Multiple-scales process
NS	Numerical solutions
RK-4	Fourth-order Runge-Kutta
VFB	Velocity feedback
CVF	Cubic velocity feedback
PD	Proportional-derivative
PDEs	Partial differential equations
ME	Modulation equations

The closed loop feedback control with the ship dynamical model is shown in Figure 1.

Figure 1.

NDF controller with the mathematical modeling of the dynamical system.

Now, we are going to apply the MSP ⁴⁷ to obtain an approximate solution of equation (1) up to the second order. In MSP, a small parameter $\epsilon$ is presented, where $0<\epsilon \ll 1$ . The solutions $\phi$ and $y$ can be presented as follows, according to MSP

$\varphi ={\displaystyle \sum }_{k=0}^1{\epsilon }^k{\varphi }_k\left({\tau }_0,{\tau }_1\right)+O({\epsilon }^2,y={\displaystyle \sum }_{k=0}^1{\epsilon }^k{y}_k\left({\tau }_0,{\tau }_1\right)+O\left({\epsilon }^2\right).$

(2)

In this context, ${\tau }_k={\epsilon }^k\tau \left(k=0,1\right)$ represents distinct time scales, where ${\tau }_0$ corresponds to the fast time scale, while ${\tau }_1$ represents the slow one. To accurately analyze the system’s behavior across these scales, it is necessary to express derivatives with respect to $\tau$ in terms of these multiple time scales. This transformation is achieved by applying the following differential operators

$\frac{d}{d\tau }=\frac{\partial }{\partial {\tau }_0}+\epsilon \frac{\partial }{\partial {\tau }_1},\frac{d^2}{d{\tau }^2}=\frac{{\partial }^2}{\partial {\tau }_0^2}+2\epsilon \frac{{\partial }^2}{\partial {\tau }_0\partial {\tau }_1}+O\left({\epsilon }^2\right).$

(3)

These operators explicitly demonstrate that higher-order terms, specifically those of $O\left({\epsilon }^2\right)$ and beyond, are disregarded due to their negligible contribution. By substituting equations (2) and (3) into (1), we derive a series of partial differential equations (PDEs), each associated with a distinct power of $\epsilon$

(a) Order ${\epsilon }^0$

$\frac{{\partial }^2{\phi }_0}{\partial {\tau }_0^2}+{\omega }_0^2{\phi }_0=0,$

(4)

$\frac{{\partial }^2{y}_0}{\partial {\tau }_0^2}+{\omega }_1^2{y}_0=0,$

(5)

(b) Order of $\epsilon$

$\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0^2}+{\omega }_0^2{\phi }_1=G_1\frac{\partial y_0}{\partial {\tau }_0}-2\frac{{\partial }^2{\phi }_0}{\partial {\tau }_0\partial {\tau }_1}-\frac{1}{2}q{\phi }_0\left(e^{i\Omega {\tau }_0}-e^{-i\Omega {\tau }_0}\right)-\left[2{\mu }_1+{\gamma }_1{\phi }_0^2+{\gamma }_2{\left(\frac{\partial {\phi }_0}{\partial {\tau }_0}\right)}^2\right]\frac{\partial {\phi }_0}{\partial {\tau }_0}-\left({\delta }_1+{\delta }_2{\phi }_0^2\right){\phi }_0^3.$

(6)

$\frac{{\partial }^2{y}_1}{\partial {\tau }_0^2}+{\omega }_1^2{y}_1=-2\frac{{\partial }^2{y}_0}{\partial {\tau }_0\partial {\tau }_1}-{\mu }_2\frac{\partial y_0}{\partial {\tau }_0}-G_2\frac{\partial {\phi }_0}{\partial {\tau }_0}.$

(7)

Consequently, the general solution of equations (4) and (5) is derived as follows

${\phi }_0=A{e}^{i{\omega }_0{\tau }_0}+\overline{A}{e}^{-i{\omega }_0{\tau }_0},$

(8)

$y_0=B{e}^{i{\omega }_1{\tau }_0}+\overline{B}{e}^{-i{\omega }_1{\tau }_0}.$

(9)

Notably, $A$ and $B$ are unknown functions, while their complex conjugates are denoted by $\overline{A}$ and $\overline{B}$ . The mentioned functions are dependent variables on ${\tau }_0$ and ${\tau }_1$ . By substituting solutions equations (8) and (9) into equations (6) and (7), the following secular terms emerge. These terms grow unbounded over time, leading to the breakdown of the perturbation expansion and rendering the solution physically invalid

$2i{\omega }_0\frac{\partial A}{\partial {\tau }_1}+2i{\omega }_0{\mu }_{1}A+3i{\omega }_0^3{A}^2\overline{A}+3{\delta }_1{A}^2\overline{A}+10{\delta }_2{A}^3{\overline{A}}^2=0,$

(10)

$2i{\omega }_1\frac{\partial B}{\partial {\tau }_1}+i{\omega }_1{\mu }_{2}B=0.$

(11)

In conclusion, we remove the previously mentioned secular terms and solve equations (6) and (7), as below

${\phi }_1=\frac{A^3{e}^{3i{\omega }_0{\tau }_0}}{8{\omega }_0^2}\left[{\delta }_1+i{\omega }_0\left({\gamma }_2{\omega }_0^2-{\gamma }_1\right)+5{\delta }_{2}A\overline{A}\right]+\frac{{\delta }_2{A}^5{e}^{5i{\omega }_0{\tau }_0}}{24{\omega }_0^2}+\frac{i{G}_1{\omega }_{1}B{e}^{i{\omega }_1{\tau }_0}}{{\omega }_0^2-{\omega }_1^2}-\frac{q\overline{A}{e}^{i\left(\Omega -{\omega }_0\right){\tau }_0}}{2\left[{\omega }_0^2-{\left(\Omega -{\omega }_0\right)}^2\right]}-\frac{qA{e}^{i\left(\Omega +{\omega }_0\right){\tau }_0}}{2\left[{\omega }_0^2-{\left(\Omega +{\omega }_0\right)}^2\right]}+CC,$

(12)

$y_1=-\frac{i{G}_2{\omega }_{0}A{e}^{i{\omega }_0{\tau }_0}}{{\omega }_1^2-{\omega }_0^2}+CC.$

(13)

The analytical solution of the ship roll motion using the MSP plays a crucial role in understanding the system’s nonlinear dynamic behavior. This perturbation technique enables the approximation of solutions in the presence of small nonlinearities and parametric excitations, providing insight into resonance conditions, stability boundaries, and long-term behavior. Unlike purely numerical methods, the MSP yields explicit expressions that reveal how system parameters influence the roll amplitude and phase, making it a powerful tool for designing control strategies and improving the safety and performance of marine vessels under various sea conditions.

Analysis and classification of resonance

This section seeks to categorize various resonance cases, examine a particular scenario, and formulate the modulation equations (ME). Resonance conditions occur when any denominator in the final solution nears zero, resulting in singularities. Consequently, the following possibilities arise.

• Sub-harmonic resonance at $\Omega =2{\omega }_0,$

• Internal resonance case at ${\omega }_1={\omega }_0.$

The system will behave in a complicated way if any of the resonance requirements listed above are satisfied. Additionally, it’s vital to remember that even when vibrations deviate from resonance, the resulting solutions remain valid. Let us now examine the issue in more detail, focusing on both cases of external and internal resonances that occur at $\Omega \approx 2{\omega }_0$ and ${\omega }_1\approx {\omega }_0$ . To achieve that, detuning parameters ${\sigma }_j\left(j=1,2\right)$ can be introduced, which quantifies how close $\Omega$ and ${\omega }_1$ are to ${\omega }_1$ and ${\omega }_1$ , respectively. Hence, we can express the following

$\Omega =2{\omega }_0+\epsilon {\sigma }_1,{\omega }_1={\omega }_0+\epsilon {\sigma }_2.$

(14)

Eliminating secular elements can yield the solvability criterion. Consequently, the subsequent equations reflect the probability of fulfilling the requirements

$2i{\omega }_0\frac{\partial A}{\partial {\tau }_1}+2i{\omega }_0{\mu }_{1}A+3i{\omega }_0^3{A}^2\overline{A}+3{\delta }_1{A}^2\overline{A}+10{\delta }_2{A}^3{\overline{A}}^2+0.5q\overline{A}{e}^{i{\sigma }_1{\tau }_1}-i{\omega }_{1}B{G}_1{e}^{i{\sigma }_2{\tau }_1}=0,$

(15)

$2i{\omega }_1\frac{\partial B}{\partial {\tau }_1}+i{\omega }_1{\mu }_{2}B+i{\omega }_0{G}_{2}A{e}^{-i{\sigma }_2{\tau }_1}=0.$

(16)

The established solvability conditions equations (15) and (16) contain two nonlinear PDEs for $A$ and $B$ both dependent on ${\tau }_1$ . Therefore, we can represent these functions in polar form as

$A\left({\tau }_1\right)=\frac{a_1\left({\tau }_1\right)}{2}{e}^{i{\Psi }_1\left({\tau }_1\right)};B\left({\tau }_1\right)=\frac{a_2\left({\tau }_1\right)}{2}{e}^{i{\Psi }_2\left({\tau }_1\right)}.$

(17)

The first-order derivative operator may be expressed in the following form as the functions $A$ and $B$ are independent of the ${\tau }_0$ scale

$\frac{\partial A}{\partial {\tau }_0}=0,\frac{\partial A}{\partial \tau }=\epsilon \frac{\partial A}{\partial {\tau }_1},\frac{\partial B}{\partial {\tau }_0}=0,\frac{\partial B}{\partial \tau }=\epsilon \frac{\partial B}{\partial {\tau }_1}.$

(18)

Based on equation (18), the following modified phase can be used to transform the solvability conditions of PDEs equations (15) and (16) into ordinary differential equations as follows

${\theta }_1\left({\tau }_1\right)={\tau }_1{\sigma }_1-2{\Psi }_1\left({\tau }_1\right),{\theta }_2\left({\tau }_1\right)={\tau }_1{\sigma }_2+{\Psi }_2\left({\tau }_1\right)-{\Psi }_1\left({\tau }_1\right).$

(19)

By inputting equation (17)–(19) into equations (15) and (16) and distinguishing the real and imaginary parts, we derive the following system from the prior analysis

$a_1\frac{d{\theta }_1}{d\tau }=a_1{\sigma }_1-\frac{a_1^3}{8{\omega }_0}\left(6{\delta }_1+5{\delta }_2{a}_1^2\right)-\frac{1}{2{\omega }_0}\left(2{\omega }_1{G}_1{a}_2\hspace{0.17em}\sin \hspace{0.17em}{\theta }_2+q{a}_1\hspace{0.17em}\cos \hspace{0.17em}{\theta }_1\right)$

$\frac{d{a}_1}{d\tau }=\frac{1}{4{\omega }_0}\left(2{\omega }_1{G}_1{a}_2\hspace{0.17em}\cos \hspace{0.17em}{\theta }_2-q{a}_1\hspace{0.17em}\sin \hspace{0.17em}{\theta }_1\right)-{\mu }_1{a}_1-\frac{a_1^3}{8{\omega }_0}\left({\gamma }_1+3{\gamma }_2{\omega }_0\right)$

(20)

$a_2\frac{d{\theta }_2}{d\tau }=\frac{1}{2{\omega }_1}\left[{\omega }_0{G}_2{a}_1\hspace{0.17em}\sin \hspace{0.17em}{\theta }_2-{\omega }_1{a}_2\left({\sigma }_1-2{\sigma }_2+\frac{d{\theta }_1}{d\tau }\right)\right]$

$\frac{d{a}_2}{d\tau }=-\frac{1}{2{\omega }_1}\left({\mu }_2{a}_2{\omega }_1+{\omega }_0{G}_2{a}_1\hspace{0.17em}\cos \hspace{0.17em}{\theta }_2\right)$

A comprehensive explanation of the importance of constructing these equations can be provided as follows: They help capture the slow evolution of amplitude and phase in nonlinear dynamical systems, provide a simplified mathematical representation of complex vibration behaviors, and facilitate an understanding of the system’s long-term dynamics and stability characteristics. They describe the gradual changes in amplitude variations and phase shifts. The phase modulations reveal the synchronization mechanism, the coupling between different system components, and the examination of stability zones. They serve as a powerful analytical tool for understanding and predicting the intricate dynamics of controlled vibrating systems.

Steady-state approach and stability analysis

This section primarily focuses on the solutions corresponding to the steady-state condition of the dynamical system. In this context, the solutions emerge at the end of the transient processes. For this analysis, we assume $\frac{d{a}_j}{d\tau }=\frac{d{\theta }_j}{d\tau }=0;\left(j=1,2\right)$ ³⁹ . Consequently, equation (20) can be employed to derive the following algebraic equations

$a_1{\sigma }_1-\frac{a_1^3}{8{\omega }_0}\left(6{\delta }_1+5{\delta }_2{a}_1^2\right)-\frac{1}{2{\omega }_0}\left(2{\omega }_1{G}_1{a}_2\hspace{0.17em}\sin {\theta }_2+q{a}_1\hspace{0.17em}\cos {\theta }_1\right)=0$

(21)

$\frac{1}{4{\omega }_0}\left(2{\omega }_1{G}_1{a}_2\hspace{0.17em}\cos {\theta }_2-q{a}_1\hspace{0.17em}\sin {\theta }_1\right)-{\mu }_1{a}_1-\frac{a_1^3}{8{\omega }_0}\left({\gamma }_1+3{\gamma }_2{\omega }_0\right)=0$

(22)

$\sin {\theta }_2=-\frac{{\omega }_1{a}_2\left({\sigma }_1-2{\sigma }_2\right)}{{\omega }_0{G}_2{a}_1}$

(23)

$\text{os}{\theta }_2=-\frac{{\mu }_2{a}_2{\omega }_1}{{\omega }_0{G}_2{a}_1}$

(24)

Substituting equations (23) and (24) into (21) and (22), respectively, the resultant equations will be

$\frac{q}{2{\omega }_0}\mathit{cos}\hspace{0.17em}{\theta }_1={\sigma }_1-\frac{a_1^2}{8{\omega }_0}\left(6{\delta }_1+5{\delta }_2{a}_1^2\right)+\frac{{\omega }_1^2{G}_1\left({\sigma }_1-2{\sigma }_2\right)}{{\omega }_0^2\kappa G_2},$

(25)

$\frac{q}{2{\omega }_0}\mathit{sin}\hspace{0.17em}{\theta }_1=-\frac{{\omega }_1^2{G}_1{\mu }_2}{{\omega }_0^2\kappa G_2}-2{\mu }_1-\frac{{\gamma }_1+3{\gamma }_2{\omega }_0}{4{\omega }_0}{a}_1^2.$

(26)

Based on the analysis of the previous system equations (23)−(26) with the removal of the modified phases ${\theta }_1$ and ${\theta }_2$ , the following frequency response equations are obtained

$r_1{a}_1^8+r_2{a}_1^6+r_3{a}_1^4+r_4{a}_1^2+r_5=0,$

(27)

$a_2=\pm \frac{a_1}{\sqrt{\kappa }},$

(28)

where

$r_1=H_6^2,r_2=2H_5{H}_6,$

$r_3={\left(H_2+H_3\right)}^2+H_5^2+2H_6\left(H_4-\frac{{\sigma }_1}{2}\right),$

$r_4=2\left(H_1+{\mu }_1\right)\left(H_2+H_3\right)+2H_5\left(H_4-\frac{{\sigma }_1}{2}\right),$

$r_5=2{\left(H_1+{\mu }_1\right)}^2+{\left(H_4-\frac{{\sigma }_1}{2}\right)}^2-\frac{q^2}{16{\omega }_0^2},$

$\kappa ={\left(\frac{{\omega }_1{\mu }_2}{{\omega }_0{G}_2}\right)}^2+{\left[\frac{{\omega }_1\left({\sigma }_1-2{\sigma }_2\right)}{{\omega }_0{G}_2}\right]}^2,$

$H_1=\frac{{\omega }_1^2{G}_1{\mu }_2}{2{\omega }_0^2{G}_2\kappa },H_2=\frac{{\gamma }_1}{8{\omega }_0},H_3=\frac{3{\gamma }_2{\omega }_0^2}{8},H_4=\frac{{\omega }_1^2{G}_1\left({\sigma }_1-2{\sigma }_2\right)}{2{\omega }_0^2{G}_2\kappa },H_5=\frac{3{\delta }_1}{8{\omega }_0},H_6=\frac{5{\delta }_2}{16{\omega }_0}.$

(29)

The system is formulated by linearizing equation (20) through the Lyapunov first (indirect) method, utilized to evaluate the stability characteristics of these solutions ⁴⁰

$\left[\begin{array}{l}{\dot{a}}_1\\ {\dot{\theta }}_1\\ {\dot{a}}_2\\ {\dot{\theta }}_2\end{array}\right]=\left[\begin{array}{l}\frac{\partial {\dot{a}}_1}{\partial a_1}\frac{\partial {\dot{a}}_1}{\partial {\theta }_1}\frac{\partial {\dot{a}}_1}{\partial a_2}\frac{\partial {\dot{a}}_1}{\partial {\theta }_2}\\ \frac{\partial {\dot{\theta }}_1}{\partial a_1}\frac{\partial {\dot{\theta }}_1}{\partial {\theta }_1}\frac{\partial {\dot{\theta }}_1}{\partial a_2}\frac{\partial {\dot{\theta }}_1}{\partial {\theta }_2}\\ \frac{\partial {\dot{a}}_2}{\partial a_1}\frac{\partial {\dot{a}}_2}{\partial {\theta }_1}\frac{\partial {\dot{a}}_2}{\partial a_2}\frac{\partial {\dot{a}}_2}{\partial {\theta }_2}\\ \frac{\partial {\dot{\theta }}_2}{\partial a_1}\frac{\partial {\dot{\theta }}_2}{\partial {\theta }_1}\frac{\partial {\dot{\theta }}_2}{\partial a_2}\frac{\partial {\dot{\theta }}_2}{\partial {\theta }_2}\end{array}\right]\left[\begin{array}{l}{a}_1\\ {\theta }_1\\ a_2\\ {\theta }_2\end{array}\right]=\left[\begin{array}{l}{J}_{11}{J}_{12}{J}_{13}{J}_{14}\\ J_{21}{J}_{22}{J}_{23}{J}_{24}\\ J_{31}{J}_{32}{J}_{33}{J}_{34}\\ J_{41}{J}_{42}{J}_{43}{J}_{44}\end{array}\right]\left[\begin{array}{l}{a}_1\\ {\theta }_1\\ a_2\\ {\theta }_2\end{array}\right],$

(30)

where

$J_{11}=-\frac{1}{8{\omega }_0}\left(2q\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_1+3{\gamma }_1{a}_1^2+9{\gamma }_2{a}_1^2{\omega }_0^3\right)-{\mu }_1,J_{12}=-\frac{q{a}_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_1}{4{\omega }_0},$

$J_{13}=\frac{{\omega }_1{G}_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_2}{2{\omega }_0},J_{13}=-\frac{{\omega }_1{G}_1{a}_2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_2}{2{\omega }_0},$

$J_{21}=-\frac{1}{16{\omega }_0{a}_1}\left(8q\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_1+36{\omega }_0{\delta }_1{a}_1^2+50{\delta }_2{a}_1^4-16{\omega }_0{\sigma }_1\right),J_{22}=\frac{q\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_1}{2{\omega }_0},$

$J_{23}=-\frac{{\omega }_1{G}_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_2}{{\omega }_0{a}_1},J_{24}=-\frac{{\omega }_1{G}_1{a}_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_2}{{\omega }_0{a}_1},J_{31}=-\frac{{\omega }_0{G}_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_2}{2{\omega }_1},$

$J_{32}=0,J_{33}=-\frac{{\mu }_2}{2},J_{34}=\frac{{\omega }_0{G}_2{a}_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_2}{2{\omega }_1},J_{41}=\frac{J_{21}}{2}+\frac{{\omega }_0{G}_2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_2}{2{\omega }_1{a}_2},J_{42}=\frac{J_{22}}{2},J_{43}=\frac{\left(2{\sigma }_2-{\sigma }_1\right)}{2a_2}+\frac{J_{23}}{2},J_{44}=\frac{J_{24}}{2}+\frac{{\omega }_0{G}_2{a}_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_2}{2{\omega }_1{a}_2}.$

(31)

Therefore, the capacity to evaluate and resolve the characteristic equation dictates the stability of the steady-state solution. The eigenvalues of the Jacobian matrix can be expressed as follows ⁴¹

$\left|\begin{array}{l}{J}_{11}-\lambda J_{12}{J}_{13}{J}_{14}\\ J_{21}{J}_{22}-\lambda J_{23}{J}_{24}\\ J_{31}{J}_{32}{J}_{33}-\lambda J_{34}\\ J_{41}{J}_{42}{J}_{43}{J}_{44}-\lambda \end{array}\right|=0,$

(32)

or

${\lambda }^4+{\Gamma }_1{\lambda }^3+{\Gamma }_2{\lambda }^2+{\Gamma }_3\lambda +{\Gamma }_4=0,$

(33)

where ${\Gamma }_1$ negative of the sum of the diagonal elements of $\left[J\right]$ , ${\Gamma }_2$ sum of the principal minors (determinants of $2\times 2$ sub-matrices on the diagonal), ${\Gamma }_3$ negative of the sum of the determinants of $3\times 3$ sub-matrices, ${\Gamma }_4$ determinant of $\left[J\right]$ . These coefficients are equal

${\Gamma }_1=-\left(J_{11}+J_{22}+J_{33}+J_{44}\right),$

${\Gamma }_2=J_{11}{J}_{22}+J_{11}{J}_{33}+J_{11}{J}_{44}+J_{22}{J}_{33}+J_{22}{J}_{44}+J_{33}{J}_{44}$

$-\left(J_{12}{J}_{21}+J_{13}{J}_{31}+J_{14}{J}_{41}+J_{23}{J}_{32}+J_{24}{J}_{42}+J_{34}{J}_{43}\right),$

${\Gamma }_3=-\left[J_{11}\hspace{0.17em}\mathit{det}\left(M_{11}\right)+J_{22}\hspace{0.17em}\mathit{det}\left(M_{22}\right)+J_{33}\hspace{0.17em}\mathit{det}\left(M_{33}\right)+J_{44}\hspace{0.17em}\mathit{det}\left(M_{44}\right)\right],{\Gamma }_4=\mathit{det}\left(J\right).$

(34)

Here

$M_{11}=\left[\begin{array}{l}{J}_{22}{J}_{23}{J}_{24}\\ J_{32}{J}_{33}{J}_{34}\\ J_{42}{J}_{43}{J}_{44}\end{array}\right],M_{22}=\left[\begin{array}{l}{J}_{11}{J}_{13}{J}_{14}\\ J_{31}{J}_{33}{J}_{34}\\ J_{41}{J}_{43}{J}_{44}\end{array}\right],M_{33}=\left[\begin{array}{l}{J}_{11}{J}_{12}{J}_{14}\\ J_{21}{J}_{22}{J}_{24}\\ J_{41}{J}_{42}{J}_{44}\end{array}\right],M_{44}=\left[\begin{array}{l}{J}_{11}{J}_{12}{J}_{13}\\ J_{21}{J}_{22}{J}_{23}\\ J_{31}{J}_{32}{J}_{33}\end{array}\right].$

(35)

The solution is asymptotically stable if and only if all real parts of the roots of equation (33) are negative, based on the Routh–Hurwitz criterion. ⁴¹ This condition is met when the coefficients comply with

${\Gamma }_1>0,{\Gamma }_3\left({\Gamma }_1{\Gamma }_2-{\Gamma }_3\right)-{\Gamma }_4{{\Gamma }_1}^2>0,{\Gamma }_1{\Gamma }_2-{\Gamma }_3>0,{\Gamma }_4>0.$

(36)

Which undergoes numerical analysis.

Results and discussion

In order to illustrate the system’s amplitude and phase performance, this section looks at both the controlled and uncontrolled scenarios, including the resonance state. We analyze the impact of certain elements on system amplitude. We also look at the relationship between the system amplitude and the detuning parameters, ${\sigma }_j\left(j=1,2\right)$ . To ascertain the system’s stability and chaotic behavior, bifurcation analysis, and Poincaré maps are also computed. Furthermore, the behavior of the functions $\varphi \left(\tau \right),y\left(t\right)$ , and their associated phase planes at resonance with and without control may be demonstrated by numerically solving system equation (20) of the previously discussed modulation using the RK-4. Consequently, the following initial conditions and data are considered

${\mu }_1=0.1,{\mu }_2=0.1,{\omega }_0=2,{\omega }_1={\omega }_0,\Omega =2{\omega }_0,q=5,{\gamma }_1=0.1,{\gamma }_2=0.25,{\delta }_1=0.3,{\delta }_2=2,a_1\left(0\right)=0,a_2\left(0\right)=0,{\theta }_1\left(0\right)=0.01,{\theta }_2\left(0\right)=0.01.$

Time response

The RK-4 method can be applied to solve the system numerically equation (1) after incorporating the NDF control, as mentioned earlier. Figures 2(a), and (b) display, respectively, the time history of the function $\varphi \left(\tau \right)$ and the phase portrait $\varphi \dot{\varphi }$ at resonance $\Omega =2{\omega }_0$ without applying the NDF control. We observe that the solution $\varphi \left(\tau \right)$ is reached 1.192, indicating a threat to the system. The system must remain operational by controlling and minimizing these harmful vibrations. On the other hand, Figure 3(a) displays the time history of the same function $\varphi \left(\tau \right)$ at resonance with applying the NDF control without delay time ( ${\tau }^{*}=0$ ), in which the amplitude of the function $\varphi \left(\tau \right)$ equal to $1.746\times 10^{-5}$ ; therefore, the amplitude is reduced by 99.9985%, where the NDF controller time response $y\left(\tau \right)$ is illustrated in Figure 3(b). For deciding the best controller that can be used for the ship model system, a comparison between different controller types on the system amplitude $\varphi \left(\tau \right)$ is shown in Figure 3(c). Obviously, the NDF controller is more effective than the positive position feedback (PPF), the cubic velocity feedback (CVFB), and the linear velocity feedback controller (LVFB). The time response $\varphi \left(\tau \right)$ is presented in Figure 3(d) applying the NDF controller with and without time delay. From this figure, we observe that there is a slight effect on the amplitude $\varphi \left(\tau \right)$ due to inserting the delay time on the NDF controller. Noting that the equations of motion for the ship model

Figure 2.

Time response solution $\varphi \left(\tau \right)$ , and the phase plane of $\varphi \dot{\varphi }$ at resonance instance and without control ( $\Omega =2{\omega }_0$ and $G_1=G_2=0$ ).

Figure 3.

Time response solution at resonance instance ( $\Omega =2{\omega }_0$ and $G_1=5,G_2=3$ ); (a) System amplitude $\varphi \left(\tau \right)$ with NDF control without delay time, (b) NDF controller amplitude $y\left(\tau \right)$ without delay time, (c) the effect of different controller types on the system amplitude $\varphi \left(\tau \right)$ , and (d) A comparison between the system amplitude $\varphi \left(\tau \right)$ with NDF control without delay time ${\tau }^{*}=0$ (blue curve) and with delay time ${\tau }^{*}=0.001$ (red curve).

• with PPF controller is given by

$\ddot{\varphi }+{\omega }_0^2\varphi +\epsilon q\varphi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t+\epsilon \left(2{\mu }_1+{\gamma }_1{\varphi }^2+{\gamma }_2{\dot{\varphi }}^2\right)\dot{\varphi }+\epsilon \left({\delta }_1+{\delta }_2{\varphi }^2\right){\varphi }^3=\epsilon G_{1}y,\ddot{y}+{\omega }_1^{2}y+\epsilon {\mu }_2\dot{y}=\epsilon G_2\varphi ,$

(37)

• with the CVFB controller is given by the following

$\ddot{\varphi }+{\omega }_0^2\varphi +\epsilon q\varphi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t+\epsilon \left(2{\mu }_1+{\gamma }_1{\varphi }^2+{\gamma }_2{\dot{\varphi }}^2\right)\dot{\varphi }+\epsilon \left({\delta }_1+{\delta }_2{\varphi }^2\right){\varphi }^3=-\epsilon G_1{\dot{\varphi }}^3,$

(38)

• with the LVFB controller is given by the following

$\ddot{\varphi }+{\omega }_0^2\varphi +\epsilon q\varphi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t+\epsilon \left(2{\mu }_1+{\gamma }_1{\varphi }^2+{\gamma }_2{\dot{\varphi }}^2\right)\dot{\varphi }+\epsilon \left({\delta }_1+{\delta }_2{\varphi }^2\right){\varphi }^3=-\epsilon G_1\dot{\varphi },$

(39)

• with the NDF controller with delay time is given by

$\ddot{\varphi }\left(\tau \right)+{\omega }_0^2\varphi \left(\tau \right)+\epsilon q\varphi \left(\tau \right)\mathit{cos}\hspace{0.17em}\Omega \tau +\epsilon \left(2{\mu }_1+{\gamma }_1{\varphi }^2\left(\tau \right)+{\gamma }_2{\dot{\varphi }}^2\left(\tau \right)\right)\dot{\varphi }\left(\tau \right)$

$+\epsilon \left({\delta }_1+{\delta }_2{\varphi }^2\left(\tau \right)\right){\varphi }^3\left(\tau \right)=\epsilon G_1\dot{y}\left(\tau -{\tau }^{*}\right),\ddot{y}\left(\tau \right)+{\omega }_1^{2}y\left(\tau \right)+\epsilon {\mu }_2\dot{y}\left(\tau \right)=\epsilon G_2\dot{\varphi }\left(\tau -{\tau }^{*}\right).$

(40)

System parameter effect on the response solution

We will demonstrate how the system parameters $q,{\mu }_1,{\gamma }_1,{\gamma }_2,{\delta }_1,{\delta }_2,$ as well as the control gains $G_1$ and $G_2$ , alter the amplitude of $\varphi \left(\tau \right)$ . It is observed that the parametric force amplitude $q$ is proportional to the amplitude $\varphi \left(\tau \right)$ as presented in Figure 4(a). Contrarily, the solution $\varphi \left(\tau \right)$ is inversely proportional to the parameters ${\mu }_1,{\gamma }_1,{\gamma }_2,G_1,G_2,{\delta }_1,$ and ${\delta }_2$ as depicted in Figures 4(b)–(d) and 5, respectively. To gain further insight into the transient and temporal behavior, a 3D surface plot is shown in Figure 6. The plot depicts the amplitude $\varphi \left(\tau \right)$ as a function of both time $\tau$ and the control parameter $G_1,G_2$ , the parametric excitation amplitude $q,$ and the damping coefficient ${\mu }_1$ where a gradual decay in oscillation amplitude is evident as $G_1,G_2,$ and ${\mu }_1$ increase, as shown in Figure 6(a), (b), (d). This decaying behavior is also reflected in the surface’s color distribution, where the colormap represents the instantaneous value of the amplitude $\varphi \left(\tau \right)$ . As can be seen, regions with high positive and negative voltages are highlighted in red and blue, respectively, while areas near zero amplitude are shown in yellow and green, as the control parameter $G_1,G_2,$ and ${\mu }_1$ increase. This visual pattern confirms the damping role of the control-damper system and emphasizes the importance of tuning $G_1,G_2,$ and ${\mu }_1$ to maintain efficient dynamic oscillations suppression. Furthermore, system amplitude $\varphi \left(\tau \right)$ is proportional to the parametric force amplitude $q$ as presented in Figure 6(c)).

Figure 4.

System parameter effects $q,{\mu }_1,{\gamma }_1,$ and ${\gamma }_2$ on the amplitude of the system $\varphi \left(\tau \right)$ .

Figure 5.

System parameter effects $G_1,G_2,{\delta }_1,$ and ${\delta }_2$ on the amplitude of the system $\varphi \left(\tau \right)$ .

Figure 6.

Three-dimensional surface plot of the amplitude of the system $\varphi \left(\tau \right)$ as a function of time $\tau$ at resonance case with (a) The control parameter $G_1$ for $q=5$ and $G_2=2,$ (b) The control parameter $G_2$ for $q=5$ and $G_1=3,$ (c) The parametric force amplitude $q$ for $G_1=3$ and $G_2=2,$ and (d) the damping coefficient ${\mu }_1$ for $G_1=3$ and $G_2=2,$ and $q=5$ .

The resonant response curves

The equilibrium bands of fixed points will be demonstrated by the resonant curves of responses that are derived from the NS of equation (21), which are displayed in Figures 7–26. Consider that solid lines indicate stable fixed points while dashed lines indicate unstable ones. It has been shown that specific stability criteria parameters, such as control system gain $G_1$ and $G_2$ greatly impact stability analysis, the damping factor ${\mu }_1,{\mu }_2,$ the nonlinear terms’ coefficients ${\gamma }_1,{\gamma }_2,{\delta }_1,{\delta }_2,$ the internal frequency ${\omega }_0,$ and the excitation amplitude of the force $q$ . Figures 7–16 show the system amplitude versus the detuning parameter ${\sigma }_1$ as the values given for the system parameters change. Figures 7–11 demonstrate that the amplitude is inversely proportional to $G_1,G_2,{\mu }_1,{\omega }_0,$ and ${\gamma }_2$ , whereas Figures 12–16 illustrates that the amplitude increases as ${\gamma }_1,q,{\mu }_2,{\delta }_1,$ and ${\delta }_2$ grows. As the settings provided for the system factors vary, the structure magnitude versus the detuning parameter ${\sigma }_2$ is displayed in Figures 17–26. While Figures 17–22 show that the amplitude increases as $G_1,G_2,{\gamma }_2,{\mu }_1,{\mu }_2,$ and ${\omega }_0$ decrease, Figure 23 shows that the amplitude is proportional to $q$ . A slight change has occurred due to the effect of the parameters ${\gamma }_1,{\delta }_1,$ and ${\delta }_2$ , as shown in Figures 24–26, respectively. Table 2 shows the stability and the instability region for the system and the controller amplitudes at various system parameters within the domain ${\sigma }_j\in \Lambda ,\Lambda =\left[-0.5,0.5\right],j=1,2$ .

Figure 7.

Resonance curves for amplitudes $a_j\left(j=1,2\right)$ versus ${\sigma }_1$ , with NDF control at various values of $G_1;{\sigma }_2=0,G_2=0.1$ .

Figure 8.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1,$ with NDF control at various values of $G_2;{\sigma }_2=0,G_1=1$ .

Figure 9.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0,$ with NDF control at various values of ${\mu }_1$ .

Figure 10.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0$ , with NDF control at various values of ${\omega }_0$ .

Figure 11.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0,$ with NDF control at various values of ${\gamma }_2$ .

Figure 12.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0,$ with NDF control at various values of ${\gamma }_1$ .

Figure 13.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0,$ with NDF control at various values of $q$ .

Figure 14.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0$ , with NDF control at various values of ${\mu }_2$ .

Figure 15.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0$ , with NDF control at various values of ${\delta }_1$ .

Figure 16.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_1;{\sigma }_2=0$ , with NDF control at various values of ${\delta }_2$ .

Figure 17.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of $G_1;{\sigma }_1=0,G_2=0.1$ .

Figure 18.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of $G_2;{\sigma }_1=0,G_1=1$ .

Figure 19.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of ${\gamma }_2;{\sigma }_1=0$ .

Figure 20.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of ${\mu }_1;{\sigma }_1=0$ .

Figure 21.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of ${\mu }_2;{\sigma }_1=0$ .

Figure 22.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of ${\omega }_0;{\sigma }_1=0$ .

Figure 23.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of $q;{\sigma }_1=0$ .

Figure 24.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of ${\gamma }_1;{\sigma }_1=0$ .

Figure 25.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of ${\delta }_1;{\sigma }_1=0$ .

Figure 26.

Resonance curves for amplitudes $a_j$ versus ${\sigma }_2$ , with NDF control at various values of ${\delta }_2;{\sigma }_1=0$ .

Table 2.

Stability and instability regions zones for the system and the controller amplitudes.

Fig. Number	Parameter used	System amplitude $a_1$ vs. ${\sigma }_j$		Controller amplitude $a_2$ vs. ${\sigma }_j$
		Stability region	Instability region	Stability region	Instability region
Figure 7	$G_1=1$	$\Lambda -\left(-0.05,0.092\right)$	$-0.05\le {\sigma }_1<0.092$	$\Lambda -\left(-0.05,0.094\right)$	$-0.05\le {\sigma }_1<0.094$
	$G_1=0.9$	$\Lambda -\left(-0.034,0.078\right)$	$-0.034\le {\sigma }_1<0.078$	$\Lambda -\left(-0.034,0.076\right)$	$-0.034\le {\sigma }_1<0.076$
	$G_1=0.8$	$\Lambda -\left(-0.01,0.054\right)$	$-0.01\le {\sigma }_1<0.054$	$\Lambda -\left(-0.012,0.054\right)$	$-0.012\le {\sigma }_1<0.054$
Figure 8	$G_2=0.105$	$\Lambda -\left(-0.05,0.1\right)$	$-0.05\le {\sigma }_1<0.1$	$\Lambda -\left(-0.052,0.1\right)$	$-0.052\le {\sigma }_1<0.1$
	$G_2=0.1$	$\Lambda -\left(-0.05,0.094\right)$	$-0.05\le {\sigma }_1<0.094$	$\Lambda -\left(-0.051,0.094\right)$	$-0.051\le {\sigma }_1<0.094$
	$G_2=0.08$	$\Lambda -\left(-0.014,0.06\right)$	$-0.014\le {\sigma }_1<0.06$	$\Lambda -\left(-0.01,0.054\right)$	$-0.01\le {\sigma }_1<0.054$
Figure 9	${\mu }_1=0.1$	$\Lambda -\left(-0.05,0.092\right)$	$-0.05\le {\sigma }_1<0.092$	$\Lambda -\left(-0.05,0.09\right)$	$-0.05\le {\sigma }_1<0.09$
	${\mu }_1=0.07$	$\Lambda -\left(-0.042,0.078\right)$	$-0.042\le {\sigma }_1<0.078$	$\Lambda -\left(-0.044,0.086\right)$	$-0.044\le {\sigma }_1<0.086$
	${\mu }_1=0.05$	$\Lambda -\left(-0.038,0.082\right)$	$-0.038\le {\sigma }_1<0.082$	$\Lambda -\left(-0.038,0.078\right)$	$-0.038\le {\sigma }_1<0.078$
Figure 10	${\omega }_0=1.8$	$\Lambda -\left(-0.024,0.066\right)$	$-0.024\le {\sigma }_1<0.066$	$\Lambda -\left(-0.02,0.066\right)$	$-0.02\le {\sigma }_1<0.066$
	${\omega }_0=2$	$\Lambda -\left(-0.052,0.092\right)$	$-0.052\le {\sigma }_1<0.092$	$\Lambda -\left(-0.052,0.094\right)$	$-0.052\le {\sigma }_1<0.094$
	${\omega }_0=2.1$	$\Lambda -\left(-0.054,0.093\right)$	$-0.054\le {\sigma }_1<0.093$	$\Lambda -\left(-0.054,0.092\right)$	$-0.054\le {\sigma }_1<0.092$
Figure 11	${\gamma }_2=0.2$	$\Lambda -\left(-0.044,0.088\right)$	$-0.044\le {\sigma }_1<0.088$	$\Lambda -\left(-0.044,0.086\right)$	$-0.044\le {\sigma }_1<0.086$
	${\gamma }_2=0.25$	$\Lambda -\left(-0.052,0.094\right)$	$-0.052\le {\sigma }_1<0.094$	$\Lambda -\left(-0.056,0.096\right)$	$-0.056\le {\sigma }_1<0.096$
	${\gamma }_2=0.3$	$\Lambda -\left(-0.054,0.095\right)$	$-0.054\le {\sigma }_1<0.095$	$\Lambda -\left(-0.042,0.084\right)$	$-0.042\le {\sigma }_1<0.084$
Figure 12	${\gamma }_1=0.05$	$\Lambda -\left(-0.046,0.09\right)$	$-0.046\le {\sigma }_1<0.09$	$\Lambda -\left(-0.045,0.088\right)$	$-0.045\le {\sigma }_1<0.088$
	${\gamma }_1=0.1$	$\Lambda -\left(-0.05,0.09\right)$	$-0.05\le {\sigma }_1<0.09$	$\Lambda -\left(-0.048,0.09\right)$	$-0.048\le {\sigma }_1<0.09$
	${\gamma }_1=0.15$	$\Lambda -\left(-0.052,0.091\right)$	$-0.052\le {\sigma }_1<0.091$	$\Lambda -\left(-0.042,0.084\right)$	$-0.042\le {\sigma }_1<0.084$
Figure 13	$q=4.8$	$\Lambda -\left(-0.058,0.098\right)$	$-0.058\le {\sigma }_1<0.098$	$\Lambda -\left(-0.058,0.1\right)$	$-0.058\le {\sigma }_1<0.1$
	$q=5$	$\Lambda -\left(-0.052,0.088\right)$	$-0.052\le {\sigma }_1<0.088$	$\Lambda -\left(-0.05,0.094\right)$	$-0.05\le {\sigma }_1<0.094$
	$q=5.2$	$\Lambda -\left(-0.044,0.086\right)$	$-0.044\le {\sigma }_1<0.086$	$\Lambda -\left(-0.044,0.086\right)$	$-0.044\le {\sigma }_1<0.086$
Figure 14	${\mu }_2=0.1$	$\Lambda -\left(-0.05,0.092\right)$	$-0.05\le {\sigma }_1<0.092$	$\Lambda -\left(-0.052,0.92\right)$	$-0.052\le {\sigma }_1<0.92$
	${\mu }_2=0.105$	$\Lambda -\left(-0.032,0.078\right)$	$-0.032\le {\sigma }_1<0.078$	$\Lambda -\left(-0.038,0.084\right)$	$-0.038\le {\sigma }_1<0.084$
	${\mu }_2=0.11$	$\Lambda -\left(-0.031,0.076\right)$	$-0.031\le {\sigma }_1<0.076$	$\Lambda -\left(-0.032,0.078\right)$	$-0.032\le {\sigma }_1<0.078$
Figure 15	${\delta }_1=0.25$	$\Lambda -\left(-0.05,0.092\right)$	$-0.05\le {\sigma }_1<0.092$	$\Lambda -\left(-0.044,0.88\right)$	$-0.044\le {\sigma }_1<0.88$
	${\delta }_1=0.3$	$\Lambda -\left(-0.044,0.09\right)$	$-0.044\le {\sigma }_1<0.09$	$\Lambda -\left(-0.046,0.09\right)$	$-0.046\le {\sigma }_1<0.09$
	${\delta }_1=0.35$	$\Lambda -\left(-0.042,0.088\right)$	$-0.042\le {\sigma }_1<0.088$	$\Lambda -\left(-0.052,0.092\right)$	$-0.052\le {\sigma }_1<0.092$
Figure 16	${\delta }_2=0.2$	$\Lambda -\left(-0.052,0.092\right)$	$-0.052\le {\sigma }_1<0.092$	$\Lambda -\left(-0.05,0.09\right)$	$-0.05\le {\sigma }_1<0.09$
	${\delta }_2=0.4$	$\Lambda -\left(-0.04,0.086\right)$	$-0.04\le {\sigma }_1<0.086$	$\Lambda -\left(-0.042,0.088\right)$	$-0.042\le {\sigma }_1<0.088$
	${\delta }_2=0.6$	$\Lambda -\left(-0.038,0.084\right)$	$-0.038\le {\sigma }_1<0.084$	$\Lambda -\left(-0.04,0.086\right)$	$-0.04\le {\sigma }_1<0.086$
Figure 17	$G_1=1$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$
	$G_1=0.9$	$\Lambda -\left(-0.022,0.03\right)$	$-0.022\le {\sigma }_2<0.03$	$\Lambda -\left(-0.022,0.03\right)$	$-0.022\le {\sigma }_2<0.03$
	$G_1=0.8$	$\Lambda -\left(-0.008,0.016\right)$	$-0.008\le {\sigma }_2<0.016$	$\Lambda -\left(-0.01,0.014\right)$	$-0.01\le {\sigma }_2<0.014$
Figure 18	$G_2=0.08$	$\Lambda -\left(-0.008,0.016\right)$	$-0.008\le {\sigma }_2<0.016$	$\Lambda -\left(-0.008,0.016\right)$	$-0.008\le {\sigma }_2<0.016$
	$G_2=0.1$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$
	$G_2=0.105$	$\Lambda -\left(-0.036,0.044\right)$	$-0.036\le {\sigma }_2<0.044$	$\Lambda -\left(-0.036,0.044\right)$	$-0.036\le {\sigma }_2<0.044$
Figure 19	${\gamma }_2=0.25$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	${\gamma }_2=0.35$	$\Lambda -\left(-0.038,0.044\right)$	$-0.038\le {\sigma }_2<0.044$	$\Lambda -\left(-0.038,0.044\right)$	$-0.038\le {\sigma }_2<0.044$
	${\gamma }_2=0.45$	$\Lambda -\left(-0.04,0.046\right)$	$-0.04\le {\sigma }_2<0.046$	$\Lambda -\left(-0.04,0.046\right)$	$-0.04\le {\sigma }_2<0.046$
Figure 20	${\mu }_1=0.1$	$\Lambda -\left(-0.036,0.044\right)$	$-0.036\le {\sigma }_2<0.044$	$\Lambda -\left(-0.038,0.044\right)$	$-0.038\le {\sigma }_2<0.044$
	${\mu }_1=0.07$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$	$\Lambda -\left(-0.03,0.04\right)$	$-0.03\le {\sigma }_2<0.04$
	${\mu }_1=0.05$	$\Lambda -\left(-0.03,0.038\right)$	$-0.03\le {\sigma }_2<0.038$	$\Lambda -\left(-0.03,0.038\right)$	$-0.03\le {\sigma }_2<0.038$
Figure 21	${\mu }_2=0.1$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$
	${\mu }_2=0.11$	$\Lambda -\left(-0.034,0.04\right)$	$-0.034\le {\sigma }_2<0.04$	$\Lambda -\left(-0.032,0.041\right)$	$-0.032\le {\sigma }_2<0.041$
	${\mu }_2=0.13$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.041\right)$	$-0.036\le {\sigma }_2<0.041$
Figure 22	${\omega }_0=1.8$	$\Lambda -\left(-0.026,0.032\right)$	$-0.026\le {\sigma }_2<0.032$	$\Lambda -\left(-0.026,0.032\right)$	$-0.026\le {\sigma }_2<0.032$
	${\omega }_0=2$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	${\omega }_0=2.1$	$\Lambda -\left(-0.04,0.048\right)$	$-0.04\le {\sigma }_2<0.048$	$\Lambda -\left(-0.04,0.048\right)$	$-0.04\le {\sigma }_2<0.048$
Figure 23	$q=4.8$	$\Lambda -\left(-0.04,0.046\right)$	$-0.04\le {\sigma }_2<0.046$	$\Lambda -\left(-0.04,0.046\right)$	$-0.04\le {\sigma }_2<0.046$
	$q=5$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	$q=5.2$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$	$\Lambda -\left(-0.032,0.04\right)$	$-0.032\le {\sigma }_2<0.04$
Figure 24	${\gamma }_1=0.1$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	${\gamma }_1=0.2$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	${\gamma }_1=0.3$	$\Lambda -\left(-0.034,0.042\right)$	$-0.034\le {\sigma }_2<0.042$	$\Lambda -\left(-0.034,0.042\right)$	$-0.034\le {\sigma }_2<0.042$
Figure 25	${\delta }_1=0.2$	$\Lambda -\left(-0.036,0.044\right)$	$-0.036\le {\sigma }_2<0.044$	$\Lambda -\left(-0.036,0.044\right)$	$-0.036\le {\sigma }_2<0.044$
	${\delta }_1=0.3$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	${\delta }_1=0.4$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
Figure 26	${\delta }_2=0.2$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	${\delta }_2=0.4$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$
	${\delta }_2=0.6$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$	$\Lambda -\left(-0.036,0.042\right)$	$-0.036\le {\sigma }_2<0.042$

System amplitude–frequency spectrum

The frequency response curve of a dynamical system demonstrates how the amplitude of the system’s response fluctuates with the excitation range’s frequency. Figure 27 displays the amplitude-frequency behavior of the uncontrolled system under different ${\mu }_1,q,{\gamma }_1,{\gamma }_2,{\delta }_1,$ and ${\delta }_2$ values. We observe that the amplitude is proportional to $q$ , and inversely proportional to ${\mu }_1,{\gamma }_1,{\gamma }_2,$ as shown in Figures 27(a)–(d), respectively. While the effect of ${\delta }_1$ and ${\delta }_2$ on the amplitude is very small, as shown in Figures 27(e) and (f). Noted that the maximum amplitude for the system occurs at the resonance instance. These results are supported by the previous results on the frequency response given in Fig. (26). A comparison between the frequency spectrum curves without and with applying NDF control is demonstrated in Figure (28) to enhance the effectiveness of the controller that used for such a system, the results of this figure is enhanced by the time response illustrated in Figures 2 and 3.

Figure 27.

Amplitude–frequency spectrum; ${\sigma }_1={\sigma }_2=G_1=G_2=0$ (a) at various values of ${\mu }_1$ , (b) at various values of $q$ , (c) at various values of ${\gamma }_1$ , (d) at various values of ${\gamma }_2$ , (e) at various values of ${\delta }_1$ , and (f) at various values of ${\delta }_2$ .

Figure 28.

Amplitude–frequency spectrum without and with NDF control.

AS with NS confirmation

The AS driven by MSP and the NS utilizing RK-4 at the resonance case ( $\Omega =2{\omega }_0,{\omega }_1={\omega }_0$ ) without control are contrasted by the curves in Figure 29(a). However, Figure 29(b) shows the comparison of the AS and the NS for the resonance case with NDF control at $G_1=3,G_2=1$ , additionally a comparison between the NS and the AS with NDF control at $G_1=1,G_2=0.1$ is depicted in Figure 29(c). There is a great deal of agreement between the numerical and analytical solutions.

Figure 29.

The variation between AS, shown by the blue dashed line, and NS, shown by the red line at; $\Omega ={\omega }_0,{\omega }_1={\omega }_0$ (a) $G_1=G_2=0$ , (b) with NDF control and $G_1=3,G_2=1$ , and (c) with NDF control and $G_1=1,G_2=0.1$ .

Bifurcation analysis

Bifurcation analysis plays a crucial role in understanding the behavior of nonlinear dynamical systems, particularly in identifying transitions between different types of motion, including periodic, quasiperiodic, and chaotic motion. By analyzing bifurcation diagrams, phase portraits, and Poincaré maps, one can gain deep insights into a system’s stability and dynamic characteristics as parameters are varied. These tools provide a graphical representation of how a system’s behavior changes with respect to a bifurcation parameter, revealing regions of stability, periodicity, and chaos.

In this study, we investigate the bifurcation behavior of a dynamical system by using the excitation amplitude $q$ as the bifurcation parameter. Specifically, we simulate the system’s behavior under two different conditions: with control and without control. This enables us t investigate how the presence of a controller affects the system’s dynamic properties.

Bifurcation analysis without control

In Figure 30, the bifurcation diagram ⁴⁸ reveals two distinct types of motion. In the region $0<q\le 0.7$ , the system exhibits periodic motion. This periodic behavior is further confirmed by simulating the phase portraits and Poincaré maps shown in Figures (a) and (b). The phase portraits clearly display the trajectory of the system and are simulated with blue curves, while the Poincaré maps highlight periodic points (marked as red points), which indicate the regular, repetitive nature of the motion within this range of $q$ .

Figure 30.

Bifurcation diagram of $\phi$ versus the bifurcation parameter $q$ in the case without control at $G_1=G_2=0$ .

Figure 31.

Phase portraits and Poincaré maps illustrating two distinct types of motion: (a), (b) periodic motion at $q=0.1$ and (c), (d) chaotic motion at $q=5.0$ . Notably, the sections of the figure labeled (b) and (d) are the magnifications of the corresponding parts (a) and (c).

As $q$ increases and enters the range $0.7<q\le 5.0$ , the system undergoes a transition to chaotic behavior. This is reflected in Fig. (30), where we observe a messy, irregular distribution of periods, indicating the onset of chaos. To further validate this chaotic motion, phase portraits and Poincaré maps are presented in Figures 31(c) and (d). In these plots, the distribution of red points appears random and scattered, a hallmark of chaotic motion, confirming the transition from periodic to chaotic dynamics as the excitation amplitude increases.

Without the control mechanism, the system exhibits a small region of periodic motion and predominantly chaotic behavior. This analysis demonstrates the significant role that the bifurcation parameter $q$ plays in influencing the system’s motion, and how the system behaves more unpredictably as $q$ exceeds a certain threshold.

Bifurcation analysis with control

Next, we repeat the bifurcation analysis, but this time with a control mechanism in place. The gain parameters are set to $G_1=G_2=0.05$ , and the results are shown in Figure 32. The bifurcation diagram reveals that the stability region for periodic motion is significantly extended compared to the uncontrolled case. In this controlled scenario, the system remains stable and exhibits periodic motion in the range $0<q\le 1.2$ . As $q$ exceeds 1.2, the system transitions to chaotic behavior in the range $0<q\le 5.0$ , similar to the behavior observed in the uncontrolled case, but with a more clearly defined boundary between periodic and chaotic regions.

Figure 32.

Bifurcation diagram of ${\phi }_1$ versus the bifurcation parameter $q$ in the presence of control at $G_1=G_2=0.05$ .

Phase portraits and Poincaré maps for the controlled system are presented in Figsure 33. These figures further demonstrate the two different types of motion. In the periodic region, the phase portraits show smooth, closed trajectories, while the Poincaré maps display discrete, evenly spaced red points, indicating periodic behavior. In the chaotic region, the red points in the Poincaré maps are scattered randomly, illustrating the system’s transition to chaotic motion.

Figure 33.

Phase portraits and Poincaré maps illustrating two distinct types of motion: (a), (b) periodic motion at $q=0.1$ and (c), (d) chaotic motion at $q=5.0$ . Notably, the sections of the figure labeled (b) and (d) are the magnifications of the corresponding parts (a) and (c).

Comparison of controlled and uncontrolled systems

The comparison between the two cases, with and without control, highlights the significant impact that the controller has on the system’s dynamics. The introduction of control not only increases the stability of the system but also extends the range of periodic motion. The controlled system exhibits a larger stability region, with periodic behavior observed for a broader range of $q$ (up to $q=1.2$ ), compared to the uncontrolled system, where periodic motion was confined to $0<q\le 0.7$ . Additionally, the chaotic behavior in the controlled system is more confined, with a clearer transition from periodic to chaotic motion as $q$ increases beyond 1.2.

These findings emphasize the importance of incorporating a control mechanism to enhance the stability of nonlinear dynamical systems. By increasing the range of periodic motion and reducing the extent of chaotic behavior, the controller helps maintain the system in a more predictable and stable state, which can be crucial for practical applications where reliability and control are essential.

Studying the bifurcation analysis of the ship’s dynamical system is essential for predicting and preventing irregular and potentially dangerous roll motions. Tools such as bifurcation diagrams, phase portraits, and Poincaré maps provide a deeper understanding of how the system responds to changes in parameters, revealing transitions from stable to chaotic motion. These analyses help identify critical thresholds and nonlinear phenomena that cannot be captured by linear models, thus aiding in the development of more effective control mechanisms and enhancing the safety and performance of ships in diverse marine environments.

Conclusion

This study discusses a dynamic model of a rolling ship subjected to a parametric force. Particularly in the resonance zone, the NDF control is employed to increase stability and reduce undesired vibrations. The AS is obtained by solving this equation using the MSP. By comparing the AS with the computed NS using the RK-4, the accuracy of the AS is confirmed. The ME are obtained by analyzing resonance instances and solvability requirements. MATLAB is used to graphically show the frequency responses and time histories of the solutions. Additionally, a visual comparison of the generated solutions’ time histories in controlled and uncontrolled scenarios is provided. Resonance curves are also used to describe steady-state solutions and stability analysis. The nonlinear dynamics of the model are further investigated using bifurcation analysis, illustrated through bifurcation diagrams, to identify the crucial factors that lead to qualitative changes in the system’s behavior. The Poincaré map is used to study periodic and chaotic oscillations, providing insight into the system’s long-term stability and the emergence of chaotic dynamics. At the worst resonance condition, identified as sub-harmonic resonance $\Omega =2{\omega }_0$ , the NDF controller was found to have decreased the amplitude by 99.9985%. The integration of the proposed NDF controller with the studied model significantly enhances system stability and mitigates chaotic behavior. Analyzing and controlling ship roll dynamics enhances safety, stability, and performance by predicting hazardous behaviors and enabling proactive roll suppression.

We hope to work beyond routine parametric in future studies, such as extending the control strategy to full 6-DOF ship models, incorporating real sea-state data, and exploring adaptive or machine-learning-based controllers for enhanced robustness under uncertain operating conditions. These additions have been included in the revised manuscript to reflect better the practical relevance and research potential of our study.