Effective statistical control strategies for complex turbulent dynamical systems

Control of complex turbulent dynamical systems involving strong nonlinearity and high degrees of internal instability is an important topic in practice. Different from traditional methods for controlling individual trajectories, controlling the statistical features of a turbulent system offers a more robust and efficient approach. Crude first-order linear response approximations were typically employed in previous works for statistical control with small initial perturbations. This paper aims to develop two new statistical control strategies for scenarios with more significant initial perturbations and stronger nonlinear responses, allowing the statistical control framework to be applied to a much wider range of problems. First, higher-order methods, incorporating the second-order terms, are developed to resolve the full control-forcing relation. The corresponding changes to recovering the forcing perturbation effectively improve the performance of the statistical control strategy. Second, a mean closure model for the mean response is developed, which is based on the explicit mean dynamics given by the underlying turbulent dynamical system. The dependence of the mean dynamics on higher-order moments is closed using linear response theory but for the response of the second-order moments to the forcing perturbation rather than the mean response directly. The performance of these methods is evaluated extensively on prototype nonlinear test models, which exhibit crucial turbulent features, including non-Gaussian statistics and regime switching with large initial perturbations. The numerical results illustrate the feasibility of different approaches due to their physical and statistical structures and provide detailed guidelines for choosing the most suitable method based on the model properties.


Introduction
Complex turbulent dynamical systems emerge throughout fields in science and technology, including in geophysics, engineering, neural science, and plasma 1 arXiv:2307.15637v1[math.DS] 28 Jul 2023 physics, to name a few [15,16,25,40,43,53,60,62,64,66].These turbulent systems are characterized by high-dimensional state spaces with substantial nonlinear energy transfers between scales [3,43,56].The control of such turbulent dynamical systems has grand importance and broad applications.The goal of control is to design an optimal course of action (the control) to drive the perturbed state of interest back to a desired final target state under minimum cost within a finite time window.For example, active and passive control strategies can reduce the aerodynamic drag of vehicles such as aircraft and cargo ships [2,9,34,55,63] and increase the efficiency of liquid and gas transport through pipelines [5,31,57].Various control strategies have also been designed for applications in industrial mixing and manufacturing [1,11,33,54].In addition, control of turbulent systems has significant implications for climate change mitigation involving large-scale models with high uncertainties [20,29,37,38].However, controlling turbulent systems, in general, has proven a formidable challenge.The central obstacles involve the development of efficient algorithms to deal with the genuinely high dimensionality and a large number of unstable modes.Linear control theory is a well-developed field [10,18,61] and linear control strategies have been applied to turbulent systems under a variety of circumstances where the system can be linearized about a fixed-point and stabilized by the control.For example, linear control with closed-loop feedback from the system has successfully been used to delay the transition from laminar to turbulent flow [6,7,21].However, linear control techniques do not scale well computationally when the dimensionality of the system becomes large [19,32].Thus low-dimensional reduced-order approximations are frequently employed through model reduction and system identification [13,23,30,59,67].Besides linear control, there have been many other innovative approaches to the control of turbulent systems.Machine learning, for example, has been used to design control laws for turbulent systems using data [12,14,19].In addition, there have been open-loop approaches where a predetermined control is applied [57,65].
Statistical control offers a fundamentally different approach compared with the traditional trajectory control methods.Statistical control aims to control certain statistical features of the underlying system.These features can be considered as each statistical moment of the critical model state and are estimated by averaging an ensemble of model states.Note that even deterministic systems can fall under the statistical control framework due to the uncertainties in the initial conditions and observations.These initial uncertainties are propagated and amplified in time as a response to the strong instability considering the turbulent nature of the system.Statistical control has several unique advantages.First, there is no need to exploit exhausting procedures to resolve the full solution of each high-dimensional chaotic trajectory of the underlying turbulent dynamics.The control strategy is achieved effectively by considering only the contributions of the leading order moments.This significantly reduces the computational cost and avoids the randomness of individual trajectories in affecting the control results.Second, although individual trajectories are turbulent, the time evolution of the statistics is deterministic and is more stable.This can be seen by noticing that the Fokker-Planck equation, which is the time evolution of the probability density function (PDF) of the state variables, is always linear despite the associated underlying dynamical system being highly nonlinear and turbulent.Therefore, the statistics are more controllable.Third, statistical control can naturally account for uncertainties and incorporate stochastic reduced-order models [49,50].It allows a large degree of freedom to design suitable strategies for efficient controls.
Statistical energy is a measure of the total statistical mean and variance of the system state across all scales.It is a natural scalar quantity to consider in the context of controlling turbulence [22,42,58].Previous works [48,50,51] have demonstrated that statistical responses in the key states can be successfully controlled using the statistical energy by making use of an energy-conservation principle that appears in numerous turbulent dynamical systems [27,42,43].In particular, exploiting symmetry in the total statistical energy dynamics avoids the inherent nonlinear structure containing instability.Consequently, this approach eliminates the need to track and control a large number of unstable modes.In its current formulation, this statistical control strategy is run in an open-loop manner without requiring online feedback from the system, allowing for the prescribed control forcing to be determined offline for efficient computation.In addition, using the scalar-valued total statistical energy as the control object circumvents the computational issues raised by high-dimensional systems.These factors point to promising applications of statistical turbulent control considering different dynamical features of the targeting turbulent systems.
The statistical control strategy aims to control the statistical energy from a perturbed state back to the target equilibrium state by exerting an external control forcing in the underlying turbulent system.From a high-level description, the strategy consists of two consecutive steps: calculating the optimal energy control and the inversion of a nonlocal control-forcing relation.In the first step, the explicit dynamics of statistical energy are derived using the aforementioned statistical energy-conservation principle.The original control of the high-dimensional system gets reduced to a linear control problem for the scalar energy, which can then be solved using the Hamilton-Jacobi-Bellman (HJB) equations directly [4,8].In the second step, a nonlocal inversion problem is solved to find the deterministic external control forcing in the underlying system, which yields the optimal control of the original turbulent system.This nonlocal inversion problem uses the coupling of the optimal energy control found in the first step with both the deterministic external forcing and the response of the statistical mean of the system to the external forcing.Direct simulation of the response of the mean would be prohibitively expensive, so a crude first-order linear response approximation was employed in previous works instead [41,46,50].Notably, the second-order feedback term in the full control-forcing relation was truncated.This simplifies the analysis and remains consistent with the firstorder approximation valid for small amplitude perturbations within the linear regime.This approach effectively controlled the statistical energy to the target equilibrium state from small initial perturbations.
This paper aims to develop new statistical control strategies for scenarios with more significant initial perturbations and stronger nonlinear responses, allowing the statistical control framework to be applied to a much wider range of problems.Note that the second-order term in the control-forcing relation, consisting of the product of the external forcing perturbation and the mean response to the forcing, is significant for large initial perturbations from the equilibrium state and thus cannot be neglected.While this term can be truncated for small perturbations, it must be included to guarantee proper performance under most large perturbations.In scenarios where the initial mean perturbation is the dominant component of the initial energy perturbation, the inclusion of the second-order term is reflected by the initial external forcing perturbation prescribed by the strategy.Even when the initial mean perturbation is small, the required strong external forcing coupled with dominant nonlinear terms to efficiently control the system usually has a correspondingly strong mean response.In such a case, the external forcing perturbation is strongly influenced by the second-order term in the control-forcing relation.
To address these difficulties, two new statistical control methods are developed in this paper.First, the higher-order methods, incorporating the secondorder term, is developed to fully resolve the control-forcing relation given the mean response.The corresponding changes to recovering the forcing perturbation effectively improve the performance of the statistical control strategy in most test cases.Second, the accuracy of the mean response used in the existing statistical control methods to the external forcing also dramatically impacts the performance of the statistical control strategy.With large perturbations, although linear response theory can provide reasonable results in some special cases [28,46,52], the assumptions justifying the use of linear response theory breaks down and the existing methods often lead to significant errors.In particular, the mean linear response is inadequate when the system is perturbed into a different dynamical regime than the equilibrium state.Due to these limitations, a mean closure model for the mean response as an alternative to the mean linear response is developed in this work.The mean closure model is based on the explicit mean dynamics given by the underlying turbulent dynamical system.The dependence of the mean dynamics on higher-order moments is closed using linear response theory but for the response of the second-order moments to the forcing perturbation rather than the mean response directly.Despite still incorporating linear response theory, the introduction of explicit dynamical information from the underlying model allows the mean closure model to better reflect the properties of the perturbed regime compared to the mean linear response, which only contains information from the equilibrium statistics.
The rest of the paper is organized as follows.Section 2 reviews the general strategies of energy-conserving turbulent dynamical systems, including the relevant assumptions and properties needed for the statistical control strategy.The statistical energy control problem is formulated in Section 3 along with the strategies for recovering the optimal forcing perturbation from the optimal energy control, namely combining the low-order or high-order methods with either the mean linear response or the mean closure model.In Section 4, the control strategies are evaluated in detail based on two prototype models.The first model is a prototypical test model that can exhibit various behaviors and dynamical regimes.The second model is the Lorenz '96 model, which is highdimensional and exhibits multiple dynamical regimes based on the magnitude of the external forcing.Section 5 discusses the results and provides guidance and suggestions for when the various strategies should be applied.Lastly, Section 6 concludes the paper and offers potential future research directions.

