# Harmonic Modeling and Control under Variable-Frequency

Maxime Grosso, Pierre Riedinger, Jamal Daafouz, Serge Pierfederici, Hicham

Janati Idrissi, Blaise Lapôtre

► To cite this version:

Maxime Grosso, Pierre Riedinger, Jamal Daafouz, Serge Pierfederici, Hicham Janati Idrissi, et al. Harmonic Modeling and Control under Variable-Frequency. 2026. hal-05535041

HAL Id: hal-05535041

https://hal.science/hal-05535041v1

Preprint submitted on 3 Mar 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Harmonic Modeling and Control under Variable-Frequency

Maxime Grosso<sup>a,c</sup>, Pierre Riedinger<sup>a</sup>, Jamal Daafouz<sup>a,d</sup>, Serge Pierfederici<sup>b</sup>,
Hicham Janati-Idrissi<sup>c</sup>, Blaise Lapôtre<sup>c</sup>

<sup>a</sup>Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France

<sup>b</sup>Université de Lorraine, CNRS, LEMTA, F-54000 Nancy, France

<sup>c</sup>Safran Electronics & Defense, F-91344 Massy, France.

<sup>d</sup>Institut Universitaire de France, F-75000 Paris, France.

## Abstract

This paper develops a harmonic-domain framework for systems with variable fundamental frequency. A variable-frequency sliding Fourier decomposition is introduced in the phase domain, together with necessary and sufficient conditions for time-domain realizability. An exact harmonic-domain differential model is derived for general nonlinear systems under variable frequency, without assumptions on the frequency variation. An explicit parameter-varying approximation is then obtained, along with a tight error bound expressed in terms of local relative frequency variation, providing a non-conservative validity criterion and clarifying the limitations of classical heuristics. A main result shows that, for linear phase-periodic systems with affine frequency dependence, stability analysis and control synthesis can be carried out without approximation and without assumptions on the frequency variation, provided the frequency evolves within a prescribed interval. As a consequence, both problems reduce to harmonic Lyapunov inequalities evaluated at the two extreme frequency values, yielding a convex LMI characterization. The framework is illustrated on a variable-speed permanent magnet synchronous motor.

Key words: Harmonic Modeling, variable-frequency, Dynamic Phasor, Parameter-Varying Systems, Toeplitz Operators.

## 1 Introduction

Periodic signals can be naturally decomposed into harmonic components, each associated with an integer multiple of a fundamental frequency. This harmonic structure provides a useful lens for analyzing and modeling systems subject to periodic excitations or exhibiting periodic steady-state behavior. Over the past decades, two main classes of approaches have emerged to address such phenomena. The first class, rooted in control theory, focuses on the synthesis of control laws ensuring asymptotic tracking or rejection of periodic signals. It includes output regulation schemes based on the internal model principle, repetitive control, and their adaptive or robust extensions [1-5]. The second class, historically developed within electrical engineering, adopts a harmonic or frequency-domain viewpoint. Seminal contributions include harmonic state-space representations [6, 7] and dynamic phasors [8-11] which have demonstrated their effectiveness in applications such as power electronics, electrical drives, and grid-connected converters [12, 13], where rich harmonic content challenges classical time-domain modeling. By representing periodic

behaviors through the evolution of harmonic components, these approaches recast time-periodic dynamics into time-invariant (LTI) models in the harmonic domain, enabling the systematic use of standard LTI tools.

Most existing harmonic modeling and control frameworks rely on the assumption of a constant fundamental frequency [14-18]. This assumption is violated in many emerging applications where the fundamental frequency varies with operating conditions or is itself a control variable. Representative examples include variable-speed motor drives, aircraft electrical systems, renewable energy generation, and autonomous power systems. In such settings, frequency variations directly affect harmonic interactions, rendering fixed-frequency harmonic models inaccurate and potentially destabilizing when used for control design. Early attempts to address variable-frequency effects introduced dynamic phasors with time-varying frequency [6, 19]. These approaches typically relied on heuristic slowly-varying assumptions, introduced nonlinear or delayed terms, and lacked rigorous guarantees regarding time-domain realizability and control-oriented model validity. More
recent contributions have clarified fundamental aspects of harmonic modeling, including the coincidence condition for time-domain realizability [18] and computational tools for harmonic Lyapunov and LMI-based analysis [20]. However, these advances remain largely restricted to the constant-frequency case, and a rigorous framework for harmonic modeling and control under time-varying frequency is still lacking.

To address this gap, we start from a sliding Fourier analysis formulated in the phase domain and introduce a variable-frequency harmonic decomposition, together with necessary and sufficient conditions for time-domain realizability. Building on this representation, we derive an exact harmonic-domain differential model that captures the true dynamics of harmonic components under frequency variation, without relying on slowly-varying frequency assumptions. This exact formulation serves as the foundation for a systematic parameter-varying approximation, together with a precise and computable error bound that quantifies its accuracy in terms of local relative frequency variation. Within the proposed approach, we further show that for a broad and practically relevant class of linear phase-periodic systems whose harmonic representation depends affinely on the instantaneous frequency, stability and state-feedback synthesis can be addressed exactly in the harmonic domain. In particular, we establish a Lyapunov stability condition that can be certified by evaluating harmonic Lyapunov inequalities at a finite number of frequency vertices, without imposing assumptions on the frequency variation. This result bridges harmonic modeling with LMI-based control design enabling stability guarantees over a prescribed frequency interval. A preliminary version of this work appeared in [21]. Here, we significantly extend that study by removing restrictive assumptions, deriving the harmonic formulation from first principles, and establishing precise conditions under which an exact and well-posed parameter-varying control framework can be obtained.

The remainder of the paper is organized as follows. Section 2 recalls the foundations of harmonic modeling under the classical constant-frequency assumption. Section 3 introduces the proposed variable-frequency sliding Fourier decomposition and establishes its fundamental properties. Section 4 derives the exact harmonic-domain dynamical model governing the evolution of variable-frequency harmonic components. Section 5 develops a parameter-varying approximation together with a quantitative and non-conservative validity criterion. Section 6 establishes exact Lyapunov stability and control synthesis results for phase-periodic systems with affine frequency dependence. Finally, Section 7 illustrates the proposed framework on a variable-speed permanent magnet synchronous motor.

**Notations:** The transpose of a matrix $A$ is denoted $A^T$ and $A^*$ denotes the complex conjugate transpose $A^* = \bar{A}^T$. The $n$-dimensional identity matrix is denoted

$Id_n$. The infinite identity matrix is denoted $\mathcal{I}$. $C^n$ denotes the space of absolutely continuous functions, $L^p([a, b], \mathbb{C}^n)$ (resp. $\ell^p(\mathbb{C}^n)$) denotes the Lebesgue spaces of $p$-integrable functions on $[a, b]$ with values in $\mathbb{C}^n$ (resp. $p$-summable sequences of $\mathbb{C}^n$) for $1 \le p \le \infty$. $L_{\text{loc}}^p$ is the set of locally $p$-integrable functions. Finally, $\otimes$ denotes the Kronecker product.

## 2 Fixed-Frequency Harmonic Modeling

Before introducing the proposed framework for variable-frequency systems, we briefly recall the fundamental principles of harmonic modeling under the classical assumption of a constant fundamental frequency [18].

**Definition 1** For a time-signal $x \in L_{\text{loc}}^2(\mathbb{R}, \mathbb{C})$ and a fixed period $T_0 = 2\pi/\omega_0$, its sliding Fourier decomposition (SFD) is the mapping $\mathcal{F}_{T_0}: x \mapsto X \in C^a(\mathbb{R}, \ell^2(\mathbb{C}^n))$ where for every $t$, the components of the sequence $X(t) = (X_k(t))_{k \in \mathbb{Z}}$, are given by:

$$X_k(t) = \frac{1}{T_0} \int_{t-T_0}^{t} x(\tau)e^{-jk\omega_0\tau} d\tau. \quad (1)$$

$X_k$ are called the $k$-th dynamic phasors. For a time vector-valued function $x \in L_{\text{loc}}^2(\mathbb{R}, \mathbb{C}^n)$ (or matrix-valued), this extends to

$$X = \mathcal{F}(x) = (\mathcal{F}(x_1), \dots, \mathcal{F}(x_n)). \quad (2)$$

Conversely, the arrival space, denoted by $\mathcal{H}_{T_0}$ and commonly referred to as the harmonic domain [18], is characterized as follows.

**Theorem 1** For a given $X \in C^a(\mathbb{R}, \ell^2(\mathbb{C}^n))$, there exists $x \in L_{\text{loc}}^2(\mathbb{R}, \mathbb{C}^n)$ such that $X = \mathcal{F}_{T_0}(x)$ if and only if

$$\dot{X}_k(t) = \dot{X}_0(t)e^{-jk\omega_0 t}, \quad \text{a.e. for all } k \in \mathbb{Z}. \quad (3)$$

This condition defines a subspace of $C^a(\mathbb{R}, \ell^2(\mathbb{C}^n))$ called the harmonic domain $\mathcal{H}_{T_0}$. The inverse mapping $\mathcal{F}_{T_0}^{-1}: \mathcal{H}_{T_0} \to L_{\text{loc}}^2$ is given for a continuous signal $x$ by:

$$x(t) = \sum_{k \in \mathbb{Z}} X_k(t) e^{jk\omega_0 t} + \frac{T_0}{2} \dot{X}_0(t). \quad (4)$$

This theorem is commonly referred to as the coincidence condition. It characterizes the compatibility between the time evolution of the harmonic coefficients and an underlying time-domain signal, and guarantees that the harmonic representation admits a unique time-domain realization. We now recall the following result [18].

**Theorem 2** Consider a general dynamical system

$$\dot{x}(t) = f(t, x(t)) \quad (5)$$
If $x$ is a Carathéodory solution of (5) then, for any $T_0 > 0$, $X = \mathcal{F}_{T_0}(x) \in \mathcal{H}_{T_0}$ and is a solution of:

$$ \dot{X}(t) = \mathcal{F}_{T_0}(f(t, x(t))) - \omega_0 \mathcal{N} X(t), \quad (6) $$

where $\mathcal{N} = \text{Id}_n \otimes \text{diag}(\mathbb{Z})$. Conversely, if $X$ is a solution of (6) that belongs to $\mathcal{H}_{T_0}$, its representative $x(t)$ is a solution of (5).

This theorem is particularly useful for $T_0$-periodic systems (i.e., $f(t, x) = f(t + T_0, x)$), since in this case the harmonic representation (6) reduces to a time-invariant system. A key advantage is that $\mathcal{F}_{T_0}(f(t, x(t)))$ can then be computed explicitly when $f$ is polynomial in $x$ with periodic coefficients, using the following property [18].

**Property 1** Let $x \in L_{\text{loc}}^2$ and let $A$ and $B$ be $L^\infty$-matrix valued functions. Then, the following properties hold:

$$ \mathcal{F}_{T_0}(Ax) = \mathcal{T}_{T_0}(A) \mathcal{F}_{T_0}(x) = AX, \quad (7a) $$

$$ \mathcal{T}_{T_0}(AB) = \mathcal{T}_{T_0}(A) \mathcal{T}_{T_0}(B) = AB, \quad (7b) $$

where $\mathcal{T}_{T_0}$ is the Toeplitz transformation. For a scalar function $a \in L^\infty$ with Fourier coefficients $(a_k)_{k \in \mathbb{Z}}$, the corresponding infinite-dimensional Toeplitz operator, bounded on $\ell^2$, is defined as:

