Abstract

Translate

A Bayes–Prony oscillating modal identification and hierarchical control strategy for low-frequency oscillation (LFO) of a ship microgrid (SM) is presented in this paper. The modal probabilistic estimation of the proposed algorithm replaces point estimation of the traditional Prony method and improves the robustness of modal identification. The hierarchical control strategy first performs modal identification by means of the batch least squares Prony (BLS-Prony) algorithm. The modal identification results are calibrated by the explanatory variance score (EVS), and the control process is transferred to recursive least squares Prony (RLS-Prony) real-time detection. The third layer of decision making transfers to Bayesian Prony (Bayes–Prony) identification in case of a loss of modality or failure of identification. The designed Bayes–Prony algorithm identifies the oscillatory modal of signals with a signal-to-noise ratio (SNR) equal to 2 dB. Compared to BLS-Prony and RLS-Prony, Bayes–Prony reduces the SNR convergence domain of the signal by 30 dB as the last layer of hierarchical control. Therefore, the third-layer decision commands are used as a scheduling reference for damping control in SM power plants. The proposed algorithms and strategies maximize the saving of computational resources while ensuring that the modal identification is effective. Finally, the correctness of the proposed algorithm and strategy is verified by the LFO waveforms of the experimental platform.

Full text

Turn on search term navigation

Translate

1. Introduction

The International Maritime Organization (IMO) is committed to supporting the United Nations Sustainable Development Goals and combating climate change by reducing greenhouse gas (GHG) emissions from shipping [1]. The main regulatory and enforcement measures taken by the International Maritime Organization since 2011 are summarized in Ref. [2] and projections to 2050 are provided. Based on the premise of the initiative, the all-electric ship (AES) represented by the integrated electric power system (IPS) has become the most promising technological path to realize the reduction in GHG emissions and green development in the shipping industry [3,4]. As a revolutionary technology in the field of advanced ship design and smart grid, the IPS not only helps the green transformation of SM [5], but more importantly, it realizes comprehensive improvement in layout and intelligence through the core idea of power integration and unified energy management.

The architecture of IPS and AES has brought enhancement to the development of SM. However, the limitations of ship size and environment result in an isolated microgrid with more severe stability challenges, for example, LFO and voltage dips caused by low damping and weak inertia [6,7,8]. The essential reason for stability problems in SM is due to the power imbalance between sources and loads at multiple time scales, whether low damping or weak inertia. In response, energy storage and converters have become the most widely used solution to the SM stability problem due to their flexible regulation of power density and energy density. For instance, the former is represented by flywheel [9], supercapacitor [10], and lithium battery [11]; the latter is represented by voltage source converter (VSC) [12] and virtual synchronous generator (VSG) [13].

Although the above energy storage devices have been widely used in microgrids, the contradiction between power density and energy density in SM cannot be reconciled. Therefore, LFO as an inherent problem of AC power systems still exists in SM and causes more serious grid stability problems. LFO belongs to the category of small disturbance stabilization of power systems and is usually analyzed by linearization methods [14,15]. In addition, the time-domain analysis method can also analyze the small disturbances in stability by using simulation or recorded wave data of the actual power system [16]. With the rise of wide-area measurement systems (WAMS) [17,18], the extraction of oscillatory modes by analyzing real-time recorded signals has also become a hot direction of attention for researchers.

In [19], Yutthagowith P. et al. proposed a fast curve-fitting technique applicable to lightning surge voltage and current reference curves. The technique is based on the least squares Prony analysis of the double integral of the recorded waveform to evaluate the waveform. The proposed technique provides a better match for the initial state of the voltage and current. In [20], Zhao. J et al. proposed a novel robust Prony method for power system oscillation pattern monitoring that is robust to synchronized signal measurement noise and outliers. It is implemented by constructing the weighted minimum total variance (WTLS) estimation problem based on the least squares Prony to eliminate the effect of measurement noise. In [21], researchers proposed a hybrid Prony model (HPM). By setting the L1 and L2 in the HPM for sparse expansion, the HPM can be converted to a separating HPM (SHPM). Similar to the Prony method, the nonlinear parameter estimation problem in SHPM can be converted into a problem of finding the roots of a polynomial. In [22], Khodaparast J. et al. examined a recursive multichannel Prony algorithm; the frequency discrimination performance of Prony is proved in the paper, but at the same time the authors present the challenge of the inability of the Prony model to converge under noisy conditions. Ultimately, the proposed algorithm is effective when the SNB is greater than 30 dB.

In [23], four typical approaches to fault diagnosis for IPS are described, including model-based, data-driven, knowledge-based, and hybrid-based approaches. Challenges in the above frameworks include difficulties in data acquisition and imbalance, limitations in sensing and onboard computational power, and reduced diagnostic reliability in the case of coupled faults. In [24], A. Mariscotti carried out a study on the problem of metrics of power quality in microgrids, proposing for pulsed loads the metrics of area, energy, and half-life in relation to the evaluation of power trajectories and their derivatives. In [25], Z. Du et al. proposed to integrate the load power parameters of the AES into the PMS. This method optimizes the energy balance to achieve power balance and frequency stability. In [26], M. Kaloev et al. proposed a methodology and framework for finding the parameters of algorithms by deep learning. The prior distribution of the parameters determines the Bayes–Prony stability of the proposed strategy, and deep learning provides a future research direction for the solution of the prior distribution.

In summary, the research and characteristics related to the study of Prony algorithms with application to power system states are summarized as follows:

(1). Since the Prony algorithm’s architecture is computationally intensive, it is usually used for offline reconstruction of recorded waveform data or historical data on long time scales.

(2). The rigidity feature of the Prony model makes it difficult for the solution to converge to the real mode when facing a noisy signal, but there are fewer studies on improving its robustness and anti-interference performance.

(3). The Prony algorithm in the literature focuses on modal identification of sparse signals, but the low inertia and weak damping characteristics of the SM make it possible to have multiple oscillatory modals in the system.

Thus, Figure 1 is a block diagram of the SM studied, and the work and contribution of this paper can be summarized as follows:

(1). A Bayes–Prony algorithm is proposed, which replaces the point estimation of the Prony algorithm with probabilistic estimation and solves the problem of the modal equations of noise-containing signals being unable to converge.

(2). The proposed strategy implements hierarchical decision making based on the sparsity and noise situation of the real signal. Among them, BLS-Prony and RLS-Prony are used to extract modals from high SNR signals, and Bayes–Prony is invoked as the last layer when their recognized modals cannot pass the calibration.

(3). The proposed strategy configures the corresponding invocation strategy based on the computational effort of the three algorithms. The PMS of the SM is invoked for Bayes computation only when RLS and BLS are disabled. The computational resources of the IPS are maximally saved.

The rest of this paper is organized as follows. Section 2 describes the solution principles and derivation of the three algorithms, Prony model, autoregressive model, BLS-Prony, and RLS-Prony. Section 3 describes the framework of the Bayes–Prony algorithm and the principle of MAP optimal estimation. In addition, the hierarchical decision control framework and process are also described in this section. In Section 4, an experimental platform for simulating SM was built, and the dynamic properties of the above three algorithms as well as the accuracy of modal extraction under different SNRs are verified, respectively. Finally, the work of this paper is summarized and conclusions are made about the research in this paper in Section 5.

2. Description of the Modal and Identification Algorithms

The study of LFO in microgrids as a small disturbance stability implies that the power oscillations will always decay and converge to a steady state with a corresponding time constant. The principle of modal identification can be summarized as follows: (1) The actual oscillating waveforms are sampled to form discrete time series, and these data contain information about the waveforms of the actual power oscillations. (2) Finding the autoregressive relationship between the discrete data will allow the corresponding modal information to be extracted from this relationship.

2.1. Model of the Oscillatory Modal and Autoregressive Prediction

Firstly, the oscillating signal in the time domain is decomposed into a superposition of multiple decaying sinusoidal components according to the Prony model. The analytical equation is as follows

(1) $y (t) = \sum_{i = 1}^{p} A_{i} \cdot \exp (σ_{i} t) \cdot \cos (ω_{i} t + ϕ_{i})$

where

i

is the

i

th modal.

A_{i}

and

ϕ_{i}

are the phase and amplitude, respectively.

σ_{i}

is the attenuation coefficient.

ω_{i}

is the angular frequency of the oscillation.

p

stands for the order of the oscillatory modes contained in

y (t)

. Equation (1) is rewritten in complex form using Euler’s Equation as follows:

(2) $\cos (ω_{i} t + ϕ_{i}) = \frac{1}{2} {\exp [j (ω_{i} t + ϕ_{i})] + \exp [- j (ω_{i} t + ϕ_{i})]}$