Background on Statistical Modeling 2.1 Statistical formulation of the turbulent systems
Turbulent dynamical systems with quadratic energy-conserving nonlinearity can be represented in the following general canonical form [43,45,49] on the state variable u ∈ R N satisfying the dynamics Above, the linear component of the operator is decomposed into two matrices: a skew-symmetric matrix representing linear dispersion effects, L, and a negative definite matrix representing dissipation effects, D. The quadratic nonlinearity is given by a bilinear operator, B(•, •), which satisfies the following energy conservation law, Here "•" denotes the Euclidean inner product.The last two terms of equation (1) represent the external forcing of the system, which is composed of the deterministic component of the forcing, F(t), and the random component, σ(t) Ẇ(t), where Ẇ is Gaussian noise.
It is useful to decompose the state u into a deterministic mean state, ū(t) = ⟨u(t)⟩, and the stochastic fluctuations about each mode where ⟨•⟩ denotes statistical expectation, and e k is the predetermined orthonormal basis.The covariance matrix of u is defined as R(t) = ⟨ZZ * ⟩ where Z = (Z 1 , . . ., Z N ) ⊺ and • * denotes the conjugate transpose.Using the above mean-fluctuation decomposition of u and equation (1) the dynamics of the mean of u can be explicitly written as The operator L u incorporates the energy transfers between modes from the linear dispersion and dissipation effects Importantly, note that the mean dynamics given in equation ( 4) are not closed due to the dependence on the covariance through the interactions with the nonlinearity.Further, even by including the next-order covariance dynamics from the equation.The covariance equation for R can be derived as Q F is the energy flux that accounts for the energy transfer from higher-order non-Gaussian statistics, and ) is positive definite and gives the energy transfer from the stochastic component of the external forcing.
The system with (6) is still not closed due to the third-order moments that appear in the energy flux term Q F .This means that the statistical mean response to external forcing cannot be fully resolved in this hierarchical approach, so in practice, various approximations and closures are needed [49].