$$ \mathcal{T}_{T_0}(a) = \begin{bmatrix} \cdots & \vdots & \vdots & \vdots \\ \cdots & a_0 & a_{-1} & a_{-2} \\ \cdots & a_1 & a_0 & a_{-1} \\ \cdots & a_2 & a_1 & a_0 \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix}. \quad (8) $$

For a matrix-valued function $A = (a_{ij})_{1 \le i \le n, 1 \le j \le m}$, the Toeplitz-block operator is $\mathcal{A} = (\mathcal{T}_{T_0}(a_{ij}))_{1 \le i \le n, 1 \le j \le m}$. Moreover, the $\ell^2$-norm of $\mathcal{A}$ satisfies $\|\mathcal{A}\|_{\ell^2} = \|\mathcal{A}\|_{L^\infty}$.

For illustration, consider an LTP system of the form

$$ \dot{x} = A(t)x + B(t)u \quad (9) $$

where $A(t)$ and $B(t)$ are $T_0$-periodic, $L^\infty$-matrix valued functions. The harmonic representation (6) reduces to an infinite-dimensional LTI system:

$$ \dot{X}(t) = (\mathcal{A} - \omega_0 \mathcal{N})X(t) + \mathcal{B}U(t). \quad (10) $$

This LTP-to-LTI transformation enables the application of classical LTI analysis and synthesis tools to periodic systems [22, 23]. The challenge is to extend this powerful and rigorous framework to systems in which the fundamental frequency, $\omega_0$, is no longer constant.

### 3 The variable-frequency SFD

We now relax the constant-frequency assumption and consider a time-varying fundamental frequency, denoted $\omega(t)$. The central idea is to replace the fixed period $T_0$ with a time-varying “pseudo-period” $T(t)$, following the approach originally proposed in [6, 19]. In the following, we distinguish between time-domain variables $t$, phase-domain variables $\phi$, and the phase function $\theta(t)$ that maps time to phase.

**Definition 2** Let $\omega \in C^0(\mathbb{R}, \mathbb{R}^+)$ be a strictly positive, continuous function. The associated phase function $\theta$ is a $C^1$ diffeomorphism defined by

$$ \theta(t) \triangleq \theta(0) + \int_{0}^{t} \omega(s) \, ds. \quad (11) $$

Its inverse, called the phase-to-time mapping, is denoted by $p(\phi) = \theta^{-1}(\phi)$ and satisfies $\frac{d}{d\phi}p(\phi) = \frac{1}{\omega(p(\phi))}$.

This diffeomorphism provides a natural way to express any time-domain signal in terms of its phase.

**Definition 3** Let $x \in L_{\text{loc}}^2(\mathbb{R}, \mathbb{C}^n)$. Its phase-domain representation is the mapping

$$ \phi \mapsto x_{\theta}(\phi) \triangleq x(p(\phi)), $$

so that for $\phi = \theta(t)$, we have: $x(t) = x_{\theta}(\theta(t))$.

Importantly, this transformation is not a mere time reparameterization: it enables the definition of a fixed-length sliding Fourier window in the phase domain, which is essential for extending harmonic analysis to variable-frequency signals. The definition of the sliding Fourier decomposition for variable-frequency signals is recalled below [6, Appendix A].

**Definition 4** Let $\theta$ be a phase function as defined in (11). The SFD of a time-domain signal $x$ is denoted by $\mathcal{F}_{\theta}$ and defined by $t \mapsto X(t) := X_{\theta}(\theta(t))$, where $X_{\theta}(\varphi) = \mathcal{F}_{2\pi}(x_{\theta})(\varphi)$ and $\mathcal{F}_{2\pi}$ is the standard SFD, introduced in Definition 1, over a window of length $2\pi$. The resulting $\tilde{X}_{\theta}(\varphi)$ is called the phase-domain phasor.

From this definition, we have the following property.

**Property 2** The $k$-th phasor of $x$ is given by

$$ X_k(t) = \frac{1}{2\pi} \int_{\theta(t)-2\pi}^{\theta(t)} x_{\theta}(\phi) e^{-jk\phi} d\phi \quad (12) $$

$$ = \frac{1}{2\pi} \int_{t-T(t)}^{t} x(\tau) e^{-jk\theta(\tau)} \omega(\tau) \, d\tau, \quad (13) $$

where the time-varying pseudo-period $T(t)$ is defined as the unique positive solution to

$$ \theta(t) - \theta(t - T(t)) = 2\pi, \quad \forall t \in \mathbb{R}. \quad (14) $$
Moreover, if $\omega$ is continuous ($\omega \in C^0$), then $T(t)$ is a $C^1$ function satisfying

$$ \dot{T}(t) = 1 - \frac{\omega(t)}{\omega(t - T(t))}. \qquad (15) $$

**PROOF.** Applying Definition 1 to $\mathcal{F}_{2\pi}(x_\theta)(\varphi)$ yields

$$ X_k(\varphi) = \frac{1}{2\pi} \int_{\varphi-2\pi}^{\varphi} x_\theta(\phi) e^{-jk\phi} d\phi, $$

noting that in (1) we have $\omega_0 = 1$ when $T_0 = 2\pi$. Equation (12) follows directly by substituting $\varphi = \theta(t)$. Equations (13) and (14) are obtained by performing the change of variables $\phi = \theta(\tau)$. Finally, differentiating (14) with respect to $t$ concludes the proof. □

**Remark 1** When $\omega$ is constant, (13) reduces to the classical fixed-frequency Fourier decomposition (1). If $\omega(t)$ is phase-periodic, the right-hand side of (15) vanishes, implying that the pseudo-period $T(t)$ is constant; denote this constant by $T_c$. Differentiating (14) then shows that $\omega(t)$ is also $T_c$-periodic. Since $T_c$ generally differs from $2\pi/\omega(t)$, the constant phasor associated with $\omega(t)$ through (13) should not be confused with the constant phasor obtained from the classical Fourier decomposition (1) with $T_0 = T_c$. These quantities are defined with respect to different referentials, namely the phase and the time domains, and therefore represent fundamentally distinct objects.

The following theorem generalizes the classical coincidence condition of Theorem 1 to the variable-frequency setting and establishes a necessary and sufficient condition for time-domain realizability.

**Theorem 3** Let $X \in C^a(\mathbb{R}, \ell^2(\mathbb{C}^n))$. There exists a signal $x \in L_{\text{loc}}^2(\mathbb{R}, \mathbb{C}^n)$ such that $X = \mathcal{F}_\theta(x)$ if and only if

$$ \dot{X}_k(t) = \dot{X}_0(t) e^{-jk\theta(t)}, \quad \text{for a.e. } t \text{ and all } k \in \mathbb{Z}. \quad (16) $$

The subspace of $C^a(\mathbb{R}, \ell^2(\mathbb{C}^n))$ satisfying this condition is called the variable-frequency harmonic domain and denoted $\mathcal{H}_\theta$. The inverse mapping $\mathcal{F}_\theta^{-1}: \mathcal{H}_\theta \to L_{\text{loc}}^2$ is given, for a continuous signal $x$, by

$$ x(t) = \sum_{k \in \mathbb{Z}} X_k(t) e^{jk\theta(t)} + \frac{\pi}{\omega(t)} \dot{X}_0(t). \quad (17) $$

**PROOF.** By definition, for every $t$ we have $X(t) := X_\theta(\theta(t))$. Setting $\varphi = \theta(t)$ gives $t = p(\varphi)$, so that $X(p(\varphi)) = X_\theta(\varphi)$, which expresses the phasor sequence as a function of the phase $\varphi$. Since $X_\theta(\varphi) = \mathcal{F}_{2\pi}(x_\theta)(\varphi)$, its phasor-derivative with respect to $\varphi$ satisfies, for every $k \in \mathbb{Z}$ (see Theorem 1):

$$ \frac{d}{d\varphi} X_{\theta,k}(\varphi) = \frac{d}{d\varphi} X_{\theta,0}(\varphi) e^{-jk\varphi}. \quad (18) $$

On the other hand, the chain rule gives, for every $k \in \mathbb{Z}$,

$$ \dot{X}_k(t) = \omega(t) \frac{d}{d\varphi} X_{\theta,k}(\varphi) \Big|_{\varphi=\theta(t)}. $$

Combining this with (18) immediately yields Eq. (16). Similarly, for a continuous phase-domain signal $x_\theta$, we have (see (4) with $T_0 = 2\pi$ and $\omega_0 = 1$):

$$ x_{\theta}(\varphi) = \sum_{k \in \mathbb{Z}} X_{\theta,k}(\varphi) e^{jk\varphi} + \pi \frac{d}{d\varphi} X_{\theta,0}(\varphi), $$

which directly yields the formula (17). □

## 4 Exact Variable-Frequency Harmonic Model

Building on the variable-frequency harmonic representation introduced in the previous section, we now derive an exact and compact harmonic-domain differential model describing the evolution of variable-frequency phasors for general nonlinear, time-varying systems. This result provides a rigorous foundation for both approximation and control synthesis developed in the sequel. Consider a general nonlinear time-varying system

$$ \dot{x}(t) = f(t, x(t)), \quad (19) $$

where $x(t) \in \mathbb{R}^n$ and $f$ satisfies the Carathéodory conditions ensuring existence of solutions. Using the phase-domain representation $x_\theta(\phi) = x(p(\phi))$, system (19) can be equivalently expressed in the phase domain as

$$ \frac{d}{d\phi} x_{\theta}(\phi) = \frac{f(p(\phi), x_{\theta}(\phi))}{\omega(p(\phi))}. \quad (20) $$

Applying Theorem 2 to the phase-domain system yields the infinite-dimensional phasor dynamics

$$ \frac{d}{d\phi} X_{\theta}(\phi) = \mathcal{F}_{2\pi} \left( \frac{f_{\theta}(\cdot, x_{\theta}(\cdot))}{\omega_{\theta}(\cdot)} \right) (\phi) - \mathcal{N} X_{\theta}(\phi), \quad (21) $$

where $f_\theta(\phi, x_\theta(\phi)) = f(p(\phi), x_\theta(\phi))$ and $\omega_\theta(\phi) = \omega(p(\phi))$. We now translate this phase-domain formulation back to the time domain, yielding an exact variable-frequency harmonic model.

**Theorem 4** The exact dynamics of the variable-frequency phasor $X(t) = \mathcal{F}_\theta(x)(t)$ associated with (19) are given by

$$ \dot{X}(t) = \mathcal{G}(\omega(t)) \mathcal{F}_{\theta}(f(t, x(t))) - \omega(t) \mathcal{N} X(t), \quad (22) $$

where $\mathcal{G}(\omega(t)) = \text{Id}_n \otimes \omega(t) \mathcal{T}_\theta^{-1}(\omega)(t)$, and $\mathcal{N} = \text{Id}_n \otimes \text{diag}_{k \in \mathbb{Z}}(jk)$.
**PROOF.** By definition, $X(t) = X_\theta(\theta(t))$. Applying the chain rule yields $\dot{X}(t) = \omega(t) \frac{d}{d\phi} X_\theta(\phi)|_{\phi=\theta(t)}$. Substituting (21) gives

$$ \dot{X}(t) = \omega(t) \mathcal{F}_{2\pi} \left( \frac{f_\theta}{\omega_\theta} \right) (\theta(t)) - \omega(t) \mathcal{N} X(t). $$

