1 Introduction↩︎

Output regulation is an essential control systems design problem [1], [2] that aims to have a system track a class of desired signals while rejecting external disturbances [1], [3], [4]. Accordingly, the desired signals and external disturbances can be lumped together as exogenous signals that are generated by an autonomous differential equation called the exosystem. Various formulations of the output regulation problems have been investigated in the past decade, such as linear systems in [5] and nonlinear systems in [3] and [4] with or without uncertainties in the exosystem. Feedforward and feedback control are widely employed generic schematics for addressing output regulation problems. In terms of feedforward control, the results in [1] showed that the output regulation of nonlinear systems can be solved by a feedforward control synthesized from specific solvable nonlinear partial differential equations called nonlinear regulator equations. The solvability of the nonlinear output regulation using feedforward control strictly relies on the perfect knowledge of the plant and the exosystem dynamics in the absence of uncertainties.

The pivotal technique employed to handle uncertainties and to solve output regulation problems in feedback control system design is the internal model principle. In this formulation, the internal model can be interpreted as an observer of the steady-state generator providing online estimates of the steady-state input. Output regulation in linear dynamical systems was shown to be achievable by solving pole assignment problems in [5]. Since the steady-state tracking error is a linear function of the exogenous signals, the resulting linear internal model has poles placed at the poles of the original exosystem. Regarding nonlinear output regulation, [6] revealed that the steady-state tracking error in a nonlinear system is a nonlinear function of the exogenous signals. As a result, the feedforward control and the linear internal model become invalid in the presence of unknown parameters and nonlinearities arising from the controlled system plant and exosystem.

Multifarious versions of the internal models for various nonlinear system dynamics and diverse exosystems have been extensively provided over the last decades, such as the canonical linear internal model [7], [8], nonlinear internal models [9], [10], and general or generic internal models [3], [11]. The canonical linear internal model has been successfully applied in solving heterogeneous nonlinear output regulation problems for known and uncertain linear exosystems with adaptive methods [4], [8], [12], [13], and even recently was employed to address the disturbance rejection problem of Euler–Lagrange systems [14]. To handle nonlinear exosystems subject to more general exosystem non-sinusoidal signal classes in the absence of uncertainties, two different nonlinear internal models were introduced by [9] and [10]. To remove various assumptions on the steady-state input, a generic internal model for addressing output regulation of minimum and non-minimum phase nonlinear systems was initially proposed in [3]. More articles and reviews related to internal models in control theory, bioengineering, and neuroscience are cited in [15].

In terms of the generic internal model, as pointed out by [15], a significant advantage of [3] is that the generic internal model does not rely on any specific expression of steady-state input as long as the steady-state generator exists. In addition, the nonlinear regulator of the generic internal model ensures robust asymptotic regulation against unstructured uncertainties, as shown in [16]. Moreover, the generic internal model can directly provide the unknown parameters arising from the exosystem, eliminating the need for the adaptive control technique and having a significant advantage compared to the canonical linear internal model [7]. In fact, the adaptive control approach faces two challenges that impede further research on the adaptive output regulation problem. Firstly, it requires the construction of a specific Lyapunov function for the nonlinear time-varying adaptive system with a positive semidefinite derivative to ensure the convergence of the partial state by using the LaSalle–Yoshizawa Theorem. This yields weaker stability, resulting in the absence of effective analysis techniques. Therefore, it only applies to a class of uncertain nonlinear systems in the parametric form. Moreover, another feature of adaptive techniques requires a known and explicit regressor determined by the controlled system structure and the exosystem, which can be found in [17] and in the nonlinear regression case [18]. Secondly, when dealing with the nonlinear time-varying adaptive system in the presence of external inputs, the derivative of the Lyapunov function constructed through the adaptive method incorporates a negative semidefinite term and external inputs. Nevertheless, the negative semidefinite term is insufficient to guarantee the boundedness property, even with a small input term (like noise input), as illustrated using a counterexample in [19].

In contrast, the nonadaptive method proposed by [20] removes the need for constructing such Lyapunov functions for the closed-loop system with positive semidefinite definite derivatives. It also did not require the studied nonlinear system to be in parametric form. As a consequence, the nonadaptive method for solving output regulation problems has received considerable attention [3], [21], [22]. Nevertheless, the generic internal model-based method for solving output regulation strictly relies on the explicit construction of a nonlinear continuous mapping function, which is assumed to exist [23]. Just as [16] pointed out, no general analytical expression is known for the nonlinear regulator construction of the generic internal model-based method. Consequently, approximation methods have been proposed using system identification least-squares techniques to select optimized parameters for the construction of nonlinear continuous mapping functions [3], [21]. It has been demonstrated that the steady-state tracking error in a nonlinear system is a nonlinear function of the exogenous signals [6]. In recent developments, [11] proposed a nonparametric learning framework for the construction of the nonlinear mapping in a general setting by establishing connections between the Generalized Sylvester Matrix Equation, the generic internal model, and a nonlinear Luenberger observer design. This framework enables the steady-state generator to be polynomial in the exogenous signal for nonlinear output regulation while relaxing restrictive and somewhat ad hoc assumptions on the exosystem, such as the stringent requirement for an even dimension and the absence of a zero eigenvalue imposed in some results.

Based on the aforementioned statements, this article aims to address the nonlinear robust output regulation for general nonlinear output feedback systems with higher relative degrees using nonparametric learning methods. These methods are in contrast with existing techniques that are based on common adaptive control methods [17], [24]–[26]. The output regulation problem for uncertain nonlinear systems with higher relative degree remains an active research area [24], [27], with some interesting results on the design of low-complexity, approximation-free, output-feedback controller to achieve output tracking with prescribed transient and steady-state performance [27]. The proposed nonparametric learning framework addresses the global robust output regulation problem via the construction of a linear generic internal model contingent on the existence of a continuous nonlinear mapping. This approach enables the transformation of the nonlinear robust output regulation problem into a robust nonadaptive stabilization design problem for an augmented system with Input-to-State Stable (ISS) dynamics. Furthermore, integrating a nonparametric learning framework ensures the viability of the nonlinear mapping, demonstrating its capability to capture intricate and nonlinear relationships without the need for the existence of a predefined model of the steady-state input behavior.

In contrast to the methodologies in [7], [24], [26] and [17], which rely on system dynamics to satisfy the linear-in-parameter (parameterized) condition and employ adaptive learning, the nonparametric learning framework directly learns the system dynamics of the steady-state input with a data-driven technique that requires the construction of a Hankel matrix using the internal model signal. Moreover, this nonparametric approach can be considered a model-free method, as it eliminates the need for a regressor, which, in adaptive methods, is typically dictated by the structure of the system and the exosystem. The proposed methodology also overcomes the construction of a specific Lyapunov function whose derivatives can only be shown to be negative semidefinite for a restricted class of systems in adaptive methods. In addition, this study also provides an alternative proof for Lemma 3 of [11] in terms of the solution for a time-varying equation. In particular, oscillatory disturbances affecting the current of a repulsive magnetic levitation system, the power angle of a virtual synchronous generator, and the heading of a marine surface vessel induce unknown steady-state chattering. This behavior is significantly more complex than the constant offset disturbances addressed in earlier studies. Overall, the nonparametric learning framework can contribute to advancing robust control strategies and learning the unknown steady-state behaviour of nonlinear systems with broader implications in many engineering applications.

The rest of this article is organized as follows. Section 2 introduces some standard assumptions and lemmas. Section 3 is devoted to the presentation of the main results, which is followed by simulation examples in Section 4 and brief conclusions in Section 5.

Notation: \(\|\cdot\|\) is the Euclidean norm. \(\emph{Id} : \mathbb{R}\rightarrow \mathbb{R}\) is an identity function. For \(X_i\in \mathbb{R}^{n_i\times m}\) with \(i=1,\dots,N\), let \(col(X_1,\dots,X_N)=[X_{1}^{\!\top},\dots ,X^{\!\top}_{N}]^{\!\top}\) and \[\textrm{diag}(X_1,\dots,X_N)=\textstyle \left[\begin{matrix} X_1 & & \\ & \ddots & \\ & & X_N\end{matrix}\right].\] A function \(\alpha: \mathbb{R}_{\geq 0}\rightarrow \mathbb{R}_{\geq 0}\) is of class \(\mathcal{K}\) if it is continuous, positive definite, and strictly increasing. The notation \(\mathcal{K}_o\) and \(\mathcal{K}_{\infty}\) identifies the subclasses of bounded and unbounded \(\mathcal{K}\) functions, respectively. For functions \(f_1(\cdot)\) and \(f_2(\cdot)\) with compatible dimensions, their composition \(f_{1}(f_2(\cdot))\) is denoted by \(f_1\circ f_2(\cdot)\). For two continuous and positive definite functions \(\kappa_1(\varsigma)\) and \(\kappa_2(\varsigma)\), \(\kappa_1\in \mathcal{O}(\kappa_2)\) means that \(\limsup_{\varsigma\rightarrow 0^{+}}\frac{\kappa_1(\varsigma)}{\kappa_2(\varsigma)}<\infty\).

2 Problem Formulation and Preliminaries↩︎

2.1 System Model and Problem Formulation↩︎

We consider a class of nonlinear control systems modelled by \[\tag{1} \begin{align} \dot{z}&=f(z,y,v,w),\tag{2}\\ \dot{x}&= A_cx+g(z,y,v,w)+B_cbu,\\ y &=C_cx, \\ e &= y - h(v,w), \end{align}\] where the state vector comprises \(z\in \mathbb{R}^{n_z}\) and \(x\in \mathbb{R}^{r}\). The \(z\)-subsystem represents the internal dynamics, while the \(x\)-subsystem represents a chain of \(r\) integrators (\(r\geq 1\)) perturbed by the nonlinear term \(g(z,y,v,w)=col(g_1,\dots, g_r)\). In this setup, \(y\in \mathbb{R}\) is the measured output, \(e\in \mathbb{R}\) is the tracking error, \(u\in \mathbb{R}\) is the control input, and \(b\) is a positive constant (unknown high-frequency gain). The vector \(w\in \mathbb{W}\subset \mathbb{R}^{n_w}\) represents uncertain parameters, where \(\mathbb{W}\) is a prescribed compact set containing the origin.

