Deep Learning-Aided Model Predictive Control of Wind Farms for AGC Considering the Dynamic Wake Effect

Though the dynamic WF model is extremely high dimensional in state space, it is observed that the dominant dynamics of spatial-temporal flow data always evolve in low-rank subspaces spanned by a few spatial coordinates [19, 27, 28]. Our goal is to infer a dynamic WF ROM in which optimal control can be readily performed. To this end, intrinsic coordinates $\textbf{z}=\phi(\textbf{x})$ that dominate the WF flow data need to be identified so that high-dimensional x can be mapped to low-dimensional vectors $\textbf{z}\in\mathbb{R}^{n_{z}}$ with $n_{z}\ll(N_{x}\times N_{y})$. At the same time, the system dynamics with control are also required to be identified in the latent space $\widehat{\textbf{z}}_{t+1}=f^{\textrm{ROM}}(\textbf{z}_{t},\textbf{u}_{t})$ so that the control problem can be solved in reduced-order space z instead of the original full state space x.

To this end, the ROM has to fulfill three properties: (1) Reconstruction. The reduced-order states $\textbf{z}=\phi(\textbf{x})$ must capture sufficient information about the full states $\textbf{x}_{t}$ (enough to enable reconstruction). (2) Future state prediction. The dynamic ROM $f^{\textrm{ROM}}$ must allow for accurate prediction of the next few reduced-order states $\textbf{z}_{t+T_{p}}$ in the whole feasible region of control ($T_{p}$ is the prediction horizon in MPC). (3) Control compatibility. The formulation of the WF ROM must be readily incorporated into optimal controller design.

Given that each grid point has local states denoting the current velocity at that point, the meshed full-state snapshot $\textbf{x}_{t}$ that describes the full WF flow field at time instant $t$ exhibits multiscale spatial physics [19]. Such data provide an opportunity for DNN to make an impact in the modeling and analysis of WF flow fields. For data that have spatial dependencies, convolutional neural networks (CNNs) have been shown to be effective in learning informative high-level features [29]. Another motivation for adopting the DNN for the WF ROM is its flexible architecture that allows for accurately fitting temporal dynamics in the desired form without resorting to hand-designed features [28]. In the following, we will introduce the proposed DNN-based WF ROM in detail.

To learn an appropriate WF ROM meeting the three requirements, in this paper, we propose a specialized DNN architecture. The proposed DNN belongs to the class of autoencoders. The CNN is utilized as the encoder and decoder layers for its power of handling spatial features to extract the intrinsic coordinates that dominate the WF flow data. A linear dynamic system in the reduced-order space is embedded into the neural network as dynamic bottleneck layers to support the optimal control implementation.

To identify the state equation of WF ROM, a sequence of time snapshots $[\textbf{x}_{1},\ldots,\textbf{x}_{T}]$ and corresponding control inputs $[\textbf{u}_{1},\ldots,\textbf{u}_{T}]$ are required as the training data set; $T$ is the size of the sequential training data set. Fig. 1 shows a sketch of the proposed dynamic autoencoder architecture. The high-dimensional wind flow states $\textbf{x}_{t}$ are compressed through the encoder network $\phi^{\textrm{en}}(\textbf{x}_{t})$, which serves as the compression mapping and produces the low-dimensional state $\textbf{z}_{t}$. We denote the encoder as $\textbf{z}_{t}=\phi^{\textrm{en}}(\textbf{x}_{t})$. To learn corresponding dynamic in the latent space, the latent state $\textbf{z}_{t}$ is fed into the designed reduced-order state equation, $\widehat{\textbf{z}}_{t+1}=f^{\textrm{ROM}}(\textbf{z}_{t},\textbf{u}_{t})$. For optimal control compatibility, we enforce $\widehat{\textbf{z}}_{t+1}=f^{\textrm{ROM}}(\textbf{z}_{t},\textbf{u}_{t})$ as a linear dynamic form based on Koopman theory [30]:

where $A\in\mathbb{R}^{n_{z}\times n_{z}}$ $B\in\mathbb{R}^{n_{z}\times 2N}$ are two global linearization matrices whose parameters are optimized during training along with the neural network parameters. Note that while the dynamics operate on the reduced-order space, it takes untransformed controls into account. The hat symbol $\wedge$ of $\widehat{\textbf{z}}_{t+1}$ denotes the predicted latent states derived from $f^{\textrm{ROM}}$, distinguished from $\textbf{z}_{t+1}$ derived directly from encoder $\phi^{\textrm{en}}(\textbf{x}_{t+1})$. Reconstruction of the full state from $\textbf{z}_{t}$ is performed by passing $\textbf{z}_{t}$ through the decoder network $\varphi^{\textrm{de}}(\textbf{z}_{t})$.

The purpose of the dynamic autoencoder is to map the high-dimensional nonlinear system (6) to a low-dimensional linear system (8). Koopman theory proves that a nonlinear system can be mapped to an infinite-dimensional linear system through the transformation of system states [30], and that the infinite-dimensional linear system can be reduced to a finite-dimensional system if the Koopman invariant subspace is found [30]. This gives the rationality of the system mapping in this paper. Besides, previous researches have observed the fact that the dominant dynamics of spatial-temporal flow data always evolve in a low-rank subspaces spanned by a few spatial coordinates [19, 27]. So a linear low dimension system is potential to capture the dominant dynamics of the high dimensional CFD model. However, there are no analytical methods to directly derive such model reduction [19]. The challenge of identifying and representing Koopman invariant subspace provides motivation for the use of deep learning methods for the powerful learning ability. [28] and [31] have shown previous success of finding Koopman invariant subspace through deep learning. In the proposed architecture, we enforce the encoder and decoder to extract reduced-order states that evolve linearly in time constrained by (8), attempting to obtain a finite-dimensional approximation of the Koopman operator. The encoder and decoder are the approximations of observable functions that span a Koopman invariant subspace.

The complete loss functions used in training the dynamic autoencoder can be divided into three parts corresponding to the three high-level requirements for the network. The reconstruction accuracy of the autoencoder is achieved using the following loss:

where $\|\cdot\|_{F}$ represents the Frobenius norm, i.e. square root of the sum of squares of all entries in a tensor, similar to the $2$-norm of a vector. The main purpose of $\mathcal{L}_{\textrm{recon}}$ is to extract the latent states from the high-dimensional data through the encoder and decoder architecture. Thus states at the same time step are utilized to formulate the reconstruction error.

The linear dynamics in WF ROM are constrained through the following loss functions. We enforce that $\widehat{\textbf{z}}_{t+m}$ predicted by the reduced-order state equation (8) to be consistent with $\textbf{z}_{t+m}$ derived from encoding $\phi^{\textrm{en}}(\textbf{x}_{t+m})$:

where $S_{\mathrm{p}}$ is a hyperparameter for how many steps to check in the loss function. $\widehat{\textbf{z}}_{t+m}$ is derived from evolving $\textbf{z}_{t}$ through the dynamic equation (8) $m$ steps forward in time.

To enable optimal control in ROM, accurate future state prediction is necessary. To this end, the state equation must capture accurate dynamic relations between control inputs and reduced-order states. This requirement is reflected through the following loss:

The overall loss function to train the dynamic autoencoder is the sum of the above three parts:

We apply $L_{2}$ regularization to the network parameters. The weights $\alpha$ and $\beta$ are hyperparameters to combine the three losses. The overall flowchart of the training is illustrated in Fig. 2. The encoder consists of the widely used CNN structure, ResNet convolutional layers [32], while the decoder inverts all operations performed by the encoder. The detailed network parameters are elaborated in Section 5.1.