Since $\mathcal{F}_{2\pi} \left( \frac{f_\theta}{\omega_\theta} \right) (\theta(t)) = \mathcal{F}_\theta \left( \frac{f(t,x(t))}{\omega(t)} \right)$, Property 1 yields

$$ \mathcal{F}_{\theta} \left( \frac{f}{\omega} \right) = \mathcal{T}_{\theta} (\mathrm{Id}_n / \omega) \mathcal{F}_{\theta} (f) = (\mathrm{Id}_n \otimes \mathcal{T}_{\theta}^{-1}(\omega)) \mathcal{F}_{\theta} (f). $$

The result follows by linearity and noting that $\omega(t)\mathcal{N}X_\theta(\theta(t)) = \omega(t)\mathcal{N}X(t)$. $\square$

When $\omega(t) = \omega_0$, $\mathcal{T}_\theta(\omega) = \omega_0\mathcal{I}$ and (22) reduces exactly to the classical fixed-frequency harmonic model. We conclude this section by relating (22) to existing exact formulations of variable-frequency dynamic phasors. In [6, Appendix A], the evolution of the phasors are given by

$$ \dot{X}(t) = x(t - T(t)\omega(t - T(t)))\dot{T}(t))e^{-\mathcal{N}\theta(t)} + \mathcal{F}_{\theta}(f(t, x(t))) + (\mathcal{T}_{\theta}(\frac{\dot{\omega}}{\omega}) - \mathcal{N}\mathcal{T}_{\theta}(\omega))X. \quad (23) $$

Although (23) is equivalent to (22) under suitable regularity assumptions, it involves explicit delays, nonlinear couplings, and dependence on $\dot{\omega}(t)$. In contrast, (22) concentrates all variable-frequency effects into the single operator $\mathcal{G}(\omega(t))$, avoids explicit delays and frequency derivatives, and requires only $C^1$ regularity of the phase $\theta(t)$, making it particularly amenable to control design.

# 5 Parameter-Varying (PV) Approximation

While the exact model of Theorem 4 provides a complete description of variable-frequency harmonic dynamics, its direct use in control synthesis is generally impractical due to the time-varying operator $\mathcal{G}(\omega(t))$, which depends non-linearly on the frequency signal $\omega$. This section derives a tractable parameter-varying approximation and establishes a precise and computable bound on the resulting modeling error.

## 5.1 An approximated PV harmonic model

Starting from the exact model (22), the dynamics can be decomposed as

$$ \dot{X}(t) = \underbrace{(\mathcal{F}_{\theta}(f(t,x)) - \omega(t)\mathcal{N}X(t))}_{\text{Nominal part}} + \mathcal{E}(t), \quad (24) $$

where the error term $\mathcal{E}(t)$ is given by

$$ \mathcal{E}(t) = (\mathrm{Id}_n \otimes \Delta_{\omega}(t)) \mathcal{F}_{\theta}(f(t, x)), \quad (25) $$

with $\Delta_\omega(t) = \omega(t)\mathcal{T}_\theta^{-1}(\omega)(t) - \mathcal{I}$. The nominal part of (24) has the structure of a parameter-varying system whereas $\mathcal{E}(t)$ captures the deviation with respect to the exact model (22). In order to quantify this deviation, we introduce the following notion.

**Definition 5** For a given time $t$ and $\tau \in [t - T(t), t]$, the relative variation of the frequency $\omega$ at time $\tau$ with respect to its value at time $t$ is defined as

$$ \delta_{\omega}(\tau, t) \triangleq \frac{\omega(t) - \omega(\tau)}{\omega(\tau)}. \quad (26) $$

The following theorem establishes a direct link between the Toeplitz operator $\Delta_\omega(t)$ and the relative frequency variation over the sliding pseudo-period, and provides a simple characterization of its induced $\ell^2$-norm.

**Theorem 5** For every $t$,

$$ \Delta_{\omega}(t) = \mathcal{T}_{\theta}(\delta_{\omega}(\cdot, t))(t), \quad (27) $$

and its induced $\ell^2$-operator norm satisfies

$$ \epsilon(t) := \|\Delta_{\omega}(t)\|_{\ell^2} = \sup_{\tau \in [t-T(t), t]} \left| \frac{\omega(t) - \omega(\tau)}{\omega(\tau)} \right|. \quad (28) $$

**PROOF.** Using the linearity of the Toeplitz operator and as $\omega(t)$ is constant with respect to the integration variable $\tau$ in the Fourier decomposition, we have:

$$ \begin{aligned} \mathcal{T}_{\theta}(\delta_{\omega}(\cdot, t))(t) &= \mathcal{T}_{\theta} \left( \frac{\omega(t)}{\omega(\cdot)} - 1 \right) (t) \\ &= \omega(t) \mathcal{T}_{\theta}^{-1}(\omega)(t) - \mathcal{I} = \Delta_{\omega}(t). \end{aligned} $$

Finally, (28) follows from the relation (see Property 1 or [24]): $\Delta_\omega(t)\|_{\ell^2} = \|\delta_\omega(\cdot, t)\|_{L^\infty}$. $\square$

The quantity $\epsilon(t)$, given by (28), measures the maximal relative frequency variation over one pseudo-period. When $\epsilon(t)$ is small, the system locally behaves as if the frequency were constant. Motivated by this observation, the following simplified model is obtained from the exact model (22) by applying the approximation

$$ \omega(t) \mathcal{T}_{\theta}^{-1}(\omega)(t) \approx \mathcal{I}. \quad (29) $$

**Definition 6** Under the condition $\epsilon(t) < \epsilon$ for some positive $\epsilon$, the approximated PV harmonic model is given by

$$ \dot{X}(t) \approx \mathcal{F}_{\theta}(f(t, x)) - \omega(t)\mathcal{N}X(t), \quad (30) $$

with $\omega(t)$ treated as a time-varying parameter.
**Theorem 6** The relative modeling error induced by the PV approximation satisfies

$$ \frac{\|\mathcal{E}(t)\|_{\ell^2}}{\|\mathcal{F}_{\theta}(f(t, x))\|_{\ell^2}} \le \epsilon(t), \qquad (31) $$

with $\epsilon(t) = 0$ in the constant-frequency case.

**PROOF.** The result follows directly from (25) and $\|\mathrm{Id}_n \otimes \Delta_{\omega}(t)\|_{\ell^2} = \|\Delta_{\omega}(t)\|_{\ell^2}$. $\square$

## 5.2 Analysis and comparison with classical heuristics

We now analyze the behavior of the error measure $\epsilon(t)$ introduced in Section 5 and clarify its relationship with the heuristic assumptions commonly used to justify PV approximations in the literature. Our objective is to show that $\epsilon(t)$ provides a precise and non-conservative quantification of approximation accuracy. Existing justifications of PV approximations of model (23) typically rely on an informal “slowly varying $\omega$” assumption [6,19], leading to:

(1) $\dot{T}(t) \approx 0,$

(2) $\omega(t) \approx \omega(t - T(t)),$

(3) $\dot{\omega}(t)/\omega(t) \approx 0.$

The precise implications of these assumptions for the modeling error are rarely made explicit. From (15), the condition $\dot{T}(t) \approx 0$ implies $\omega(t) \approx \omega(t - T(t))$, which in turn corresponds to $\delta_{\omega}(t - T(t), t) \approx 0$. This reasoning, however, only constrains the frequency mismatch at the boundary of the pseudo-period. It does not control the frequency variations occurring inside the interval $[t - T(t), t]$. In contrast, the proposed quantity $\epsilon(t)$ uniformly bounds the relative frequency variation over the entire pseudo-period. As such, $\epsilon(t) \approx 0$ constitutes a non-conservative condition for the validity of the PV approximation.

To make this distinction explicit, consider a linearly increasing frequency profile $\omega(t) = at + \omega_0$, with $a = \dot{\omega} > 0$. In this case, the error admits the closed-form expression

$$ \epsilon(t) = \left(1 - \frac{4\pi a}{\omega(t)^2}\right)^{-\frac{1}{2}} - 1 \approx \frac{2\pi\omega\dot{\omega}(t)}{\omega(t)^2}. \qquad (32) $$

This expression shows that $\epsilon(t)$ decreases monotonically with $\omega(t)$. Hence, for a constant acceleration, the PV approximation becomes increasingly accurate at higher frequencies.

For a prescribed admissible modeling error $\epsilon$, solving $\epsilon(t) = \epsilon$ in (32) yields the maximum allowable frequency rate

$$ a_{\max} = \frac{\omega(t)^2}{4\pi} \left(1 - \frac{1}{(1 + \epsilon)^2}\right). \qquad (33) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. $\epsilon(t)$ provides an accurate measure of the modeling error, whereas the heuristic $|\dot{\omega}/\omega|$ and $\dot{\omega}$ are poor indicators.</figcaption></figure>

This relation makes explicit how the admissible rate of frequency variation scales with $\omega(t)$, a dependence that is not captured by classical heuristics. Enforcing $\epsilon(t) = \epsilon$ for all $t$ leads to the differential equation

$$ \dot{\omega}(t) = K\omega(t)^2, \quad K = \frac{1}{4\pi} \left(1 - \frac{1}{(1 + \epsilon)^2}\right), \qquad (34) $$

whose solution $\omega(t) = (\omega_0^{-1} - Kt)^{-1}$ exhibits finite-time blow-up. In this case, $\dot{\omega}(t)/\omega(t) = K\omega(t)$, which grows unboundedly with $\omega(t)$, even though $\epsilon(t)$ is held constant by construction. This example demonstrates that the heuristic condition $\dot{\omega}/\omega \approx 0$ does not reliably characterize model validity. Figure 1 provides simulations with more complex frequency profiles to confirm this observation: both $|\dot{\omega}|$ and $|\omega/\omega|$ may take large values while $\epsilon(t)$ decreases, indicating improved approximation accuracy.

Overall, these results demonstrate that $\epsilon(t)$ provides a rigorous and non-conservative characterization of the small local relative frequency variation assumption underlying the PV harmonic framework. Unlike classical heuristics, it directly quantifies the modeling error induced by frequency variations over the pseudo-period.

# 6 Affine Linear phase-periodic Systems

This section shows that, for a practically relevant class of phase-periodic systems, variable-frequency effects can be handled exactly in the harmonic domain. In contrast to the approximation-based approach of Section 5, we focus on linear phase-periodic systems whose harmonic representation depends affinely on the instantaneous frequency. Such a structure naturally arises in applications such as electrical drives, aircraft power systems, and modulated converters. For this class, we establish an exact Lyapunov stability condition for the autonomous dynamics, which then serves as the foundation for state-feedback controller synthesis via an appropriate system augmentation.
## 6.1 Exact Lyapunov stability condition

We consider linear systems whose coefficients are periodic in the phase and depend affinely on the instantaneous frequency. These systems are described by the following affinely frequency-modulated phase-periodic (AFM-LPP) model.

**Definition 7** AFM-LPP systems are described by

$$ \dot{x}(t) = (A_0(\theta(t)) + \omega(t)A_1(\theta(t)))x(t) + (B_0(\theta(t)) + \omega(t)B_1(\theta(t)))u(t), \quad (35) $$

where $A_i(\cdot)$ and $B_i(\cdot)$ are $L^\infty$ matrix-valued functions that are $2\pi$-periodic in $\theta$, and $\theta$ denotes a phase function as defined in Definition 2.

