Data-Enabled Predictive Control for Grid-Connected Power Converters

TL;DR

This paper introduces a data-driven control approach for grid-connected power converters, combining DeePC and ARMA-based MPC to improve stability and performance without relying on explicit system models.

Contribution

It presents a novel combination of DeePC and ARMA-based MPC methods, addressing scalability issues and demonstrating effectiveness in stabilizing power converters.

Findings

DeePC eliminates undesired oscillations in power converters.

ARMA-based MPC improves scalability for high-order systems.

The combined approach stabilizes a voltage source converter in a power system.

Abstract

We apply a novel data-enabled predictive control (DeePC) algorithm in grid-connected power converters to perform safe and optimal control. Rather than a model, the DeePC algorithm solely needs input/output data measured from the unknown system to predict future trajectories. We show that the DeePC can eliminate undesired oscillations in a grid-connected power converter and stabilize an unstable system. However, the DeePC algorithm may suffer from poor scalability when applied in high-order systems. To this end, we present a finite-horizon output-based model predictive control (MPC) for grid-connected power converters, which uses an N-step auto-regressive-moving-average (ARMA) model for system representation. The ARMA model is identified via an N-step prediction error method (PEM) in a recursive way. We investigate the connection between the DeePC and the concatenated PEM-MPC method, and then analytically and numerically compare their closed-loop performance. Moreover, the PEM-MPC is applied in a voltage source converter based HVDC station which is connected to a two-area power system so as to eliminate low-frequency oscillations. All of our results are illustrated with high-fidelity, nonlinear, and noisy simulations.

Figures

Tables

Table 1. TABLE I: Parameters of the Power Converter System

Base Values for Per-unit Calculation
Voltage base value: $U_b=380V$	Power base value: $S_b=50kVA$
Frequency base value: $f_b=50Hz$
Parameters of the Power Part (per-unit values)
Converter-side inductor: $L_F=0.05$	LCL capacitor: $C_F=0.05$
Grid-side inductor: $L_g=0.35$	Grid-side resistor: $R_g=0.02$
Parameters of the Control Part
PI gains of the PLL: $103.1\left(\mathrm{rad}/s\right),5311.5\left(\mathrm{rad}/s\right)$
PI gains of the current control loop: $0.3(p.u.),10(p.u.)$
Control frequency for the digital signal processor: $10kHz$
Parameters of the DeePC
Length of the initial trajectories: $T_{\mathrm{ini}}=40$
Length of the prediction horizon: $N=30$
Length of the data to construct the Hankel matrix: $T=500$
Sampling frequency for the DeePC: $1kHz$
Control cost: $R=I$	Output cost: $Q=400\times I$
Reference vector : $r={\mathrm{𝟏}}_N\otimes \mathrm{col}(1,0,1)$
Constraint set : $𝒰=\{u:-2\times {\mathrm{𝟏}}_{mN}\le u\le 2\times {\mathrm{𝟏}}_{mN}\}$
Constraint set : $𝒴=\{u:r-2\times {\mathrm{𝟏}}_{pN}\le y\le r+2\times {\mathrm{𝟏}}_{pN}\}$
Regularization parameter: ${\lambda }_g=10$

Equations

$$\mathscr{H}_{L}(u):=\left[{\begin{array}[]{*{20}{c}}{{u_{1}}}&{{u_{2}}}&\cdots&{{u_{T-L+1}}}\\ {{u_{2}}}&{{u_{3}}}&\cdots&{{u_{T-L+2}}}\\ \vdots&\vdots&\ddots&\vdots\\ {{u_{L}}}&{{u_{L+1}}}&\cdots&{{u_{T}}}\end{array}}\right]$$