A sketch of the data flow in the training procedure is shown in Fig. 2. In practice, $\textbf{z}_{t}$ is evolved through the linear dynamic equation (8) $S_{p}$ steps forward in time to calculate $\mathcal{L}_{\textrm{pre}}$ and $\mathcal{L}_{\textrm{lin}}$. To learn a reduced order system for optimal control application, the three parts of the loss function are integral and necessary. Minimizing $\mathcal{L}_{\textrm{recon}}$ enforces the mapping is invertible so that the dimension reduction is valid, which is the typical loss function for a standard autoencoder. Minimizing $\mathcal{L}_{\textrm{pre}}$ enforces that the derived dynamic system can accurately simulate the time evolution of original CFD model. Minimizing $\mathcal{L}_{\textrm{lin}}$ enforces the latent state uniqueness. Without enforcing the consistency between $\widehat{\textbf{z}}_{t+m}$ and $\textbf{z}_{t+m}$, (8) may produce a latent state from which reconstruction of x is possible but that is not a valid encoding (i.e. the corresponding latent state derived by the encoder).

In the proposed dynamic autoencoder, the parameters of the matrix A and B are trained with the encoder and decoder parameters by minimizing the loss function (12). The state reconstruction and system dynamics are simultaneously considered in the training procedure. Because the latent states for purely reconstruction can not cover important system dynamic information [19] [28], such integration is necessary for the model reduction of dynamic systems.

To control the WF dynamic power production based on WF ROM, the power output equation in the reduced-order space, $\textbf{P}_{t}=f^{\textrm{pow}}(\textbf{z}_{t},\textbf{u}_{t})$, is required to be identified in the form that can be incorporated into optimal control. The reduced-order states $[\textbf{z}_{1},\ldots,\textbf{z}_{T}]$ derived from the trained encoder, control inputs $[\textbf{u}_{1},\ldots,\textbf{u}_{T}]$ and corresponding wind farm power productions $[\textbf{P}_{1},\ldots,\textbf{P}_{T}]$ are used as the training data set.

Note that although this output relation is straightforward because no dynamics exist, it needs to fulfill the control compatibility requirement. Directly applying the nonlinear output equation in control implementation will cause great difficulty for optimization [20]. We circumvent this problem by imposing a locally linear form on this output equation in the training process:

where $C_{t}\in\mathbb{R}^{N\times n_{z}}$, $D_{t}\in\mathbb{R}^{N\times 2N}$ are local Jacobians, and $o_{t}\in\mathbb{R}^{N}$ is the offset. These matrices $C_{t},D_{t}$ and $o_{t}$ formulate the local linearization of the nonlinear output equation and are determined at each time step as functions of the current reduced-order states and control inputs.

Therefore, instead of directly learning the output equation, a separate neural network is designed to learn the linearization matrices, as shown in Fig. 3. This neural network is a fully-connected one and maps from the current reduced-order states and control inputs $(\overline{\textbf{z}}_{t},\overline{\textbf{u}}_{t})$ to the corresponding local linearization matrices. We denote this neural network as $[C_{t},D_{t},o_{t}]=\Psi(\overline{\textbf{z}}_{t},\overline{\textbf{u}}_{t})$. To ensure the accuracy of the local linearization, in the training process, both $(\overline{\textbf{z}}_{t},\overline{\textbf{u}}_{t})$ and its neighborhood point $(\textbf{z}_{t},\textbf{u}_{t})$ are required. To this end, we adds a random noise of small amplitude to the original $(\textbf{z}_{t},\textbf{u}_{t})$ in the training data set to generate $(\overline{\textbf{z}}_{t},\overline{\textbf{u}}_{t})$. $(\overline{\textbf{z}}_{t},\overline{\textbf{u}}_{t})=(\textbf{z}_{t},\textbf{u}_{t})+\mathcal{N}(0,\epsilon)$. Then, the loss function for training $\Psi$ is as follows:

After training, the local linearization matrices can be directly derived through $\Psi$. In this way, in the optimal control algorithm, the locally linear form (13) can be readily formulated based on the derived linearization matrices.

Through training the proposed two networks, the reduced-order state equation, Eq. (8), and the power output equation, Eq. (13), can be finally derived and compose the WF ROM.