Statistical Energy and Response to External Forcing
One quantity of primary interest is the total statistical energy of the system, defined as a combination of the energy in the mean and total covariance Here, ū is the mean vector of u, and tr(R) is the trace of the covariance matrix.
In addition to the energy conservation principle given in equation ( 2), the following assumptions, detailed in [42], are needed to formulate the dynamics of the statistical energy E, that is, and e i • [B(e j , e i ) + B(e i , e j )] = 0 for any i, j.
The above identities characterize the general symmetry in the system that the self-interactions and the closed interactions between pairs of modes vanish under the quadratic nonlinearity.To simplify the notation used in this discussion, we also assume uniform damping D = −dI with d > 0. Under these assumptions and using equations ( 4) and ( 6), the total statistical energy satisfies Statistical energy generally decays to the equilibrium state exponentially.Notably, the dynamics depends only directly on the external forcing and first-order mean state, not the covariance or higher-order moments.This allows for controlling the response to external forcing from determining the total statistical energy and solely considering the external forcing and the response of the mean to the forcing.
To determine the response of the statistical energy to perturbations of the deterministic external forcing, denote the statistical energy of the system under the equilibrium distribution as E eq and denote the energy perturbation as E ′ (t) = E(t) − E eq .The equilibrium energy satisfies dE eq /dt = 0, so, using equation ( 11), the equilibrium energy can be explicitly computed using only first-order mean state by Using the above equation and further denoting the deterministic forcing perturbation as δF = F − F eq and the corresponding mean perturbation as δ ū = ū − ūeq , the energy perturbation, E ′ , satisfies Again, for simplicity, assuming that the stochastic component of the external forcing is not perturbed and does not contribute to the energy perturbation dynamics.It is useful to further decompose the response of the energy perturbation for each mode, i.e.
where κ k is the kth component of δF.Likewise, F eq,k , ūeq,k , and δ ūk are the kth components of F eq , ūeq , and δ ū respectively.To emphasize the central role of the forcing perturbation in each component, the mean response, δ ūk , dependence on the forcing perturbation, κ = (κ 1 , . . ., κ n ) ⊺ , is noted explicitly in the above equation.

Methods on Statistical Control
This section describes the control strategies for high-dimensional turbulent systems, including large-amplitude perturbations.The method is generally split into two consecutive steps: i) the recovery of the optimal energy control for the total statistical energy, and ii) the attribution of the forcing contribution for each detailed spectral mode.Especially, high-order accuracy is achieved by considering different ways to approximate the mean responses and high-order feedback in the control.We illustrate the general idea in the diagram in Figure 1.

Optimal Control of the Perturbed Energy
As the first step of the statistical control strategy, the objective is to drive the total statistical energy from a perturbed state back to a target equilibrium state with a minimized cost through prescribing the deterministic external forcing perturbation in each mode, κ k .The direct optimal control of the energy dynamics is achieved by controlling the much simpler scalar equation independent of the full dimensionality of the system.In the second step, the external forcing perturbation that yields this optimal control will be recovered by attributing the contribution to the total energy from each individual mode.
In this and future sections, the energy perturbation will be denoted as E to simplify the notation.Define the energy control problem as where the objective is to control the energy perturbation, E, back to zero over the time interval [0, T ] using the controls C k , representing the total contribution to the energy perturbation from each mode.The cost functional of this control problem is proposed as where α k gives the relative weights between each control and the total energy perturbation.k T is the cost coefficient for the final energy perturbation from the equilibrium state at time T .From the above setups, the control of total energy becomes a standard linear control problem with a quadratic cost, and so the Hamilton-Jacobi-Bellman (HJB) equation [4,8] can be applied to find the optimal control C * k .In addition, the total statistical energy is a scalar quantity, so the associated energy control problem is tractable even when the underlying system is high-dimensional.
Solving the HJB equations [51] leads to the following Riccati equation which is solved backward in time.The quantity K is a factor of the value function that appears in the HJB equations.The full details of this application of the HJB equations can be found in [51].Using the solution of K, the optimal response of the energy perturbation, E * , is given by the forward equation Finally, the optimal control in each mode, C * k , can be calculated as

Inversion of the Control-Forcing Relation
The goal is to find the forcing perturbation in each mode κ k that yields the optional control discovered from the energy perturbation.For simplicity in notation, the optimal control of each mode found in the previous section will be denoted by C k , where the superscript " * " is dropped.Using equation ( 13) yields the following relation between the energy control and the forcing perturbation in each mode: The mean perturbation response, δ ūk , depends on the forcing perturbation κ.
Inverting this relation for κ k involves solving an ODE system which couples the optimal control, C k , the forcing perturbation, κ k , and the mean response, δ ūk .

The Low-Order Method
In previous works, small forcing perturbation and mean state response are assumed.Thus, the second-order term, κ k • δ ūk , is omitted in equation ( 22) so that the relation only accounts for the dominant leading-order response of the energy to the forcing perturbation [48,51].This assumption is justified by using leading order approximations for the mean response for small initial perturbations, which introduce O(δ 2 ) errors in the same order as the second-order perturbation term.In doing so, the inversion relation is linearized, leading to simplified analytical solutions.The ODE resulting from this linearized relation, referred to in this paper as the "low-order method", is given by Here dC k /dt can be explicitly calculated using equations ( 17), (19), and ( 21) The term δ ūk (0) denotes the mean perturbation in the initial perturbed state.Solving this ODE system involves approximating the response of the mean, dū k /dt, to the forcing perturbation.Strategies for computing the mean responses are detailed in sections 3.3.