When $A_1 \equiv 0$ and $B_1 \equiv 0$, the model reduces to a classical linear phase-periodic system. The AFM-LPP structure therefore strictly generalizes the fixed-frequency case while preserving a form that admits exact harmonic-domain analysis.

Applying the variable-frequency SFD to (35) and invoking Theorem 4 yields the exact harmonic representation

$$ \dot{X}(t) = \mathcal{G}(\omega(t))(A_0 X(t) + B_0 U(t)) + \omega(t)(A_1 X(t) + B_1 U(t) - \mathcal{N}X(t)), \quad (36) $$

where the Toeplitz block operators $A_i = \mathcal{T}_\theta(A_i)$ and $B_i = \mathcal{T}_\theta(B_i)$, $i = 0, 1$, are constant and bounded on $\ell^2$.

We first focus on the autonomous dynamics

$$ \dot{x}(t) = (A_0(\theta(t)) + \omega(t)A_1(\theta(t)))x(t), \quad (37) $$

with associated harmonic model

$$ \dot{X}(t) = (\mathcal{G}(\omega(t))A_0 + \omega(t)(A_1 - \mathcal{N}))X(t). \quad (38) $$

The following theorem provides an exact Lyapunov stability condition for these dynamics. Its key feature is that stability under arbitrary admissible frequency trajectories can be certified by enforcing Lyapunov inequalities at a finite number of frequency vertices only. This is made possible by two structural properties induced by affine frequency dependence: convexity with respect to $\omega$, and commutativity of Toeplitz operators generated by scalar phase functions.

**Theorem 7** Consider the constant Toeplitz block operators $A_0$ and $A_1$, bounded on $\ell^2$. Assume that there exists a Toeplitz block operator $\mathcal{P} = \mathcal{P}^* > 0$, bounded on $\ell^2$, such that

$$ \mathcal{P}(A_0 + \omega_i(A_1 - \mathcal{N})) + (A_0 + \omega_i(A_1 - \mathcal{N}))^* \mathcal{P} < 0, \quad (39) $$

for all $\omega_i \in \{\omega_{\min}, \omega_{\max}\}$. Let $\omega(t)$ be any time-varying frequency of the form

$$ \omega(t) = \alpha(t) \omega_{\min} + (1 - \alpha(t)) \omega_{\max}, \quad \alpha(t) \in (0, 1). $$

Then the harmonic Lyapunov inequality

$$ \mathcal{P}(\mathcal{G}(\omega(t))A_0 + \omega(t)(A_1 - \mathcal{N})) + (\cdot)^* < 0 \quad (40) $$

holds for all $t$. Moreover, since $\mathcal{P}$, $A_0$, and $A_1$ admit $2\pi$-periodic and $L^\infty$ representatives, (40) is equivalent to the following differential Lyapunov inequality, which holds almost everywhere:

$$ \dot{P}(\theta) + P(\theta)(A_0(\theta) + \omega(t)A_1(\theta)) + (\cdot)^\top < 0, \quad (41) $$

with $P(\theta) = P(\theta)^\top > 0$ and $\dot{\theta}(t) = \omega(t)$.

**PROOF.** Consider the Toeplitz transformation:

$$ \mathcal{T}_\theta(\omega) = \mathcal{T}_\theta(\alpha) \omega_{\min} + (\mathcal{I} - \mathcal{T}_\theta(\alpha)) \omega_{\max}. $$

Recall that the product of two positive definite operators remains positive definite if these operators commute. Since $\alpha(t) \in (0, 1)$, both operators $\mathcal{T}_\theta(\alpha)$ and $\mathcal{I} - \mathcal{T}_\theta(\alpha)$ are positive definite. Moreover, as $\alpha$ is a scalar function, $\mathcal{T}_\theta(\alpha)$ commutes with any Toeplitz operator and $\mathcal{T}_\theta(\alpha) = \mathcal{T}_\theta(\alpha)^*$. Using these properties, we multiply (39) on the left by $\text{Id}_n \otimes \mathcal{T}_\theta(\alpha)$ for $\omega_i = \omega_{\min}$, and by $\text{Id}_n \otimes (\mathcal{I} - \mathcal{T}_\theta(\alpha))$ for $\omega_i = \omega_{\max}$. Summing the resulting inequalities preserves negative definiteness and yields

$$ \mathcal{P}(A_0 + (\text{Id}_n \otimes \mathcal{T}_\theta(\omega))(A_1 - \mathcal{N})) + (\cdot)^* < 0. \quad (42) $$

Multiplying (42) on the left by

$$ \mathcal{G}(\omega(t)) = \omega(t)(\text{Id}_n \otimes \mathcal{T}_\theta(\omega))^{-1}, $$

which is positive definite, self-adjoint and commutes with the involved operators, directly gives (40). Finally, since $A_0$, $A_1$, and $\mathcal{P}$ admit $2\pi$-periodic representatives, multiplying (42) on the left by $(\text{Id}_n \otimes \mathcal{T}_\theta(\omega))^{-1}$ leads to

$$ \mathcal{P}\left[\mathcal{T}_{2\pi}\left(\frac{A_0(\phi)}{\omega(p(\phi))} + A_1(\phi)\right) - \mathcal{N}\right] + (\cdot)^* < 0. $$

This Lyapunov inequality admits an equivalent differential formulation in the phase domain [20]:

$$ P(\phi)\left(\frac{A_0(\phi)}{\omega(p(\phi))} + A_1(\phi)\right) + (\cdot)^\top + \frac{dP(\phi)}{d\phi} < 0. $$

The inequality (41) follows by multiplication by $\omega(t)$ and application of the chain rule $\dot{P} = \omega \frac{dP}{d\phi}$. $\square$
**Corollary 1** Let $\mathcal{P}$ and $\mathcal{P}(\theta)$ satisfy the conditions of Theorem 7. Then the Lyapunov functions $\mathcal{V}(X(t)) = X(t)^*\mathcal{P}X(t)$ and $v(t, x) = x(t)^\top\mathcal{P}(\theta(t))x(t)$ guarantee asymptotic stability of the systems (38) and (37), respectively, for any frequency function $\omega(t)$ taking values in $[\omega_{\min}, \omega_{\max}]$.

## 6.2 Harmonic control design

We now extend the stability result of Theorem 7 to control design. To this end, we introduce an augmented system representation that incorporates the controller dynamics while preserving the AFM-LPP structure. As a result, the vertex-based Lyapunov conditions of Theorem 7 remain applicable to the closed-loop system. The guiding idea is to incorporate integral action and harmonic cancellation systematically while remaining fully compatible with the Toeplitz operator structure. To this end, we introduce a forwarding variable $z$, whose dynamics implement an internal-model principle:

$$ \dot{z}(t) = \omega(t)(J(\theta(t))z(t) + L(\theta(t))C(\theta(t))x(t)) $$

where $J, L, C$ are given $2\pi$-periodic phase functions. In general, the matrix $C$ defines the regulated outputs, while $J$ and $L$ are design choices selected to satisfy specific performance objectives. For instance, choosing $J = 0$ and $L = I$ yields a pure integral action, ensuring steady-state tracking of the desired outputs. In the application section, we illustrate how appropriate selections of $J, L$, and $C$ can be used to selectively target and attenuate specific harmonic components. In harmonic variables, we have

$$ \dot{Z}(t) = \omega(t) ((\mathcal{J} - \mathcal{N})Z(t) + \mathcal{L}\mathcal{C}X(t)) \quad (43) $$

where the operators $\mathcal{J} = \mathcal{T}_0(J)$, $\mathcal{L} = \mathcal{T}_0(L)$, $\mathcal{C} = \mathcal{T}_0(C)$ are assumed constant and bounded on $\ell^2$. Stacking $\tilde{X} = (X, Z)$ gives the augmented harmonic system

$$ \dot{\tilde{X}}(t) = (\tilde{\mathcal{A}}(\omega(t)) - \omega(t)\mathcal{N})\tilde{X}(t) + \tilde{\mathcal{B}}(\omega(t))U(t), \quad (44) $$

with

$$ \tilde{\mathcal{A}}(\omega) = \begin{bmatrix} \mathcal{G}(\omega)\mathcal{A}_0 + \omega\mathcal{A}_1 & 0 \\ \omega\mathcal{L}\mathcal{C} & \omega\mathcal{J} \end{bmatrix}, \tilde{\mathcal{B}}(\omega) = \begin{bmatrix} \mathcal{G}(\omega)\mathcal{B}_0 + \omega\mathcal{B}_1 \\ 0 \end{bmatrix}. $$

Let $U(t) = -\mathcal{K}\tilde{X}(t)$ where $\mathcal{K}$ is a state-feedback operator to be designed. Applying Theorem 7 to the closed-loop operator $\tilde{\mathcal{A}}_{cl}(\omega) = (\tilde{\mathcal{A}}(\omega) - \omega\mathcal{N}) - \tilde{\mathcal{B}}(\omega)\mathcal{K}$ yields, at each frequency vertex, a quadratic Lyapunov inequality in $\mathcal{P}$ and $\mathcal{K}$. Introducing the standard change of variables $\mathcal{Y} = \mathcal{KS}$, where $\mathcal{Y}$ and $\mathcal{S} = \mathcal{P}^{-1}$ are bounded Toeplitz block operators, yields the following set of linear matrix inequalities, which are convex in the decision variables:

$$ (\tilde{\mathcal{A}}_{\omega_i} - \omega_i \mathcal{N}) \mathcal{S} + \mathcal{S}(\tilde{\mathcal{A}}_{\omega_i} - \omega_i \mathcal{N})^* - \tilde{\mathcal{B}}_{\omega_i} \mathcal{Y} - \mathcal{Y}^* \tilde{\mathcal{B}}_{\omega_i}^* < 0, \quad (45) $$

and where $\tilde{\mathcal{A}}_{\omega_i} = \tilde{\mathcal{A}}(\omega_i)$, $\tilde{\mathcal{B}}_{\omega_i} = \tilde{\mathcal{B}}(\omega_i)$ with $\omega_i \in \{\omega_{\min}, \omega_{\max}\}$ and $\mathcal{G}(\omega_i) = \mathcal{I}$. The LMIs (45) are therefore a direct convex reformulation of the Lyapunov conditions of Theorem 7 applied to the augmented closed-loop AFM-LPP system. The feasibility of these LMIs guarantees closed-loop stability for all admissible time-varying frequencies and provides a quadratic Lyapunov function given by $\tilde{\mathcal{X}}^*\mathcal{P}\tilde{\mathcal{X}}$. The resulting design yields a bounded Toeplitz state-feedback gain $\mathcal{K} = \mathcal{Y}\mathcal{S}^{-1} = [\mathcal{K}_x, \mathcal{K}_z]$, admitting a time domain phase-periodic representation,

$$ K(\theta(t)) = \sum_{h \in \mathbb{Z}} K_h e^{j h \theta(t)}. \quad (46) $$

By introducing Toeplitz block weighting operators $\mathcal{R}$ and $\mathcal{Q}$, the LMIs (45) can be interpreted in an LQR-like framework as classically done in the finite dimension LTI setting [25]. The resulting controller is a guaranteed-cost state-feedback law that ensures closed-loop stability under time-varying frequency and is obtained by solving the following optimization problem:

$$ \min_{\mathcal{S}, \mathcal{Y}} \mathrm{Tr}_0(\mathcal{M}), \quad \text{s.t.} \quad (47) $$