The functions \(h(\cdot)\), \(f(\cdot)\), and \(g_i(\cdot)\) are sufficiently smooth and vanish at the origin, satisfying \(h(0,w)=0\), \(f(0,0,0,w)=0\), and \(g_i(0,0,0,w)=0\) for all \(w\in \mathbb{W}\). The exogenous signal \(v(t)\in \mathbb{R}^{n_v}\) represents the reference input and/or disturbance, generated by the exosystem: \[\begin{align} \label{eqn:32exosystem32system} \dot{v}=S(\sigma)v, \end{align}\tag{3}\] where \(\sigma \in \mathbb{S}\subset \mathbb{R}^{n_\sigma}\) represents the exosystem uncertainty, and \(S(\sigma)\) is a matrix depending on \(\sigma\). The matrices \(A_c\in \mathbb{R}^{r\times r}\), \(B_c\in \mathbb{R}^r\), and \(C_c\in \mathbb{R}^{1\times r}\) form a canonical representation: \[\begin{align} A_c=\left[\begin{matrix}\boldsymbol{0}& I_{r-1}\\ 0& \boldsymbol{0} \end{matrix}\right]\!, \quad C_c=col(1, \boldsymbol{0}_{r-1})^\top, \quad B_c=col(\boldsymbol{0}_{r-1},1), \end{align}\] where \(I_{r-1}\) is the identity matrix and \(\boldsymbol{0}_{r-1}\) is the zero vector of dimension \(r-1\). The system notation and description is summarized in Table 1.

The nonlinear robust output regulation problem is formulated as follows:

2.2 Assumptions and Coordinate Transformation↩︎

We first state the standard assumptions regarding the exosystem and the regulator equations.

To facilitate the control design for the system with relative degree \(r \geq 2\), we employ the input-driven filter as in [29], \[\begin{align} \label{input-filter} \dot{\hat{x}}=A\hat{x}+B_cu, \end{align}\tag{4}\] where \(A=A_c-\lambda C_c\) is Hurwitz with \(\lambda=col(\lambda_1,\dots, \lambda_r)\). Applying the coordinate transformation \(\tilde{x}_i=b^{-1}x_i-\hat{x}_i\) (\(i=1,\dots,r\)) yields the transformed system: \[\label{nonlinear-systems-trans} \begin{align} \dot{z}&=f(z,y,v,w),\\ \dot{\tilde{x}}&= A\tilde{x}+b^{-1}(\lambda y+g(z,y,v,w)),\\ \dot{y} &=b\hat{x}_2+b\tilde{x}_2+g_1(z,y,v,w), \\ \dot{\hat{x}}_i&= \hat{x}_{i+1}-\lambda_i\hat{x}_1,\quad i=2,\dots,r-1,\\ \dot{\hat{x}}_r&=u-\lambda_r\hat{x}_1. \end{align}\tag{5}\]

Ideally, the aim is for the filter state \(\hat{x}(t)\) to track a steady-state profile induced by the ideal control \(\boldsymbol{u}(v,w,\sigma)\).

The following lemma establishes the structure of the steady-state generator for the filter states, which is crucial for the internal model design.

Since \(S(\sigma)\) has eigenvalues on the imaginary axis (Assumption 1), the generator matrix \(\Phi_u(\sigma)\) for the polynomial basis \(\tau_u(v)\) also shares this spectral property. Given that \(A\) is Hurwitz, the spectra of \(A\) and \(\Phi_u(\sigma)\) are disjoint. Thus, the Sylvester equation \(P_u \Phi_u - A P_u = B_c \Gamma_u\) has a unique solution \(P_u\). Differentiating \(\boldsymbol{\hat{x}} = P_u \tau_u\) yields \[\begin{align} \dot{\boldsymbol{\hat{x}}} &= P_u \Phi_u \tau_u\\ &= (A P_u + B_c \Gamma_u)\tau_u \\ &= A \boldsymbol{\hat{x}} + B_c \boldsymbol{u}. \end{align}\] This completes the proof. \(\Box\)

Defining the error \(\boldsymbol{E}(\mu)=b^{-1}\boldsymbol{x}(\mu)-\boldsymbol{\hat{x}}(\mu)\), the solution to the regulator equations for the composite system is \[\{\boldsymbol{z}(\mu), \boldsymbol{E}(\mu), \boldsymbol{y}(\mu), \boldsymbol{\hat{x}}(\mu), \boldsymbol{u}(\mu)\}.\]

2.3 Generic Internal Model Design↩︎

Under Assumptions 1 and 4, there exists a positive integer \(n\) such that the steady-state signal \(\boldsymbol{\hat{x}}_2(\mu)\) satisfies the autonomous differential equation: \[\begin{align} \label{aode-explicit} \frac{d^{n}\boldsymbol{\hat{x}}_2}{dt^{n}}+a_{n}(\sigma)\frac{d^{n-1}\boldsymbol{\hat{x}}_2}{dt^{n-1}}+\dots+a_{1}(\sigma)\boldsymbol{\hat{x}}_2=0, \end{align}\tag{6}\] where \(a_{i}(\sigma) \in \mathbb{R}\). By defining the state vector \(\boldsymbol{\xi} = col(\boldsymbol{\hat{x}}_2, \dot{\boldsymbol{\hat{x}}}_2, \dots, \boldsymbol{\hat{x}}_2^{(n-1)})\), 6 can be represented in state-space form, referred to as the steady-state generator, \[\label{stagerator} \begin{align} \dot{\boldsymbol{\xi}} &= \Phi(a(\sigma)) \boldsymbol{\xi}, \\ \boldsymbol{\hat{x}}_2 &= \Gamma \boldsymbol{\xi}, \end{align}\tag{7}\] where the companion matrix \(\Phi(a(\sigma))\) and the output vector \(\Gamma\) are defined as \[\begin{align} \Phi(a(\sigma)) &= \left[ \begin{array}{c|c} \boldsymbol{0}_{(n-1)\times 1} & I_{n-1} \\ \hline -a_{1}(\sigma) & -a_{2}(\sigma),\dots,-a_{n}(\sigma) \\ \end{array} \right],\\ \Gamma &= \left[ 1, 0, \dots, 0 \right]_{1\times n}. \end{align}\]

To design an internal model that reproduces \(\boldsymbol{\hat{x}}_2\), introduce a dynamic system with state \(\eta \in \mathbb{R}^{2n}\), \[\begin{align} \label{explicit-mas1} \dot{\eta} = M\eta + N\hat{x}_2, \end{align}\tag{8}\] where the pair \((M, N)\) is chosen in the controllable canonical form: \[\label{MNINter} \begin{align} M&= \left[ \begin{array}{c|c} \boldsymbol{0}_{(2n-1)\times 1} & I_{2n-1} \\ \hline -m_{1} & -m_{2},\dots,-m_{2n} \\ \end{array} \right], \\ N&= \left[ 0, \dots, 0, 1 \right]_{1\times 2n}^{\top}. \end{align}\tag{9}\] The coefficients \(m_j\) are selected such that the matrix \(M\) is Hurwitz. Let \(P_M(\lambda) = \lambda^{2n} + \sum_{j=1}^{2n} m_j \lambda^{j-1}\) denote the characteristic polynomial of \(M\). Since the eigenvalues of \(\Phi(a)\) are distinct with zero real parts (by Assumption 4) and \(M\) is Hurwitz, their spectra are disjoint. This condition ensures nonsingularity of the matrix polynomial function: \[\Xi(a) = \underbrace{\Phi(a)^{2n}+\sum\nolimits_{j=1}^{2n}m_{j}\Phi(a)^{j-1}}_{P_M(\Phi(a))} \in \mathbb{R}^{n \times n }.\]

The Generalized Sylvester Equation↩︎

To immerse the generator dynamics 7 into the internal model 8 , we seek a transformation matrix \(Q \in \mathbb{R}^{2n \times n}\) satisfying the Generalized Sylvester Matrix Equation (see [11]): \[\begin{align} M Q - Q \Phi(a(\sigma)) = -N\Gamma. \label{MNGAMMAPhi} \end{align}\tag{10}\] Due to the specific canonical structure of \((M, N)\), the unique solution \(Q\) admits an explicit form. Let \(Q\) be partitioned as \(Q = col(Q_1, \dots, Q_{2n})\) where \(Q_j \in \mathbb{R}^{1 \times n}\). Using the commutativity property \(\Phi(a)\Xi(a)^{-1}=\Xi(a)^{-1}\Phi(a)\) and the structural identity \(col(\Gamma, \Gamma \Phi(a),\dots, \Gamma\Phi(a)^{n-1} )=I_n\), the upper block of \(Q\) can be explicitly derived as \[\begin{align} col& (Q_1(a),\dots,Q_n(a))\nonumber\\ =&\;col(\Gamma \Xi(a)^{-1},\Gamma \Phi(a)\Xi(a)^{-1},\dots, \Gamma \Phi(a)^{n-1}\Xi(a)^{-1})\nonumber\\ =&\;\underbrace{col(\Gamma,\Gamma \Phi(a),\dots, \Gamma \Phi(a)^{n-1})}_{I_n} \Xi(a)^{-1}, \label{XIQA-explicit} \end{align}\tag{11}\] with the general term given by \(Q_{j}(a)=\Gamma \Xi(a)^{-1}\Phi(a)^{j-1}\).