The High-Order Method
The second-order term in the control forcing relation was truncated in the loworder method.However, large errors might be introduced due to this omission when the perturbation amplitude grows large.Still, the same analysis can be Step Step 1 is the calculation the optimal control, C k , for each mode.First, a Riccati equation, equation (17), is solved.This is used to calculate the optimal energy response using equation (19).The optimal control is then calculated using equation (21).
Step 2 consists of finding the forcing perturbation, κ k , which yields the optimal control in each mode.Inverting the control-forcing relation involves solving coupled equations for the forcing and the mean response to that forcing.There are two choices for the forcing equations: the low order equations, equation (23), and the high order equations, equation (26).For the mean response there are two strategies: a linear response for the mean, equation (28), and the mean dynamics with a closure for the higher-order moments, equation (35).Choosing one strategy from each category yields four strategies total.done by including this additional term.The ODE resulting from inverting equation (22), including all contributing terms, is given by where dC k /dt is given by equation ( 25) and δ ūk (0) is the initial mean perturbation in the perturbed state.Compared to equation (23), equation ( 26) contains one more perturbation term that accounts for the higher-order contributions to the energy response.For large perturbations from the target equilibrium state, the high-order feedback term becomes necessary to accurately recover the true forcing forms.The inversion of the full control-forcing relation will be referred to in this paper as the "high-order method" in which the response of the energy to the forcing and mean perturbations is fully resolved.On the other hand, the extra terms in the denominator may lead to additional numerical complications, but the improved performance is generally worth the additional sophistication.
We will discuss the performance in the detailed numerical tests in Section 4.

Different Strategies to Recover Mean Responses
The inversion of the control-forcing relation given by equations ( 23) and ( 26) requires the response of the mean to the forcing perturbation.Unfortunately, a direct simulation for the mean responses would be prohibitively expensive considering the extremely high dimensional problem, so approximations for the mean response are developed instead.In this subsection, we develop two approaches to efficiently estimate the mean responses without directly solving the full equations.The first strategy uses linear response theory as a convenient way to recover the leading-order mean response based on the Fluctuation-Dissipation Theorem (FDT) [35,39].The second approximation provides a more accurate high-order approximation based on solving the explicit mean dynamical equation with a high-order closure.

Mean Linear Response to Forcing
The linear mean response to the forcing perturbation κ is computed by where the linear response operator is given by and G ℓ (u) = −p eq div u (e ℓ p eq (u)), (30) where p eq (u) is the equilibrium probability density of u.The basis vector e ℓ corresponds with the ℓth mode.The forcing perturbation δF p,ℓ is the constant external forcing perturbation corresponding to the initial perturbed state.This is contrasted with κ, which is the time-dependent forcing perturbation based on the control that takes over at time t = 0.The linear response operator only depends on the PDF of the target equilibrium state.This explicit formula provides a great computation reduction for convenient calculation of the mean responses but only in the leading order.Then, the linear response approximation is incorporated into the low-order and high-order methods given in equations ( 23) and (26).The equations of the mean response are given by with initial condition While the linear response includes contributions from all modes, typically, only the contribution from a few modes is relevant, allowing the rest to be truncated to further save the computation cost.As further notice, the linear response operator can still be difficult to calculate.In this paper, we adopt the quasi-Gaussian approximation [24,35,44] where the equilibrium probability density, p eq (u), is approximated by a Gaussian distribution.In this case is used in equation (29).Notice this approximation is quasi-Gaussian since the higher-order moments will still be involved due to the probability expectation in (29) based on the true solution rather than only the Gaussian closure.There are many other strategies for approximating the linear response operator [50].While the quasi-Gaussian approximation is more than sufficient for our purposes, other strategies can be investigated.The previous works on the statistical control strategy have used an exponential fit based on an information theory criterion to approximate the linear response operator This is justified with a quasi-Gaussian approximation and a diagonal covariance matrix where the linear response operator reduces to an autocorrelation function of u, which frequently have exponential and oscillatory structures.While this exponential fit is not utilized in this paper, it is noted for its potential application as a useful approximation in practical settings and for the theoretical convenience of having an explicit form of the linear response operator.

Mean Dynamical Equation Closure for Response to Forcing
While directly using linear response theory for the mean response provides a useful approximation, as in section 3.3.1, the linear response is proved inadequate in many cases with large perturbations.This paper develops another strategy to incorporate the mean response based directly on the mean dynamical equation explicitly given in equation ( 4) subject to the forcing perturbation κ.For this method, the dependence of the mean dynamics on the secondorder moments through the nonlinearity is closed using a linear response for the covariance.Compared to the linear response of the mean, which only incorporates statistical information from the target equilibrium distribution, the use of the mean dynamical equation incorporates crucial additional dynamical information.This can be particularly useful when the system is perturbed into a different dynamical regime.
The mean dynamical equation in the mean closure method is given by ( 4) This equation still requires the solution of the second-order covariances R inversely dependent on the form of the forcing perturbation κ.A closure model is constructed using linear response theory for the response of the covariance, R ij , to the external forcing where the linear response operator for the covariance is given by and G ℓ (u) = −p eq div u (e ℓ p eq ).
The higher-order closure of the mean equation enables a better characterization of the mean responses respecting its explicit nonlinear dynamics.The errors from the linear response approximation then appear in the second-order covariances rather than the first-order mean.The linear response in (37) will require the computation of lagged third moments, adding more non-Gaussian information into the approximation.As in Section 3.3.1,a quasi-Gaussian approximation is used for the linear response operator to efficiently compute the response operators.