$$ \begin{pmatrix} \Xi(\omega_i) & * & * \\ \mathcal{R}^{\frac{1}{2}}\mathcal{Y} & -\mathcal{I} & * \\ \mathcal{Q}^{\frac{1}{2}}\mathcal{S} & 0 & -\mathcal{I} \end{pmatrix} \le 0 \quad \begin{pmatrix} \mathcal{M} & \mathcal{I} \\ \mathcal{I} & \mathcal{S} \end{pmatrix} > 0 $$

with

$$ \Xi(\omega_i) = (\tilde{\mathcal{A}}_{\omega_i} - \omega_i \mathcal{N})\mathcal{S} + \mathcal{S}(\tilde{\mathcal{A}}_{\omega_i} - \omega_i \mathcal{N})^* - \tilde{\mathcal{B}}_{\omega_i} \mathcal{Y} - \mathcal{Y}^* \tilde{\mathcal{B}}_{\omega_i}^* $$

where $\mathrm{Tr}_0$ is the average-trace operator defined in [20], and $\mathcal{M}$ is a Toeplitz block operator bounded on $\ell^2$.

**Remark 2** The stability and synthesis results established in this section are exact at the infinite-dimensional harmonic level and do not rely on frequency-rate bounds or operator approximation. In practice, numerical implementation requires finite-dimensional truncations of the Toeplitz operators, for which feasibility-preserving convergence properties are well known and do not alter the theoretical guarantees [20]. The open-source PhasorArray Toolbox [26, 27] automates the construction of truncated Toeplitz operators and provides routines for harmonic Lyapunov, Riccati, and LMI problems, making implementation straightforward.

## 6.3 Classical harmonic controllers as special cases

The harmonic-domain controller structure introduced above is designed to preserve the AFM-LPP property of the closed-loop system and thereby admit an exact Lyapunov-based analysis under variable-frequency conditions. In this section, it is shown that this formulation
encompasses several classical harmonic control strategies, including proportional-resonant [28, 29], repetitive controllers [30] and modulation-based control strategies [31] as special cases. This observation establishes a direct connection between the proposed approach and well-known designs, and provides further insight into its practical relevance. For instance, the ideal repetitive controller given by $C(s) = \frac{e^{-sT_0}}{1-e^{-sT_0}}$ possesses infinite gain at integer multiples of the fundamental frequency. Truncating its Fourier expansion yields a bank of resonant oscillators, each associated with a specific harmonic. Indeed, with $H$ the number of harmonics retained one has $C(s) \approx \sum_{h=1}^{H} \frac{K_h s}{s^2+(h\omega_0)^2}$. Each harmonic therefore introduces a resonant subsystem of the form:

$$ \dot{z}_h = \begin{bmatrix} 0 & h\omega_0 \\ -h\omega_0 & 0 \end{bmatrix} z_h + \begin{bmatrix} 1 \\ 0 \end{bmatrix} e. \quad (48) $$

Note that, under the rotating-frame change of variables $w_h = R(h\theta)z_h$, this formulation is equivalent to [31]

$$ \dot{w}_h = \begin{bmatrix} \cos(h\theta) \\ \sin(h\theta) \end{bmatrix} e, \quad (49) $$

which corresponds to the integral structure used in modulation-based control schemes. When $\omega(t) = \omega_0(t)$ is time-varying and the system dynamics are phase-periodic, the situation differs, as the targeted harmonics are defined in the phase domain rather than in time, leading to index-based harmonic resonance matrix:

$$ \frac{\partial}{\partial \varphi} z_{\theta h}(\varphi) = \begin{bmatrix} 0 & h \\ -h & 0 \end{bmatrix} z_{\theta}(\varphi) + \begin{bmatrix} 1 \\ 0 \end{bmatrix} e_{\theta}(\varphi). \quad (50) $$

Consequently, within the proposed framework, the corresponding phase-domain realization is:

$$ \dot{Z}_h(t) = \omega(t) \left( \begin{bmatrix} -\mathcal{N} & h\mathcal{I} \\ -h\mathcal{I} & -\mathcal{N} \end{bmatrix} Z_h(t) + \begin{bmatrix} \mathcal{I} \\ 0 \end{bmatrix} E(t) \right), \quad (51) $$

which is directly compatible with the forwarding structure previously introduced. Stacking these oscillators for all selected harmonics produces a harmonic predictor that can be tuned via the same methodology used for the main controller. This viewpoint clarifies that PR and repetitive controllers are not separate design philosophies, but rather structured instances of the general internal-model construction enabled by the AFM-LPP harmonic representation. Moreover, it extends their applicability to systems with time-varying frequency.

# 7 Application to Variable-Speed Motor Drives

To illustrate the relevance of the proposed AFM-LPP framework, we apply it to a permanent magnet synchronous motor (PMSM) operating under variable-speed conditions. We show that the harmonic representation of the PMSM naturally fits into the AFM-LPP class introduced in Section 6, which enables the direct application of the exact stability and synthesis results developed in that section. This allows the construction of harmonic controllers ensuring tracking, harmonic mitigation, and closed-loop stability over a prescribed speed range. Throughout the remainder of this section, we assume that rotor speed $\omega_m(t)$ is measurable.

## 7.1 PMSM as a phase-periodic system

Consider a surface-mounted PMSM with $p$ pole pairs. Let the state vector be $x = (i_{abc}, \omega_m)$, where $i_{abc} = (i_a, i_b, i_c)^\top$ denotes the stator phase currents and $\omega_m$ the rotor mechanical speed. The control input is the applied stator voltage $u = v_{abc}$, and the external disturbance is the load torque $w = \Gamma_L$. The PMSM dynamics can then be written in state-space form as

$$ \dot{x}(t) = A(\theta(t))x(t) + B_u u(t) + B_w w(t), \quad (52) $$

where the system matrices are given by

$$ A(\theta) = \begin{bmatrix} -\frac{r}{L}I_3 & \frac{p\psi_f}{L}\Phi_{abc}(p\theta) \\ -\frac{p\psi_f}{J}\Phi_{abc}(p\theta)^\top & -\frac{B_f}{J} \end{bmatrix}, \\ B_u = \begin{bmatrix} \frac{1}{L}I_3 \\ 0_{1 \times 3} \end{bmatrix}, \quad B_w = \begin{bmatrix} 0_{3 \times 1} \\ -\frac{1}{J} \end{bmatrix}. $$

The machine parameters are the stator phase inductance $L$, phase resistance $r$, rotor inertia $J$, viscous friction coefficient $B_f$, and permanent magnet flux linkage $\psi_f$. The vector $\Phi_{abc}(p\theta)$ represents the back-EMF waveform, which is $2\pi$-periodic in the electrical angle $p\theta$. For an ideal surface-mounted PMSM, this waveform is purely sinusoidal and given by

$$ \Phi_{abc}(p\theta) = \begin{bmatrix} \sin(p\theta) \\ \sin(p\theta - 2\pi/3) \\ \sin(p\theta + 2\pi/3) \end{bmatrix} \quad (53) $$

The state matrix $A(\theta)$ depends periodically on the electrical angle through $\Phi_{abc}(p\theta)$, which makes the PMSM dynamics phase-periodic. Consequently, the PMSM provides a canonical example of an AFM-LPP system and is therefore directly compatible with the harmonic modeling and control framework developed in Section 6. Applying the variable-frequency Fourier decomposition $\mathcal{F}_\theta$
to the state-space model (52), where the phase variable $\theta$ satisfies $\dot{\theta} = \omega$ with $\omega$ a time-varying parameter, leads to the following AFM-LPP harmonic representation which matches the exact harmonic PMSM model when $\omega = \omega_m$:

$$ \dot{X}(t) = (\mathcal{G}(\omega(t))\mathcal{A} - \omega(t)\mathcal{N})X(t) + \mathcal{G}(\omega(t))(\mathcal{B}_uU(t) + \mathcal{B}_wW(t)) \quad (54) $$

where $X = \mathcal{F}_{\theta}(x) = (I_{abc}, \Omega)$, $U = V_{abc} = \mathcal{F}_{\theta}(u)$ and $W = \mathcal{F}_{\theta}(w)$ and where the constant Toeplitz-block operators $\mathcal{A}, \mathcal{B}_u$, and $\mathcal{B}_w$ are given by:

$$ \mathcal{A} = \begin{bmatrix} -\frac{r}{L}\mathcal{I}_3 & \frac{\psi_f p}{L}\mathcal{T}_{\theta}(\Phi_{abc}) \\ -\frac{p\psi_f}{J}\mathcal{T}_{\theta}(\Phi_{abc})^* & -\frac{B_f}{J}\mathcal{I}_1 \end{bmatrix}, \quad (55a) $$

$$ \mathcal{B}_u = \begin{bmatrix} \frac{1}{L}\mathcal{I}_3 \\ 0_{1 \times 3} \end{bmatrix}, \quad \mathcal{B}_w = \begin{bmatrix} 0_{3 \times 1} \\ -\frac{1}{J}\mathcal{I}_1 \end{bmatrix} \quad (55b) $$

where $\mathcal{I}_n = \mathcal{T}_{\theta}(Id_n)$, and $\mathcal{T}_{\theta}(\Phi_{abc})$ represents the back-EMF harmonics. For the ideal sinusoidal back-EMF (53), the Toeplitz operators $\mathcal{T}_{\theta}(\Phi_{abc})$ exhibit a sparse structure. For example, since $\Phi_a(\theta) = \sin(p\theta) = \frac{1}{2j}(e^{jp\theta} - e^{-jp\theta})$, the Toeplitz matrix $\mathcal{T}_{\theta}(\Phi_a)$ has only two nonzero diagonals, at offsets $\pm p$, and each of their diagonal elements are equal to $\pm \frac{j}{2}$. The harmonic model (54) is affine in the frequency $\omega$ and Toeplitz structured, and therefore satisfies the AFM-LPP assumptions required by the analysis and synthesis results of Section 6.2.

## 7.2 Control objectives and synthesis

The control design aims to enforce phase-periodic steady-state behavior and harmonic mitigation, while preserving stability under variable-frequency operation. In particular, the following objectives are considered:

(1) Speed regulation: accurate tracking of the mechanical speed $x_4(t) = \omega_m(t)$ to a given reference $\omega_{m,0}^{\text{ref}}$ in the presence of load torque disturbances $\Gamma_L$.

(2) Power factor optimization: regulation of the d-axis current to zero (or to a prescribed reference) in order to maximize efficiency.

(3) Harmonic mitigation: rejection of phase-periodic current harmonics induced by torque ripple at $2\omega_m$, which generates current harmonics in $i_{abc}$ of order $p \pm 2k$, $k = 1, 2, \dots$.

These objectives are addressed by enforcing zero steady-state values of selected phase-domain output components using integral action in the harmonic domain. Speed regulation is achieved by integrating the error between the DC component $\Omega_0(t)$ of the mechanical speed phasor and the reference $\omega_{m,0}^{\text{ref}}$.

Power factor optimization is achieved by regulating the DC component of the d-axis current $i_d$ in the rotating dq frame where $i_{dq} = T_p(\theta)(i_{abc})$ and $v_{dq} = T_p(\theta)(v_{abc})$ and where the generalized Park transformation [32] is provided by (here, for $k = p$):