Finally, defining the steady-state internal model coordinate \(\boldsymbol{\eta}^{\star} = Q \boldsymbol{\xi}\), we construct the Hankel matrix \(\Theta(\boldsymbol{\eta}^{\star})\) as in [30]: \[\begin{align} \Theta (\boldsymbol{\eta}^{\star}) \equiv \left[\begin{matrix}\eta^{\star}_{1} &\dots&\eta^{\star}_{n}\\ \vdots&\ddots&\vdots\\ \eta^{\star}_{n} &\dots&\eta^{\star}_{2n -1} \end{matrix}\right] \in \mathbb{R}^{n \times n}. \end{align}\] As shown in Lemma 3 of [11], substituting the solution of 10 into the integral form of the internal model ensures that \(\boldsymbol{\eta}^{\star}\) satisfies the target dynamics \(\dot{\boldsymbol{\eta}}^{\star} = M\boldsymbol{\eta}^{\star} + N\boldsymbol{\hat{x}}_2\), completing the design.

2.4 Error Dynamics↩︎

Perform coordinate and input transformations on the composite systems 3 , 5 , and 8 to give \[\begin{align} \bar{z}&= z-\boldsymbol{z},& \bar{x}= \tilde{x}-\boldsymbol{E}, \\ \bar{\eta}&= \eta-\boldsymbol{\eta}^{\star}-Nb^{-1}e,& e=y-\boldsymbol{x}_1, \end{align}\] which yields an error system in the form: \[\label{Main-sys1}\begin{align} \dot{\bar{z}} &= \bar{f}(\bar{z},e,\mu),\\ \dot{\bar{x}}&=A\bar{x}+ b^{-1}\!\big[\bar{g}(\bar{z}, e, \mu)+\lambda e\big],\\ \dot{\bar{\eta}}&=M\bar{\eta} -N\!\left(\bar{x}_2 -b^{-1}e+b^{-1}\bar{g}_1(\bar{z}, e, \mu)\right)\!,\\ \dot{e}&=b(\hat{x}_2-\boldsymbol{\hat{x}}_2)+b\bar{x}_2 +\bar{g}_1(\bar{z}, e, \mu) ,\\ \dot{\hat{x}}_i &= \hat{x}_{i+1}-\lambda_i\hat{x}_1,\quad i=2,\dots,r-1,\\ \dot{\hat{x}}_r &= u-\lambda_r\hat{x}_1, \end{align}\tag{12}\] where \(\mu=col(\sigma,v,w)\), \[\begin{align} \bar{f}(\bar{z},e,\mu) &=f(\bar{z}+\boldsymbol{z},e+\boldsymbol{x}_1,\mu)-f(\boldsymbol{z},\boldsymbol{x}_1,\mu),\\ \bar{g}(\bar{z}, e, \mu) &= g(\bar{z}+\boldsymbol{z},e+\boldsymbol{x}_1,\mu)-g(\boldsymbol{z},\boldsymbol{x}_1,\mu). \end{align}\] It can be verified that, for all \(\mu\in \mathbb{V}\times\mathbb{W} \times\mathbb{S}\), \(\bar{f}(0,0,\mu)=0\) and \(\bar{g}(0, 0, \mu)=0\), Problem 1 can be solved if a control law can be found to stabilize the system 12 .

Let \(\bar{x}_c=col(\bar{\eta}, \bar{x})\) and \(\bar{G}_c(\bar{z}, e, \mu)=b^{-1}col\big(Ne-N\bar{g}_1(\bar{z}, e, \mu), \bar{g}(\bar{z}, e, \mu)+\lambda e\big)\); the system 12 can be rewritten into the form \[\tag{13}\begin{align} \dot{\bar{z}} &=\bar{f}(\bar{z},e,\mu),\\ \dot{\bar{x}}_c &= \underbrace{\left[\begin{matrix}M & -NC_cA_c\\ \boldsymbol{0}&A \end{matrix}\right]}_{M_c}\bar{x}_c+\bar{G}_c(\bar{z}, e, \mu) ,\tag{14}\\ \dot{e} &=b(\hat{x}_2-\chi(\boldsymbol{\eta}^*))+b\bar{x}_2 +\bar{g}_1(\bar{z}, e, \mu) ,\\ \dot{\hat{x}}_i &= \hat{x}_{i+1}-\lambda_i\hat{x}_1, \quad i=2,\dots,r-1,\\ \dot{\hat{x}}_r &=u-\lambda_r \hat{x}_1. \end{align}\] It can be verified that, for all \(\mu\in \mathbb{V}\times\mathbb{W} \times\mathbb{S}\), \(\bar{G}_c(0, 0, \mu)=\boldsymbol{0}\) and the matrix \(M_c\) is Hurwitz. Hence, the \((\bar{z}, \bar{x}_c)\)-subsystem in system 13 is in a similar form as the system (8) of [11]. As a result, the \((\bar{z}, \bar{x}_c)\)-subsystem in system 13 , under Assumptions 1, 2, and 3, admits the following properties (see Properties 1 and 2 in [11]):

3 Main Results↩︎

3.1 Nonparametric Learning in Robust Output Regulation↩︎

We now provide an alternative proof of [11] Lemma 3 in terms of the time-varying equation \[\begin{align} \label{a-explicit0} \Theta (\boldsymbol{\eta} ^{\star} )a +col(\boldsymbol{\eta}_{n +1}^{\star},\dots,\boldsymbol{\eta} ^{\star}_{2n})=\boldsymbol{0}, \end{align}\tag{15}\] where \(\boldsymbol{\eta} ^{\star}=\underbrace{Q \boldsymbol{\xi}}_{\theta}\).

From \(Q = col(Q_1, \dots, Q_{2n})\) and \(\theta =Q \boldsymbol{\xi}\), \(Q_{j}(a)=\Gamma \Xi(a)^{-1}\Phi(a )^{j-1} \in \mathbb{R}^{1\times n }\), \(1\leq j\leq 2n\), the real Hankel matrix is given by \[\begin{align} \Theta (\theta)=& \left[\begin{matrix}Q_1 \boldsymbol{\xi}&Q_2 \boldsymbol{\xi}&\cdots&Q_n \boldsymbol{\xi}\\ Q_2 \boldsymbol{\xi}&Q_3 \boldsymbol{\xi}&\cdots&Q_{n+1} \boldsymbol{\xi}\\ \vdots&\vdots&\ddots&\vdots\\ Q_{n}\boldsymbol{\xi} &Q_{n +1}\boldsymbol{\xi}&\cdots&Q_{2n -1}\boldsymbol{\xi} \end{matrix}\right]\nonumber\\ =& \left[\begin{matrix}Q_1 \boldsymbol{\xi}&Q_1\Phi(a) \boldsymbol{\xi}&\cdots&Q_1 \Phi(a)^{n -1}\boldsymbol{\xi}\\ Q_2 \boldsymbol{\xi}&Q_2\Phi(a) \boldsymbol{\xi}&\cdots&Q_{2} \Phi(a)^{n -1}\boldsymbol{\xi}\\ \vdots&\vdots&\ddots&\vdots\\ Q_{n}\boldsymbol{\xi} &Q_{n}\Phi(a)\boldsymbol{\xi}&\cdots&Q_{n}\Phi(a)^{n -1}\boldsymbol{\xi} \end{matrix}\right]\nonumber\\ =&\;\underbrace{col(Q_1,\dots,Q_n)}_{\Xi(a)^{-1}} \underbrace{\left[\begin{matrix}\boldsymbol{\xi} &\Phi(a) \boldsymbol{\xi} &\dots & \Phi(a)^{n -1} \boldsymbol{\xi} \end{matrix}\right]}_{\Pi}.\nonumber \end{align}\] where the columns of the Krylov matrix \(\Pi\) form the Krylov subspace. Under Assumption 1, the matrix \(\Phi(a)\) is diagonalizable with distinct eigenvalues \(\lambda_1=\imath \hat{\omega}_{1},\dots,\lambda_n=\imath \hat{\omega}_{n}\). Moreover, the matrix \(\Phi(a)\) is in companion form. Therefore, from [31] and [32], there exists a diagonalizable matrix \(\Lambda\) with \(\Lambda=\textrm{diag}(\lambda_1,\dots,\lambda_n)\) and a nonsingular Vandermonde matrix \[P_{\Lambda}=\left[\begin{matrix}1 & 1 & \dots &1\\ \lambda_1 & \lambda_2 & \cdots& \lambda_n\\ \vdots &\vdots& \ddots & \vdots \\ \lambda_1^{n-1} & \lambda_2^{n-1} & \cdots& \lambda_n^{n-1}\\ \end{matrix}\right]\] such that \(\Phi(a)= P_{\Lambda}\Lambda P_{\Lambda}^{-1}\). As a result, from 7 , let \(\nu(t)=P_{\Lambda}^{-1}\boldsymbol{\xi}(t)\), which results in \[\nu(t)=\underbrace{\textrm{col}(e^{\lambda_1 t}\nu_1(0),\dots, e^{\lambda_n t}\nu_n(0))}_{e^{\Lambda t}\nu(0)}.\] Hence, the time-varying matrix \(\Pi(t)\) admits the form \[\begin{align} \Pi(t)=&\;\left[\begin{matrix}\boldsymbol{\xi}(t) &\Phi(a) \boldsymbol{\xi}(t) &\cdots & \Phi(a)^{n-1} \boldsymbol{\xi}(t) \end{matrix}\right]\\ =&\;P_{\Lambda} \left[\begin{matrix}\nu(t) &\Lambda \nu(t) &\cdots & \Lambda^{n -1}\nu(t) \end{matrix}\right]\\ =&\;P_{\Lambda} \textrm{diag}(e^{\lambda_1 t}\nu_1(0),\dots, e^{\lambda_n t}\nu_n(0))P_{\Lambda}^{\top}. \end{align}\] It is noted from \(\boldsymbol{\xi}= col\!\left(\boldsymbol{\hat{x}}_2,\frac{d\boldsymbol{\hat{x}}_2}{dt},\dots,\frac{d^{n-1}\boldsymbol{\hat{x}}_2}{dt^{n-1}}\right)\!\) and ?? that \[\begin{align} \nu(0)&= P_{\Lambda}^{-1}\boldsymbol{\xi}(0)\\ =&\; P_{\Lambda}^{-1} col\Big(\sum\limits_{i=1}^{n}C_{j}(v(0), w,\sigma),\dots, \sum\limits_{j=1}^{n}C_{j}(v(0), w,\sigma)\lambda_j^{n-1}\Big)\\ =&\; \underbrace{P_{\Lambda}^{-1}P_{\Lambda}}_{I_n}col(C_{1}(v(0), w,\sigma),\dots, C_{n}(v(0), w,\sigma)). \end{align}\] From Assumption 5, for any \(v(0)\in \mathbb{V}\), \(w\in \mathbb{W}\), and \(\sigma\in \mathbb{S}\), \(C_{i}(v(0), w,\sigma)\neq 0\) results in \(\nu_i(0)\neq 0\), for \(i=1,\dots,n\). Hence, the matrix \(\Pi(t)\) is nonsingular due to the fact that \(\nu_i(0)\neq 0\), \(\forall i in\in\{i,\dots,n\}\). Therefore, \(\Theta (\theta)\) is nonsingular. Equation 11 admits the solution \[\begin{align} \label{a-explicit} a =\underbrace{-\Theta(\theta)^{-1}\textrm{col}(\theta_{n+1},\dots,\theta_{2n})}_{\check{a} (\theta)}. \end{align}\tag{16}\] \(\Box\)