Numerical Results
We examine the performance of the methods developed in Section 3 using detailed numerical tests.Combining the high-order or low-order methods with the mean equation closure or mean linear response gives a total of four approaches to compare.These methods concern different aspects of component approximations of the statistical control strategy, and the choice to implement each one can be made accordingly based on the specific problem.These four strategies are evaluated on two complex nonlinear models exhibiting various dynamical and statistical behaviors.The first test model is a prototypical triad nonlinear model [49] focusing on a generic coupling between three modes of a turbulent system.It can exhibit a wide variety of nonlinear and non-Gaussian behaviors.Two regimes of this model are considered, including a nearly Gaussian regime with nonlinear energy transfers between modes and a non-Gaussian regime exhibiting an energy cascade representing the transition to turbulence.The second test model is the classic Lorenz '96 model [36] with 40 dimensions which shows multiple dynamical regimes depending on the magnitude of the external forcing.Large perturbations inducing regime switching will be the primary consideration.These test models will illustrate the differences between the statistical control strategies.
The experimental setup is as follows.First, the external forcing perturbation is calculated offline using different statistical control strategies.The optimal energy control is calculated using equations ( 17), (19), and (21).The forcing perturbation which yields this optimal control is found using either the highorder method given in equation (26) or the low-order method given in equation (23).The mean response to the forcing perturbation is approximated by the mean equation closure model described in equations (35) and (36) or by the mean linear response in equations ( 28) and (31).Second, the external forcing perturbation is applied to a Monte Carlo simulation of the underlying dynamical system.An initial ensemble of model trajectories of size M = 1 × 10 4 is drawn from the initial distribution.The perturbed initial state is created by a deterministic forcing perturbation.The deterministic component of the external forcing F eq is perturbed by a constant forcing amplitude δF p to drive the ensemble is forced into a perturbed state away from the original statistical equilibrium.
At the time t = 0, the control strategy takes over.Thus the previous forcing perturbation is replaced by the statistical control forcing, κ.The response of the statistical energy to the forcing is calculated from the ensemble and is tracked as the system is controlled back to the equilibrium state.The energy responses from each strategy are compared to the theoretically optimal energy response in equation (19).The uncontrolled case where no control forcing perturbation is applied, in which the energy naturally decays back to the original equilibrium state, is also used as a point of comparison.Other quantities, such as the mean response, variance response, and empirical control, are also compared for tracking the performance.

A Prototype Nonlinear Triad Model
The first test model is a prototypical 3-dimensional model equipped with a quadratic energy-conserving nonlinearity, which is referred to in this paper as the triad model.The triad model represents a generic nonlinear coupling between three variables universal in turbulent flows [17,47,49].Such a triad interacting structure would emerge as the bare truncation of three identified modes from a high-dimensional turbulent model.The triad model can also generate a wide variety of nonlinear and non-Gaussian behaviors due to the dominant role of the nonlinear coupling term.This sets a desirable first test model to evaluate the skills of the different proposed approaches considering the high-order contributions in the mean state and the energy equation.Despite the nonlinearity, the triad model is analytically tractable in terms of the equilibrium statistics when the linear parts have certain special structures [47], making it an appropriate test model used to have an in-depth study of the different features of the four proposed strategies.
The state variables of the triad model are represented by u = (u 1 , u 2 , u 3 ) ⊺ with governing differential equations: In addition, the quadratic coupling coefficients satisfy which ensures the general energy-conserving property (2) of the turbulent dynamical system framework.Figure 2 shows two typical regimes of the triad model used to evaluate the strategies: one has near-Gaussian statistics, while the other is highly non-Gaussian.

Control on a Near-Gaussian Regime
The first test regime for the triad model is a near-Gaussian regime which nonetheless contains strong nonlinear energy transfers between modes to reach the equipartition of energy.Sample trajectories and equilibrium distributions for this regime are pictured in Regime I in Figure 2. Nearly Gaussian and weakly non-Gaussian features are common in practice, such as in fully turbulent flow with strong mixing.The damping coefficients for this regime are perturbed by δF p,3 = −4 until time t = 0.While only F 3 is perturbed initially, all modes are used to control the system back to the equilibrium state.
The control of the system from the perturbed state back to the equilibrium state under different strategies is shown in Figure 3.All the methods show faster convergence than the no-control scenario to efficiently return the unperturbed equilibrium state.Notably, the high-order method achieves the most accurate near-optimal performance, particularly the high-order method with a mean equation closure.In contrast, the low-order methods overshoot the response incurring relatively higher costs.This is the first confirmation of the crucial role of the higher-order correction in the energy equation when nonlinearity is dominant, as in the triad system.There is a stark difference in the control forcing perturbations between the high-order and low-order methods which can be seen by κ 2 and κ 3 in panel (d) of Figure 3.The initial forcing perturbation for κ 2 differs entirely from the low-order and high-order methods.This is because the corresponding initial mean perturbation, δ ū2 , pictured in panel (e), is relatively large; thus, the initial contribution of the second-order κ 2 • δ ū2 term to the control-forcing relation is quite significant.Indeed, the low-order method overcompensates for lacking this high-order term with a large initial forcing perturbation.In contrast, the initial forcing perturbation for the high-order method is significantly smaller.For the third mode, κ 3 , the initial forcing perturbation is comparable between the high-order and low-order methods due to the initial mean perturbation δ ū3 (0) being relatively small.However, the high-order method produces a stronger forcing perturbation in κ 3 shortly after the initial time.This is explained by the observation that in all methods, the mean perturbation response, δ ū3 , quickly takes on a large value.So the second-order κ 3 • δ ū3 term in the control-forcing relation is non-negligible.Because of this, only the high-order method can account for the contribution of the second-order term to the energy response.One can also see a difference between the initial forcing perturbation for κ 1 between the mean linear response and mean equation closure methods.This is an example where the mean linear response does not accurately produce the initial mean perturbation, while the mean closure model directly incorporates the initial mean perturbation.