$$ \begin{aligned} i_{dq}^{k} &= T_{k}(\theta) i_{abc} \\ &= \frac{2}{3} \left[ \begin{array}{l} \cos(k\theta) \quad \cos(k\theta - \frac{2\pi}{3}) \quad \cos(k\theta + \frac{2\pi}{3}) \\ -\sin(k\theta) \quad -\sin(k\theta - \frac{2\pi}{3}) \quad -\sin(k\theta + \frac{2\pi}{3}) \end{array} \right] i_{abc}. \end{aligned} \quad (56) $$

This corresponds to enforcing $I_{d,0} = 0$, where $I_d = \mathcal{F}_{\theta}(i_d)$ is the harmonic representation of the d-axis current.

Harmonic mitigation is implemented by exploiting generalized Park transformations $T_k(\theta)$ for $k \neq p$ which map the k-th harmonic components of the three-phase currents into DC signals in the rotating frame. Integral action on these quantities enforces rejection of the corresponding harmonics in phase-periodic steady state. The resulting control output operator is defined as

$$ C(\theta) = \begin{bmatrix} 0 & 1 \\ [1 \ 0]T_p(\theta) & 0 \\ \text{stack}_{k=0,2,6,8}(T_k(\theta) \ 0) \end{bmatrix}, \quad (57) $$

which selects the relevant speed, d-axis, and harmonic components to be regulated.

In the harmonic-domain framework, tracking of phase-periodic operating trajectories requires shifting the origin to the harmonic equilibrium associated with the reference speed and disturbance. The phase-periodic steady state dynamics can be characterized from the definition of an equilibrium in phase domain:

**Definition 8** In phase-domain, an equilibrium $(X^{\text{ref}}, U^{\text{ref}})$ for a phase-periodic disturbance $W$ satisfies the nonlinear equation:

$$ 0 = (\mathcal{A} - \mathcal{T}_{\theta}(\omega_m^{\text{ref}})\mathcal{N})X^{\text{ref}} + \mathcal{B}_uU^{\text{ref}} + \mathcal{B}_wW, \quad (58) $$

where the phase frequency $\omega_m^{\text{ref}}$ corresponds to the speed rotor reference.

Assuming a known phase-periodic disturbance $W$ and the reference speed trajectory $\omega_m^{\text{ref}}(t)$ takes values in the interval $[\omega_{\min}, \omega_{\max}]$, the full derivation of the equilibrium $(X^{\text{ref}}, U^{\text{ref}}) = \mathcal{F}_{\theta}(x^{\text{ref}}, u^{\text{ref}})$ is given in Appendix. Note that, when the phase function $\omega$ differs from the reference $\omega_m^{\text{ref}}$, the corresponding equilibrium quantities evolve according to the harmonic dynamics associated with the actual phase trajectory

$$ \dot{X}^{\text{ref}}(t) = (\mathcal{G}(\omega(t))\mathcal{A} - \omega(t)\mathcal{N})X^{\text{ref}}(t) + \mathcal{G}(\omega(t))(\mathcal{B}_uU^{\text{ref}}(t) + \mathcal{B}_wW). \quad (59) $$
We now rewrite the system around the phase-periodic equilibrium in order to synthesize a stabilizing controller for the deviation dynamics. Let the state and input errors be given by $E = X - X^{\text{ref}}$ and $\tilde{U} = U - U^{\text{ref}}$, where $(X^{\text{ref}}, U^{\text{ref}})$ satisfies the harmonic equilibrium equation (59). The resulting error dynamics are given by

$$ \dot{E}(t) = (\mathcal{G}(\omega(t))A - \omega(t)\mathcal{N})E(t) + \mathcal{G}(\omega(t))B_u\tilde{U}(t) $$

Using the control output as defined by (57), the harmonic forwarding dynamics are introduced according to (43) by choosing $\mathcal{J} = 0$ and $\mathcal{L} = \mathcal{I}$, yielding

$$ \dot{Z}(t) = \omega(t)(-\mathcal{N}Z(t) + \mathcal{C}E(t)). \quad (60) $$

The corresponding harmonic-augmented system then reads with $\tilde{E} = (E, Z)$

$$ \tilde{E}(t) = (\tilde{\mathcal{A}}(\omega(t)) - \omega(t)\mathcal{N})\tilde{E}(t) + \tilde{\mathcal{B}}_v(\omega(t))\tilde{U}(t) \quad (61) $$

where

$$ \tilde{\mathcal{A}}(\omega) = \begin{bmatrix} \mathcal{G}(\omega)A & 0 \\ \omega\mathcal{C} & 0 \end{bmatrix}, \quad \tilde{\mathcal{B}}_v(\omega) = \begin{bmatrix} \mathcal{G}(\omega)B_u \\ 0 \end{bmatrix}. \quad (62) $$

The augmented deviation system (61) is an AFM-LPP system. The synthesis results of Section 6 therefore yield a Toeplitz-block feedback gain $\mathcal{K} = [\mathcal{K}_x \ \mathcal{K}_z]$ obtained by solving the harmonic synthesis problem (47) and using the reconstruction procedure (17). The corresponding time-domain control law is

$$ u(t) = K_x(\theta(t))(x(t) - x^{\text{ref}}(t)) + K_z(\theta(t))z(t) + u^{\text{ref}}(t). \quad (63) $$

This controller guarantees asymptotic convergence $E(t) = X(t) - X^{\text{ref}}(t) \to 0$ provided that $\omega \in [\omega_{\min}, \omega_{\max}]$. However, for the PMSM, $\omega$ is not a free variable but is constrained by the mechanical speed, i.e., $\omega = \omega_m$ and $\theta = \theta_m$. Therefore, the parameter $\omega$ coincides with a state variable, and the stability guarantee holds as long as the state trajectory remains inside the admissible domain $[\omega_{\min}, \omega_{\max}]$. Proposition 1, given and proved in the appendix, provides a sufficient condition ensuring this property. This result ensures that the frequency admissibility condition required by the AFM-LPP stability theorem is preserved along closed-loop trajectories. It is obtained by determining the largest quadratic Lyapunov level set contained inside the admissible domain. The admissible level is uniquely determined by the minimum distance between the DC component $\omega_{m,0}^{\text{ref}}$ of the rotor speed reference and the bounds $\omega_{\min}$ and $\omega_{\max}$. More specifically, let the Lyapunov function $\mathcal{V}(\tilde{E}) = \tilde{E}^* \mathcal{P} \tilde{E}$, where $\tilde{E}$ is the augmented state defined in (61) and $\mathcal{P} = S^{-1}$, where $S$ is the solution of (47). Let

$S$ be the representative of $S$ in time domain and denote by $s_{44}$ the fourth diagonal entry of $S$ and by $S_{44,0}$ the zero-order phasor of $S_{44} = \mathcal{F}_{\theta}(s_{44})$. For a given step reference $\omega_{m,0}^{\text{ref}}(t_0)$, if the initial augmented state $\tilde{E}(t_0)$ satisfies

$$ \mathcal{V}(\tilde{E}(t_0)) < L_{\max}(\omega_{m,0}^{\text{ref}}(t_0)), \quad (64) $$

where

$$ L_{\max}(\omega_{m,0}^{\text{ref}}) = \frac{\delta^2(\omega_{m,0}^{\text{ref}})}{S_{44,0}}, $$

with

$$ \delta(\omega_{m,0}^{\text{ref}}) = \min(|\omega_{\min} - \omega_{m,0}^{\text{ref}}|, |\omega_{\max} - \omega_{m,0}^{\text{ref}}|), $$

then the closed-loop trajectory satisfies

$$ E(t) \to 0 \quad \text{and} \quad \omega_m(t) \in [\omega_{\min}, \omega_{\max}] $$

for all $t \ge t_0$. Moreover, any subsequent step change of the reference $\omega_{m,0}^{\text{ref}}(t_k)$ $k = 1, 2, \dots$ preserving (64) is admissible and does not affect stability.

### 7.3 Simulation results

To validate the proposed harmonic control framework, numerical simulations are performed on the PMSM operating under variable-speed conditions and phase-periodic load torque disturbances. The harmonic synthesis problem (47) is solved using the weighting matrices $Q = \mathcal{I}$ and $R = 100\mathcal{I}$. The operating speed range is chosen as $[\omega_{\min}, \omega_{\max}] = [10, 200]$ rad/s (approximately 100–2000 rpm), representative of low-to-medium speed operation. Using a 10th-order truncation scheme, the computation time required to obtain a sufficiently accurate solution is 114 s. The PMSM parameters used in simulation are reported in Table 1.

<table><caption>Table 1
PMSM parameters</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Stator resistance</td><td>r</td><td>0.5 Ω</td></tr><tr><td>Phase inductance</td><td>L</td><td>1.5 mH</td></tr><tr><td>Flux linkage</td><td>ψ<sub>f</sub></td><td>0.14 Wb</td></tr><tr><td>Rotor inertia</td><td>J</td><td>0.03 kg m<sup>2</sup></td></tr><tr><td>Viscous friction coefficient</td><td>B<sub>f</sub></td><td>0.02 N m s/rad</td></tr><tr><td>Number of pole pairs</td><td>p</td><td>4</td></tr></tbody></table>

Figure 2 shows the closed-loop time-domain responses obtained with the synthesized harmonic controller, including the phase currents $i_{abc}$, dq-currents, rotor speed $\omega_m$, load torque $\Gamma_L$ and the applied control voltages $v_{abc}$. The corresponding harmonic spectra are displayed in Figure 3. For comparison, Figures 4 and 5 report the system response obtained when harmonic mitigation is not
<figure><img src="image_3.png" /><figcaption>Figure 2. Phase currents $i_{abc}$, current $i_{dq}$, rotor speed $\omega_m$ (blue), $\omega_m^{ref}$ (yellow) and $\omega_{m,0}^{ref}$ (red), load torque $\Gamma_L$ and the applied control $u = v_{abc}$.</figcaption></figure>

included in the control objectives, while keeping identi-
cal speed references and torque disturbances. The sim-
ulation results highlight the following key properties:

*   *Speed regulation:* The rotor speed $\omega_m$ accurately tracks its time-varying reference $\omega_m^{ref}$, defined in (66c) (Appendix) and computed from $\omega_{m,0}^{ref}$, over the entire operating range, as illustrated in the zoom of Figure 2. This supports the mechanical filtering assumption underlying its construction.

* *d-axis current regulation:* the mean value of the direct-axis current $i_d$ is maintained close to zero, ensuring power factor optimization.

*   *Harmonic mitigation:* when harmonic rejection is enabled, the steady-state phase currents $i_{abc}$ exhibit no harmonic components other than the fundamental p-th harmonic, independently of the reference speed, as shown in Figure 3. In particular, this implies $i_d = 0$ and $i_q = \text{const}$ in phase-periodic steady state.

Figure 2 includes two zoomed views: one at quasi-
constant rotor speed and one under time-varying speed.
These plots highlight that phase-periodicity with respect
to the electrical angle $\theta$ should not be confused with
periodicity in the time domain, as the latter is gener-
ally non-sinusoidal under variable-speed operation (see
the zoomed area around $t = 8$ s). Finally, the level-set
guarantee is illustrated in Figure 6, which compares the

<figure><img src="image_4.png" /><figcaption>Figure 3. Magnitude of harmonics illustrating targeted suppression in $i_a$ and $i_d$. Top to bottom: $I_{a,k}$, $I_{d,k}$, $\Omega_{m,k}$, $\Gamma_k$ and $V_{a,k}$, $k = 0, 1, \dots, 8$.</figcaption></figure>