From Lemma 3 in [11], the existence of a nonlinear mapping \(\chi\left(\eta, \check{a} (\eta )\right)\) strictly relies on the solution of a time-varying equation, \[\Theta(\eta)\check{a} (\eta)+\textrm{col}(\eta_{n +1},\dots,\eta_{2n})=\boldsymbol{0}.\] It is noted that \(\Theta(\eta(t))\) is not always invertible over \(t\geq 0\), and there may be time instants where the inverse of \(\Theta(\eta(t))\) may not be well-defined. From Assumptions 1 and 4, it follows that \(\boldsymbol{\eta} ^{\star}\) and \(a\) belong to some compact set \(\mathbb{D}\). For the composite system 1 , as shown in Fig. 1, we propose the regulator \[\tag{17} \begin{align} \dot{\hat{a}} &=- k_a \Theta(\eta)^{\!\top}\left[\Theta (\eta)\hat{a} +\textrm{col}(\eta_{n +1},\dots{},\eta_{2n})\right],\\ u &= \alpha_{s,r}(\epsilon_1,\epsilon_2,\dots,\epsilon_r,k^*,\eta,\hat{a}),\tag{18} \end{align}\] where \(k_a\) is positive constant, \(\eta\) is generated in 8 , \(\hat{a}\) is the estimate of the unknown parameter vector \(a\), \(\rho (\cdot)\geq 1\) is a positive smooth function, \(\epsilon_1=e\),\[\begin{align} \label{alpha-function-s} \epsilon_{i+1}=&\;\hat{x}_{i+1}-\alpha_{s,i}(\epsilon_1,\dots,\epsilon_i,k^*,\eta,\hat{a},\hat{x}_1),\nonumber\\ \alpha_{s,1}( \epsilon_1,k^*,\eta,\hat{a})=&-k^*\rho(\epsilon_1)\epsilon_1+\chi_s(\eta,\hat{a}),\nonumber\\ \alpha_{s,2}(\epsilon_1,\epsilon_2,k^*,\eta,\hat{a},\hat{x}_1)=&-b\epsilon_1-\epsilon_2+\lambda_{2}\hat{x}_1+\frac{\partial \alpha_{s,1} }{\partial \eta}\dot{\eta}\nonumber\\ &+b\frac{\partial \alpha_{s,1}}{\partial \epsilon_1}(\epsilon_2-{k}^*\rho(\epsilon_1)\epsilon_1)\nonumber\\ &-\frac{1}{2}\epsilon_2\!\left(\frac{\partial \alpha_{s,1}}{\partial \epsilon_1}\right)^{\!2}+\frac{\partial \alpha_{s,1} }{\partial \hat{a}}\dot{\hat{a}}\nonumber\\ &+\frac{\partial \alpha_{s,1} }{\partial k^*}\dot{k}^*,\nonumber\\ \alpha_{s,i}(\epsilon_1,\dots,\epsilon_i,k^*,\eta,\hat{a},\hat{x}_1)=&-\epsilon_{i-1}-\epsilon_i+\lambda_{i}\hat{x}_{1}+\frac{\partial \alpha_{s,i-1} }{\partial \eta}\dot{\eta}\nonumber\\ &+\frac{\partial \alpha_{s,i-1} }{\partial \hat{x}_1}\dot{\hat{x}}_1+ \sum_{j=2}^{i-1}\frac{\partial \alpha_{s,i-1} }{\partial \epsilon_{j}}\dot{\epsilon}_j\nonumber\\ &+b\frac{\partial \alpha_{s,i-1}}{\partial \epsilon_1}(\epsilon_2- {k}^*\rho(\epsilon_1)\epsilon_1)\nonumber\\ &-\frac{1}{2}\epsilon_i\!\left(\frac{\partial \alpha_{s,i-1}}{\partial \epsilon_1}\right)^{\!2}\nonumber\\ &+\frac{\partial \alpha_{s,i} }{\partial \hat{a}}\dot{\hat{a}}+\frac{\partial \alpha_{s,1} }{\partial k^*}\dot{k}^*,\nonumber\\ & \quad\quad\quad\quad\quad i=3,\dots,r, \end{align}\tag{19}\] where \(\dot{k}^*\) will be zero when \(k^*\) is a constant, \(\hat{a}\) is generated in 17 , \(\hat{x}_{1}, \dots, \hat{x}_{r}\) and \(\eta\) are generated in 4 and 8 , respectively. The smooth function \(\chi_s(\eta,\hat{a})\) is given by \[\begin{align} \label{chisatu} \chi_{s}(\eta, \hat{a})=\chi(\eta, \hat{a})\Psi(\delta+1-\|col(\eta, \hat{a})\|^2), \end{align}\tag{20}\] where \[\chi(\eta,\hat{a})\equiv\Gamma \Xi (\hat{a})\textrm{col}(\eta_{1},\dots{},\eta_{n}),\] with \(\Psi(\varsigma)=\frac{\psi(\varsigma)}{\psi(\varsigma)+\psi(1-\varsigma)}\), \(\delta=\max\limits_{(\eta,\hat{a})\in \mathbb{D}}\|col(\eta, \hat{a})\|^2\) and \[\psi(\varsigma)=\left\{\begin{array}{cc}e^{-1/\varsigma} & \textrm{for}\;\varsigma >0,\\ 0 & \textrm{for}\;\varsigma \leq0. \end{array}\right.\] Performing the coordinate/input transformations \[\begin{align} &\bar{\eta}_e =\bar{\eta}+b^{-1}N\epsilon_1,\quad\bar{a}=\hat{a}-a,\quad\hat{x}_2=\epsilon_2+\alpha_{s,1},\\ &\bar{\chi}_s(\bar{\eta}_e,\bar{a},\mu)= \chi_s(\bar{\eta}_e+\boldsymbol{\eta}^*,\bar{a}+a)-\chi(\boldsymbol{\eta}^*,a), \end{align}\] leads to the augmented system: \[\tag{21}\begin{align} \dot{\bar{z}} =&\;\bar{f}(\bar{z},\epsilon_1,\mu),\\ \dot{\bar{x}}_c= &\;M_c\bar{x}_c+\bar{G}_c(\bar{z}, \epsilon_1, \mu), \\ \dot{e}_1=&\;b(\epsilon_2-k^*\rho(\epsilon_1)\epsilon_1)+b\bar{\chi}_s(\bar{\eta}_e,\bar{a},\mu)\nonumber\\ &\;+b\bar{x}_2 +\bar{g}_1(\bar{z}, \epsilon_1, \mu), \\ \dot{\hat{x}}_i=&\;\hat{x}_{i+1}-\lambda_i\hat{x}_1,\;\;i=2,\dots,r-1,\;\;\\ \dot{\hat{x}}_r=&\;u-\lambda_r\hat{x}_1,\\ \dot{\bar{a}}=&-k_a \Theta(\boldsymbol{\eta}^{\star})^{\!\top}\Theta (\boldsymbol{\eta}^{\star} )\bar{a} - k_1 \bar{O}(\bar{\eta}_e,\bar{a}),\tag{22} \end{align}\] where \[\begin{align} \bar{O}(\bar{\eta}_e,\bar{a})=&\; \Theta(\boldsymbol{\eta}^{\star})^{\!\top}\Theta (\bar{\eta}_e)\bar{a}+\Theta (\bar{\eta}_e)^{\!\top}\Theta (\boldsymbol{\eta}^{\star} )\bar{a}\\ &+\Theta(\bar{\eta}_e)^{\!\top} \Theta(\bar{\eta}_e)a +\Theta(\bar{\eta}_e)^{\!\top} \Theta(\bar{\eta}_e)\bar{a}\\ &+\Theta(\boldsymbol{\eta}^{\star})^{\!\top} \Theta(\bar{\eta}_e){a} + \Theta(\boldsymbol{\eta}^{\star})^{\!\top} \textrm{col}(\bar{\eta}_{e,n +1},\dots{},\bar{\eta}_{e,2n})\\ &+\Theta(\bar{\eta}_e)^{\!\top}\textrm{col}(\bar{\eta}_{e,n +1},\dots{},\bar{\eta}_{e,2n}). \end{align}\] Moreover, the \(\bar{a}\)-subsystem 22 has a form similar to system (27d) in [11]. As a result, the \(\bar{a}\)-subsystem 22 , under Assumptions 1, 2, and 3, admits the following lemma (see Lemma 4 in [11]).

From Property 1, the changing supply rate technique [28] can be applied to show that, given any smooth function \(\Delta_{Z}(\bar{Z})>0\), there exists a continuous function \(V_{1}(\bar{Z})\) satisfying \[\underline{\alpha}_{1}\big(\big\|\bar{Z} \big\|^2\big)\leq V_{1}\big( \bar{Z} \big)\leq\overline{\alpha}_{1}\big(\big\|\bar{Z} \big\|^2\big)\] for some class \(\mathcal{K}_{\infty}\) functions \(\underline{\alpha}_{1}(\cdot)\) and \(\overline{\alpha}_{1}(\cdot)\) such that, for all \(\mu\in \Sigma\), along the trajectories of the \(Z\) subsystem, \[\dot{V}_{1} \leq-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2+ \hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\epsilon_1^2,\] where \(\hat{\gamma}^*\) is known positive constant and \(\hat{\gamma} \left(\cdot\right)\geq 1\) is a known smooth positive definite function.

Define the Lyapunov function \(U_1(\bar{Z}, \epsilon_1)=V_{1}\big( \bar{Z} \big)+ \epsilon_1^2\). Then, the time derivative of \(U_1\equiv U_1(\bar{Z}, \epsilon_1)\) along the trajectory of \(\epsilon_1\)-subsystem with \(\hat{x}_2=\epsilon_2+\alpha_1\) and \(\eta=\bar{\eta}+\boldsymbol{\eta}^{\star}+Nb^{-1}\epsilon_1\) leads to \[\begin{align} \label{U1-derivative} \dot{U}_1(\bar{Z}, \epsilon_1)=&\;\dot{V}_{1}\big( \bar{Z} \big)+ 2\epsilon_1\dot{\epsilon}_1\nonumber\\ \leq &-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2+ \hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\epsilon_1^2+2\epsilon_1\bar{g}_1(\bar{z}, \epsilon_1, \mu) \nonumber\\ &+2b\epsilon_1(\underbrace{\epsilon_2+\alpha_1( \epsilon_1,k^*,\eta)}_{\hat{x}_2}-\chi(\boldsymbol{\eta}^*))+2b\epsilon_1\bar{x}_2\nonumber\\ \leq &-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2+ \hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\epsilon_1^2+2\epsilon_1\bar{g}_1(\bar{z}, \epsilon_1, \mu)\nonumber\\ &+2b\epsilon_1(\epsilon_2+\underbrace{\alpha_1( \epsilon_1,k^*,\eta)}_{-k^*\rho(\epsilon_1)\epsilon_1+\chi(\eta)}-\chi(\boldsymbol{\eta}^*))+2b\epsilon_1\bar{x}_2 \nonumber\\ \leq &-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2-\big(2b k^*\rho(\epsilon_1)-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2\nonumber\\ &+2b\epsilon_1(\epsilon_2-\bar{\chi}(\bar{\eta},\epsilon_1,\mu))\nonumber\\ &+2b\epsilon_1\bar{x}_2 +2\epsilon_1\bar{g}_1(\bar{z}, \epsilon_1, \mu)\nonumber\\ \leq &-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2-\big(2b k^*\rho(\epsilon_1)-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2\nonumber\\ &+2b\epsilon_1\epsilon_2-2b\epsilon_1\bar{\chi}(\bar{\eta},\epsilon_1,\mu)\nonumber\\ &+2b\epsilon_1\bar{x}_2 +2\epsilon_1\bar{g}_1(\bar{z}, \epsilon_1, \mu)\nonumber\\ \leq &-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2-\big(2b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2\nonumber\\ &+2b\epsilon_1\epsilon_2+\Delta_1( \epsilon_1,\bar{Z},\mu) \end{align}\tag{24}\] where \[\begin{align} \Delta_1( \epsilon_1,\bar{Z},\mu)&=b^2\bar{x}_2^2 +\bar{g}_1(\bar{z}, \epsilon_1, \mu)^2+b^2\bar{\chi}(\bar{\eta}, \epsilon_1, \mu)^2,\\ \bar{\chi}(\bar{\eta}, \epsilon_1, \mu) &\equiv \chi(\bar{\eta}+\boldsymbol{\eta}^*+Nb^{-1} \epsilon_1)-\chi(\boldsymbol{\eta}^*). \end{align}\] Now let \(U_2(\bar{Z}, \epsilon_1, \epsilon_2)=U_1(\bar{Z}, \epsilon_1)+\epsilon_2^2\). The time derivative of \(U_2\equiv U_2(\bar{Z}, \epsilon_1, \epsilon_2)\) along the trajectory of \(\epsilon_2\)-subsystem with \(\hat{x}_3=\epsilon_{3}+\alpha_2\) is given by \[\begin{align} \dot{U}_2 \leq &\;\dot{U}_1+2 \epsilon_2\dot{\epsilon}_2\nonumber\\ \leq&-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2-\big(b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2\nonumber\\ &+2b\epsilon_1\epsilon_2+\Delta_1( \epsilon_1,\bar{Z},\eta)+2 \epsilon_2(\epsilon_{3}+\alpha_2-\lambda_2\hat{x}_1-\dot{\alpha}_1)\nonumber\\ \leq&-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2-\big(b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2\nonumber\\ &+2b\epsilon_1\epsilon_2+\Delta_1( \epsilon_1,\bar{Z},\eta)\nonumber\\ &+2 \epsilon_2\Big(\epsilon_{3}+\alpha_2-\lambda_2\hat{x}_1-\frac{\partial \alpha_1}{\partial \epsilon_1}\dot{\epsilon}_1-\frac{\partial \alpha_1 }{\partial \eta}\dot{\eta}-\frac{\partial \alpha_1 }{\partial k^*}\dot{k}^*\Big)\nonumber\\ \leq&-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2-\big(b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2\nonumber\\ &+2b\epsilon_1\epsilon_2+\Delta_1( \epsilon_1,\bar{Z},\eta)+2 \epsilon_2\epsilon_{3}\nonumber\\ &+2 \epsilon_2\Big(\alpha_2-\lambda_2\hat{x}_1-\frac{\partial \alpha_1 }{\partial \eta}\dot{\eta}-b\frac{\partial \alpha_1}{\partial \epsilon_1}(\underbrace{\epsilon_2+\alpha_1}_{\hat{x}_2}-\chi(\eta))\nonumber\\ &-\frac{\partial \alpha_1}{\partial \epsilon_1}\big[b\bar{\chi}(\bar{\eta},\epsilon_1,\mu)+b\bar{x}_2 +\bar{g}_1(\bar{z}, \epsilon_1, \mu)\big]-\frac{\partial \alpha_1 }{\partial k^*}\dot{k}^*\Big)\nonumber\\ \leq&-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2-\big(b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2\nonumber\\ &+\Delta_1( \epsilon_1,\bar{Z},\eta)+2 \epsilon_2\epsilon_{3}\nonumber\\ &+2 \epsilon_2\Big(\alpha_2-\epsilon_2+ \mathop{\vtop{\m@th\ialign{##\crcr \hfil\displaystyle{\epsilon_2+b\epsilon_1-\lambda_2\hat{x}_1-\frac{\partial \alpha_1 }{\partial \eta}\dot{\eta}}\hfil\crcr \noalign{\kern 3\p@\nointerlineskip} \m@th \setbox\z@\braceld \bracelu\leaders\vrule \@height\ht\z@ \@depth\z@\hfill \kern\p@\vrule \@width\p@\kern\p@\vrule \@width\p@\kern\p@\vrule \@width\p@ \crcr\noalign{\kern 3\p@}}}}\limits\nonumber\\ &\mathop{\vtop{\m@th\ialign{##\crcr \hfil\displaystyle{{}-b\frac{\partial \alpha_1}{\partial \epsilon_1}(\epsilon_2\underbrace{-k^*\rho(\epsilon_1)\epsilon_1+\chi(\eta)}_{\alpha_1}-\chi(\eta))-\frac{\partial \alpha_1 }{\partial k^*}\dot{k}^*{}}\hfil\crcr \noalign{\kern 3\p@\nointerlineskip} \m@th \setbox\z@\braceld\vrule \@width\p@\kern\p@\vrule \@width\p@\kern\p@\vrule \@width\p@\kern\p@ \leaders\vrule \@height\ht\z@ \@depth\z@\hfill\bracerd \braceld\leaders\vrule \@height\ht\z@ \@depth\z@\hfill \kern\p@\vrule \@width\p@\kern\p@\vrule \@width\p@\kern\p@\vrule \@width\p@ \crcr\noalign{\kern 3\p@}}}}\limits_{-\alpha_2 }\nonumber\\ &\mathop{\vtop{\m@th\ialign{##\crcr \hfil\displaystyle{{}+\frac{1}{2}\epsilon_2\Big(\frac{\partial \alpha_1}{\partial \epsilon_1}\Big)^2}\hfil\crcr \noalign{\kern 3\p@\nointerlineskip} \m@th \setbox\z@\braceld\vrule \@width\p@\kern\p@\vrule \@width\p@\kern\p@\vrule \@width\p@\kern\p@ \leaders\vrule \@height\ht\z@ \@depth\z@\hfill\braceru\crcr\noalign{\kern 3\p@}}}}\limits\Big)+\epsilon^2_2\nonumber\\ &+\underbrace{\big[b^2\bar{\chi}(\bar{\eta},\epsilon_1,\mu)^2+b^2\bar{x}_2^2 +\bar{g}_1(\bar{z}, \epsilon_1, \mu)^2\big]}_{\Delta_1( \epsilon_1,\bar{Z},\eta)}\nonumber\\ \leq&-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2+2\Delta_1( \epsilon_1,\bar{Z},\mu)+2 \epsilon_2\epsilon_{3}\nonumber\\ &-\big(2b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2-\epsilon_2^2.\nonumber \end{align}\] Now let \(U_i(\bar{Z}, \epsilon_1, \dots,\epsilon_i)=U_{i-1}(\bar{Z}, \epsilon_1,\dots,\epsilon_{i-1})+\epsilon_i^2\). The time derivative of \(U_i\equiv U_i(\bar{Z}, \epsilon_1,\dots,e_{i})\) along the trajectory of \(\epsilon_i\)-subsystem with \(\hat{x}_{i+1}=\epsilon_{i+1}+\alpha_i\) is given by \[\begin{align} \dot{U}_i\leq&-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2+i\Delta_1( \epsilon_1,\bar{Z},\mu)+2 \epsilon_i\epsilon_{i+1}\nonumber\\ &-\big(2b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2-\sum\nolimits_{j=2}^{i}\epsilon_{j}^2.\nonumber \end{align}\] Finally, at \(i=r\) results in \[\begin{align} \label{U1-derivative-r} U_r(\bar{Z}, \epsilon_1, \dots,\epsilon_r)=V_{1}\big( \bar{Z} \big) +\sum\nolimits_{i=1}^{r}\epsilon_i^2. \end{align}\tag{25}\] The time derivative of \(U_r\equiv U_r(\bar{Z}, \epsilon_1, \dots,\epsilon_r)\) along the trajectory of systems 19 and 21 with \({e}_{r+1}=0\) is given by \[\begin{align} \label{U1-derivative-sr} \dot{U}_r\leq&-\Delta_{Z}(\bar{Z} )\big\|\bar{Z} \big\|^2+r\Delta_{a}( \epsilon_1,\bar{Z},\bar{a},\mu)\nonumber\\ &-\big(2b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\big)\epsilon_1^2-\sum\nolimits_{j=2}^{r}\epsilon_{j}^2, \end{align}\tag{26}\] where \[\Delta_a( \epsilon_1,\bar{Z},\bar{a},\mu)=b^2\bar{x}_2^2 +\bar{g}_1^2(\bar{z}, \epsilon_1, \mu)+b^2\bar{\chi}_s^2(\bar{\eta}_e,\bar{a},\mu).\] From Properties 2 and 4, there are positive smooth functions \(\gamma_{g0}(\cdot)\) and \(\gamma_{g1}(\cdot)\), positive constants \(\phi_0\), \(\phi_1\), and \(\phi_2\) such that \[\begin{align} \Delta_a( \epsilon_1,\bar{Z},\bar{a},\mu)\leq & \;(\gamma_{g0}(\bar{Z})+\phi_1)\|\bar{Z}\|^2\nonumber\\ &+\epsilon_1^2(\gamma_{g1}(\epsilon_1)+\phi_0)+\phi_{2}\alpha_{\bar{a}}(V_{\bar{a}}). \end{align}\] Now, define a Lyapunov function by \[\begin{align} \label{U-zer-a} U(\bar{Z}, \epsilon_1, \dots,\epsilon_r, \bar{a})=U_r(\bar{Z}, \epsilon_1, \dots,\epsilon_r)+\phi_{\bar{a}} V_{\bar{a}}(\bar{a}), \end{align}\tag{27}\] where the positive constant \(\phi_{\bar{a}}\) is to be specified, and the integral Input-to-State Stable Lyapunov functions \(V_{\bar{a}}(\bar{a})\) is given in Property 3 of Lemma 3. The time derivative of \(U\equiv U(\bar{Z}, \epsilon_1, \dots,\epsilon_r, \bar{a})\) along the trajectories of 21 with the control input ?? satisfies \[\begin{align} \label{Uderivative} \dot{U}=&\;\dot{U}_r+\phi_{\bar{a}} \dot{V}_{\bar{a}}(\bar{a})\nonumber\\ \leq & -\big(\Delta_{Z}(\bar{Z} )-r\gamma_{g0}(\bar{Z})-r\phi_1-\phi_{\bar{a}}c_{ae}\big)\big\|\bar{Z} \big\|^2\nonumber\\ &-\big(2b k^*\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)-r\gamma_{g1}(\epsilon_1)-r\phi_0\nonumber\\ &-\phi_{\bar{a}}c_{ae}\big)\epsilon_1^2-\sum\nolimits_{j=2}^{r}\epsilon_{j}^2-(\phi_{\bar{a}}-r\phi_{2})\alpha_{\bar{a}}(V_{\bar{a}}). \end{align}\tag{28}\] Let the parameter and the smooth functions be \[\begin{align} \label{rho-inquality-2} \phi_{\bar{a}}& \geq r\phi_{2}+1,\notag\\ \Delta_{Z}(\bar{Z} ) &\geq r\gamma_{g0}(\bar{Z})+r\phi_1+\phi_{\bar{a}}c_{ae}+1,\notag\\ \rho(\epsilon_1) &\geq \max\{\gamma_{g1}(\epsilon_1), \hat{\gamma} \left(\epsilon_1\right), 1\},\notag\\ k^* &\geq {(3+ \hat{\gamma}^*+ r+ r\phi_0+\phi_{\bar{a}}c_{ae}+1)}/({2b}). \end{align}\tag{29}\] Hence, 28 yields the inequality \[\begin{align} \label{strlyafun}\dot{U}\leq -\big\|\bar{Z}\big\|^2-\sum\nolimits_{j=1}^{r}\epsilon_{j}^2-\alpha_{\bar{a}}(V_{\bar{a}}). \end{align}\tag{30}\] Finally, because \(U(\bar{Z}, \epsilon_1, \dots,\epsilon_r, \bar{a})\) is positive definite and radially unbounded and satisfies a strict Lyapunov function satisfying inequality 30 , it follows that the closed-loop system is uniformly asymptotically stable for all \(col(v,w,\sigma)\in \mathbb{V}\times \mathbb{W}\times \mathbb{S}\). This completes the proof. \(\Box\)

From Theorem 1, we can also use the adaptive method to estimate the gain \(k^*\). As a result, Theorem 1 can admit the following corollary.

Define the Lyapunov function \[\begin{align} V_{r}(\bar{Z}, \epsilon_1,\dots, \epsilon_r, \bar{a},\tilde{k})=\underbrace{\overbrace{V_{1}\big( \bar{Z} \big) +\sum\nolimits_{i=1}^{r}\epsilon_i^2}^{U_r(\bar{Z}, \epsilon_1, \dots,\epsilon_r)}+\phi_{\bar{a}} V_{\bar{a}}(\bar{a})}_{U(\bar{Z}, \epsilon_1, \dots,\epsilon_r,\bar{a})}+b\tilde{k}^2 \end{align}\] for the system 21 , where \(\tilde{k}=\hat{k}-k^*\) and \(U_r(\bar{Z}, \epsilon_1, \dots,\epsilon_r)\) are the same as in 25 and \(U(\bar{Z}, \epsilon_1, \dots,\epsilon_r,\bar{a})\) is the same as that in 27 . It can be verified that \(V_{r}(\bar{Z}, \epsilon_1,\dots, \epsilon_r, \bar{a},\tilde{k})\) is globally positive definite and radially unbounded. Then, the time derivative of \(V_{r}(\bar{Z}, \epsilon_1,\dots, \epsilon_r, \bar{a}, \tilde{k})\) along the trajectories of 21 with the control input ?? satisfies \[\begin{align} \label{Vrderivative} \dot{V}_{r}=&\;\underbrace{\dot{U}_r+\phi_{\bar{a}} \dot{V}_{\bar{a}}(\bar{a})}_{\dot{U}} +2b\tilde{k}\dot{\tilde{k}}\nonumber\\ \leq & -\big(\Delta_{Z}(\bar{Z} )-r\gamma_{g0}(\bar{Z})-r\phi_1-\phi_{\bar{a}}c_{ae}\big)\big\|\bar{Z} \big\|^2\nonumber\\ &-\big(2b (\underbrace{\hat{k}-k^*}_{\tilde{k}}+k^*)\rho(\epsilon_1)-3-\hat{\gamma}^* \hat{\gamma} \left(\epsilon_1\right)\nonumber\\ &-r\gamma_{g1}(\epsilon_1)-r\phi_0-\phi_{\bar{a}}c_{ae}\big)\epsilon_1^2-\sum\nolimits_{j=2}^{r}\epsilon_{j}^2\nonumber\\ &-(\phi_{\bar{a}}-r\phi_{2})\alpha_{\bar{a}}(V_{\bar{a}})+2b\tilde{k}\rho(\epsilon_1)\epsilon_1^2 \end{align}\tag{31}\] Let the parameter and the smooth functions be \[\begin{align} \label{rho-inquality-Vr} \phi_{\bar{a}}& \geq r\phi_{2}+1,\notag\\ \Delta_{Z}(\bar{Z} ) &\geq r\gamma_{g0}(\bar{Z})+r\phi_1+\phi_{\bar{a}}c_{ae}+1,\notag\\ \rho(\epsilon_1) &\geq \max\{\gamma_{g1}(\epsilon_1), \hat{\gamma} \left(\epsilon_1\right), 1\},\notag\\ k^* &\geq {(3+ \hat{\gamma}^*+ r+ r\phi_0+\phi_{\bar{a}}c_{ae}+1)}/({2b}). \end{align}\tag{32}\] Hence, 31 yields the inequality \[\begin{align} \label{V-strlyafun} \dot{V}_{r}\leq -\big\|\bar{Z}\big\|^2 - \sum_{j=1}^{r}\epsilon_{j}^2 - \alpha_{\bar{a}}(V_{\bar{a}}). \end{align}\tag{33}\] Since \(\alpha_{\bar{a}}(\cdot)\) is a class \(\mathcal{K}\) function, the term on the right-hand side is negative semi-definite. By invoking the LaSalle–Yoshizawa Theorem [33], the boundedness of \(V_r(t)\) implies that the signals \(\bar{Z}\), \(\epsilon_1\), \(\dots\), \(\epsilon_r\), \(\bar{a}\), and \(\tilde{k}\) are globally uniformly bounded. Furthermore, it follows that \[\lim_{t \to \infty} \Bigg( \big\|\bar{Z}(t)\big\|^2 + \sum_{j=1}^{r}\epsilon_{j}^2(t) + \alpha_{\bar{a}}(V_{\bar{a}}(t)) \Bigg) = 0,\] which implies \(\lim_{t \to \infty} \|\bar{Z}(t)\| = 0\), \(\lim_{t \to \infty} \epsilon_j(t) = 0\) for \(j=1,\dots,r\), and \(\lim_{t \to \infty} \alpha_{\bar{a}}(V_{\bar{a}}(t)) = 0\). Therefore, the controller ?? solves Problem 1. \(\Box\)

4 Practical Examples↩︎

4.1 Regulation of a Virtual Synchronous Generator (VSG)↩︎

To demonstrate the application of the proposed framework in modern power systems, consider a grid-connected inverter controlled as a virtual synchronous generator (VSG). The VSG control strategy forces the power electronics interface to emulate the electromechanical dynamics of a traditional synchronous machine, thereby providing inertia and damping support to the grid.

The core dynamics are governed by the virtual swing equation [34]: \[\label{vsg-dynamics} \begin{align} \dot{x}_1 &= x_2, \\ \dot{x}_2 &= \frac{1}{J}(u - D x_2 - P_e(x_1) + d(t)), \\ y &= x_1, \quad e = y - y_r, \end{align}\tag{34}\] where \(x_1 = \delta\) is the power angle (phase difference between inverter and grid), \(x_2 = \dot{\delta}\) is the frequency deviation from the nominal grid frequency, the control input \(u\) represents the mechanical power reference (\(P_\textrm{ref}\)), and \(J\) and \(D\) are the virtual inertia and virtual damping coefficients, respectively, which are programmable parameters in the inverter’s microcontroller.

The nonlinearity arises from the electrical power transfer equation \(P_e(x_1)\). For a connection to an infinite grid with voltage \(V_g\) through a line impedance \(X\), this relationship is sinusoidal: \[\begin{align} P_e(x_1) = \frac{E_g V_g}{X} \sin(x_1 + \phi), \end{align}\] where \(E_g\) is the inverter output voltage magnitude. The parameters \(V_g\), \(X\), and \(\phi\) are determined by the physical grid conditions and are typically unknown and time-varying (represented by the disturbance \(d(t)\)).

Differentiating the output \(y=x_1\) twice yields \[\begin{align} \dot{y} &= x_2, \\ \ddot{y} &= \frac{1}{J}\left( u - D x_2 - \frac{E_g V_g}{X}\sin(x_1+\phi) + d(t) \right). \end{align}\] The control input \(u\) appears in the second derivative, indicating a relative degree of \(r=2\). The control objective is to maintain grid synchronization (regulate \(x_1\) to a desired load angle \(y_r\) or track a frequency reference) despite fluctuations in grid voltage \(V_g\) and impedance \(X\). The term \(\sin(x_1)\) represents a strong trigonometric nonlinearity that the internal model must learn to ensure stable power delivery.

In the numerical simulation, the physical and virtual parameters of the VSG system are configured as follows. The virtual inertia and damping coefficient are set to \(J=0.1\) and \(D=5\), respectively. The electrical parameters characterizing the grid interconnection are chosen as

• Inverter output voltage magnitude: \(E_g = 10\) V;

• Nominal grid voltage magnitude: \(V_g = 10\) V;

• Line impedance: \(X = 10\,\Omega\);

• Initial phase offset: \(\phi = 0\).

Based on these parameters, the theoretical maximum power transfer capability is \(P_{\max} = \frac{E_g V_g}{X} = 10\). While fixed values of \(V_g\) and \(X\) are used to simulate the plant dynamics, these parameters are treated as unknown and potentially time-varying by the proposed nonparametric controller. The controller must learn the unknown sinusoidal relationship \(P_e(x_1) = 10 \sin(x_1)\) online to achieve the tracking objective.

To validate the regulation capability under standard dispatch conditions, the control objective is to regulate the power angle to a constant operating point, \[\begin{align} y_r(t) = 0.5 \, \text{rad}, \end{align}\] which simulates a fixed power dispatch command from the transmission system operator. Although the reference is constant, the steady-state control input \(u_{ss}\) required to maintain this angle is unknown a priori due to the uncertain grid impedance \(X\) and disturbances. The internal model must learn the inverse of the nonlinear power transmission characteristic to generate the correct mechanical power reference. The interaction dynamics between the grid and the VSG, as defined in 34 , combined with unknown disturbances, result in an uncertain steady-state behavior. This makes deriving an explicit solution mathematically intractable, particularly since only the system output is available for feedback. To address this, assume that the steady-state input can be described by a linear internal model of the form: \[\begin{align} \frac{d^{7}\hat{u}}{dt^{7}} + a_7 \frac{d^{6}\hat{u}}{dt^{6}} + \dots + a_2 \frac{d\hat{u}}{dt} + a_1 \hat{u} = 0, \end{align}\] where \(\boldsymbol{a} = col(a_1, a_2, \dots, a_7)\) is a vector of unknown constant coefficients that the controller must identify or adapt to.

To rigorously test the robustness of the proposed nonparametric controller against unknown periodic disturbances, the disturbance \(d(t)\) is designed as a composite signal: \[\begin{align} d(t) = 0.1 \sin(2\pi t). \end{align}\] The internal model is expected to learn and compensate for this unknown frequency component (\(1\,\text{Hz}\)) to achieve asymptotic regulation.

For the control law ?? , we can choose \(\rho(e)=1+e^2\) based on 29 to make the inequality 28 hold \(k_a=20\) is any positive number in 17 , \(\lambda_1=2\) and \(\lambda_2=2\) are chosen to make \(A\) in 4 to be Hurwitz. \(m_1=1\), \(m_2=9.5144\), \(m_3=44.7616\), \(m_4=137.7619\), \(m_5=309.4184\), \(m_6=535.9283\), \(m_7=737.6421\), \(m_8=819.2345\), \(m_9=737.6421\), \(m_{10}=535.9283\), \(m_{11}=309.4184\), \(m_{12}=137.7619\), \(m_{13}=44.7616\), and \(m_{14}=9.5144\) are chosen to make \(M\) in 8 defined in 9 to be Hurwitz. The simulation starts with the initial conditions \(x(0)=col(0.1,0)\), \(\eta(0)=\boldsymbol{0}_{14}\), \(\hat{k}(0)=2\), and \(\hat{a}(0)=\boldsymbol{0}_7\).

The results in Figs. 3 (a)–4 confirm that the VSG successfully synchronizes with the grid, with the internal model reconstructing the unknown grid interaction dynamics \(P_e(x_1)\). In Fig. 4, the solid black line represents the true, unknown grid power transfer characteristic \(P_e(\delta) = P_{\max}\sin(\delta)\). The colored scatter points represent the instantaneous estimate \(\hat{P}_e(t)\) generated by the nonparametric internal model against the actual power angle \(\delta(t)\). The convergence of the trajectory (from blue to red) onto the theoretical manifold demonstrates that the internal model has successfully learned the unknown sinusoidal nonlinearity profile without any a priori structural knowledge.

4.2 Heading Control of a Surface Vessel↩︎

Consider the robust course-keeping control problem for a surface vessel. The yaw dynamics are often described by the nonlinear Norrbin model [35], [36], which captures the essential maneuvering characteristics, including the nonlinear damping effects: \[\label{eq:ship95dynamics} \begin{align} \dot{x}_1 &= x_2, \\ \dot{x}_2 &= -\frac{1}{T} x_2 - \frac{\alpha}{T} x_2^3 + \frac{K}{T} u + d(t), \\ y &= x_1, \quad e = y - y_r, \end{align}\tag{35}\] where \(x_1 = \psi\) is the heading angle, \(x_2 = r\) is the yaw rate, and \(u = \delta\) is the rudder angle. \(T\) and \(K\) are the time constant and gain indices, respectively. The parameter \(\alpha\) represents the coefficient of the nonlinear cubic damping, which is crucial for describing the vessel’s turning behavior at high speeds. \(d(t)\) denotes the unknown environmental disturbances induced by wind, waves, and ocean currents.

Differentiating the output \(y=x_1\) twice yields \[\begin{align} \dot{y} &= x_2, \\ \ddot{y} &= -\frac{1}{T} x_2 - \frac{\alpha}{T} x_2^3 + \frac{K}{T} u + d(t). \end{align}\] The control input \(u\) appears in the second derivative, confirming a relative degree of \(r=2\). The nonlinearity \(-\frac{\alpha}{T} x_2^3\) depends on the unmeasured state \(x_2\) (assuming only heading is measured or to test the observer) and the unknown parameter \(\alpha\). The disturbance \(d(t)\) acts as a matched uncertainty. The control objective is to steer the vessel to a desired heading \(y_r\) while rejecting the wave-induced yaw moments.

In the numerical simulation, the physical parameters of the USV are configured based on a standard prototype model:

• Gain index: \(K = 0.5\) s\(^{-1}\);

• Time constant: \(T = 3\) s;

• Nonlinear damping coefficient: \(\alpha = 1\) s\(^2\).

While fixed values are used for the plant simulation, these parameters are treated as unknown by the controller. The controller must learn the inverse of the nonlinear steering dynamics online.

To validate the regulation capability under dynamic sea states, the control objective is to track a constant heading reference: \[\begin{align} y_r(t) =\frac{\pi}{4}. \end{align}\] Although the reference trajectory is defined, the steady-state control input \(u_{ss}\) required to maintain this course is unknown a priori due to the parametric uncertainties and external disturbances from the ocean environment. The interaction dynamics between the vessel and the fluid environment, combined with the disturbance \(d(t)\), result in an uncertain steady-state behavior. We assume the steady-state input can be described by a linear internal model of the form: \[\begin{align} \frac{d^{5}\hat{u}}{dt^{5}} + a_5 \frac{d^{4}\hat{u}}{dt^{4}} + a_4 \frac{d^{3}\hat{u}}{dt^{3}} + a_2 \frac{d\hat{u}}{dt} + a_1 \hat{u} = 0, \end{align}\] where \(\boldsymbol{a} = col(a_1, \dots, a_5)\) is a vector of unknown constant coefficients.

To rigorously test robustness, the disturbance \(d(t)\) is designed as a sinusoidal signal simulating regular wave encounters: \[\begin{align} d(t) = 0.5 \sin(t). \end{align}\] The internal model is expected to compensate for this unknown frequency component (\(1\,\text{rad/s}\)) to achieve asymptotic tracking.

For the control law ?? , we choose \(\rho(e)=1+e^2\) to satisfy the inequality condition. The stabilization parameter is set to \(k_a=15\) in 17 . The filter parameters \(\lambda_1=1.5\) and \(\lambda_2=1.5\) are selected to ensure that the matrix \(A\) in 4 is Hurwitz. To construct the internal model dynamics, the parameters \(m_1=10\), \(m_2=55.11\), \(m_3=151.35\), \(m_4=270.59\), \(m_5=346.41\), \(m_6=329.72\), \(m_7=234.84\), \(m_8=122.69\), \(m_9=44.52\) and \(m_{10}=9.95\) are chosen to make the matrix \(M\) in 9 Hurwitz. The simulation starts with initial conditions \(x(0)=col(1,0)\), \(\eta(0)=\boldsymbol{0}_{12}\), \(\hat{k}(0)=2\), and \(\hat{a}(0)=\boldsymbol{0}_6\).

The results in Figs. 5–4 confirm that the USV successfully tracks the desired heading reference. Fig. 4 specifically illustrates the X-Y plane trajectory, where the vessel maintains a constant surge speed \(U=10\) m/s while adjusting its heading \(\psi\) according to the control law \(\dot{X}=U\cos(\psi), \dot{Y}=U\sin(\psi)\), demonstrating effective path following despite the presence of nonlinear damping and wave disturbances.

4.3 Regulation of a Repulsive Magnetic Levitation System↩︎

To evaluate the proposed strategy for an unstable nonlinear system, consider a repulsive magnetic levitation system, in which a permanent magnet is levitated above an electromagnet coil. The objective is to regulate the levitated height \(y\) to a desired trajectory \(y_r\). The dynamics are governed by the balance between the electromagnetic repulsion force and gravity, modelled as \[\label{maglev-dynamics} \begin{align} \dot{x}_1 &= x_2,& \\ \dot{x}_2 &= \frac{C}{m_g} \frac{x_3}{x_1^2} - g, &\\ \dot{x}_3 &= -\frac{R}{L} x_3 - \frac{K_b}{L} x_2 + \frac{1}{L} u + d(t),& \\ y &=x_1, & e = y - y_r, \end{align}\tag{36}\] where \(x_1\) denotes the vertical distance (air gap) between the magnet and the coil, \(x_2\) is the vertical velocity, and \(x_3\) is the current in the coil. The control input \(u\) is the voltage applied to the coil. The system parameters include the mass of the floating magnet \(m_g\), gravity \(g\), electromagnetic force constant \(C\), coil resistance \(R\), inductance \(L\), and the back-electromotive force coefficient \(K_b\). An external disturbance \(d(t)\) represents voltage fluctuations or unmodelled electrical dynamics.

Unlike the attractive suspension model, the repulsive magnetic force is modelled as \(F_m = C x_3 / x_1^2\), assuming the current \(x_3\) is positive. Differentiating the output \(y=x_1\) successively until the input appears gives the input-output dynamics: \[\label{eq:maglev95standard} y^{(3)} = f(x,t) + \beta(x_1) u.\tag{37}\] The system has a relative degree of \(r=3\). The state-dependent high-frequency gain is derived as \(\beta(x_1) = \frac{C}{m_g L x_1^2}\), and the nonlinear drift dynamics \(f(x,t)\) is given by \[f(x,t) = \frac{C}{m_g x_1^2} \left( -\frac{R}{L} x_3 - \frac{K_b}{L} x_2 - \frac{2 x_2 x_3}{x_1} + d(t) \right).\] Both \(f(x,t)\) and \(\beta(x_1)\) contain strong nonlinearities coupled with the unmeasured states \(x_2\) and \(x_3\). Specifically, the term \(-2 x_2 x_3 / x_1\) introduces a geometric nonlinearity related to the change of magnetic force with position.

Since only the position output \(y\) is available, the proposed input-driven filter is used to recover the state information. The controller aims to maintain the magnet at a stable levitation height despite the inherent instability of the open-loop dynamics. The complex dynamics of the repulsive magnetic levitation system 36 and the unknown exosystem 3 result in an unknown steady-state behavior, making it challenging and impossible to derive an explicit solution, especially considering that only the output is available. Therefore, assume that the system \[\begin{align} \frac{d^{5}\boldsymbol{\hat{u}}}{dt^{5}}+a_1 \boldsymbol{\hat{u}}+a_2 \frac{d\boldsymbol{\hat{u}}}{dt}+a_3 \frac{d^{2}\boldsymbol{\hat{u}}}{dt^{2}}+a_4\frac{d^{3}\boldsymbol{\hat{u}}}{dt^{3}}+a_5\frac{d^{4}\boldsymbol{\hat{u}}}{dt^{4}}=0, \end{align}\] can describe the steady-state input, where \(a=col(a_1,a_2,a_3,a_4,a_5)\) is the unknown constant vector.

In the simulation, the physical parameters are chosen as \(m_g=0.1\) kg, \(g=9.8\) m/s\(^2\), \(C=1\times 10^{-4}\), \(R=2\,\Omega\), \(L=0.05\) H, and \(K_b = 0.01\) V\(\cdot\)s/m. The reference trajectory is set to \(y_r(t) = 0.5\) m (levitating at a height of \(50~cm\)). The disturbance is set as \(d(t) = 2\sin(0.5t)\). For the control law ?? , we can choose \(\rho(e)=20+20e^2\) based on 29 to make the inequality 28 negative definite, \(k_a=2\) is any positive number in 17 , \(\lambda_1=4.5\), \(\lambda_2=6.5\) and \(\lambda_3=3\) are chosen to make \(A\) in 4 to be Hurwitz. \(m_1=1\), \(m_2=6.9552\), \(m_3=23.6871\), \(m_4=51.8675\), \(m_5=80.7073\), \(m_6=93.1412\), \(m_7=80.7073\), \(m_8=51.8675\), \(m_9=23.6871\), and \(m_{10}=6.9552\) are chosen to make \(M\) in 8 defined in 9 to be Hurwitz. The initial position is \(x_1(0)=0.2\) m. The initial velocity is \(x_2(0)=0\) m/s. The initial current in the coil is \(x_3(0)=-10\). The initial conditions for the controller are set to \(\eta(0)=\boldsymbol{0}_{10}\), \(\hat{a}(0)=\boldsymbol{0}_5\), and \(\hat{k}(0)=0.5\).

The simulation results are shown in Figs. 9 (a)–9 (c). Specifically, Fig. 9 (a) demonstrates that the tracking error converges to zero, indicating that the proposed method successfully drives the magnet to the target height. Fig. 9 (b) illustrates the desired vertical distance (air gap) between the magnet and the coil, showing that the magnet is successfully levitated to the target height of 0.5 m. Finally, Fig. 9 (c) illustrates the evolution of the parameter estimates \(\hat{a}(t)\), confirming the adaptation capability of the learning algorithm.

The repulsive magnetic levitation example is representative of a broader class of electromechanical systems with non-affine nonlinear coupling between mechanical states and electrical dynamics. While its state-space model does not strictly fall into the high-order normal form 1 assumed throughout the paper, the proposed controller ?? still achieves the prescribed regulation/tracking objective in this benchmark. This observation suggests that the design may possess a degree of robustness to structural mismatches beyond the nominal model class. A rigorous extension of the analysis to cover such general nonlinear electromechanical systems—where the “next-state affine” structure is not satisfied—will be pursued in future work.

5 Conclusion↩︎

We have presented a nonparametric learning solution for the robust output regulation of nonlinear systems via output feedback. This framework effectively converts the regulation task into a robust stabilization problem for systems with integral Input-to-State Stable (iISS) inverse dynamics. Unlike traditional adaptive methods, our approach eliminates the need for explicit regressor construction and strict Lyapunov function design for parameter learning. The approach is illustrated in three numerical examples, involving a repulsive magnetic levitation system, a virtual synchronous generator, and a marine surface vessel induce unknown steady-state chattering, showing convergence of the parameter estimation error and of the tracking error to zero. It is worth noting that the repulsive magnetic levitation system in Subsection 4.3 is not directly covered by the normal form 1 , due to its non-affine structural characteristics and complex state-input coupling. Nevertheless, the proposed control scheme remains effective and meets the objective, demonstrating its broad applicability to generalized nonlinear systems. Future work will focus on establishing a unified framework that explicitly accommodates such general electromechanical structures and provides corresponding stability and performance guarantees, as well as on validating the proposed methodology in experiments, including the control of impinging-jet mixers for the manufacturing of solid lipid nanoparticle.

References↩︎

---

1. This research was supported by the U.S. Food and Drug Administration under the FDA BAA-22-00123 program, Award Number 75F40122C00200.
   Shimin Wang and Richard D. Braatz are with Massachusetts Institute of Technology, Cambridge, MA 02142, USA (bellewsm@mit.edu, braatz@mit.edu).
   Martin Guay is with Queen’s University, Kingston, ON K7L 3N6, Canada (martin.guay@queensu.ca).
   Corresponding author: Richard Braatz↩︎