Control on a highly non-Gaussian Regime
The second test regime, shown in Regime II of Figure 2, has highly non-Gaussian statistics and intermittency.It features an energy cascade from u 1 to u 2 and u 3 reminiscent of the transition to turbulence.The damping coefficients are and the linear dispersion coefficients are L 1 = 0.03, L 2 = 0.02, and L 3 = −0.01.The unperturbed deterministic external forcing is given by F 1 = F 2 = F 3 = 2 and the stochastic external forcing coefficients are σ 1 = 2 and σ 2 = σ 3 = 1.The perturbed state is achieved through constant deterministic forcing perturbations δF p,1 = δF p,2 = δF p,3 = 2 until time t = 0. Figure 4 shows the results of applying the control strategies to Regime II of the triad model.We focus on the performance of the dominant mode u 1 .The other two modes u 2 u 3 have qualitatively similar performance and are omitted for a cleaner representation.Similar to the near-Gaussian regime, there is a strong forcing perturbation in κ 1 for the high-order method after the initial time.The initial forcing using different methods is comparable due to the relatively small initial mean perturbation δ ū1 (0) but the mean response shortly after the initial time requires the high-order method to capture the subsequent response by the energy.In addition, this example illustrates how the mean equation closure model can produce more accurate forcing perturbations under a strongly nonlinear non-Gaussian regime even when the initial perturbation is similar.As expected, stronger non-Gaussianity requires a more accurate calibration of the mean responses taking into account the higher-order statistics.The high-order equation closure gains a more accurate estimation of the mean state, thus leading to the most accurate result.For low-order methods, lacking the higher-order correction term often leads to larger errors.We suspect that the agreement in the low-order linear response approach comes as an accidental cancellation of errors.

A High-Dimensional Model with Multiple Regimes
The Lorenz '96 model is a standard test model which mimics geophysical waves and exhibits phenomena such as mid-latitude baroclinic instability [36].The shows the autocorrelation functions (ACFs) for each regime.Note that the ACF for the F = 5 regime exhibits long-term oscillatory behavior while the ACF of the F = 8 regime decays very fast.Panel (c) shows the equilibrium distribution for each regime.The F = 5 regime is highly non-Gaussian while the F = 8 regime is nearly Gaussian.
model is defined in a 40-dimension vector state by Here, the variables are indexed periodically, e.g., u 41 = u 1 .The external forcing F is the same for each mode, and so the equilibrium statistics of the system are spatially invariant.Note that quadratic nonlinearity satisfies the energy conservation law given in equation ( 2), which can be shown by the symmetry of the nonlinearity in equation (42).
A key property of the Lorenz '96 model is that it exhibits a variety of dynamical and statistical regimes by altering the value of F .Multiple dynamical regimes are typical of complex turbulent systems and represent a classic obstacle to effective control.Figure 5 exhibits two dynamical regimes corresponding to F = 5 and F = 8.The F = 5 regime is weakly chaotic.It has highly non-Gaussian statistics and a long decorrelation time.Meanwhile, the regime corresponding to F = 8 features near-Gaussian statistics and strongly chaotic dynamics with a correspondingly short decorrelation time.Previous results [48] have shown the statistical control strategy to be effective at controlling small perturbations in the Lorenz '96 model back to the equilibrium state.
In the current experiment, we show the efficacy of the strategies on a large perturbation which drive the system into a different dynamical regime.This leads to a much more challenging problem since the model state goes through a statistical transition between two distinctive regimes.The linear response estimation is no longer valid since the model moves far beyond the linear and near-Gaussian regime.The higher-order corrections become necessary to guarantee effective control performance.In this case, the F eq = F = 5 regime is taken as the equilibrium state and δF p = 3 so that the perturbed state is in the F = 8 regime.Figure 6 shows the control strategies for this large perturbation.In this case, the high-order method with the mean dynamical equation closure shows the most effective strategy.In fact, it shows that combining both the high-order and the mean closure methods is essential to achieve good performance.This is a typical example to confirm the necessity of including high-order corrections when nonlinear and non-Gaussian features become dominant.

Further Discussions
The statistical control strategies extend the effective control methods beyond the small perturbation scenario and demonstrate promise for applications to a broader range of turbulent situations.To summarize, they enjoy several attractive features.Using the total energy as the object of control, there is no need to track and control a large dimension of instabilities due to the energyconservation principle.Thus, the computational cost is significantly reduced and independent of the dimensionality of the system.Further, the control can be determined entirely offline and only requires statistical information about the target equilibrium state, which is usually available from history observation data in many realistic applications.show the response of the energy perturbation to the control forcing for the low-order strategies and high-order strategies respectively.The energy perturbation is normalized by the dimension of the system.Panels (c)-(f) show the controls, forcing, mean response, and variance response for each mode.Note the system is translationally invariant, so the corresponding values for each mode are the same.
Here, we discuss several key features in the new high-order control strategies based on the observations from the numerical experiments.