theoretical maximum admissible Lyapunov level $L_{\max}$ with the closed-loop Lyapunov trajectory. The results show that the Lyapunov function remains strictly below the admissible bound during reference step changes, confirming closed-loop stability is preserved.

**Conclusion**

This paper establishes a unified and rigorous framework for harmonic modeling and control of systems under variable-frequency. An exact harmonic-domain formulation is derived, providing a principled basis for analysis and systematic approximation, together with a quantitative validity criterion that characterizes the modeling error induced by frequency variations. For linear phase-periodic systems with affine dependence on the instantaneous frequency, stability and state-feedback synthesis are shown to be certified exactly by enforcing Lyapunov inequalities at a finite number of frequency vertices. This result bridges harmonic modeling and LMI-based control design, enabling tractable certification over continuous frequency intervals. The approach is illustrated through harmonic mitigation in variable-speed electrical drives, demonstrating both practical relevance and effectiveness for systems operating under nonstationary frequency conditions. Beyond electrical applications, the proposed framework provides a general harmonic-domain methodology for stability analysis and control of phase-periodic systems with time-varying frequency.
<figure><img src="image_5.png" /><figcaption>Figure 4. Simulation results without harmonic mitigation: the phase currents $i_{abc}$, the dq currents $i_{dq}$, the rotor speed $\omega_m$ and the applied control $v_{abc}$, for the same load torque and reference speed profile as in the previous case.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6. Admissible level $L_{max}$ (red) and Lyapunov function $V(\tilde{E})$ (blue) for the reference $\omega_{m,0}^{ref}$ of Fig. 2.</figcaption></figure>

# 8 Appendix

This appendix provides technical details supporting the application section. First, the phase-periodic equilibrium of the PMSM is derived in the harmonic domain. Then, the maximum admissible Lyapunov level ensuring frequency admissibility is characterized.

<figure><img src="image_7.png" /><figcaption>Figure 5. Magnitude of harmonics without harmonic mitigation objective. Top to bottom: $I_{a,k}$, $I_{d,k}$, $\Omega_{m,k}$ and $V_{a,k}$, $k = 0, 1, \dots, 8$.</figcaption></figure>

## 8.1 PMSM Phase-periodic equilibrium

Given a known phase-periodic disturbance $W$, the harmonic equilibrium ($X^{ref}$, $U^{ref}$) = $\mathcal{F}_{\theta}(x^{ref}, u^{ref})$ is obtained from the dq-frame PMSM model derived using the Park transform $T_p$ defined in (56). The corresponding phase-periodic harmonic model is given by

$$ 
\begin{aligned}
L\dot{I}_d &= \mathcal{G}(\omega(t))(-rI_d + V_d) + \omega(t)L(pI_q - NI_d) \\
L\dot{I}_q &= \mathcal{G}(\omega(t))(-rI_q + V_q) - \omega(t)(p\psi\mathcal{I} + pLI_d + LN\mathcal{I}_q) \\
J\dot{\Omega} &= \mathcal{G}(\omega(t))\left(\frac{3}{2}p\psi I_q - W - b_f\Omega\right) - \omega(t)JN\Omega
\end{aligned}
$$

where $I_{dq} = \mathcal{F}_{\theta}(i_{dq})$ and $V_{dq} = \mathcal{F}_{\theta}(v_{dq})$ denote the harmonic representations of the dq currents and voltages.

Since the physical PMSM satisfies $\omega = \omega_m$, and under phase-periodic steady-state operation (i.e., when $\Omega = \mathcal{F}_{\theta}(\omega_m)$ is constant), the equilibrium equation is obtained by left-multiplying by $\mathcal{G}^{-1}(\omega)$:

$$ 
\begin{aligned}
0 &= \frac{1}{L}(-rI_d + V_d) + \mathcal{T}_{\theta}(\omega_m)(pI_q - NI_d) \\
0 &= \frac{1}{L}(-rI_q + V_q) - \mathcal{T}_{\theta}(\omega_m)\left(\frac{p\psi}{L}\mathcal{I} + pI_d + NI_q\right) \\
0 &= \frac{1}{J}\left(\frac{3}{2}p\psi I_q - W - b_f\Omega\right) - \mathcal{T}_{\theta}(\omega_m)N\Omega
\end{aligned}
$$
Under the selected integral actions, $\Omega_0 = \omega_{m,0}^{\text{ref}}$, $I_{d,k} = 0$ for all $k \in \mathbb{Z}$, and $I_{q,k} = 0$ for all $k \neq 0$, the k-th harmonic component of the last equation with $k \neq 0$ yields the nonlinear harmonic recursion

$$ \Omega_k = H(jk\omega_{m,0}^{\text{ref}}) \left( \frac{1}{J} W_k + \sum_{p=-\infty, p \neq k}^{+\infty} \Omega_{k-p} \Omega_p(jp) \right) \quad (65) $$

where

$$ H(jk\omega) = \frac{-1}{\frac{b_f}{J} + jk\omega}. $$

Equation (65) defines a nonlinear harmonic fixed-point relation. The equilibrium can therefore be computed using a fixed-point iteration initialized from the fundamental component. For $k=0$, using the symmetry property $\Omega_{-p} = \Omega_p^*$, the classical steady-state relation is recovered:

$$ I_{q,0} = \frac{2}{3p\psi} (W_0 + b_f \omega_{m,0}^{\text{ref}}), \quad \Omega_0 = \omega_{m,0}^{\text{ref}}. $$

In particular, if $W_k = 0$ for all $k \neq 0$, then $\Omega_k = 0$ for all $k \neq 0$. The recursion (65) reveals that higher-order mechanical harmonics arise from nonlinear harmonic coupling. However, due to the low-pass filtering effect of the mechanical dynamics, higher-order components are strongly attenuated, which justifies truncating the recursion in practice. As an illustration, consider a disturbance containing only a second harmonic component $W_2$. Starting from the initialization $\Omega_k = 0$ for all $k \neq 0$, the first iteration yields

$$ \Omega_2 = H(2j\omega_{m,0}^{\text{ref}}) \frac{1}{J} W_2 \text{ and } \Omega_{-2} = \bar{\Omega}_2 $$

Higher-order components are generated through harmonic coupling, for example, at second iteration appears:

$$ \Omega_4 = H(4j\omega_{m,0}^{\text{ref}})\Omega_2^2(2j). $$

Because of mechanical filtering, these higher-order terms are typically negligible, leading to the approximation $\Omega_{\pm 4} \approx 0$. Under this approximation, the diffusion effect vanishes and the harmonic equilibrium is given by

$$ \begin{aligned} I_{q,0}^{\text{ref}} &= \frac{2}{3p\psi_f} (W_0 + B_f \omega_{m,0}^{\text{ref}}), \\ \Omega_0^{\text{ref}} &= \omega_{m,0}^{\text{ref}}, \\ \Omega_2^{\text{ref}} &= \frac{1}{J} H(2j\omega_{m,0}^{\text{ref}}) W_2, \\ \Omega_k^{\text{ref}} &= 0, \quad k \neq 0, \pm 2, \end{aligned} $$

These harmonic-domain equilibrium quantities uniquely define the corresponding phase-periodic time-domain

reference trajectories:

$$ i_d^{\text{ref}}(t) = 0, \quad (66a) $$

$$ i_q^{\text{ref}}(t) = \frac{2}{3p\psi_f} (W_0 + B_f \omega_{m,0}^{\text{ref}}), \quad (66b) $$

$$ \omega_m^{\text{ref}}(t) = \omega_{m,0}^{\text{ref}} + 2|\Omega_2| \cos(2\theta(t) - \varphi), \quad (66c) $$

where $\varphi = \angle \Omega_2$. The control reference $U^{\text{ref}}$ is then obtained from the electrical equilibrium equations:

$$ V_d^{\text{ref}} = -\mathcal{T}_{\theta}(\omega_m^{\text{ref}}) pL I_q^{\text{ref}}, \quad V_q^{\text{ref}} = rI_q^{\text{ref}} + \mathcal{T}_{\theta}(\omega_m^{\text{ref}}) p\psi_f, $$

or, equivalently, in the time domain,

$$ \begin{aligned} v_d^{\text{ref}}(t) &= -\omega_m^{\text{ref}}(t) pL i_q^{\text{ref}}(t), \\ v_q^{\text{ref}}(t) &= r i_q^{\text{ref}}(t) + \omega_m^{\text{ref}}(t) p\psi_f. \end{aligned} $$

Finally, the corresponding phase-periodic equilibrium $(x^{\text{ref}}, u^{\text{ref}})$ is obtained via the inverse Park transformation, $i_{abc}^{\text{ref}} = T_p^{-1}(i_{dq}^{\text{ref}})$ and $v_{abc}^{\text{ref}} = T_p^{-1}(v_{dq}^{\text{ref}})$.

## 8.2 Maximum Admissible Lyapunov Level

The following proposition provides the largest Lyapunov level set ensuring that the rotor speed trajectory remains inside the admissible interval $[\omega_{\min}, \omega_{\max}]$. Geometrically, this corresponds to determining the largest ellipsoidal level set of the Lyapunov function that is fully contained inside the admissible speed domain. Since the Lyapunov function is quadratic and radially unbounded, the maximal admissible level is necessarily attained when the speed trajectory reaches one of the admissible boundaries. The proof therefore reduces to characterizing the Lyapunov level associated with these boundary conditions.

**Proposition 1** Consider, at time $t_0$, a speed reference $\omega_{m,0}^{\text{ref}}(t_0)$ and the corresponding equilibrium pair $(X^{\text{ref}}, U^{\text{ref}})$. Define the Lyapunov function

$$ V(\tilde{E}) = \tilde{E}^* P \tilde{E}, $$

where $\tilde{E}$ is the augmented state defined in (61), and $P = S^{-1}$ with $S$ is the solution of problem (47). If the initial condition satisfies

$$ V(\tilde{E}(t_0)) < L_{\max}(\omega_{m,0}^{\text{ref}}), \quad (68) $$

where $L_{\max}$ denotes the largest Lyapunov level set for which the rotor speed component $x_4 = \omega_m$ remains inside the admissible interval $\omega_m(t) \in [\omega_{\min}, \omega_{\max}]$, then the error state satisfies $E(t) \to 0$ as $t \to \infty$. Furthermore, any subsequent step change in the reference $\omega_{m,0}^{\text{ref}}(t_k)$, $k = 1, 2, \dots$, that preserves condition (68) is admissible
<figure><img src="image_8.png" /><figcaption>Figure 7. Projection of the maximum level set on the subspace $\Omega \setminus \Omega_0 \times \Omega_0$</figcaption></figure>

and does not affect stability. The admissible level $L_{\max}$ is given by

$$L_{\max}(\omega_{m,0}^{\text{ref}}) = \frac{\delta^2(\omega_{m,0}^{\text{ref}})}{S_{44,0}}, \quad (69)$$

where $\delta(\omega_{m,0}^{\text{ref}}) = \min(|\omega_{\min} - \omega_{m,0}^{\text{ref}}|, |\omega_{\max} - \omega_{m,0}^{\text{ref}}|)$, and $S_{44,0}$ denotes the zero-order phasor of $S_{44} = \mathcal{F}_{\theta}(s_{44})$, with $s_{44}$ being the fourth diagonal entry of the time-domain realization $S$.