A layout of a $3\times 3$ WF is considered and simulated with WFSim. The typical NREL 5 MW Type III WT is assumed to be installed.

The rotor diameter D $=126$ m [37]. WTs have 5D spacing in the streamwise direction and 3D in the spanwise direction. The simulation parameters are set as follows [16, 14]. The air density $\rho=1.2\textrm{kg}/\textrm{m}^{3}$. The upper and lower bound of the thrust coefficient are set as $0.1\leq C_{T}^{\prime}\leq 2$. The value for $C_{T}^{\prime\textrm{max}}=2$ corresponds to the Betz-optimal value. $C_{T}^{\prime\textrm{min}}=0.1$ indicating that we do not allow turbines to shut down completely. The turbine electrical dynamics whose time scale is in millisecond are not considered in this study, which is common practice in WF level researches [5]. The ramping limit $\triangle C_{T}^{\prime\textrm{max}}$ is $0.2/\textrm{s}$ such that turbines can not de- and uprate instantaneously. The control limits for yaw angles are as follows. $-25^{\textrm{o}}\leq\gamma\leq 25^{\textrm{o}}$. $\triangle\gamma^{\textrm{max}}=0.3^{\textrm{o}}/\textrm{s}$. In the simulation, the simulated wind field is $2520\times 1560\textrm{m}^{2}$, with a grid of $100\times 55$ cells ($N_{x}\times N_{y}$). The corresponding cell size of the discretization is $25.2m\times 28m$. The sample period is 1 s. The incoming wind speed in the experiment is $10\textrm{m}/\textrm{s}$. The wake-loss-heavy situation is studied when the wind direction aligns with the row of turbines.

Analogous to frequency sweeps in system identification [31], we excite the system with sinusoidal control inputs with random frequencies during simulation to generate data. Every 500 s of simulation, the control input frequencies of all turbines are altered and randomly selected. Time snapshots of the current $\textbf{x}_{t}$, $\textbf{u}_{t}$ and $\textbf{P}_{t}$ are stored every one second of simulation. In total, the training set contains $10^{5}$ sequential snapshots of the dynamic WF simulations. To test the generalization ability of the trained models, the validation data set contains 5000 samples generated in the same way.

For the dynamic autoencoder, the WF full-state dimension is $100\times 55\times 2$, i.e. 11000. The reduced-order space dimension $n_{z}$ is set to 20. The architecture of the encoder is shown in Table 1. The widely used CNN structure, residual convolutional blocks [32], with ReLU activations construct the encoder. The decoder inverts all operations performed by the encoder. For the loss function, $S_{p}=50$. $\alpha=300$ and $\beta=1/S_{p}$ are set to balance the three losses so that they are at about the same order of magnitude.

We compare the prediction accuracy of the proposed DNN method with the mainstream classic model reduction method, dynamic mode decomposition with control (DMDc) method. The DMDc method has been recently used to construct the WF ROM by [27, 38]. DMDc tries to devise a low-rank linear system that best approximates the system dynamics. DMD is based on singular value decomposition (SVD) to derive the interpretable modes that characterize spatiotemporal flow data. The dimension reduction is done by a linear projection to the modes provided by SVD. To test the prediction accuracy, WF ROMs derived from the DMDc and DNN approaches are used to predict the future state and power production trajectories. Then, the predicted trajectories are tested against the WFSim simulation results.

300 tests are conducted for each method with a power prediction of 400 time steps forward in each test. Fig. 5 shows the relative error against the full-order WFSim simulation both in the training data set and testing data set. Solid lines represent the mean prediction error across 300 tests, while the shaded regions correspond to one standard deviation about the mean. Both DMDc and the proposed DNN obtain a good accuracy for reconstructing the already observed time evolution in the training data set. However, the DMDc performance degrades considerably in predicting the test data. This occurs because the linear nature of DMDc limits its generalization ability. Besides, SVD encounters difficulty in handling multiscale spatial features [19]. In fact, DMDc can be viewed as a simplified autoencoder with a single unactivated fully connected layer as the encoder. In contrast, the proposed DNN is able to generate stable predictions even in the test data. Its relative error of power prediction remains less than $0.1$. The CNN inherently handles spatial features and the deep autoencoder provides nonlinear coordinations.