The High-Order Correction
The control-forcing relation, which encodes the energy response to the external deterministic forcing perturbation, has a second-order perturbation term, δF • δ ū, which accounts for the higher-order contributions of the deterministic forcing perturbation to the energy response.This term consists of the product of the forcing perturbation and the mean perturbation in response to the forcing.Under the circumstances with small perturbations from the equilibrium state, both the mean perturbation and the forcing perturbation are small, so this term can be truncated without compromising the accuracy of the energy response.This method, where only the leading-order contributions to the energy response are considered, is the low-order method.However, for most large perturbations, the second-order term becomes large and significantly affects the energy response.The high-order approach incorporates this second-order term in the inversion of the control-forcing relation, fully resolving the energy response given the forcing perturbation and mean response.We explore several circumstances where the second-order perturbation term significantly impacts the energy response; thus, adding the high-order method can yield significant improvements over the low-order method.
When the initial mean perturbation is large, the initial value of the external forcing perturbation is greatly affected by the presence of the second-order term.This can be seen in equations ( 23) and (26) where the initial condition for the high-order method includes an extra term for the initial mean perturbation δ ū(0).The effect of the initial mean perturbation on the resulting forcing perturbation can also be seen in κ 2 of Figure 3 in the control of Regime I in Section 4.1.However, a large initial mean perturbation is unnecessary, and the second-order term can still have a significant effect even when the initial mean perturbation is relatively small if there is still a large initial energy perturbation due to the variance.Because the initial energy perturbation is large, the relative balance of the mean and variance in the total energy perturbation can shift over time, resulting in a potentially large mean perturbation after the initial time.In this case, the second-order term significantly impacts the evolution of the forcing perturbation even with the same initial conditions.Several examples of this phenomenon can be seen in the numerical tests in Sections 4.1 and 4.2, especially in Figure 6, where a drastic phase transition is shown.
As a further comment, whether to use the low-order or high-order methods can be made independently for each mode.For example, in a multiscale system, the effect of the second-order term may be small relative to the total energy response for small-scale modes and only be significant for larger-scale modes.In this case, the high-order method could be applied to only a subset of largescale modes, while the low-order method is used for the rest, simplifying the dynamics in those small-scale modes without compromising performance.

The Mean Closure Equation
The response of statistical energy to the external forcing depends directly on the mean response to the forcing.It is indirectly linked to the higher-order moments through the mean dynamical equation.This property is critical to formulating the statistical control strategy, allowing for attributing an external forcing perturbation to the optimal control by solving the control-forcing relation.Therefore, accurately approximating the mean response to external forcing is vital to the success of the statistical control strategy.Linear response theory effectively approximates the mean response for small perturbations, and it can perform well for larger perturbations in some cases when non-Gaussian statistics is not so important.However, its skill degenerates when the system is largely perturbed to a different dynamical regime where the linear response operator, based solely on the unperturbed dynamics, can provide very little information for the future perturbed state.
Using a mean closure equation for the mean response, which directly incorporates dynamics from the model, is expected to show improved performance when the initial perturbation spans multiple dynamical regimes.A mean dynamical closure equation is based on the explicit mean dynamics given in equation ( 4), in which the dependence on higher-order moments is closed using a suitable approximation.The closure considered in this paper utilizes a linear response for the higher-order contribution of the covariance described in equations ( 35) and (36).While this mean closure equation still relies on the linear response for the covariance, the mean dynamics still provide more information about the perturbed regime than the mean linear response.
The initial forcing perturbation is affected by the choice of mean response.The mean linear response cannot directly use the initial mean perturbation and instead must use the initial mean perturbation predicted by the linear response to a constant forcing perturbation.This is necessary to guarantee the convergence of the forcing perturbation to zero in the linear response case.The mean closure model, however, can utilize the initial mean perturbation directly.This is illustrated by κ 1 in Figure 3 in Regime I of Section 4.1 where the initial mean response differs between the linear response and mean closure methods.In this case, the mean equation closure achieves better performance.Even when the initial mean perturbation is similar to the initial mean perturbation predicted by the linear response, the mean closure model can provide more accurate dynamics in many cases.In Section 4.1, the mean equation closure method performs better among all cases, especially in Regime II.In Section 4.2, the Lorenz '96 model is perturbed from a non-Gaussian regime to a near-Gaussian turbulent regime.The mean equation closure again performs better than the linear response in this case with multiple dynamical regimes.

Convergence to the Equilibrium State
The total statistical energy bounds the total mean and variance.Ideally, one hopes that the efficient control of the equilibrium energy will also achieve the shows the forcing perturbations in each mode.Note that the forcing perturbation for the high-order method does not converge to zero.Panel (c) shows the equilibrium distribution, the perturbed distribution, and the alternative distribution achieved by the high-order method that yields the same statistical energy.efficient control of the mean and variance back to the equilibrium state.While this appears to be the case in most applications, this is not mathematically guaranteed by the statistical control strategy.Figure 7 illustrates an example where the optimal energy response is achieved through the high-order mean closure strategy, but the external forcing perturbation converges to a constant non-zero state.Essentially the system converges to a different equilibrium with a different constant external forcing but the same equilibrium energy.
The cost functional given in equation ( 16) only penalizes the strength of the direct energy control C k rather than the external forcing perturbations κ k that yield that control.In addition, the control-forcing relation given in Equation ( 22) admits multiple solutions in the limit in both the low-order and high-order formulations.One natural fix for such an issue is incorporating additional terms into the cost function, for example, explicitly excluding the mean state.