**PROOF.** Consider a reference $\omega_{m,0}^{\text{ref}} \in [\omega_{\min}, \omega_{\max}]$ and the associated equilibrium $(X^{\text{ref}}, U^{\text{ref}})$. Assume that the corresponding rotor speed reference $\omega_m^{\text{ref}}$ remains inside the admissible interval $[\omega_{\min}, \omega_{\max}]$. The objective is to determine the largest Lyapunov level set $L_{\max}$ for which the rotor speed component $x_4 = \omega_m$ remains inside this interval.

Since $\omega(t) > 0$ if and only if $\mathcal{T}_{\theta}(\omega)(t) > 0$, admissibility of the rotor speed is equivalent to the harmonic constraints

$$\omega_{\min} \mathcal{I} \le \mathcal{T}_{\theta}(\omega_m)(t) \le \omega_{\max} \mathcal{I}.$$

Because the Lyapunov function is radially unbounded, the maximal admissible level set is necessarily attained on the boundary of these constraints. Considering, for instance, the lower bound, there exists $X \in \ell^2$, $X \neq 0$, such that

$$(\mathcal{T}_{\theta}(\omega_m) - \omega_{\min} \mathcal{I})X = 0.$$

In the time domain, this condition implies $(\omega_m(t) - \omega_{\min})x(t) = 0$ almost everywhere on $[t-T(t), t]$. Since $x$ is not identically zero and $\omega_m$ is continuous, it follows that $\omega_m(t) = \omega_{\min}$ on this interval. Therefore, the maximal level set is obtained by considering boundary trajectories satisfying $\omega_m(t) = \omega_{\min}$ or $\omega_m(t) = \omega_{\max}$.

Because the Lyapunov function defines convex elliptical level sets centered at the equilibrium, the active constraint is necessarily the boundary closest to the equilibrium speed $\Omega^{\text{ref}}$ (see Fig. 7).

Therefore, the admissible level is determined by

$$\delta(\omega_{m,0}^{\text{ref}}) = \min(|\omega_{\min} - \omega_{m,0}^{\text{ref}}|, |\omega_{\max} - \omega_{m,0}^{\text{ref}}|).$$

Without loss of generality, assume that the minimum is attained at $\omega_{\min}$.

Let $\Omega = \mathcal{F}_{\theta}(\omega_m)$. Since $\omega_m$ is constant on the constraint boundary, we set $\Omega_k = 0$ for all $k \neq 0$ in the augmented state $\tilde{E}$. The computation of the maximal level set then reduces to the quadratic optimization problem

$$\max_{X,Z} \mathcal{V}(\tilde{E})$$

subject to the constraint $C_0(\Omega_0) = \Omega_0 - \omega_{\min} \ge 0$. Using the KKT optimality conditions, there exists a Lagrange multiplier $\alpha$ such that

$$\tilde{E} = \frac{1}{2} \alpha \mathcal{S} \nabla_{[X,Z]} C_0(\Omega_0),$$

where $\mathcal{S} = \mathcal{P}^{-1}$ and $\nabla_{[X,Z]} C_0(\Omega_{m,0})$ denotes the gradient of the constraint with respect to $(X, Z)$. Considering the zero-order phasor associated with $\Omega$ (i.e., the central row of the fourth Toeplitz block), we obtain

$$\Omega_0 - \omega_{m,0}^{\text{ref}} = \frac{1}{2} \alpha (\mathcal{S} \nabla_{[X,Z]} C_0(\Omega_0))_{4,0}.$$

Since $\nabla_{[X,Z]} C_0(\Omega_0)$ is zero everywhere except for the component associated with $\Omega_0$, which is equal to 1, the Toeplitz structure of $\mathcal{S}$ yields

$$(\mathcal{S} \nabla_{[X,Z]} C_0(\Omega_0))_4 = S_{44} = \mathcal{F}_{\theta}(s_{44}).$$

Therefore,

$$\Omega_0 - \omega_{m,0}^{\text{ref}} = \frac{1}{2} \alpha S_{44,0}.$$

Enforcing the boundary condition $\Omega_0 = \omega_{\min}$ gives

$$\alpha = 2 \frac{\omega_{\min} - \omega_{m,0}^{\text{ref}}}{S_{44,0}}.$$

Substituting the corresponding expression of $\tilde{E}$ into the Lyapunov function yields (69), which concludes the proof. □

## References

[1] H. Muramatsu and S. Katsura. An adaptive periodic disturbance observer for periodic disturbance suppression. *IEEE Trans. on Industrial Informatics*, 14(10), 2018.

[2] Q. Han, L. Wang, L. Marconi, Z. Liu, and H. Su. Adaptive control of nonlinear systems with sinusoidal uncertainties via internal model principle. In *IFAC PapersOnLine*, 2023.

[3] D. Astolfi, L. Praly, and L. Marconi. Harmonic internal models for structurally robust periodic output regulation. *Systems & Control Letters*, 161:105154, 2022.
[4] A. Cecilia, D. Astolfi, M. Bin, and R. Costa-Castelló. Canceling output disturbances in observer design through internal model filters. Automatica, 162:111529, 2024.

[5] C.H. Chung and M.S. Chen. A robust adaptive feedforward control in repetitive control design for linear systems. Automatica, 48(1):183–190, 2012.

[6] S.R. Sanders, J.M. Noworolski, X.Z. Liu, and G.C. Verghese. Generalized Averaging Method for Power Conversion Circuits. IEEE Trans. on Power Electronics, 6(2), 1991.

[7] X. Kwon, J. and Wang, F. Blaabjerg, and C.L. Bak. Frequency-Domain Modeling and Simulation of DC Power Electronic Systems Using Harmonic State Space Method. IEEE Trans. on Power Electronics, 32(2):1044–1055, 2017.

[8] P. Mattavelli, G.C. Verghese, and A.M. Stankovic. Phasor dynamics of thyristor-controlled series capacitor systems. IEEE Trans. on Power Systems, 12(3):1259–1267, 1997.

[9] A.M. Stankovic, B.C. Lesieutre, and T. Aydin. Modeling and analysis of single-phase induction machines with dynamic phasors. IEEE Trans. on Power Systems, 14(1):9–14, February 1999.

[10] Y. Levron and J. Belikov. Modeling power networks using dynamic phasors in the dq0 reference frame. Electric Power Systems Research, 144:233–242, March 2017.

[11] E. Karami, G.B. Gharehpetian, M. Madrigal, and J.J. Chavez. Dynamic phasor-based analysis of unbalanced three-phase systems in presence of harmonic distortion. IEEE Trans. on Power Systems, 33(6):6642–6654, 2018.

[12] Y. Cai, Y. He, H. Du, and J. Liu. Harmonic state-space modeling and system characteristic analysis of grid-connected inverter parallel-operation system considering asynchronous carriers. IEEE Trans. on Power Electronics, 39(7):8645–8666, July 2024.

[13] R. Kong, S. Sahoo, S. Ou, X. Meng, G. Gao, and F. Blaabjerg. Data-light oscillation mode identification for fast stability assessment of grid-tied converters. IEEE Trans. on Power Electronics, 40(8):11170–11181, August 2025.

[14] J. B. Kwon. Harmonic State Space (HSS) Modeling for Power Electronic Based Power Systems. PhD thesis, Aalborg Universitetsforlag, 2017.

[15] J.E. Ormrod. Harmonic State Space Modelling of Voltage Source Converters. PhD thesis, University of Canterbury, Christchurch, New Zealand, 2013.

[16] Y. Cai, Y. He, H. Zhang, H. Zhou, and J. Liu. Research on harmonic state-space modeling and calculation analysis of low-switching-frequency grid-connected inverter considering the impact of digitization. IEEE Trans. on Power Electronics, 38(1):1003–1021, 2023.

[17] J. Lyu, X. Zhang, X. Cai, and M. Molinas. Harmonic State-Space Based Small-Signal Impedance Modeling of a Modular Multilevel Converter With Consideration of Internal Harmonic Dynamics. IEEE Trans. on Power Electronics, 34(3):2134–2148, 2019.

[18] N. Blin, P. Riedinger, J. Daafouz, L. Grimaud, and P. Feyel. Necessary and Sufficient Conditions for Harmonic Control in Continuous Time. IEEE Trans. on Automatic Control, 67(8):4013–4028, 2022.

[19] T. Yang, S. Bozhko, J.M. Le-Peuvedic, G. Asher, and C.I. Hill. Dynamic Phasor Modeling of Multi-Generator Variable Frequency Electrical Power Systems. IEEE Trans. on Power Systems, 31(1):563–571, 2016.

[20] F. Vernerey, P. Riedinger, and J. Daafouz. A TBLMI Framework for Harmonic Robust Control. IEEE Trans. on Automatic Control, 70(4):2537–2551, 2025.

[21] M. Grosso, P. Riedinger, J. Daafouz, S. Pierfederici, H. Janati-Idrissi, and B. Lapôtre. Frequency-varying harmonic domain control for PMSMs with current harmonic mitigation. In IEEE Energy Conversion Conference Congress and Exposition (ECCE), 2025.

[22] M. Grosso, P. Riedinger, J. Daafouz, S. Pierfederici, H. Janati-Idrissi, and B. Lapôtre. Control of three-phase PWM rectifier using multiple frame dq transform and harmonic modeling. In 50th Annual Conference of the IEEE Industrial Electronics Society, IECON, 2024.

[23] M. Grosso, P. Riedinger, J. Daafouz, S. Pierfederici, H. Janati-Idrissi, and B. Lapôtre. Harmonic Control of Three-Phase AC/DC Converter With Time-Domain Guarantees. IEEE Trans. on Control Systems Technology, 33(3), 2025.

[24] A. Böttcher and B. Silbermann. Introduction to Large Truncated Toeplitz Matrices. Universitext. Springer New York, New York, NY, 1999.

[25] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan. Linear Matrix Inequalities in System and Control Theory. Society for Industrial and Applied Mathematics, 1994.

[26] M. Grosso, P. Riedinger, and J. Daafouz. The PhasorArray Toolbox. In *openly available on GitHub at github.com/mxmGrss/phasorArray_Toolbox*, 2025.

[27] M. Grosso, P. Riedinger, and J. Daafouz. The phasorarray toolbox for harmonic analysis and control design. In *submitted to European Control Conference, 2026 Island*, 2026.

[28] R. Teodorescu, F. Blaabjerg, M. Liserre, and P.C. Loh. Proportional-resonant controllers and filters for grid-connected voltage-source converters. *IEE Proceedings-Electric Power Applications*, 153(5):750–762, 2006.

[29] M. Liserre, R. Teodorescu, and F. Blaabjerg. Multiple harmonics control for three-phase grid converter systems with the use of PI-RES current controller in a rotating frame. *IEEE Trans. on Power Electronics*, 21(3):836–841, May 2006.

[30] G. Weiss and M. Häfele. Repetitive control of MIMO systems using $H_\infty$ design. Automatica, 35(7):1185–1199, 1999.

[31] M. Bodson, A. Sacks, and P. Khosla. Harmonic generation in adaptive feedforward cancellation schemes. IEEE Trans. on Automatic Control, 39(9):1939–1944, 1994.

[32] Mitra Mirhosseini. Modified park transformation of three-phase signals for harmonic extraction with frequency-adaptive capability. *International Journal of Electrical Power & Energy Systems*, 150:109069, 2023.