After substituting (2) into (1), (1) can be simplified as follows:

(3) $y (t) = \sum_{i = 1}^{p} {B_{i} \cdot \exp [(σ_{i} + j ω_{i}) t] + {\bar{B}}_{i} \cdot \exp [(σ_{i} - j ω_{i}) t]}$

where

B_{i}

and

{\bar{B}}_{i}

are the positive and negative order components of the modes, respectively.

y [n]

is a discretized time series formed after

y (t)

is sampled at a fixed time interval, and

y [n]

is represented as follows:

(4) $y [n] = \sum_{i = 1}^{p} (B_{i} z_{i}^{n} + {\bar{B}}_{i} {\bar{z}}_{i}^{n}), (n = 0, 1, 2, \dots, N)$

where

n

is the time index and represents the

n

th sampling point.

z_{i}

and

{\bar{z}}_{i}

are conjugate poles, and their expressions are as follows:

(5) $z_{i} = \exp [σ_{i} \cdot Δ t] \cdot \exp (j ω_{i} \cdot Δ t) \cdot {\bar{z}}_{i} = \exp [σ_{i} \cdot Δ t] \cdot \exp (- j ω_{i} \cdot Δ t)$

Substituting (5) into (4), the Prony model of the discrete time series of oscillatory modes is obtained.

(6) $y [n] = 2 \cdot Re [\sum_{i = 1}^{p} (B_{i} z_{i}^{n})], (n = 0, 1, 2, \dots, N)$

Secondly, an autoregressive (AR) model is constructed to describe the functional relationship in $y [n]$ , because $y [n]$ contains the causal relationship of $y (t)$ . The AR model is described as follows:

(7) $x [n] = - a_{1} x [n - 1] - a_{2} x [n - 2] - a_{3} x [n - 3] \Rightarrow x [n] = - \sum_{k = 1}^{m} (a_{k} x [n - k]) + e [n]$

where

x [n]

is the measured value at the current moment and also the value that needs to be predicted based on the AR model.

m

is the order of the AR model.

a_{k}

stands for the autoregressive coefficient.

x [n - k]

is the signal value at the

k

th moment in the past.

In order to prove the equivalence of (6) and (7), a zeroed polynomial is constructed to convert (6) into a difference equation.

(8) $\prod_{i = 1}^{p} (1 - z_{i} \cdot λ_{i}) (1 - {\bar{z}}_{i} \cdot λ_{i}) \cdot y [n] = 0$

where

Δ_{i}

and

{\bar{Δ}}_{i}

are zeroed polynomials, and

Δ_{i} = 1 - z_{i} \cdot λ_{i}

{\bar{Δ}}_{i} = 1 - {\bar{z}}_{i} \cdot λ_{i}

λ_{i}

is a delay operator and

λ_{i} = y [n - i] / y [n]

. Expanding Equation (8) as follows,

(9) $y [n] = - b_{1} y [n - 1] - b_{2} y [n - 2] \dots - b_{2 p - 1} y [n - 2 p + 1] - b_{2 p} y [n - 2 p]$

Rewriting the above equation into compact form gives

(10) $y [n] = - \sum_{τ = 1}^{2 p} b_{τ} (y [n - τ])$

Equation (10) proves that the Prony model shows AR properties itself after constructing the difference equation. Comparing Equations (7) and (10), both are equivalent while ignoring the error

e [n]

. Note that the number of orders in the AR model and the number of modes in the Prony model satisfy the relation

m = 2 p

. The modal resolution problem now transforms into a problem of solving for the autoregressive coefficients and orders of the AR model.

2.2. BLS-Prony Method Extraction of Oscillatory Modal

Now, the core problem of this paper becomes the solution of Equation (10). To ensure that the time indices in Equation (10) are all positive, let $\min (n) = m$ . As a result, the number of equations constructed is $N - m + 1$ . However, the real sample data $N$ is much larger than the oscillating modal order $m$ , so $N - m + 1 ≫ m$ is necessarily satisfied. In other words, Equation (10) must be over-determined as follows:

(11) $[\begin{matrix} y [m] \\ y [m + 1] \\ ⋮ \\ y [N] \end{matrix}] = - [\begin{matrix} y [m - 1] & y [m - 2] & \dots & y [0] \\ y [m] & y [m - 1] & \dots & y [1] \\ ⋮ & ⋮ & ⋱ & ⋮ \\ y [N - 1] & y [N - 2] & \dots & y [N - m] \end{matrix}] [\begin{matrix} b_{1} \\ b_{2} \\ ⋮ \\ b_{m} \end{matrix}] \Rightarrow y = - Y b$

where

Y

is the regression vector matrix.

b

is the vector of coefficients of the AR.

y

is the observation vector. The BLS is constructed from Equation (11) and solved for

b

that satisfies the minimum residual sum.

(12) $\min J (b) = | | y + Y b | |^{2}$

where

J (b)

is the objective function of BLS. Expanding Equation (12), we solve for its gradient as follows:

(13) $\{\begin{cases} J (b) = y^{T} y + 2 y^{T} Y b + b^{T} Y^{T} Y b \\ \nabla_{b} J (b) = 2 Y^{T} Y b + 2 Y^{T} y \end{cases}$

Let

\nabla_{b} J (b) = 0

to obtain the formula for

b

\nabla_{b} J (b)

(14) $b = - {(Y^{T} Y)}^{- 1} Y^{T} y$

Now, reviewing the problem posed at the beginning of this section, we have finished solving B and will next determine the order of the AR model.

Step 1: Construct the Hankel matrix of $y [n]$ and decompose it into singular values.

(15) $H = [\begin{matrix} y [0] & y [1] & \dots & y [L] \\ y [1] & y [2] & \dots & y [L + 1] \\ ⋮ & ⋮ & ⋱ & ⋮ \\ y [M] & y [M + 1] & \dots & y [N] \end{matrix}] = U Σ V^{T} = \sum_{i = 1}^{r} σ_{i} u_{i} v_{i}^{T}$

where

H

is the Hankel matrix of

y [n]

U

is an orthogonal matrix, which contains the left singular vectors.

U

is also an orthogonal matrix, which contains the right singular vectors.

r

is the rank of

H

Σ

is a diagonal matrix and its main diagonal elements are not equal to zero.

σ_{i}

u_{i}

, and

v_{i}

represent singular values, left singular unit vectors, and right singular unit vectors, respectively.

$σ_{i}$ in descending order can be expressed as

$σ_{1} \geq σ_{2} \dots \geq σ_{r} \geq σ_{r + 1} = \dots = σ_{\min (M + 1, L + 1)} = 0 .$

The first m singular values from the above equation are truncated and these singular values correspond to the main features. The truncated approximation matrix is given by

(16) $H_{m} = U_{m} Σ_{m} V_{m}^{T} = \sum_{i = 1}^{m} σ_{i} u_{i} v_{i}^{T}$

Step 2: Combining Equations (15) and (16), the remaining equation after the interception can be expressed as follows:

(17) $H - H_{m} = \sum_{i = 1}^{r} σ_{i} u_{i} v_{i}^{T} - \sum_{i = 1}^{m} σ_{i} u_{i} v_{i}^{T} = \sum_{i = m + 1}^{r} σ_{i} u_{i} v_{i}^{T}$

We express the above residuals in terms of energy as follows:

(18) $E_{H - H_{m}} = | | H - H_{m} | |_{F}^{2} = \sum_{i = m + 1}^{r} σ_{i}^{2} | | u_{i} v_{i}^{T} | |_{F}^{2}$

Similarly, according to the above equation, the energies of A and B can be introduced separately as

(19) $E_{H} = \sum_{i = 1}^{r} σ_{i}^{2}, E_{H_{m}} = \sum_{i = 1}^{m} σ_{i}^{2}$

The singular values correspond to the features of the matrix, and the square of the singular values reflects the generalized energy of the matrix. In order to ensure that the main features of the matrix remain unchanged before and after truncation, the energy of the residual matrix needs to be less than 5%. Therefore, the order-reliant energy constraints are designed as

(20) $ρ = \frac{E_{H_{m}}}{E_{H}} = \sum_{i = 1}^{m} σ_{i}^{2} / \sum_{i = 1}^{r} σ_{i}^{2} \geq 95 %$

Now, the solution of the AR model has been completed. The oscillatory modals will be solved next.

Step 3: Expand Equation (10) as follows.

(21) $y [n] + b_{1} y [n - 1] + \dots + b_{m} y [n - m] = 0$

Substituting Equation (4) as a generalized solution into Equation (21),

(22) $\{\begin{cases} \sum_{i = 1}^{p} (B_{i} z_{i}^{n}) + b_{1} \sum_{i = 1}^{p} (B_{i} z_{i}^{n - 1}) + b_{2} \sum_{i = 1}^{p} (B_{i} z_{i}^{n - 2}) + \dots + b_{m} \sum_{i = 1}^{p} (B_{i} z_{i}^{n - m}) = 0 \\ \sum_{i = 1}^{p} ({\bar{B}}_{i} {\bar{z}}_{i}^{n}) + b_{1} \sum_{i = 1}^{p} ({\bar{B}}_{i} {\bar{z}}_{i}^{n - 1}) + b_{2} \sum_{i = 1}^{p} ({\bar{B}}_{i} {\bar{z}}_{i}^{n - 2}) + \dots + b_{m} \sum_{i = 1}^{p} ({\bar{B}}_{i} {\bar{z}}_{i}^{n - m}) = 0 \end{cases}$

For the above equation to be satisfied in

\forall n

, the factorization of the equation must be equal to zero, such that

(23) $\{\begin{cases} z_{i}^{m} + b_{1} z_{i}^{m - 1} + b_{2} z_{i}^{m - 2} + \dots + b_{m} = 0 \\ {\bar{z}}_{i}^{m} + b_{1} {\bar{z}}_{i}^{m - 1} + b_{2} {\bar{z}}_{i}^{m - 2} + \dots + b_{m} = 0 \end{cases}, (i = 1, 2, \dots, p)$

The above equation shows that the modal information parameters

z_{i}

and

{\bar{z}}_{i}

are exactly the roots of the pole equation. We ignore the repeated conjugate results, since

+ ω_{i}

and

- ω_{i}

represent the same physical frequency. Ultimately, mapping the conjugate pole

(z_{i}, {\bar{z}}_{i})

to the modal parameters according to Equation (23), the resolved modes are as follows:

(24) $\{\begin{cases} σ_{i} = Re [\ln z_{i} / Δ t] \\ ω_{i} = Im [\ln z_{i} / Δ t] / 2 π \end{cases}$

In addition, the phases and amplitudes of the modes are included in B and B. A system of equations is constructed from Equation (4), and the information contained therein is resolved.

(25) $[\begin{matrix} y [0] \\ ⋮ \\ y [N] \end{matrix}] = [\begin{matrix} \exp (α_{1} t_{0}) & \dots & \exp (α_{P} t_{0}) \\ ⋮ & ⋱ & ⋮ \\ \exp (α_{1} t_{N}) & \dots & \exp (α_{P} t_{N}) \end{matrix}] [\begin{matrix} B_{1} \\ ⋮ \\ B_{p} \end{matrix}] + [\begin{matrix} \exp (- α_{1} t_{0}) & \dots & \exp (- α_{P} t_{0}) \\ ⋮ & ⋱ & ⋮ \\ \exp (- α_{1} t_{N}) & \dots & \exp (- α_{P} t_{N}) \end{matrix}] [\begin{matrix} {\bar{B}}_{1} \\ ⋮ \\ {\bar{B}}_{p} \end{matrix}]$

The solutions to the above equations are mapped to amplitude and phase as follows:

(26) $\{\begin{cases} A_{i} = 2 | B_{i} | \\ ϕ_{i} = | \ln (2 B_{i} / A_{i}) | \end{cases}$

Equations (24) and (26) are the results of the analysis of the oscillatory modes based on BLS-Prony.

2.3. RLS-Prony Method Extraction of Oscillatory Modal

The two previous sections establish the Prony and AR model, and for the problem above we propose the BLS-Prony algorithm. However, this algorithm asks for the precondition that there are discrete sampling samples, and the identification process is based on N historical sampling values that have been obtained. This implies that the recognition process is not only offline, but also requires the storage of all or batches of historical data. As a result, this places a higher demand on the computational resources of the chip. Therefore, in this section, online identification for dynamic updating of model parameters is realized by the recursive least squares method.

First, rewrite Equation (7) in the form of a matrix as follows:

(27) $y [n] = - Y {(n)}^{T} a$

where

Y (n) = {(y [n - 1], y [n - 2], \dots, y [n - n])}^{T}

a = {(a_{1}, a_{2}, \dots, a_{n})}^{T}

y [n]

is the value of the signal observed at the current moment,

Y (n)

represents the regression vector, and

a

is the vector of coefficients to be obtained. Similarly, to ensure that the indexing of historical data is valid, the minimum number of indexes is equal to zero and the condition

n \geq m

is satisfied. To construct a weighted least squares function, and to find

a

such that the value of the objective function is minimized,

(28) $J (a) = \sum_{i = m}^{n} λ^{n - i} {(y [i] + Y {(i)}^{T} a)}^{2}$

where

J (a)

is the constructed weighted least squares objective function.

λ

is the weighting factor, and

0 < λ < 1

. Note that Equations (12) and (28) are equivalent when

λ = 1

. Then, let

\nabla_{a} J (a)

be equal to zero and solve for

a

(29) $\nabla_{a} J (a) = 2 \sum_{i = m}^{n} λ^{n - i} Y (i) (y [i] + Y {(i)}^{T} a) = 0$

Equation (29) can be further simplified and organized as

$\sum_{i = m}^{n} λ^{n - i} Y (i) Y {(i)}^{T} a = - \sum_{i = m}^{n} λ^{n - i} Y (i) y [i] \Rightarrow R (n) a = r (n)$

where

R (n)

is the autocorrelation matrix,

r (n)

is the autocorrelation matrix, and

R (n)

and

r (n)

encode the modal information of

y [n]

. Thus, RLS achieves parameter estimation by recursively updating

R (n)

and

r (n)

. Their recursive equations are as follows.

(30) $\begin{array}{l} R (n) = \sum_{i = n} λ^{n - i} Y (i) Y {(i)}^{T} + λ \sum_{i = m}^{n - 1} λ^{(n - 1) - i} Y (i) Y {(i)}^{T} = Y (n) Y {(n)}^{T} + λ R (n - 1) \\ r (n) = - \sum_{i = n} λ^{n - i} Y (i) y [i] + λ (- \sum_{i = m}^{n - 1} λ^{(n - 1) - i} Y (i) y [i]) = - Y (n) y [n] + λ r (n - 1) \end{array}$

If Equation (30) is substituted into Equation (20), the inverse matrix of

R (n)

has to be solved directly. The computational effort of direct inversion of a matrix is enormous. Therefore, the direct inverse is avoided by transforming the Sherman–Morrison formula to its low-rank correction. The computational effort of direct inversion of a matrix is enormous. Consequently, the direct inverse is avoided by transforming the Sherman–Morrison formula to its low-rank correction.

(31) $R^{- 1} (n) = \frac{1}{λ} R {(n - 1)}^{- 1} - \frac{\frac{1}{λ} R {(n - 1)}^{- 1} Y (n) Y {(n)}^{T} \frac{1}{λ} R {(n - 1)}^{- 1}}{1 + Y {(n)}^{T} \frac{1}{λ} R {(n - 1)}^{- 1} Y (n)}$

Let

R {(n)}^{- 1} = P (n) \Rightarrow R {(n - 1)}^{- 1} = P (n - 1)

. Rearranging Equation (10) yields the following form:

(32) $P (n) = \frac{1}{λ} (P (n - 1) - K (n) Y {(n)}^{T} P (n - 1))$

where

K (n)

is the Kalman gain vector, such that

$K (n) = \frac{P (n - 1) Y (n)}{λ + Y {(n)}^{T} P (n - 1) Y (n)}$

After substituting Equation (32) into Equation (29), the final expression for A is obtained by organizing and simplifying as

(33) $\begin{array}{l} a (n) = \frac{(P (n - 1) - K (n) Y {(n)}^{T} P (n - 1)) (λ r (n - 1) - Y (n) y [n])}{λ} \\ = a (n - 1) - K (n) e [n] \end{array}$

In the previous section, BLS was used to construct the Hankel for SVD decomposition based on the batch historical data, whose singular value decay points correspond to the order of the AR. RLS cannot directly construct the Hankel matrix for estimation due to the dynamic iteration of the data. Therefore, the SVD decomposition is performed through $P {(n)}^{- 1}$ , and the steep decay point of the singular values is judged to be the order of the AR. The process of order determination is as follows.

(34) $P {(n)}^{- 1} = λ P {(n - 1)}^{- 1} + Y (n) Y {(n)}^{T} \approx \sum_{i = 1}^{n} λ^{n - i} Y (i) Y {(i)}^{T}$

In Equation (2),

P {(n)}^{- 1}

is the information matrix containing all of the information of the Hankel. Therefore, we are looking for a sharp decay point of the

P {(n)}^{- 1}

singular value.

(35) $P {(n)}^{- 1} = U Σ V^{T}$

where

Σ = d i a g [\begin{matrix} σ_{1} & σ_{2} & \dots & σ_{m_\max} \end{matrix}]

. The decomposed eigenvalues are selected as the smallest according to the principle of energy accumulation so that the lost energy needs to be less than 5%.

(36) $ρ = \sum_{i = 1}^{m} σ_{i} / \sum_{i = 1}^{m_\max} σ_{i} \geq 95 %$

Finally, the additional threshold parameters A and B are used to dynamically adjust the order.

(37) $\{\begin{cases} σ_{m + 1} \geq δ σ_{1}, m = m + 1 \\ σ_{m} \leq γ σ_{1}, m = m - 1 \end{cases}$

Equation (37) requires an increase in order when the

m + 1

singular value is detected to be significantly larger than the maximum singular value

δ

times. Conversely, the order needs to be decreased.

3. Bayes–Prony Method for Probabilistic Identification

In the previous section, BLS-Prony and RLS-Prony are constructed to realize the identification of modal parameters, respectively. Both algorithms require the same prerequisites for an effective implementation: (1) the signal must be the sum of pure complex exponentials; (2) the frequency and attenuation factor are constant over the sampling time; (3) the components are independent and do not interact with each other. In other words, the above conditions make these two algorithms less robust and resistant to interference. In order to enhance the robustness and anti-interference of the above algorithms, the probabilistic identification feature of the Bayes–Prony algorithm is designed to enhance the performance of the algorithm.

3.1. Bayes Model and Posterior Probability Distribution Design

$y (t)$ introduces a noise term to Equation (1) and the new modeling is as follows:

(38) $y (t) = \sum_{i = 1}^{p} A_{i} \cdot \exp (σ_{i} t) \cdot \cos (2 π f_{i} t + ϕ_{i}) + ε (t)$

where

ε (t)

is independent identically distributed Gaussian noise. We define the expression for parameter

θ

in the system as follows:

$θ = [\begin{matrix} f_{1}, & \dots, & f_{m}, & σ_{1}, & \dots, & σ_{m}, & A_{1}, & \dots, & A_{m}, & ϕ_{1}, & \dots, & ϕ_{m}, & α^{2} \end{matrix}]$

According to Bayes’ theorem, the posterior probability is proportional to the product of the likelihood function and the prior probability, such that

(39) $P (θ | y) \propto P (y | θ) \cdot P (θ)$

where

y = [\begin{matrix} y [0], & \dots, & y [N] \end{matrix}]

;

P (θ | y)

is the likelihood function of the noise, which is expressed as follows:

(40) $\begin{array}{l} P (y | θ) = \prod_{n = 0}^{N} \frac{1}{\sqrt{2 π α^{2}}} \exp [- \frac{{(y [n] - \overset{⌢}{y} [n; θ])}^{2}}{2 α^{2}}] \\ = P (y | θ) = {(\frac{1}{2 π α^{2}})}^{(N + 1) / 2} \exp [- \frac{1}{2 α^{2}} \sum_{n = 0}^{N} {(y [n] - \overset{⌢}{y} [n; θ])}^{2}] \end{array}$

where

\overset{⌢}{y} [n; θ]

is the reconfiguration signal. The prior distributions of the remaining parameters are each designed as follows:

$P (f) = \frac{1}{\sqrt{2 π α_{f}^{2}}} \exp [- \frac{{(f - μ_{f})}^{2}}{2 α_{f}^{2}}],$

$P (σ) = \frac{1}{\sqrt{2 π α_{σ}^{2}}} \exp [- \frac{{(σ - μ_{σ})}^{2}}{2 α_{σ}^{2}}], σ < 0$

$P (A) = \frac{λ}{2} \exp (- λ | A |),$

$P (ϕ) = \frac{1}{2 π}, ϕ \in [- π, π],$

$P (α^{2}) = \frac{β_{ε}^{α_{ε}}}{Γ (α_{ε})} {(α^{2})}^{- α_{ε} - 1} \exp (- \frac{β_{ε}}{α^{2}})$

where

P (f)

and

P (σ)

are Gaussian prior distribution for

f

and

σ

;

P (A)

is the Laplace-distributed prior probability of the magnitude;

P (ϕ)

is the continuously uniformly distributed prior probability of the phase; and

P (α^{2})

is the inverse Gamma distributed prior probability of the noise variance. The joint prior probability distribution for each of the above parameters is obtained by combining the prior probability distributions.

(41) $\begin{array}{l} P (θ | y) \propto P (y | θ) \cdot P (θ) \\ \propto {{(\frac{1}{2 π α^{2}})}^{\frac{N + 1}{2}} \exp [- \frac{1}{2 α^{2}} \sum_{n = 0}^{N} {(y [n] - \overset{⌢}{y} [n; θ])}^{2}]} \times [\prod_{n = 0}^{N} P (f) \cdot P (σ) \cdot P (A) \cdot P (ϕ)] \cdot P (α^{2}) \\ \propto {(\frac{1}{2 π α^{2}})}^{\frac{N + 1}{2}} \times \exp [- \frac{1}{2 α^{2}} \sum_{n = 0}^{N} {(y [n] - \overset{⌢}{y} [n; θ])}^{2}] \times [\prod_{n = 0}^{N} {\exp [- \frac{{(f - μ_{f})}^{2}}{2 α_{f}^{2}}] \\ \times \exp [- \frac{{(σ - μ_{σ})}^{2}}{2 α_{σ}^{2}}] \times \exp (- λ | A |) \times C_{1}]} \times {(α^{2})}^{- α_{ε} - 1} \exp (- \frac{β_{ε}}{α^{2}}) \end{array}$

The final form of the posterior distribution is obtained by simplifying and organizing the above equation.

(42) $\begin{matrix} P (θ | y) \propto \exp \{- \frac{N + 1}{2} \log (α^{2}) - \frac{1}{2 α^{2}} \sum_{n = 0}^{N} {(y [n] - \overset{⌢}{y} [n; θ])}^{2} \\ - \sum_{n = 0}^{N} [\frac{{(f - μ_{f})}^{2}}{2 α_{f}^{2}} + \frac{{(σ - μ_{σ})}^{2}}{2 α_{σ}^{2}} + λ | A |] - (α_{ε} + 1) \log (α^{2}) - \frac{β_{ε}}{α^{2}}\} \end{matrix}$

The prior distribution in Equation (41) is a key piece of the Bayesian approach in the above analysis. The essential reason for this is that Bayes’ theorem on the fact that the posterior is proportional to the prior, while the prior is proportional to the likelihood function. Therefore, the modal frequency prior has zero probability in the defined region, and then its posterior probability remains zero despite the fact that the data exhibits strong evidence in that region. The discussion of prior distributions will have to be further developed on the basis of this study.

3.2. Maximum Posteriori Probability Estimation and Optimization for Solution Design

Equation (42) transforms the point estimation of the Prony model into a probabilistic or the parameters. Then the maximum point of the posterior probability distribution will be searched for, which is the probability parameter with the highest confidence level.

(43) ${\hat{θ}}_{M A P} = \arg \max_{θ} [P (θ | y)] \Leftrightarrow {\hat{θ}}_{M A P} = \arg \min_{θ} {- \log [P (θ | y)]}$

where

{\hat{θ}}_{M A P}

is the maximum estimated parameter of

P (θ | y)

. Defining the objective function,

J (θ)

is of the following form

(44) $\begin{array}{l} J (θ) = - \log [P (θ | y)] = \frac{N + 1}{2} \log (α^{2}) + \frac{1}{2 α^{2}} \sum_{n = 0}^{N} {(y [n] - \overset{⌢}{y} [n; θ])}^{2} \\ + \sum_{n = 0}^{N} [\frac{{(f - μ_{f})}^{2}}{2 α_{f}^{2}} + \frac{{(σ - μ_{σ})}^{2}}{2 α_{σ}^{2}} + λ | A |] + (α_{ε} + 1) \log (α^{2}) + \frac{β_{ε}}{α^{2}} \\ ≜ I + I I + I I I + I V \\ ≜ J_{s} (θ) + J_{n s} (θ) \end{array}$

$I = 1 / 2 α^{2} \cdot \sum_{n = 0}^{N} {(y [n] - \overset{⌢}{y} [n; θ])}^{2}$

$I I = \sum_{n = 0}^{N} [{(f - μ_{f})}^{2} / 2 α_{f}^{2} + {(σ - μ_{σ})}^{2} / 2 α_{σ}^{2}]$

$I I I = \sum_{n = 0}^{N} (λ | A |)$

$\frac{N + 1}{2} \log (α^{2}) + (α_{ε} + 1) \log (α^{2}) + \frac{β_{ε}}{α^{2}}$

In the above equation, I contains trigonometric functions that cause the Hessian matrix to be non-positive definite and make A non-convex in nature. As a result, the gradient descent converges to a local extreme, corresponding to some secondary mode and omitting the main. In addition, III is the L1 regularization term, and it cannot be used directly with gradient descent due to nontrivial zeros. Combining the above properties, $J (θ)$ can be categorized into smooth and non-smooth terms, with the smooth term $J_{s} (θ) = I + I I + I V$ and the non-smooth term $J_{n s} (θ) = I I I$ .

For the former, $J_{s} (θ)$ can directly compute its gradient due to its full-domain differentiability, such that

(45) $\nabla J_{s} (θ) = [\begin{matrix} \frac{\partial J_{s}}{\partial A} & \frac{\partial J_{s}}{\partial σ} & \frac{\partial J_{s}}{\partial f} & \frac{\partial J_{s}}{\partial ϕ} & \frac{\partial J_{s}}{\partial α^{2}} \end{matrix}]$

where the parameter update rule can be expressed as

(46) $θ [t + 1] = θ [t] - η \nabla J_{s} (θ [t])$

η

is the learning efficiency, expressed as the rate of convergence of

θ

For the latter, a proximal gradient descent can be made for A since its zeros in $J_{n s} (θ)$ are not trivial, such that

(47) $A [t^{'}] = A [t] - η \frac{\partial \nabla J_{s} (θ [t])}{\partial A}, A [t + 1] = S_{η λ} (A [t^{'}])$

where

A [t^{'}]

is the pre-updated value.

S_{η λ}

is the soft threshold operator, and its expression is

S_{η λ} (x) = s i g n (x) \cdot \max [(| x | - η λ), 0]

In addition, for either $J_{s} (θ)$ or $J_{n s} (θ)$ , there are parametric constraints on this objective function, e.g., $σ > 0$ , $f > 0$ , and $α^{2} > 0$ , depending on the physical conditions constructed by the Prony model. In the above optimization solving process, if the constraints are not considered, the final result may converge to meaningless intervals. Therefore, the meaningless intervals are projected to the feasible domain during the solution process. The updating formula for the above parameters can be modified as follows

(48) $θ [t + 1] (\cdot) = P_{C} {θ [t] (\cdot) - η \nabla J_{s} (θ [t] (\cdot))}, (\cdot) = f, σ, α^{2}$

where

P_{C}

is the projection operator of the constraint set

C = {x : x \leq 0}

3.3. Hierarchical Control Strategy Design

In the previous section, the modal information of A was extracted by the Bayes–Prony method. Then, the extracted modal information was judged for system hierarchical decision making. The time complexity of the three modal identification algorithms mentioned above needs to be analyzed before starting hierarchical control design.

BLS-Prony takes the N data points collected and finds the best-fitting parameters all at once. According to (14), the total computation of BLS is $o (m N^{2} + N^{3})$ . RLS-Prony avoids inverting the matrix directly and updating the covariance matrix recursively. Therefore, its complexity is at most $o (N^{2})$ . The time complexity of Bayes–Prony is similar to that of BLS-Prony, but the objective function is more complex to solve, with the complexity associated with computing the likelihood function and gradient at each iteration. In summary, the former algorithms can be recognized at the local control level of the energy storage system, and when the local recognition calibration fails, it is necessary to transfer to the PMS for Bayes–Prony identification. Based on the above control theory, if the proposed strategy is used in a large SM, the necessary condition is that the ESS can directly obtain the generator speed. This is an important guarantee for subsequent damping control.

The control chip is selected as TMS32-DSP28379. The main frequency of the DSP is 200 MHz. The number of sampling points $N$ is equal to 50. Therefore, the total execution time for the DSP to realize an operation of $o (N^{3})$ is equal to 31.2 ms. In addition, process flows such as data handling, loop execution, and ADC data reading result in a final execution time of approximately 50 ms. The frequency range of the LPO is equal to 0.2–2 Hz, so the 500 ms period ensures the accuracy of local identification. PMS and local communication delay occur when Bayesian control is activated. The communication delay can be ignored relative to the time scale of the LFO.

Next, the design process of the proposed hierarchical control strategy is as follows.

(49) $e [n] = y [n] - \overset{⌢}{y} (t) = y [n] - \sum_{m = 1}^{M} {\overset{⌢}{h}}_{m} {\overset{⌢}{z}}_{m}^{n} ≜ e = y - \overset{⌢}{s}$

where

{\overset{⌢}{h}}_{m} = {\overset{⌢}{A}}_{i} \cdot \exp (j {\overset{⌢}{ϕ}}_{i}) / 2

;

z_{m} = \exp [({\overset{⌢}{σ}}_{i} + j {\overset{⌢}{ω}}_{i}) \cdot Δ t]

;

e = {[e [0], e [1], \dots, e [N]]}^{T}

y = {[y [0], y [1], \dots, y [N]]}^{T}

;

\overset{⌢}{s} = {[\overset{⌢}{s} [0], \overset{⌢}{s} [1], \dots, \overset{⌢}{s} [N]]}^{T}

Next, the statistical properties of the residuals are obtained, containing the residual mean, residual standard deviation, and residual variance, and the above statistics are calculated as follows:

$\bar{e} = \frac{1}{N + 1} \sum_{n = 0}^{N} e [n], {\overset{⌢}{σ}}_{e} = \sqrt{\frac{1}{N} \sum_{n = 0}^{N} {(e [n] - \bar{e})}^{2}}, {\overset{⌢}{σ}}_{e}^{2} = \frac{1}{N} \sum_{n = 0}^{N} {(e [n] - \bar{e})}^{2}$

where

\bar{e}

is the overall mean of the residuals;

{\overset{⌢}{σ}}_{e}

represents the sample standard deviation of the residuals; and

{\overset{⌢}{σ}}_{e}^{2}

represents the sample variance of the residuals. The above three statistical variables have been able to describe the accuracy of

y [n]

; in order to be able to further quantify the confidence of

y [n]

and use the results as a basis for judging the next level of decision making, it is also necessary to introduce an energy benchmark for the system. In summary, the following definitions are made:

$P_{y} = {\overset{⌢}{σ}}_{y}^{2} = \frac{1}{N} \sum_{n = 0}^{N} {(y [n] - \bar{y})}^{2}, \bar{y} = \frac{1}{N + 1} \sum_{n = 0}^{N} y [n]$

where

{\overset{⌢}{σ}}_{y}^{2}

is the sample variance representing

y [n]

;

P_{y}

is the signal energy of

y [n]

, and

\bar{y}

is the overall mean of

y [n]

. Assuming that the above statistical parameters are accessible, the explained variance score (EVS) is introduced as a measure of the accuracy of the fit after the Prony/AR model reconstruction. The EVS is computed as follows:

(50) $E V S = 1 - \frac{P_{e}}{P_{y}} = 1 - \frac{{\overset{⌢}{σ}}_{e}^{2}}{{\overset{⌢}{σ}}_{y}^{2}}$

According to Equation (50), EVS has a range of values $[0, 1]$ which is a direct measure of the proportion of the original signal fluctuations that can be captured by the reconstructed model. An EVS equal to 1 indicates that the residual power is minimized and the reconstructed model is able to explain all of the variance of the signal.

Figure 2 shows the modal identification and hierarchical control framework. The algorithm and strategy are triggered when the speed of the generator deviates.

Layer 1: BLS-Prony first performs signal acquisition and then performs modal order determination when N satisfies the set conditions. Note that there are two sets of raw sample data in this process, the sampled data preprocessed by filtering and noise reduction is used to set the modal order, and the preprocessed data is used to parse the modal information. The reason for solving for modal information with preprocessed samples is that the noise reduction and filtering aspects introduce lags and amplitude variations in the signal. These variations can cause Prony to fail to converge to the true modal state. The credibility of the modal information is verified by EVS. If the modal information is credible, the parameters are configured to the RLS-Prony layer for real-time monitoring; if the EVS calibration fails, it directly skips the second layer to enter the third layer of Bayes–Prony identification. The reason for skipping the RLS-Prony layer is that both BLS and RLS are solved on the basis of the Prony model, and thus both have unity for sample data that do not converge.

Layer 2: Modalities that are successfully calibrated at the first layer are continuously monitored. If the modal information passes the EVS calibration, the modal information will be fed to the control system output. If the modal is lost alive during the process and fails the EVS calibration, it will jump to the Bayes–Prony layer.

Layer 3: Both RLS and BLS need to jump to this layer for probability estimation when they cannot pass the EVS calibration. During the operation of this layer, the prior distribution parameters of the model as well as the model order are configured by the BLS layer when the wide-area identification is first triggered. If the recognition result passes the EVS calibration, the modal information is used as the hierarchical control output signal and the recognition process is exited at the same time. If the recognition result fails to pass the EVS calibration, the status is reported to the upper-level power management system.

Now, the theoretical derivation of modal identification and hierarchical decision control of low-frequency oscillations is completed, and the final hierarchical decision control signal will be sent to the management of the power station for damping control. The energy storage system (ESS) will start when the power management system (PMS) of the SM receives the LFO event trigger command. The ESS starts to detect the rotational speed of the SG and exchanges power from the grid through charging and discharging. The value of the electrical energy exchanged between the pure energy system and the SM is equal to $D Δ ω$ , where $D$ is the damping coefficient and $Δ ω$ is the change in rotational speed of the SG.

4. Result and Discussion

In this section, the dynamic performance of the proposed discrimination algorithm with hierarchical decision control is validated. Firstly, the interference resistance as well as the robustness of the recognition algorithm is verified in simulation mode. Secondly, identification and decision making are carried out based on real low-frequency oscillation data sampled from the experimental platform. Figure 3 shows the experimental platform built according to the description of Figure 1.

4.1. Dynamic Response Analysis of BLS-Prony

Figure 4 shows the Prony model order identification process for a single modality. The samples in the figure have a sampling frequency of $f_{s} = 100 Hz$ , and the sampling interval is $Δ T = 0.01 s$ . The dimensions of the Hankel matrix are set to $125 \times 376$ . Singular values and the order criterion are obtained after SVD decomposition of the Hankel matrix. Figure 4a,b are the normalized singular values and $σ_{3}$ has converged to zero, so the order of the AR model is equal to 2. Figure 4c,d show the cumulative energy function used to truncate the small singular values and select the smallest $m$ for which the energy share reaches a threshold greater than 95%. Figure 5 shows the comparison between the unimodal Prony reconstruction and the input signal with the characteristic pole distribution. The setup conditions for the algorithm are as follows:

$y (t) = A \exp (σ t) \cdot \cos (2 π f t + ϕ) \Rightarrow y (t) = 1.5 \exp (- 0.5 t) \cdot \cos (2 π 0.5 t + π / 3)$

y (t)

is a given unimodal input signal which is identified by BLS-Prony as

$m = 2, z_{1, 2} = 0.9945 \pm j 0.0313$

${\overset{⌢}{A}}_{1} = 1.5, {\overset{⌢}{σ}}_{1} = - 0.5, {\overset{⌢}{f}}_{1} = 0.5, {\overset{⌢}{ϕ}}_{1} = 60.0000, ζ_{1} = 0.1572$

The BLS-Prony reconstruction is proved to be able to be recognized without error based on the fitting of the two curves in figure. In addition, the fact that the pole is located near the unit circle proves that the set system is in a critically damped state, and the recognized damping ratio also proves this feature.

Figure 6 shows the Prony model order identification process for three modes. The sampled data in the figure to construct the Hankel matrix, and the SVD decomposition of the criterion $σ_{7}$ converge to zero, proving that the order of the AR model is equal to 6. Figure 7 shows the comparison between the three modal Prony reconstructions and the input signal with the characteristic pole distribution. The setup conditions of the algorithm are as follows:

$y (t) = y_{1} + y_{2} + y_{3} \Rightarrow \{\begin{cases} y_{1} (t) = 1.5 \exp (- 0.5 t) \cdot \cos (2 π \times 0.5 t + π / 3) \\ y_{2} (t) = 1.5 \exp (- 0.2 t) \cdot \cos (2 π \times 2 t + π / 6) \\ y_{3} (t) = 1.5 \exp (- 0.5 t) \cdot \cos (2 π \times 3 t + 0) \end{cases}$

y (t)

is a given trimodal input signal, which is identified by BLS-Prony as

$m = 2, z_{1, 2} = 0.9774 \pm j 0.1864 \Rightarrow {\overset{⌢}{A}}_{1} = 1.5, {\overset{⌢}{σ}}_{1} = - 0.5, {\overset{⌢}{f}}_{1} = 0.5, {\overset{⌢}{ϕ}}_{1} = 60.0089, ζ_{1} = 0.1572$

$m = 4, z_{3, 4} = 0.9901 \pm j 0.1251 \Rightarrow {\overset{⌢}{A}}_{2} = 1.5, {\overset{⌢}{σ}}_{2} = - 0.2, {\overset{⌢}{f}}_{2} = 2, {\overset{⌢}{ϕ}}_{2} = 29.9920, ζ_{2} = 0.0159$

$m = 6, z_{5, 6} = 0.9945 \pm j 0.0313 \Rightarrow {\overset{⌢}{A}}_{3} = 1.5, {\overset{⌢}{σ}}_{3} = - 0.5, {\overset{⌢}{f}}_{3} = 3.0, {\overset{⌢}{ϕ}}_{3} = 0.0015, ζ_{3} = 0.0265$

Compared to the unimodal identification results, there are errors in the multimodal identification, and the absolute and normalized errors for each modality are as follows, respectively.

$Δ ε_{i} = | ϕ_{i} - {\overset{⌢}{ϕ}}_{i} | \Rightarrow Δ ε_{1} = 0.0089, Δ ε_{2} = 0.0080, Δ ε_{3} = 0.0015$

$ε_{i} % = Δ ε_{i} / 360^{°} \Rightarrow ε_{1} % = 0.0025 %, ε_{2} % = 0.0022 %, ε_{3} % = 0.0004 %$

In summary, the identification error proves that the observations of the BLS-Prony model are reliable, and combined with Figure 3 to Figure 6, it can be seen that both the normalized singular value and the cumulative energy function are able to accurately determine the oscillatory modes in the input signal under ideal conditions.

4.2. Dynamic Response Analysis of RLS-Prony

Figure 8 shows the oscillatory mode identification based on the RLS-Prony method, in which the input signal contains three modes, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8a show the original $y (t)$ and the reconstructed waveforms $y [t]$ of the RLS method, and Figure 8b shows the error waveforms $y (t) - y [t]$ , in which the error reaches the order of $10^{- 4}$ . Table 1 shows the solution process and iteration process for Figure 8. The initial length of the samples in the table is the modal order, and each time the iterative process does not converge, the length of the samples is increased in intervals of two. When the sample length is equal to 12, the equation converges. The converged roots and modal identification results are as follows.

$\begin{array}{l} z_{1} ({\bar{z}}_{1}) = 0.9774 \pm j 0.1865 \Rightarrow {\overset{⌢}{f}}_{1} = 0.5043 Hz, {\overset{⌢}{σ}}_{1} = - 0.5137, {\overset{⌢}{A}}_{1} = 1.5231 \\ z_{2} ({\bar{z}}_{2}) = 0.9901 \pm j 0.1251 \Rightarrow {\overset{⌢}{f}}_{2} = 2.0005 Hz, {\overset{⌢}{σ}}_{2} = - 0.2002, {\overset{⌢}{A}}_{2} = 1.5012 \\ z_{3} ({\bar{z}}_{3}) = 0.9944 \pm j 0.0315 \Rightarrow {\overset{⌢}{f}}_{3} = 3.0004 Hz, {\overset{⌢}{σ}}_{3} = - 0.5003, {\overset{⌢}{A}}_{3} = 1.4998 \end{array}$

Figure 9 shows the identification of oscillatory modes at different noise levels based on the RLS method, in which the input signals of Figure 9a,b contain two modes with an injected noise equal to 100 dB; and the input signals of Figure 9c,d contain three modes with an injected noise equal to 180 dB. The two SNR parameters above have been searched for through Monte Carlo experiments to be able to stabilize the extraction and resolution of modal information multiple times. As a result, although the SNR of the input signal is highly demanding, its error is an order of magnitude less when compared to Figure 8. There are two main reasons for the above problems. First, the rigidity of the Prony model makes it numerically sensitive, and the equations are not solved after the injection of noise, which leads to the dispersion of the amount of error in the iterative process. Second, the sample length is too low, which makes the rigidity problem more prominent. If the length is made equal to 100, it is possible to later reduce the required SNB to 120 dB.

In summary, regarding the ideal Prony model, although both BLS-Prony and RLS-Prony are able to fully recognize the modes in the signal, the recognition performance of both types of algorithms degrades when noise is injected into the signal. In comparison, BLS-Prony is able to achieve discrimination within the sensor accuracy, while RLS needs to increase the sample length to compensate for the rigidity of the equations.

4.3. Dynamic Response Analysis of Bayes–Prony

Figure 10 verifies the correctness of the maximum a posteriori probability distribution estimation and numerical solution procedure. Figure 10a shows the negative logarithmic a posteriori probability waveforms during the iterative process, with multiple peaks during the parameter optimization process which are mainly due to overshooting and tuning due to the large learning rate during the optimization process. Figure 10b–d show the learning paths for the frequency, phase, and attenuation coefficient, respectively. All parameters in the figures eventually converge to their true values. Finally, the sampled data is reconstructed according to the parameters optimized by MAP in Figure 10f, and the noise injected by the sampled data is equal to 2 dB. Note that in a 2 dB noise environment, the main features of the signal are largely submerged. MAP optimization is able to improve the discrimination accuracy by 120 dB compared to the rigid Prony model.

Figure 11 verifies the modal reconstruction performance of BLS-Prony, RLS-Prony, and Bayes–Prony under different signal-to-noise ratio conditions. In Figure 11a, the SNR is set to 200 dB, which also means an ideal input signal. The result in the figure shows that BLS-Prony and RLS-Prony can converge to the real signal perfectly. It is also the non-relaxation property of the Prony model itself that makes its AR model unique. In Figure 11b, SNR is set to 50 dB, and it can be seen that RLS-Prony firstly shows convergence error. This is due to the fact that the RLS iteration has less data and its solution cannot converge. Continuing to reduce the SNR until 20 dB, the phenomenon in Figure 11c is similar to that in Figure 11b. The SNR in Figure 11d is set to 2 dB, implying that the signal basically contains a lot of noise. At this point, only Bayes–Prony is able to converge to the true modal. In summary, BLS-Prony has the highest stability in the absence of extreme distortion of the signal. However, Bayes–Prony is able to maximize modal identification with increasing SNR.

In summary, the MAP optimization can extend the rigid point estimation of Prony to uncertain probability estimation, which greatly reduces the signal-to-noise ratio of the signal and makes it more adaptable to the real ship power grid environment. Of course, in the numerical solution process, in order to improve the optimization path of the parameters, we need adaptive design of the learning rate to eliminate the overshoot and regulate the identification path and design of gradient clipping to prevent the gradient explosion caused by data anomalies; concerning path optimization in the process of MAP solution, this paper will not provide discussion. At this point, the modal extraction of the sampled signal is completed.

5. Conclusions

A Bayes–Prony oscillation pattern recognition and hierarchical control strategy for LFO of an SM is proposed in this paper. The former replaces the point estimation of the traditional Prony algorithm with the probabilistic estimation idea, which improves the anti-noise performance of the traditional Prony algorithm. The latter serves as the decision control of whether the system is undergoing low-frequency oscillations, which determines whether the grid is put into damping control and the prerequisite for its own stabilization. Finally, the correctness and effectiveness of the proposed algorithms are verified by analyzing the time-domain characteristics of the proposed algorithms.

This paper describes the Prony/AR model applicable to the analysis of LFO. Since the Prony model itself matches the autoregressive model, the parameters in the Prony model are solved by least squares. Although both the BLS-Prony and RLS-Prony algorithms achieve the correct identification of ideal signals, this rigid solution cannot be applied to real microgrids for real signals containing noise. Therefore, the Bayes–Prony algorithm is proposed, which improves the anti-noise performance of the traditional algorithms. In addition, a hierarchical decision control strategy is designed based on the three algorithms, which provides observations to higher-authority controllers of the microgrid as a determination of the current state of the microgrid.

Although the proposed strategy can determine the current state of the microgrid with high confidence, the stability of Bayes–Prony, as the last layer of hierarchical control, is a key part of the successful realization of the whole strategy. Therefore, improving the identification performance of Bayes–Prony using, for example, weakly informative a priori or a priori sensitivity analysis will become a key research direction for engineering the algorithm in the future.

Author Contributions

Conceptualization, methodology, investigation, formal analysis, validation, visualization, writing—original draft, data curation, Y.D.; supervision, project administration, funding acquisition, resources, writing—review and editing, validation, K.Z.; methodology, software, investigation, J.D.; formal analysis, resources, L.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Dataset available on request from the authors. The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1 Block diagram of the ship microgrid studied in this paper.

Figure 2 Modal identification and hierarchical control framework.

Figure 3 Experimental test platform.

Figure 4 Results based on the BLS discrimination response of a single oscillatory mode. (a) All singular values after normalization. (b) First 6 orders of singularity. (c) The cumulative energy function of a matrix. (d) $ρ \geq 95 %$ .

Figure 5 BLS-Prony reconstructed waveforms in single oscillatory modes. (a) Reconstructed waveform. (b) Distribution of characteristic poles.

Figure 6 BLS-Prony recognition response based on three oscillatory modes. (a) All singular values after normalization. (b) First 12 orders of singularity. (c) The cumulative energy function of a matrix. (d) $ρ \geq 95 %$ .

Figure 7 BLS-Prony reconstructed waveforms in single oscillatory modes. (a) Reconstructed waveform. (b) Distribution of characteristic poles.

Figure 8 Waveforms based on the RLS-Prony response of three oscillatory modes. (a) Original and reconstructed waveforms. (b) Error waveform.

Figure 9 Modal identification of RLS for multimodal signals with different noise levels. (a) Bimodal reconstructed waveform. (b) Bimodal error waveform. (c) Trimodal reconstructed waveform. (d) Trimodal error waveform.

Figure 10 Flow of the maximum a posteriori probability distribution estimation and numerical solution procedure. (a) Posterior probability distribution. (b) Amplitude optimization. (c) Attenuation coefficient optimization. (d) Frequency optimization. (e) Gradient variation. (f) Time-domain dynamic response.

Figure 11 Modal reconstruction performance for different SNR conditions. (a) SNR = 200 dB. (b) SNR = 50 dB. (c) SNR = 20 dB. (d) SNR = 2 dB.

Table 1

The solution process and iteration process of Figure 8.

N	Number of Iterations	$a = [\begin{matrix} a_{1} & a_{2} & a_{3} & a_{4} & a_{5} & a_{6} \end{matrix}]$						$e [n]$
N	Number of Iterations	$a_{1}$	$a_{2}$	$a_{3}$	$a_{4}$	$a_{5}$	$a_{6}$	$e [n]$
8	1	−2.52	1.17	1.19	−0.27	−1.11	0.54	3.61
	2	−3.50	3.92	−0.44	−2.08	1.35	−0.23	1.40
	3	−4.10	5.74	−1.77	−3.03	2.96	−0.79	0.52
	4	−5.24	11.35	−12.86	7.99	−2.54	0.31	2.07
	5	−0.14	0.72	−1.42	1.41	−0.70	0.14	NaN
10	1	−2.61	1.29	1.27	−0.33	−1.26	0.64	0.027
	2	−3.50	3.91	−0.44	−2.08	1.34	−0.23	1.31
	3	−4.10	5.74	−1.77	−3.03	2.96	−0.79	0.52
	4	−5.24	11.35	−12.86	7.99	−2.54	0.31	2.07
	5	−0.14	0.72	−1.42	1.41	−0.70	0.14	NaN
12	1	−5.92	14.67	19.45	14.55	−5.82	0.97	0.001

References

1. Li, W.; Zhao, H.; Zhu, J.; Yang, T. A Novel Reactive Power Sharing Control Strategy for Shipboard Microgrids Based on Deep Reinforcement Learning. J. Mar. Sci. Eng.; 2025; 13, 718. [DOI: https://dx.doi.org/10.3390/jmse13040718]

2. Ma, Y.; Wang, Z.; Liu, H.; Tang, H.; Ji, Y.; Han, F. Efficient and sustainable power propulsion for all-electric ships: An integrated methanol-fueled SOFC-sCO₂ system. Renew. Energy; 2024; 230, 120822. [DOI: https://dx.doi.org/10.1016/j.renene.2024.120822]

3. International Maritime Organization (IMO). Cutting GHG Emissions. Available online: https://www.imo.org/zh/mediace-ntre/hottopics/pages/cutting-ghg-emissions.aspx (accessed on 11 November 2025).

4. Gong, Q.; Xu, J.; Ye, J.; Feng, H.; Shen, A. Nonlinear Model Predictive Control for Premixed Turbocharged Natural Gas Engine. IEEE/ASME Trans. Mechatron.; 2022; 27, pp. 3694-3705. [DOI: https://dx.doi.org/10.1109/TMECH.2021.3130910]

5. Wen, J.; Lu, J.; Zhang, S.; Liu, R.; Spataru, C.; Weng, Y.; Lv, X. Intelligent control for rapidity and security of all-electric ships gas turbine under complex mutation load using optimized neural network. Appl. Therm. Eng.; 2024; 248, 123120. [DOI: https://dx.doi.org/10.1016/j.applthermaleng.2024.123120]

6. Liu, S.; Wang, Y.; Liu, Q.; Panchal, S.; Zhao, J.; Fowler, M.; Fraser, R.; Yuan, J. Thermal equalization design for the battery energy storage system (BESS) of a fully electric ship. Energy; 2024; 312, 133611. [DOI: https://dx.doi.org/10.1016/j.energy.2024.133611]

7. Yu, Y.; Guan, Y.; Kang, W.; Vasquez, J.C.; Guerrero, J.M. A Generic Inertial-Response Preserved Active-Damping Control for Grid-Forming Inverters Emulating Synchronous Machines. IEEE Trans. Ind. Electron.; 2025; 72, pp. 5897-5905. [DOI: https://dx.doi.org/10.1109/TIE.2024.3481992]

8. Wei, C.-G.; Shangguan, X.-C.; Yang, Y.-H.; He, Y.; Zhang, C.-K. Delay-Dependent H_∞Robust Frequency Control in Microgrids: Coordination of Secondary Frequency Control and Virtual Inertia Control. IEEE Trans. Ind. Inform.; 2025; 21, pp. 6455-6465. [DOI: https://dx.doi.org/10.1109/TII.2025.3567267]

9. Liu, R.; Zhang, H.; Xue, C.; Li, Y.R. Distributed Low-Frequency Oscillation Damping in Low-Voltage Islanded Multi-Bus Microgrids With Virtual Synchronous Generators. IEEE Trans. Smart Grid; 2024; 15, pp. 5331-5345. [DOI: https://dx.doi.org/10.1109/TSG.2024.3417827]

10. Maroufi, S.M.; Karrari, S.; Rajashekaraiah, K.; De Carne, G. Power Management of Hybrid Flywheel-Battery Energy Storage Systems Considering the State of Charge and Power Ramp Rate. IEEE Trans. Power Electron.; 2025; 40, pp. 9944-9956. [DOI: https://dx.doi.org/10.1109/TPEL.2025.3546013]

11. Sun, C.; Ali, S.Q.; Joos, G.; Bouffard, F. Design of Hybrid-Storage-Based Virtual Synchronous Machine With Energy Recovery Control Considering Energy Consumed in Inertial and Damping Support. IEEE Trans. Power Electron.; 2022; 37, pp. 2648-2666. [DOI: https://dx.doi.org/10.1109/TPEL.2021.3111482]

12. Vykhodtsev, A.V.; Jang, D.; Wang, Q.; Rosehart, W.; Zareipour, H. Physics-Aware Degradation Model of Lithium-ion Battery Energy Storage for Techno-Economic Studies in Power Systems. IEEE Trans. Sustain. Energy; 2024; 15, pp. 316-327. [DOI: https://dx.doi.org/10.1109/TSTE.2023.3285453]

13. Ma, Q.; Chen, L.; Li, L.; Min, Y.; Gong, Y.; Liang, K. Effect of Grid-Following VSC on Terminal Frequency. IEEE Trans. Power Syst.; 2023; 38, pp. 1775-1778. [DOI: https://dx.doi.org/10.1109/TPWRS.2022.3233761]

14. Chen, S.; Sun, Y.; Han, H.; Fu, S.; Luo, S.; Shi, G. A Modified VSG Control Scheme with Virtual Resistance to Enhance Both Small-Signal Stability and Transient Synchronization Stability. IEEE Trans. Power Electron.; 2023; 38, pp. 6005-6014. [DOI: https://dx.doi.org/10.1109/TPEL.2023.3243025]

15. Hasan, K.N.; Preece, R. Influence of Stochastic Dependence on Small-Disturbance Stability and Ranking Uncertainties. IEEE Trans. Power Syst.; 2018; 33, pp. 3227-3235. [DOI: https://dx.doi.org/10.1109/TPWRS.2017.2779887]

16. Du, W.; Wang, Y.; Wang, H.F.; Ren, B.; Xiao, X. Small-Disturbance Stability Limit of a Grid-Connected Wind Farm With PMSGs in the Timescale of DC Voltage Dynamics. IEEE Trans. Power Syst.; 2021; 36, pp. 2366-2379. [DOI: https://dx.doi.org/10.1109/TPWRS.2020.3031351]

17. Sanchez-Gasca, J.; Chow, J. Computation of Power System Low-Order Models from Time Domain Simulation Using HANKEL Matrix. IEEE Trans. Power Syst.; 1997; 12, pp. 1461-1467. [DOI: https://dx.doi.org/10.1109/59.627842]

18. Lara, J.D.; Henriquez-Auba, R.; Ramasubramanian, D.; Dhople, S.; Callaway, D.S.; Sanders, S. Revisiting Power Systems Time-Domain Simulation Methods and Models. IEEE Trans. Power Syst.; 2014; 39, pp. 2421-2437. [DOI: https://dx.doi.org/10.1109/TPWRS.2023.3303291]

19. Yutthagowith, P.; Pattanadech, N. Improved Least-Square Prony Analysis Technique for Parameter Evaluation of Lightning ImpulseVoltage and Current. IEEE Trans. Power Deliv.; 2016; 31, pp. 271-277. [DOI: https://dx.doi.org/10.1109/TPWRD.2015.2448640]

20. Zhao, J.; Zhang, G. A Robust Prony Method Against Synchrophasor Measurement Noise and Outliers. IEEE Trans. Power Syst.; 2017; 32, pp. 2484-2486. [DOI: https://dx.doi.org/10.1109/TPWRS.2016.2612887]

21. Wang, P.; Li, X. Hybrid Sparse Expansion and Separable Hybrid Prony Method. IEEE Trans. Signal Process.; 2021; 69, pp. 3033-3043. [DOI: https://dx.doi.org/10.1109/TSP.2021.3078108]

22. Khodaparast, J.; Fosso, O.B.; Molinas, M. Recursive Multi-Channel Prony for PMU. IEEE Trans. Power Deliv.; 2024; 39, pp. 693-705. [DOI: https://dx.doi.org/10.1109/TPWRD.2023.3335999]

23. Huang, K.; Li, W.; Cao, S.; Gao, F.; Li, R.; Xu, W.; Lin, B. Recent advances in fault diagnosis of ship integrated power systems: A review. Ocean Eng.; 2026; 343, 123141. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.123141]

24. Mariscotti, A. Power Quality Phenomena, Standards, and Proposed Metrics for DC Grids. Energies; 2021; 14, 6453. [DOI: https://dx.doi.org/10.3390/en14206453]

25. Du, Z.; Yin, H.; Zhang, X.; Hu, H.; Liu, T.; Hou, M.; Giannelos, S.; Strbac, G. Decarbonisation of Data Centre Networks through Computing Power Migration. Proceedings of the 2025 IEEE 5th International Conference on Computer Communication and Artificial Intelligence (CCAI); Haikou, China, 23–25 May 2025.

26. Kaloev, M.; Krastev, G. Tailored Learning Rates for Reinforcement Learning: A Visual Exploration and Guideline Formulation. Proceedings of the 2023 7th International Symposium on Innovative Approaches in Smart Technologies (ISAS); Istanbul, Turkiye, 23–25 November 2023.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Bayesian Prony Modal Identification and Hierarchical Control Strategy for Low-Frequency Oscillation of Ship Microgrid

Content area

Abstract

Full text

1. Introduction

2. Description of the Modal and Identification Algorithms

2.1. Model of the Oscillatory Modal and Autoregressive Prediction

2.2. BLS-Prony Method Extraction of Oscillatory Modal

2.3. RLS-Prony Method Extraction of Oscillatory Modal

3. Bayes–Prony Method for Probabilistic Identification

3.1. Bayes Model and Posterior Probability Distribution Design

3.2. Maximum Posteriori Probability Estimation and Optimization for Solution Design

3.3. Hierarchical Control Strategy Design

4. Result and Discussion

4.1. Dynamic Response Analysis of BLS-Prony

4.2. Dynamic Response Analysis of RLS-Prony

4.3. Dynamic Response Analysis of Bayes–Prony

5. Conclusions