Conclusion
An efficient method of controlling the complex turbulent system with energy conserving nonlinearity is achieved through control of the total statistical energy from a perturbed state back to equilibrium without controlling the large number of multiscale and potentially unstable modes.This paper proposed new statistical control strategies by extending previous works [48,51], which had been restricted to scenarios with small perturbations from the equilibrium state.Incorporating the high-order term in the control-forcing relation accounts for the second-order contribution of the perturbations to the statistical energy, allowing the strategy to account for the response of the energy more accurately to large amplitude external forcing perturbations.Additionally, introducing a mean dynamical closure model allows the statistical control strategy to better account for large perturbations that drive the system into different dynamical regimes whose dynamics cannot be adequately reflected directly by a mean linear response approximation.These strategies allow for the practical application of the statistical control strategy to a wider variety of perturbations and regimes than previously possible.
The field of statistical control theory remains relatively underexplored, providing many promising research directions.The results presented in this paper could be further refined by developing more sophisticated methods for incorporating the mean and covariance dynamics into the mean response.Besides, other designs for mean closure models could be considered, several of which are described in [49].It may also be possible to incorporate other more suitable statistical functionals into the control strategy in addition to the statistical energy.This would allow, for example, the direct control of the statistical mean, giving finer control over the response of the system.Lastly, while the current statistical control strategy is conducted in an open-loop fashion, determining the control offline, a natural extension would be incorporating feedback from the system into a closed-loop statistical control strategy.This could be accomplished, for example, by combining measurements of the actual mean response of the system into the inversion of the control-forcing relation.The computational advantages of the statistical control strategy would be very advantageous in such a closed-loop formulation, which requires real-time incorporation of the model feedback.

1 :Figure 1 :
Figure 1: Schematic diagram of the statistical control strategy.Step 1 is the calculation the optimal control, C k , for each mode.First, a Riccati equation, equation(17), is solved.This is used to calculate the optimal energy response using equation(19).The optimal control is then calculated using equation(21).Step 2 consists of finding the forcing perturbation, κ k , which yields the optimal control in each mode.Inverting the control-forcing relation involves solving coupled equations for the forcing and the mean response to that forcing.There are two choices for the forcing equations: the low order equations, equation(23), and the high order equations, equation(26).For the mean response there are two strategies: a linear response for the mean, equation(28), and the mean dynamics with a closure for the higher-order moments, equation(35).Choosing one strategy from each category yields four strategies total.

d 1 = d 2 = d 3 = 1 . 3 = − 1 . 1 Figure 2 : 3 Figure 3 :
Figure 2: The dynamics and equilibrium distributions of the prototype triad model under two different regimes.Panels (a) and (c) show sample trajectories of each regime of the model.Panels (b) and (d) show the equilibrium marginal distributions of the state variables as well as their pairwise joint distributions in each regime.The nonlinearity in the model produces non-Gaussian distributions in each regime.In particular Regime II exhibits intermittency and highly non-Gaussian statistics.

1 Figure 4 :
Figure 4: Example of controlling a highly non-Gaussian regime in the prototypical triad model.Panel (a) shows the response of the energy to the forcing perturbation for all strategies.Panel (b) shows the control for each strategy for the u 1 mode.Panel (c) shows the forcing perturbation in the u 1 mode.Panel (d) shows the mean response in the u 1 mode.

8 FFigure 5 :
Figure 5: Sample trajectories and distributions of the 40-dimensional Lorenz 96 model for both the F = 5 (weakly chaotic; highly non-Gaussian) and F = 8 (strongly chaotic; nearly Gaussian) regimes.Panel (a) shows sample trajectories for one sample of each regime in the form of the Hovmoller diagram.Panel (b)shows the autocorrelation functions (ACFs) for each regime.Note that the ACF for the F = 5 regime exhibits long-term oscillatory behavior while the ACF of the F = 8 regime decays very fast.Panel (c) shows the equilibrium distribution for each regime.The F = 5 regime is highly non-Gaussian while the F = 8 regime is nearly Gaussian.

Figure 6 :
Figure6: The control of the Lorenz 96 model from the perturbed state of F = 8 back to the equilibrium state of F = 5.Note this is a large perturbation into a regime with very different dynamics and statistics from the equilibrium.Panels (a) and (b) show the response of the energy perturbation to the control forcing for the low-order strategies and high-order strategies respectively.The energy perturbation is normalized by the dimension of the system.Panels (c)-(f) show the controls, forcing, mean response, and variance response for each mode.Note the system is translationally invariant, so the corresponding values for each mode are the same.

3 Figure 7 :
Figure 7: An example where the optimal energy response is achieved, but the system is forced to a different equilibrium state.The triad model has parametersd 1 = d 2 = d 3 = 1, L 1 = L 2 = L 3 = 0, B 1 = 1, B 2 = −0.6,B 3 = −0.4,F 1 = F 2 = F 3 = 0.5, and σ 1 = σ 2 = σ 3 = 0.5.The perturbed state has F 3 = −1.Panel (a)shows the energy response to the low-order and high-order control strategies.Panel (b) shows the forcing perturbations in each mode.Note that the forcing perturbation for the high-order method does not converge to zero.Panel (c) shows the equilibrium distribution, the perturbed distribution, and the alternative distribution achieved by the high-order method that yields the same statistical energy.