The field visualization of streamwise velocity at a snapshot is shown in Fig. 6. The top figure is simulated by the NS-based WFSim. The middle one is predicted and reconstructed by the trained dynamic autoencoder. The bottom is the difference between the two images. In the entire training data set, the mean absolute reconstruction error of the whole $100\times 55$ wind velocity field is $0.3\textrm{m}/\textrm{s}$ with the average relative error of $3.26\%$ and maximal relative error of $9.02\%$. In the test data set, the mean reconstruction error is $0.69\textrm{m}/\textrm{s}$ with the average relative error of $6.41\%$. The results show the errors are much less than the simulated ground truth.

The widely used proportional distribution (ProD) in current studies is illustrated as the benchmark [12, 13]. At each control step, each turbine estimates its available power capabilities based on the measured wind velocity and sends it to the WF controller. Then, the total power demand is distributed proportionally to the available power capabilities of each turbine [12]. No wake interaction is considered. The greedy power $P_{\textrm{greedy}}$ is defined as the time-averaged WF power production with local greedy settings, where $C_{T}^{\prime}=2$ and $\gamma=0^{\textrm{o}}$ for all turbines in the WF [39]. $P_{\textrm{greedy}}$ is considered as the maximal trackable power by previous WF control schemes that do not consider the wake effect, such as ProD [12, 13].

The commonly used power reference signal for standard AGC qualification, the RegD signal coming from the Pennsylvania, New Jersey, Maryland (PJM) electric market [40], is applied to test the control performance, as shown in Fig. 7. Two different scenarios are evaluated. For the top figure case, the time-varying power reference is set to $P^{\textrm{ref}}_{t}=P_{\textrm{greedy}}(0.7+0.3*n_{t}^{\textrm{AGC}})$. $n_{t}^{\textrm{AGC}}$ is the RegD test signal that is normalized to $\pm 1$. Therefore, this power reference will never exceed $P_{\textrm{greedy}}$. For the bottom figure case, the power reference is $P^{\textrm{ref}}_{t}=P_{\textrm{greedy}}(0.9+0.6*n_{t}^{\textrm{AGC}})$. For a period, more power is demanded than $P_{\textrm{greedy}}$.

As shown in Fig. 7, if the power reference is lower than $P_{\textrm{greedy}}$, ProD can provide great tracking performance without considering the wake interaction. However, ProD fails when the power reference is higher. Interestingly, as shown in Fig. 7 (b), from $t=320$ s to $t=370$ s, the WF power produces more than $P_{\textrm{greedy}}$ with ProD control, which occurs because wakes of upstream turbines are not yet fully developed. When the wake arrives at downstream turbines, the available wind power decreases and the power production converges to $P_{\textrm{greedy}}$ from $t=370$ s to $t=600$ s.

In contrast, power tracking can be ensured in both cases by the proposed deep learning-aided MPC. The power regulation capacity is increased compared to ProD. This occurs because the proposed MPC can take into consideration the dynamic wake interactions among WTs in a time horizon $T_{p}$. If $T_{p}$ is long enough to cover the wake propagation time within the given WF, MPC can coordinate each turbine operation dynamically to increase the power production capacity while maintaining the power tracking performance. For this case, the prediction horizon $T_{p}=250$s, the receding horizon $T_{a}=30$s. $\varepsilon=0.02$.