$$\left\{\begin{array}[]{l}{x_{t+1}}=A{x_{t}}+B{u_{t}}\\ {y_{t}}=C{x_{t}}+D{u_{t}}\end{array}\right.\,,$$

$${\mathcal{O}}_{\ell }(A,C):=\mathrm{col}(C,CA,...,C{A}^{\ell -1})$$

$$\left[{\begin{array}[]{*{20}{c}}{{U_{P}}}\\ {{U_{f}}}\end{array}}\right]:=\mathscr{H}_{T_{\rm ini}+N}(u^{\rm{d}})\,,\;\;\left[{\begin{array}[]{*{20}{c}}{{Y_{P}}}\\ {{Y_{f}}}\end{array}}\right]:=\mathscr{H}_{T_{\rm ini}+N}(y^{\rm{d}})\,,$$

$$\left[{\begin{array}[]{*{20}{c}}{{U_{P}}}\\ {{Y_{P}}}\\ {{U_{f}}}\\ {{Y_{f}}}\end{array}}\right]g=\left[{\begin{array}[]{*{20}{c}}{{u_{\rm ini}}}\\ {{y_{\rm ini}}}\\ u\\ y\end{array}}\right]\,.$$

$$\begin{array}[]{l}\mathop{{\rm{min}}}\limits_{g,u\in\mathcal{U},y\in\mathcal{Y}}\;\;{\left\|u\right\|_{R}^{2}}+{\left\|{y-r}\right\|_{Q}^{2}}+{\lambda_{g}}{\left\|g\right\|_{2}^{2}}\\ s.t.\;\;\left[{\begin{array}[]{*{20}{c}}{{U_{P}}}\\ {{Y_{P}}}\\ {{U_{f}}}\\ {{Y_{f}}}\end{array}}\right]g=\left[{\begin{array}[]{*{20}{c}}{{u_{\rm ini}}}\\ {{y_{\rm ini}}}\\ u\\ y\end{array}}\right]\,,\end{array}$$

$$\begin{array}[]{l}\mathop{\min}\limits_{u\in\mathcal{U},y\in\mathcal{Y}}\;\;{\left\|u\right\|_{R}^{2}}+{\left\|{y-r}\right\|_{Q}^{2}}\\ s.t.\;\;y=K\left[{\begin{array}[]{*{20}{c}}{{u_{\rm ini}}}\\ {{y_{\rm ini}}}\\ u\end{array}}\right]\,,\end{array}$$

$$y_i=K_i\phi ,$$

$$\underset{K_i}{\min }\sum _{j=1}^{N_{\mathrm{trj}}}{(y_{i(j)}-K_i{\phi }_{(j)})}^2,$$

$$\left\{{\begin{array}[]{*{20}{l}}{{K_{i(j)}}={K_{i(j-1)}}+L_{(j)}^{T}\left({{y_{i(j)}}-K_{i(j-1)}{\varphi_{(j)}}}\right)}\\ {{L_{\left(j\right)}}=\frac{{{P_{\left({j-1}\right)}}{\varphi_{(j)}}}}{{1+\varphi_{(j)}^{T}{P_{\left({j-1}\right)}}{\varphi_{(j)}}}}}\\ {{P_{\left(j\right)}}={P_{\left({j-1}\right)}}-\frac{{{P_{\left({j-1}\right)}}{\varphi_{(j)}}\varphi_{(j)}^{T}{P_{\left({j-1}\right)}}}}{{1+\varphi_{(j)}^{T}{P_{\left({j-1}\right)}}{\varphi_{(j)}}}}}\end{array}}\right.\,,$$

$$K_i=[y_{i(1)}\cdots y_{i(N_{\mathrm{trj}})}][{\phi }_{(1)}\cdots {\phi }_{(N_{\mathrm{trj}})}{]}^{+},$$

$$K=[y_{(1)}\cdots y_{(N_{\mathrm{trj}})}][{\phi }_{(1)}\cdots {\phi }_{(N_{\mathrm{trj}})}{]}^{+}.$$

$$K={Y_{f}}{\left[{\begin{array}[]{*{20}{c}}{{U_{P}}}\\ {{Y_{P}}}\\ {{U_{f}}}\end{array}}\right]^{+}}\,.$$

$$g={\left[{\begin{array}[]{*{20}{c}}{{U_{P}}}\\ {{Y_{P}}}\\ {{U_{f}}}\end{array}}\right]^{+}}\left[{\begin{array}[]{*{20}{c}}{{u_{\rm ini}}}\\ {{y_{\rm ini}}}\\ u\end{array}}\right]\,,$$

$$\begin{array}[]{l}\mathop{\min}\limits_{g}\;\;{\left\|g\right\|_{2}^{2}}\\ s.t.\;\;\left[{\begin{array}[]{*{20}{c}}{{U_{P}}}\\ {{Y_{P}}}\\ {{U_{f}}}\end{array}}\right]g=\left[{\begin{array}[]{*{20}{c}}{{u_{\rm ini}}}\\ {{y_{\rm ini}}}\\ u\end{array}}\right]\,.\end{array}$$

$$\begin{array}[]{l}\mathop{\min}\limits_{u^{\prime}\in\mathcal{U},y^{\prime}\in\mathcal{Y}}\;\;{\left\|u^{\prime}\right\|_{R}^{2}}+{\left\|{y^{\prime}-r}\right\|_{Q}^{2}}\\ s.t.\;\;u^{\prime}={U_{f}}g^{\prime\star},\;y^{\prime}={Y_{f}}g^{\prime\star},\\ {\rm{where}}\;\;g^{\prime\star}=\mathop{\arg\min}\limits_{g^{\prime}}\;\;{\left\|g^{\prime}\right\|_{2}^{2}}\\ \;\;\;\;\;\;\;\;\;\;\;s.t.\;\;\left[{\begin{array}[]{*{20}{c}}{{U_{P}}}\\ {{Y_{P}}}\\ {{U_{f}}}\end{array}}\right]g^{\prime}=\left[{\begin{array}[]{*{20}{c}}{{u_{\rm ini}}}\\ {{y_{\rm ini}}}\\ u^{\prime}\end{array}}\right]\,.\end{array}$$

$$C(u,y,g)={\Vert u\Vert }_R^2+{\Vert y-r\Vert }_Q^2+{\lambda }_g{\Vert g\Vert }_2^2$$

$$C(u^{\star },y^{\star },g^{\star })\le C(u^{\prime \star },y^{\prime \star },g^{\prime \star }).$$

Full text

I Introduction

The penetration of power-electronic devices in modern power systems is ever-increasing due to the development of renewable energy, microgrids, high-voltage direct-current (HVDC) transmission systems, etc. [1, 2]. This tendency is posing great challenges to power system operations because the dynamics of power converters are substantially different from synchronous generators (SGs). For example, SGs have large rotational rotors which physically determine the output frequencies, while power converters consist of static semiconductor apparatus and have high controllability.

Conventionally, the control structure of power converters is designed according to engineering experience and the corresponding control gain tuning is based on iterative trial-and-error methods. Also, lots of effort has been put into the modeling of power converters, which provides insights into the system dynamics and criteria for control gain tuning [3, 4, 5]. However, these approaches heavily rely on rich engineering experience and lack systematicness. In addition, the control structure generally assumes a stiff power grid and may present poor robustness against variable grid conditions. For example, the most widely-used control structure, which consists of a phase-locked loop (PLL) and a current control loop, can become unstable when the power converter is connected to a weak grid with high grid impedance (or equivalently, low short-circuit ratio) [6, 7, 8].

Even though offline design and analysis (based on a nominal model) can be conducted to determine an optimal control parameter set, optimal performance can rarely be achieved during online operation because (i) the real parameters of the power converter (e.g., capacitance and inductances of the LCL filter) are hard to obtain due to different operation conditions and manufacturing inaccuracy; (ii) sometimes the underlying algorithms for the converter are designed by another manufacturer and are not obtainable, i.e., some part of the converter system is unknown; (iii) the power grid is generally an unknown system from power converter side which significantly affects the dynamic performance; and (iv) the offline design generally employs a constant power grid model (which in most cases is assumed to be an infinite bus) for the power converter, yet the real power grid is variable.

Normally, these problems are handled using robust or adaptive methods [9, 8]. However, these methods are still model-based, result in complex controllers, and suffer from scalability problems for large and uncertain (or even partially unknown) models – especially, in grid-connected applications. Inspired by recent advances in machine learning and artificial intelligence, recent control approaches entirely circumvent such model-based solutions in favor of data-driven approaches [10, 11, 12].

In this paper, we use a novel Data-enabled Predictive Control (DeePC) algorithm to compute optimal and safe control policies for grid-connected power converters, which uses real-time feedback to drive the unknown system along a desired (i.e., optimal and constrained) trajectory [13]. The DeePC algorithm presented in [13] relies on behavioural system approach [14, 15, 16, 17]. Instead of using a parametric model for system representation (e.g., state space matrices obtained from system identification), the approach in [14, 15, 16, 17] describes the input/output behaviour of the system through the subspace of the signal space in which trajectories live. This signal space of trajectories is spanned by the columns of a data Hankel matrix which results in a non-parametric and data-centric perspective on dynamical control systems.

The DeePC approach presented here relies on input/output data samples from the converter-internal signals and terminal signals towards the power grid (whose model is unknown from the perspective of the converter’s controller), and successfully eliminates undesired oscillations by solving the optimal regulation problem in a receding horizon manner – with the input/output constraints incorporated and in absence of any system model. However, when applied in large-scale systems, e.g., in the case of power transmision oscillation damping [18, 19], the optimal regulation problem in DeePC may suffer from poor scalability due to its high dimension.

To this end, we use a finite-horizon output-based model predictive control (MPC) for grid-connected power converters, wherein the unknown system is represented by an $N$-step auto-regressive-moving-average (ARMA) model and identified via least-square $N$-step prediction error method (PEM). The PEM can be solved in a recursive way which enables an iterative calculation and possible online implementation. We will show that this concatenated PEM-MPC method is scalable for large-scale unknown systems, and analytically discuss how it is related to DeePC. Namely, DeePC provably outperforms the PEM-MPC method in terms of the cost in the optimal regulation problem (see Lemma III.2) although this performance gap can be made smaller by appropriate regularizations of the optimal control and system identification problems. We propose to use the PEM-MPC in a grid-connected voltage source converter (VSC) based HVDC station to attenuate the power system oscillations which are caused by the interactions among multiple synchronous generators [20]. All of our results are illustrated with high-fidelity nonlinear simulations.

The remainder of this paper is organized as follows: in Section II we provide an overview for the DeePC approach and apply it to a grid-connected power converter. Section III presents the concatenated PEM-MPC and we discuss how it is related to DeePC. In Section IV we apply the PEM-MPC in a two-area power system which contains one VSC-HVDC station. We conclude the paper in Section V.

II Data-Enabled Predictive Control

II-A Preliminaries and Notation

For an unknown discrete-time LTI system that has $m$ inputs and $p$ outputs, we denote by $u_{i,t}\in\mathbb{R}$ the $i\rm{th}$ input of the system at time $t\in\mathbb{Z}_{\geq 0}$ and $y_{i,t}\in\mathbb{R}$ the $i\rm{th}$ output at time $t\in\mathbb{Z}_{\geq 0}$, where $\mathbb{Z}_{\geq 0}$ is the discrete-time axis. The input vector of the system is denoted by $u_{t}={\rm{col}}(u_{1,t},...,u_{m,t})\in\mathbb{R}^{m}$, and the output vector is denoted by $y_{t}={\rm{col}}(y_{1,t},...,y_{p,t})\in\mathbb{R}^{p}$, where ${\rm{col}}(a_{1},...,a_{i}):=[a_{1}^{\top}\;\cdots\;a_{i}^{\top}]^{\top}$. Let $u={\rm{col}}(u_{1},u_{2},...)$ and $y={\rm{col}}(y_{1},y_{2},...)$ be the input and output trajectories, respectively, whose dimensions can be inferred from the context.

Let $L,T\in\mathbb{Z}_{\geq 0}$ and $T\geq L$. The trajectory $u\in\mathbb{R}^{mT}$ is persistently exciting of order L if the Hankel matrix

[TABLE]

is of full row rank, i.e., the signal $u$ is sufficiently long and sufficiently rich.

Consider the following $n$-order discrete-time LTI system (minimal representation):

[TABLE]

where $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, $C\in\mathbb{R}^{p\times n}$, $D\in\mathbb{R}^{p\times m}$, and $x_{t}$ is the state of the system at time $t\in\mathbb{Z}_{\geq 0}$.

The lag of the system in (2) is defined by the smallest integer $\ell\in\mathbb{Z}_{\geq 0}$ so that the observability matrix

[TABLE]

has rank $n$.

Let $T_{\rm ini},N\in\mathbb{Z}_{\geq 0}$ such that $T\geq(m+1)(T_{\rm ini}+N+n)-1$. Consider an input trajectory $u^{\rm{d}}$ and an output trajectory $y^{\rm{d}}$ (both are of length $T\in\mathbb{Z}_{\geq 0}$, i.e., $u^{\rm{d}}\in\mathbb{R}^{mT}$ and $y^{\rm{d}}\in\mathbb{R}^{pT}$) measured from the $n$-order unknown system (2) such that $u^{\rm{d}}$ is persistently exciting of order $T_{\rm ini}+N+n$. Here we use the superscript d to indicate that these two trajectories are data sets measured from the unknown system. We use $u^{\rm{d}}$ and $y^{\rm{d}}$ to construct the Hankel matrices $\mathscr{H}_{T_{\rm ini}+N}(u^{\rm{d}})$ and $\mathscr{H}_{T_{\rm ini}+N}(y^{\rm{d}})$, which are further partitioned into two parts as

[TABLE]

where $U_{P}\in\mathbb{R}^{mT_{\rm ini}\times(T-T_{\rm ini}-N+1)}$, $U_{f}\in\mathbb{R}^{mN\times(T-T_{\rm ini}-N+1)}$, $Y_{P}\in\mathbb{R}^{pT_{\rm ini}\times(T-T_{\rm ini}-N+1)}$ and $Y_{f}\in\mathbb{R}^{pN\times(T-T_{\rm ini}-N+1)}$.

According to the behavioral system theory [15], ${\rm{col}}(u_{\rm ini},y_{\rm ini},u,y)$ is a trajectory of (2) if and only if there exist $g\in\mathbb{R}^{T-T_{\rm ini}-N+1}$ such that

[TABLE]

The trajectory ${\rm{col}}(u_{\rm ini},y_{\rm ini})$ can be thought of as an initial condition for the trajectory ${\rm{col}}(u_{\rm ini},y_{\rm ini},u,y)$ and ${\rm{col}}(u,y)$ as a future trajectory departing from this initial condition. If $T_{\rm ini}\geq\ell$, the future output trajectory $y$ is uniquely determined through (4) for every given input trajectory $u$.

II-B Review of DeePC

Instead of learning a parametric system representation through system identification, the DeePC attempts to learn the system’s behaviour and computes optimal control inputs using past data measured from the unknown system. Moreover, input/output constraints can be conveniently incorporated to ensure safety, described as follows.

After using the input/output trajectory ${\rm{col}}(u^{\rm{d}},y^{\rm{d}})$ ($u^{\rm{d}}\in\mathbb{R}^{mT}$ and $y^{\rm{d}}\in\mathbb{R}^{pT}$) to construct the Hankel matrices in (3), DeePC solves the following optimization problem at every sampling time to get the optimal future control inputs

[TABLE]

where $\mathcal{U}\subseteq\mathbb{R}^{mN}$ and $\mathcal{Y}\subseteq\mathbb{R}^{pN}$ are the input and output constraint sets, $R\in\mathbb{R}^{mN\times mN}$ is the control cost matrix (positive definite), $Q\in\mathbb{R}^{pN\times pN}$ is the output cost matrix (positive semidefinite), $\lambda_{g}\in\mathbb{R}_{\geq 0}$ is the regularization parameter, $r\in\mathbb{R}^{pN}$ is the reference vector for the output signals, $N$ is the prediction horizon, and ${\rm{col}}(u_{\rm ini},y_{\rm ini})$ consists of the most recent input/output trajectory of (2). The norm ${\left\|a\right\|_{X}^{2}}$ of the vector $a$ denotes the quadratic form $a^{T}Xa$, and the norm ${\left\|a\right\|_{2}}$ denotes $\sqrt{{a^{T}}a}$.

We note that a two-norm penalty on $g$ is included in the cost function as a regularization term to avoid overfitting. In fact, in the case when stochastic disturbances affect the output measurements, a two-norm regularization on $g$ coincides with two-norm robustness with respect to the noise affecting the output measurements [21]. The DeePC involves solving the optimization problem (DeePC) in a receding horizon manner [13], that is, after calculating the optimal control input sequence $u^{\star}$, we apply $(u_{t},...,u_{t+s})=(u_{0}^{\star},...,u_{s}^{\star})$ to the system for some $s\leq N-1$ time steps, update ${\rm{col}}(u_{\rm ini},y_{\rm ini})$ to the most recent input/output measurements and then set $t$ to $t+s+1$ for the DeePC algorithm.

II-C Application to a Grid-Connected Power Converter

Fig.1 shows the one-line diagram of a three-phase power converter which is connected to an ac power grid via an LCL. The system is nonlinear (the nonlinearity comes from the coordinate transformation), and the order of the system is $n=10$. Here, the ac power grid is modeled as an infinite bus with fixed voltage magnitude $1.0\rm{(p.u.)}$ and frequency $50\rm{Hz}$. The base values for per-unit calculation and the system parameters are given in Table I, wherein the symbol ${\bf{1}}_{n}\in\mathbb{R}^{n}$ denotes the column vector with all the entries being $1$, and ${\bf{0}}_{n}\in\mathbb{R}^{n}$ denotes the column vector with all the entries being [math], $I$ is the identity matrix whose dimension can be inferred from the text, and the Kronecker product of $A$ and $B$ is denoted by $A\otimes B$. The control part of the converter consists of a synchronous reference frame PLL, a current control loop and coordinate transformation blocks [4, 22].

The current control loop contains two proportional-integral (PI) regulators to make the converter-side d-axis and q-axis current components (i.e., $I_{d}$ and $I_{q}$ as shown in Fig.1) track their references $I_{d}^{ref}$ and $I_{q}^{ref}$, which enables fast current limiting under faults and harmonic suppression of the converter-side current.

The PLL is used for grid-synchronization. The input of the PLL is the q-axis voltage component (i.e., $V_{q}$) of the LCL’s capacitor, and the output is the angle reference $\theta$ for coordinate transformation. The $q$-axis voltage component $V_{q}$ is controlled to be zero in steady state via the PI regulator in the PLL, and therefore the voltage vector is aligned with the d-axis, that is, the d-axis voltage component (i.e., $V_{d}$) is the magnitude of the capacitor’s voltage vector in steady state.

When the converter is connected to a strong grid that features low grid-side impedance, it will present anticipated dynamic performance. However, one significant challenge for the converter’s operation is that generally the converter does not have any information about the properties of the ac power grid. If the converter is connected to a weak grid that has high grid-side impedance, e.g., $L_{g}=0.35\rm{p.u.}$ as used in this section, the PLL will have significant interaction with the current control loop as well as the grid impedance, which may result in instabilities of the system [23, 8]. This kind of small-signal instability features the oscillations of the PLL’s output and current/voltage signals, which should be eliminated to ensure the safe operation of the power grid.

In this section, we use the DeePC to perform optimal control for the converter that is connected to a weak ac grid, which can attenuate the oscillations in the converter and stabilize the system by penalizing the tracking errors in the cost function. The power converter together with the power grid is a black-box system from the view of the DeePC. The DeePC provides control inputs to this black-box system and measures its input/output trajectories with sampling frequency $1\rm{kHz}$ (i.e., sampling time $1\rm{ms}$).

We choose $I_{d}^{ref}$ and $I_{q}^{ref}$ to be the control inputs considering that the current control loop has high response speed and allows fast current limiting. The measured outputs from the black-box system are $V_{d}$, $V_{q}$ and $I_{d}$, which will be controlled to track their references via the DeePC. These output signals contain measurement noise (white noise with noise power: $5.0\times 10^{-6}$). We note that $I_{q}$, which corresponds to the reactive current, is not chosen as the measured output for the DeePC because $I_{q}$ is uniquely determined by the power flow constraint when $V_{d}$, $V_{q}$ and $I_{d}$ equal their reference values in steady state. Since the DeePC has no information about the black-box system, $T_{\rm ini}$ is chosen to be sufficiently large ($T_{\rm ini}=40$) to meet $T_{\rm ini}\geq\ell$.

Fig.2 plots the time-domain responses of the power converter. When the DeePC is activated at $t=1.0\rm{s}$, the voltage/current oscillations are effectively eliminated. By comparison, the system is unstable without the DeePC (where $I_{d}^{ref}$ and $I_{q}^{ref}$ are respectively set as $1.0{\rm{p.u.}}$ and [math]), and there are voltage/current oscillations whose magnitudes are amplified over time. These simulation results verify the effectiveness of using DeePC to perform optimal predictive control for a grid-connected power converter with unknown grid conditions (e.g., grid-side impedance and grid properties) and unknown converter model and inner controls.

III From DeePC to MPC

III-A Scalability of DeePC

The DeePC approach provides a safe and optimal solution to the regulation problem by solely using the measured data from the unknown system, and thus allows to stabilize an unknown and unstable system, as illustrated in the previous case study. However, when applied in high-order systems, the DeePC may suffer from poor scalability because the dimension of the decision variable $g$ (which is $T-T_{\rm ini}-N+1$) depends on the length of data $T$ to construct the Hankel matrix. In other words, the optimization problem in (DeePC) is of high dimension when choosing a long sequence of data to possibly eliminate the impacts of measurement noise.

To remove $g$ from the constraint and ensure the scalability of the optimization problem in (DeePC), we consider a subsequent system identification and model predictive control whose decision variables are $u$ and $y$.

III-B Model Predictive Control

In what follows, we consider the following finite-horizon output-based MPC problem

[TABLE]

where the decision variables $u$ and $y$ are the control inputs and measurement outputs over the prediction horizon, and $K\in{\mathbb{R}^{pN\times(m{T_{\rm ini}}+p{T_{\rm ini}}+mN)}}$ is the $N$-step transition matrix predicting how future outputs of the system are determined by the initial input/output data and the future inputs. Note that the optimization problem (MPC) is solved in a receding horizon manner resulting in an online feedback control.

Since the decision variable $g$ in (DeePC) does not appear in (MPC), solving (MPC) has much less computational burden than solving (DeePC), especially when the Hankel matrix has high column number (i.e., $g$ has high dimension). On the other hand, (MPC) depends on an explicit predictive model given by the transition matrix $K$ in the equality constraint.

III-C Prediction Error System Identification

Observe that the predictive model in (MPC) is an $N$-step ARMA model for the discrete-time LTI system mapping past inputs and outputs ${\rm col}(u_{\rm ini},y_{\rm ini})$ as well as future inputs $u$ to future outputs $y$. In particular, for the $i\rm{th}$ ($i\in\{1,...,pN\}$) element of $y$, we have

[TABLE]

where $K_{i}$ is the $i\rm{th}$ row of $K$ and $\varphi={\rm{col}}(u_{\rm ini},y_{\rm ini},u)$.

Given past measurements of $y_{i}$ and $\varphi$, the transition matrix $K_{i}$ can be computed offline through system identification. In the absence of measurement noise, $K_{i}$ can be computed exactly with $m{T_{\rm ini}}+p{T_{\rm ini}}+mN$ linearly independent measurements of $\varphi$ and associated $y_{i}$. However, measurement noise will significantly affect the accuracy of this approach. A standard solution to remedy this problem and to eliminate the effects of the noise is to use a larger data set and apply a least-square $N$-step PEM minimizing

[TABLE]

where $N_{\rm trj}>m{T_{\rm ini}}+p{T_{\rm ini}}+mN$ is the number of measured trajectories, and $y_{i(j)}$ and $\varphi_{(j)}$ belong to the $j\rm{th}$ trajectory. Indeed, the subsequent combination of PEM and MPC is a standard approach to model-based control that has proved itself in many applications throughout academia and industry [24, 25].

As an alternative to the batch optimization approach (PEM) combining all the measured trajectories to solve for $K$ in one step, $K$ can be obtained by adopting the recursive least-square algorithm [26]

[TABLE]

where $K_{i(j-1)}$ is the least-square solution for $K_{i}$ after combining the previous $j-1$ trajectories, and $K_{i(j-1)}$ is updated to $K_{i(j)}$ using the data (i.e., $y_{i(j)}$ and $\varphi_{(j)}$) in the $j\rm{th}$ trajectory. The matrices $L_{(j)}\in\mathbb{R}^{(m{T_{\rm ini}}+p{T_{\rm ini}}+mN)\times 1}$ and $P_{(j)}\in\mathbb{R}^{(m{T_{\rm ini}}+p{T_{\rm ini}}+mN)\times(m{T_{\rm ini}}+p{T_{\rm ini}}+mN)}$ are two intermediate variables updated through the $j\rm{th}$ trajectory in the recursive algorithm.

The recursive algorithm (6) is not only more scalable for large data sets, but it could possibly also be applied online to adapt the model $K$ used in (MPC) to a variable environment. Moreover, the size of $K$ is solely related to $T_{\rm ini}$ and $N$ and independent of the size of the data set, which leads to a smaller optimization problem size than (DeePC) when applied to high-order systems, considering that the dimension of the decision variable $g$ in (DeePC) depends on the length of the data $T$ to construct the Hankel matrix.

III-D Relation of DeePC, MPC, and PEM

In the following, we will relate the subsequent system identification through (PEM) and model-based control through (MPC) to the data-driven (DeePC) strategy. In a first step, observe that the least-square solution for (PEM) can be expressed in closed form as

[TABLE]

where the Moore-Penrose pseudoinverse of the matrix $X$ is denoted by $X^{+}$.

The transition matrix $K$ can then be obtained as

[TABLE]

Since every column in the constructed Hankel matrix (3) is one trajectory measured from the unknown system, $K$ can also be calculated by using the Hankel matrix (3) by

[TABLE]

By combining (9) with $y=K\varphi$ and $y=Y_{f}g$ one obtains

[TABLE]

which equals the solution of the optimization problem

[TABLE]

Observe that the optimization problem (LN) seeks the least-norm solution to the equality constraint. Thus, we refer to it as the least-norm (LN) problem. As derived above, for Hankel matrix data, the solution $g^{\prime\star}$ of (LN) is related to the solution of the least-squares (PEM) problem $K^{\star}$ by $y=Y_{f}g^{\prime\star}=K^{\star}\varphi$. It can also be derived using subspace identification methods [17]. We summarize this observation below.

Lemma III.1** (PEM and LN)**

Consider the least-norm problem (LN) with Hankel matrix ${\rm{col}}(U_{P},Y_{P},U_{f})$. Consider the least-square $N$-step prediction error method optimization problem (PEM), and assume that its data $y_{i(j)}$ and $\varphi_{(j)}$ is arranged into the Hankel matrix ${\rm{col}}(U_{P},Y_{P},U_{f},Y_{f})$. Then the solution $K^{\star}$ of (PEM) and the solution $g^{\prime\star}$ of (LN) are related by the equation $K^{\star}\varphi=Y_{f}g^{\prime\star}$.

Observe that if the (LN) solution $g^{\prime\star}$ is used as predictive model (equality constraint) of the (MPC) problem by setting $y=Y_{f}g^{\prime\star}$ (where $g^{\prime\star}$ is a function of $u$), then the subsequent concatenation of the (MPC) and the (LN) (or equivalently (PEM) by Lemma III.1) optimization problems reads as follows:

[TABLE]

Note that $g^{\prime\star}$ depends on the decision variable $u^{\prime}$ in (PEM-MPC). The inner problem of (PEM-MPC) is the system identification step through (LN) (or equivalently (PEM) by Lemma III.1). The outer problem on the other hand is identical to the (MPC) problem. Let

[TABLE]

be the combined cost taking into account both system performance (i.e., ${\left\|{{u}}\right\|_{R}^{2}}+{\left\|{{y}-r}\right\|_{Q}^{2}}$) and the complexity of $g$ (i.e., ${\lambda_{g}}{\left\|{{g}}\right\|_{2}^{2}}$). This is in fact the true cost we wish to minimize when noise is present in the system as large entries of $g$ result in overfitting of the noisy trajectories in the Hankel matrix in (4) [13]. Since this cost appears in (DeePC), we can compare directly the performance of (PEM-MPC) and (DeePC) with respect to this cost. We present the comparison below.

Lemma III.2** (PEM-MPC and DeePC)**

Consider the optimal solution $(u^{\prime\star},y^{\prime\star},g^{\prime\star})$ of the concatenated (PEM-MPC) problem and the optimal solution $(u^{\star},y^{\star},g^{\star})$ of the (DeePC) problem. It holds that

[TABLE]

That is, (DeePC) achieves a cost less or equal to the cost of the concatenated (PEM-MPC) with respect to (11).

Proof:

Observe that any feasible point $(u^{\prime},y^{\prime},g^{\prime\star})$ of (PEM-MPC) is also a feasible point of (DeePC). Since $g^{\prime\star}$ is a particular $g$ that satisfies the constraints of (DeePC) then the feasible set of (PEM-MPC) is a subset of (DeePC). Hence, it holds that $C(u^{\star},y^{\star},g^{\star})\leq C(u^{\prime\star},y^{\prime\star},g^{\prime\star})$. ∎

In other words, the DeePC presents better performance than the MPC formulated in (MPC) if $K$ is obtained by (PEM). On the other hand, the MPC problem in (MPC) has the advantage that it doesn’t contain the decision variable $g$ and therefore has lower computational burden, which make it scalable to high-order systems. Furthermore, a recursive algorithm can be used in PEM-MPC to obtain the $N$-step transition matrix $K$ online, which enables the implementation for digital signal processors. In the next section we numerically observe that the gap between (DeePC) and (PEM-MPC) presented in Lemma III.2 can actually be made smaller by an appropriate choice of regularization $\lambda_{g}$.

III-E Comparative Case Studies

Fig.3 displays the time-domain responses of the power converter when the PEM-MPC is used, with the other settings remaining the same as those in Fig.1. Before $t=0\rm{s}$, $I_{d}^{ref}$ and $I_{q}^{ref}$ were respectively set as $1.0{\rm{p.u.}}+\tau_{1}$ and $\tau_{2}$ for $30\rm{s}$, where $\tau_{1}$ and $\tau_{2}$ are two different white noise signals (noise power: $1.0\times 10^{-4}$). During this process, the power converter is critically stable with a grid-side inductor $L_{g}=0.34\rm{p.u.}$, and the other parameters are the same as those in Table I.

By using recursive least-square algorithm, the $N$-step transition matrix $K$ was calculated iteratively with measurements spaced by $20\rm{ms}$, that is, $K$ has been obtained during the converter’s operation by recursively solving the PEM problem in (PEM) with $1500$ trajectories before $t=0\rm{s}$.

At $t=0.2\rm{s}$, the PEM-MPC is activated, which effectively attenuates the voltage and current oscillations and stabilizes the system. Then, we test the robustness of the PEM-MPC by changing the grid-side inductance $L_{g}$ from $0.34\rm{p.u.}$ to $0.35\rm{p.u.}$ at $t=0.7\rm{s}$, and to $0.5\rm{p.u.}$ at $1.0\rm{s}$. It can be seen that the system is stable under these two disturbances, and the current and voltage signals track the references with fast dynamics, that is, the PEM-MPC shows robustness in terms of parameter changes in the black-box system. Here $K$ is not updated during the two disturbances to test how an inaccurate model in $K$ affects the performance of the PEM-MPC.

For comparison, Fig.3 also plots the time-domain responses of the power converter when the DeePC is applied. By choosing $T=500$, the DeePC also effectively eliminates the voltage and current oscillations when activated at $t=0.2\rm{s}$, and presents robustness when $L_{g}$ is changed from $0.34\rm{p.u.}$ to $0.35\rm{p.u.}$ at $t=0.7\rm{s}$, and to $0.5\rm{p.u.}$ at $1.0\rm{s}$. Note the Hankel matrix is obtained from the system with $L_{g}=0.34\rm{p.u.}$ and is not updated during the disturbances.

When choosing $T=330$ in the DeePC, the oscillations are obviously eliminated as well after the DeePC is activated. However, the voltage and current signals are oscillating after $L_{g}$ is changed from $0.35\rm{p.u.}$ to $0.5\rm{p.u.}$ at $1.0\rm{s}$. This is because the trajectory ${\rm{col}}(u^{\rm{d}},y^{\rm{d}})$ in this case contains less information of the system than that with $T=500$, and thus the prediction is more strongly affected by the measurement noise when solving the DeePC. Choosing a larger $T$ is a convenient way to resolve this problem, yet it results in a higher dimension of $g$ and thus higher computational burden in solving (DeePC), which may not be scalable to high-$\ell$ systems as discussed before.

Fig.4 plots the relationship between $T$ and the time-domain cost (from $t=0.2\rm{s}$ to $t=1.4\rm{s}$) to illustrate how different values of $T$ affect the system performance, where the time-domain cost is $\left\|{{u_{\rm time}}}\right\|_{R}^{2}+\left\|{{y_{\rm time}}-{r_{\rm time}}}\right\|_{Q}^{2}$ ($u_{\rm time}$ and $y_{\rm time}$ are the input and output trajectories measured from the system, and $r_{\rm time}$ is the reference vector for the outputs). It is shown that the time-domain cost sharply decreases with the increase of $T$ from $320$ to $400$, and remain nearly the same (or slightly increases) if further increasing $T$.

For comparison, the converter’s responses without the DeePC or the PEM-MPC are plotted by the green lines in Fig.3. In this case, the grid-side inductance $L_{g}$ is also changed from $0.34\rm{p.u.}$ to $0.35\rm{p.u.}$ at $t=0.7\rm{s}$, and to $0.5\rm{p.u.}$ at $1.0\rm{s}$. It can be seen that the system is unstable with the voltage and current signals keeping oscillating, which endangers the power system operation.

Fig.5 displays how $\lambda_{g}$ affects the optimization cost of the system (setting $T=500$). The optimization costs of DeePC and PEM-MPC are respectively obtained by solving (DeePC) and (PEM-MPC) at $t=0.2\rm{s}$, and Fig.5 shows that DeePC outperforms PEM-MPC under any $\lambda_{g}$, which is consistent with Lemma III.2. In addition, Fig.5 also plots how $\lambda_{g}$ affects the time-domain cost (from $t=0.2\rm{s}$ to $t=1.4\rm{s}$) of the DeePC, and shows that the time-domain cost is dramatically reduced when increasing $\lambda_{g}$ from $0.01$ to $20$, and it slightly increases if further increasing $\lambda_{g}$.

IV Application of PEM-MPC in Large-Scale Systems

To illustrate the effectiveness of the PEM-MPC in large-scale systems, we now provide a detailed simulation study based on a nonlinear model of the three-phase two-area test system (the order of the system is $n=90$) in Fig.6. This test system consists of four SGs and one VSC-HVDC station. The VSC-HVDC station is a large-capacity power converter with an LCL filter, which in our case, employs the control scheme given in Fig.1. The base values for this test system are: $f_{\rm{b}}=50\rm{Hz}$, $S_{\rm{b}}=350\rm{MVA}$ and $U_{\rm{b}}=220\rm{kV}$ (line to line). In the VSC-HVDC station, the grid-side inductor of the LCL filter is chosen as $L_{g}=0.05{\rm{p.u.}}$, and the other parameters are the same as those in Table I unless otherwise specified. The parameters of the SGs and the power grid are chosen as the same as those in [19].

This test system has one pair of weakly-damped inter-area modes and presents low-frequency oscillations, as caused by the fast exciters in the SGs as well as long transmission lines [20]. The PEM-MPC enables us to use the high controllability of the converter to attenuate the inter-area oscillations. We still choose $I_{d}^{ref}$ and $I_{q}^{ref}$ to be control inputs provided by the PEM-MPC, and $V_{d}$ and $V_{q}$ are two of the three measured outputs from the converter. Another output signal is the inter tie-line active power $P_{Tie-Line}$ transferred from Area 1 to Area 2, as shown in Fig.6. All the outputs contain measurement noise (white noise with noise power $1.0\times 5.0^{-6}$). We note that DeePC may not be suitable to be applied to such a high-order system because (DeePC) contains a high-dimension decision variable $g$ and thus cannot be solved in real time.

Fig.7 plots the time-domain responses of the inter tie-line active power when different $T_{\rm ini}$ and $N$ are adopted. The sampling time for PEM-MPC is $1\rm{ms}$. Note that before $t=0\rm{s}$, white noise signals (noise power: $1.0\times 10^{-4}$) were injected into the test system through $I_{d}^{ref}$ and $I_{q}^{ref}$, which lasted for 1 minute. During this process, the $N$-step transition matrix $K$ was calculated iteratively with measurements spaced by $20\rm{ms}$, that is, $K$ was obtained by solving the PEM problem with 3000 trajectories but in the following cases, $K$ is not adapted online when the PEM-MPC is performed.

Firstly, since the test system is of high order (with $90$ state variables) and is a black-box system to the PEM-MPC, we choose a sufficiently large $T_{\rm ini}$ ($T_{\rm ini}=200$) to meet $T_{\rm ini}\geq\ell$, and the prediction horizon is chosen as $N=80$. It can be seen that the inter-area oscillation is well eliminated after the PEM-MPC is activated at $t=10\rm{s}$. Then, considering that the low-frequency oscillations are caused by the interactions among the SGs’ rotors and thus can possibly be represented by an equivalent low-order system, $T_{\rm ini}$ and $N$ are chosen smaller to test the performance of the PEM-MPC. Fig.7 shows that when $(T_{\rm ini},N)$ are chosen to be $(10,10)$ and even $(5,10)$, the PEM-MPC can still effectively attenuate the inter tie-line power oscillations. Also observe that with the decrease of $T_{\rm ini}$ and $N$, some slight oscillations still exist after the PEM-MPC is activated, which is caused by the model mismatching, i.e., a smaller size of $K$ may not be able to accurately represent the full-order model.

To further illustrate how the choice of $T_{\rm ini}$ affects the overall performance, Fig.8 plots the variations of the time-domain cost with different values of $T_{\rm ini}$ (from $t=10\rm{s}$ to $t=14\rm{s}$). Fig.8 shows that the time-domain cost is significantly reduced with the increase of $T_{\rm ini}$ from $5$ to $50$, but remains nearly the same if further increasing $T_{\rm ini}$, which indicates the test system is well represented with $T_{\rm ini}=50$.

V Conclusion

We applied DeePC in grid-connected power converters to eliminate undesired oscillations caused by weak grid conditions. The DeePC has no information about the power converter or the power grid, and solely uses input/output data measured from the unknown system. We showed that the DeePC can stabilize an unstable system by performing optimal and constrained receding-horizon control. Moreover, the DeePC presents more robustness against changes in the unknown system with a longer input/output data sequence to construct the Hankel matrix, but meanwhile results in higher computational burden and poor scalability especially in high-order systems. For this reason, we presented a concatenated PEM-MPC method as an alternative model-based, optimal, and constrained receding-horizon control, wherein the PEM can be implemented in a recursive manner. We discussed the connection between DeePC and PEM-MPC, and formally showed that the DeePC outperforms PEM-MPC with regards to the cost in the optimization problem. We applied the PEM-MPC in a VSC-HVDC station which is connected to a high-order power grid that contains four SGs, and showed that the PEM-MPC effectively eliminates the low-frequency oscillations caused by the interactions among multiple SGs.