Fig. 8 and 10 further depict the coordinated operation of different turbines for the Fig. 7 (b) case. The yaw angles from upstream to downstream turbines gradually swing to $25^{\textrm{o}}$, $16.3^{\textrm{o}}$ and $0^{\textrm{o}}$ by the proposed MPC. As shown in Fig. 8, the upstream turbines purposely yaw to deflect the wake so that it will not at all or partially overlap the downwind turbines and more power can be harvested. The detailed wind speed profile is shown in Fig. 9. The instantaneous wind speed at $t=600s$ on the dashed line of $y=400\textrm{m}$ in Fig. 8 is plotted. By coordinating WF operation, the wind speed faced by downstream turbines in the MPC case is higher compared to the ProD case. Fig. 10 additionally depicts the individual turbine power. The power production of individual WTs is adjusted dynamically by the WF controller to track the power reference. In both cases, the upstream turbine produces the most power since the wind speed in front of these turbines is the highest. Furthermore, it is proven that more power can be harvested by the downstream WTs in the MPC case.

Through the proposed WF ROM and MPC framework, the trackable power signal range increases to the maximal power extraction of the WF when it operates at its global optimal point considering wake interactions. For the considered 9-turbine WF, the range of trackable AGC signals is improved by $50\%$ compared with traditional ProD control.

With the reduction of system state dimensions, the WF ROM is computational efficient for control applications. On a 3.0 GHz Intel Core with i9 processor, for the 9-turbine WF case, the averaged computation time per simulation time step of WFSim is 167 ms. On the other hand, for the WF ROM with 20 states, the computation time for one time step is only 0.452 ms, which is around $0.3\%$ of the time taken by WFSim. For optimally controlling the WF ROM with the prediction horizon $T_{p}=250s$, each optimization iteration in Algorithm 1 takes approximately 0.48 s with the CVXOPT toolbox in Python. This speed is orders of magnitude faster than optimal control directly applied to the full-order system, witch takes approximately 133 s per iteration [18].

A larger WF is further studied with 25 turbines that have 6D spacing in the streamwise direction and 4D in the spanwise direction, as shown in Fig. 11. To our best knowledge, this is the largest WF case in WF control studies considering the dynamic wake effects. For even larger wind farms, wind farm partitioning methods can be exploited to split the wind farm into smaller parts [41], and the proposed methods can be applied on individual ones. Furthermore, the performance of the deep learning-aided MPC is evaluated with a varying incoming wind speed. The incoming wind speed is shown in Fig. 12, which is a period of realistic measurement data in Denmark and is applied to the west boundary of $x=0\textrm{m}$ as inflow boundary condition. The file of case settings and the wind speed measurement data are also provided online²²2The data has been uploaded to https://github.com/kkxchen/WF-deep-learning-aided-MPC-supplementary-material for easy access .

In this simulation, the simulated wind field is $3500\times 2416\textrm{m}^{2}$, with a grid of $100\times 82$ cells ($N_{x}\times N_{y}$). The corresponding cell size of the discretization is $35m\times 30m$. The position of first turbine at the left bottom is (250m,200m). A flow field is first generated with the uniform wind speeds of 10 m/s in the longitudinal direction. Then, for the considered topology, the flow is propagated 400 seconds in advance with the greedy control of all turbines so that the wakes are fully developed. The flow field obtained after the initialization is utilized as the initial flow field for the simulation results presented in this work. In this experiment, the network architecture of the encoder is shown in Table 2.

The control performance of the proposed method under the varying inflow case is shown in Fig. 13. Through the proposed framework, the complex flow dynamics can be considered in the optimal control efficiently. Compared to the greedy operation, the trackable power signal range is improved by $68\%$ in this case. Besides, the feedback mechanism of MPC makes the proposed control method robust to the incoming wind variation. At every $T_{a}$ time, the current wind field full states are observed and fed back to the controller. The ROM prediction error for the wind field introduced by incoming speed turbulence is thus corrected.

For the 25-turbine WF case, the averaged computation time per simulation time step of WFSim is 210 ms. For the derived WF ROM with 50 states, the computation time per time step is only 1.21 ms. To optimally control the WF ROM with the proposed deep learning-aided MPC, the averaged computational cost is only 1.3 s per iteration in Algorithm 1.