Nonlinear optimization problems (NOPs) are widely employed in the atmospheric and oceanic sciences; for instance, variational data assimilation, and sensitivity analysis of the model to initial and parameter conditions ( Ishikawa et al., 2001; Mu et al., 2004; Robert et al., 2006). Consequently, optimization approaches to solve NOPs play an important role in this field, and contribute to advancing progress in the predictability of atmospheric or oceanic systems, as well as improving the forecasting skill of numerical models ( Duan et al., 2014; Mu and Duan, 2003; Zheng et al., 2012b).

Some special approaches related to NOPs have been applied in atmospheric and oceanic studies. For example, the conditional nonlinear optimal perturbation (CNOP) approach is used to investigate predictability and uncertainty ( Mu and Duan, 2003; Mu et al., 2007; Yu et al., 2011), and variational data assimilation (VDA) is employed to find the optimal initial field in order to reduce initial errors in numerical models ( Ishikawa et al., 2001; Robert et al., 2006). As well as predictability problems related to initial errors, those related to parameter errors have also been investigated. For example, optimal estimates of model parameters can be obtained by fitting model results to observations ( Tziperman and Thacker, 1989; Wang et al., 2005). In addition, an optimal forcing vector approach is applied to correct model errors for improving the model simulation ability ( Duan et al., 2014). The above approaches have the potential to improve the forecasting skill and reduce uncertainty in numerical simulations. Traditional optimization algorithms, which require the derivatives of the cost functions, are most commonly employed during optimization processes of these approaches. However, when the derivatives do not exist or cannot be calculated, it is difficult to search for the optimal values using these traditional optimization algorithms, which limits their broad application. Consequently, derivative-free optimization approaches are receiving increasing attention in solving NOPs. They have been used in VDA ( Fang et al., 2009; Gratton et al., 2014; Zheng et al., 2012b) and for optimizing model parameters ( Severijns and Hazeleger, 2005; Wang and Huo, 2013), the results of which have demonstrated their efficiency. The indication is that derivative-free optimization algorithms are powerful and important tools for future research, and thus it is necessary to further explore their application, for instance in calculating CNOPs.

In fact, the CNOP, a nonlinear generalization of linear singular vectors (LSVs), is the optimal value of a nonlinear optimization problem. The CNOP approach was first proposed to find the optimal initial perturbation (CNOP-I)of a nonlinear system ( Mu et al., 2003). Its application was then extended to identify the combined optimal mode of the initial perturbation and model parameter perturbation (CNOP-I and CNOP-P) ( Mu et al., 2010). In short, a CNOP is a kind of initial perturbation or parameter perturbation with maximum nonlinear evolution at the prediction time and satisfying a given constraint. It has been applied to investigate the predictability of ENSO (including the optimal precursor, the effects of nonlinearity on error growth, and the importance of model parameter error to the spring predictability barrier) ( Duan and Mu, 2006; Duan et al., 2004; Mu and Duan, 2003; Mu and Wang, 2007; Yu et al., 2011), the nonlinear growth of instability of the oceanic thermohaline circulation ( Mu et al., 2004), the adaptive observations of typhoons ( Mu et al., 2009), the stability of a grass ecosystem ( Mu and Wang, 2007; Sun and Mu, 2009), and ensemble prediction ( Mu and Jiang, 2008). All these applications illustrate that the CNOP method is a useful tool for predictability and stability studies associated with the atmospheric and oceanic sciences. Because of its explicit physical meaning, the CNOP approach is expected to play an increasing role in weather and climate studies in the future.

Generally, during the maximization process to calculate CNOPs, the derivative-based traditional optimization algorithms, such as the sequential quadratic programming algorithm ( SQP, Barclay et al., 1998), the spectral projected gradient algorithm ( SPG2, Birgin et al., 2000), and the limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm ( L-BFGS, Liu and Nocedal, 1989), are applied and the adjoint system is always used to determine the gradients of the cost functions. However, not all numerical models have available adjoint systems, which limits the broad application of the CNOP approach. In addition, if the cost functions are non-differentiable, these derivative-based optimization approaches cannot be used. Consequently, it is necessary to explore a derivative-free optimization algorithm for calculating CNOPs. For example, Zheng et al. (2012a) applied a genetic algorithm (GA) to get CNOPs and discussed the predictability problems related to “on-off” switches. They stated that the GA was more powerful in determining the optimal values than the traditional algorithm. Recently, a derivative-free global optimization approach, the Differential Evolution (DE) algorithm, has been employed to obtain the CNOP-Ps satisfying box constraint conditions by using a dynamic global vegetation model with intermediate complexity ( Sun and Mu, 2012, 2013). However, this algorithm was just employed to obtain the CNOP in a lower dimensional optimization problem related to terrestrial ecosystem processes. Furthermore, it has not yet been used to deal with a ball constraint for a finite-dimensional optimization problem, which is important in the CNOP approach.

In this study, a projection skill, which is also employed in some traditional optimization algorithms for solving optimization problems with ball constraint conditions, such as the SPG2 algorithm, was explanted to the DE algorithm for calculating the CNOP of the finite-dimensional NOP with a ball constraint condition. The CNOP can be associated with the initial or parameter perturbations; but here, the CNOP related to initial perturbations was computed. To test the validity of the DE algorithm and apply the DE algorithm in a discussion of the predictability of the atmosphere and ocean, a comparison was made using the DE algorithm and the traditional algorithm. Specifically, the DE and SPG2 algorithms were both employed to obtain the CNOPs of the NOPs with ball constraint conditions by using a simple model of Burgers‘equation. The CNOPs obtained by these two algorithms were then compared to identify the applications of the DE algorithm in calculating the CNOPs satisfying ball constraint conditions.

The remainder of the paper is organized as follows. The CNOP approach is introduced in section 2. Brief descriptions of the two algorithms, a simple model of Burgers’ equation, and the experiment design are also provided in this section. The preliminary experimental results are presented in section 3, and a summary and discussion of the key findings is provided in section 4.

As demonstrated in Figs. 1a and 1b, the CNOPs calculated by the SPG2 and DE algorithms had almost identical spatial distributions, expect for small differences in magnitude (see Fig. 1c) when T= 30 s and n= 101. Although differences existed, the nonlinear time evolutions (in terms of their norm squares ) of both CNOPs computed by the two algorithms were in good agreement (see Fig. 2a). By increasing the number of spatial grids from 101 to 201 and 301, we obtained similar results (see Figs. 1d-f and Fig. 2b; Figs. 1g-i and Fig. 2c). This demonstrates that by using a projection skill in the DE algorithm it is possible to calculate CNOPs satisfying ball constraints.

Figure 1

	Figure Option ViewDownloadNew Window
	Figure 1 Spatial distributions of CNOPs (m s-1) calculated by the (a, d, g) SPG2 algorithm and (b, e, h) DE algorithm, and the (c, f, i) CNOP differences (m s-1) (DE minus SPG2) when T= 30 s and n= 101, 201, 301, respectively (DE: the Differential Evolution algorithm; SPG2: the Spectral Projected Gradient algorithm).

Figure 2

	Figure Option ViewDownloadNew Window
	Figure 2 Nonlinear evolutions of the CNOPs calculated by the SPG2 and DE algorithms when T= 30 s and n = (a) 101, (b) 201, and (c) 301.

In order to study the dependence of the DE algorithmon the length of the prediction time in calculating CNOPs, we prolonged the prediction time T to 40 s and 50 s. Here, we only show the CNOPs obtained when n= 101 and T= 40 s, 50 s (see Fig. 3). Under other conditions (i.e., n= 201, 301), the results were similar. Little difference was found between the CNOPs obtained by the SPG2 and DE algorithms, and their nonlinear time evolutions were almost identical (data not shown). These findings illustrate that the DE algorithm could still be used to calculate CNOPs with increasing prediction time.

Table 1 summarizes the optimal values of the cost function in each experiment. At the same grid size and prediction time, the cost function values corresponding to the CNOPs calculated by the different algorithms were nearly the same. We also computed the correlation coefficients between the CNOPs obtained by the DE and SPG2 algorithms, all of which were greater than 0.99.

Table 1

Table 1

Table 1 Summary of the optimal values of the cost function in the 18 numerical experiments. Algorithm T = 30 s T = 40 s T = 50 s n = 101 SPG2 DE -1.2351E-05 -1.2351E-05 -7.1930E-04 -7.1930E-04 -4.2255E-02 -4.2255E-02 n = 201 SPG2 DE -1.0376E-05 -1.0376E-05 -5.7139E-04 -5.7139E-04 -3.1545E-02 -3.1539E-02 n = 301 SPG2 DE -9. 8909E-06 -9.8895E-06 -5.3727E-04 -5.3713E-04 -2.9234E-02 -2.9231E-02

Table 1 Summary of the optimal values of the cost function in the 18 numerical experiments.

In conclusion, the derivative-free DE algorithm is able to calculate the CNOPs satisfying ball constraint conditions by being explanted to a projection skill. At the same spatial grid size and prediction time, the CNOPs obtained by the DE algorithm were nearly the same as those obtained by the SPG2 algorithm in terms of their spatial distributions and nonlinear evolutions. The DE algorithm is a good and reliable choice for obtaining CNOPs in complex models without available adjoint systems, or under the condition that the cost functions are non-differentiable.

Figure 3

	Figure Option ViewDownloadNew Window
	Figure 3 Spatial distributions of CNOPs (m s-1) calculated by the (a, d) SPG2 algorithm and (b, e) DE algorithm, and the (c, f) CNOP differences (m s-1) (DE minus SPG2) when n= 101 and T= 40 s, 50 s, respectively.

Many studies have indicated that the CNOP approach is a useful tool in atmospheric and oceanic studies. However, for complex numerical models without available adjoint systems, or under the condition that the cost functions are non-differentiable, it is difficult or impossible to obtain the CNOPs by using traditional gradient-based optimization methods.

To overcome this difficulty, we introduce the DE algorithm, a derivative-free method, together with a ball constraint condition and a projection skill to calculate the CNOPs by using a simple model of Burgers‘ equation. Through comparisons between the CNOPs obtained by the DE and SPG2 algorithms, we found that for the same NOP, the CNOP calculated by the DE algorithm was almost identical to that calculated by the SPG2 algorithm in terms of their spatial distributions and nonlinear time evolutions. Consequently, the DE algorithm has the potential to play an increasingly important role in calculating the CNOPs in future studies. The numerical results also imply that this algorithm can be applied to obtain the optimal values of optimization problems, especially for non-differentiable, nonlinear, and multimodal optimization problems.

However, the DE algorithm is far more computationally expensive than the SPG2 algorithm. Take the dimension of spatial grids n= 101 and prediction time T= 30 s as an example. During the optimization process, the SPG2 algorithm took less than 100 iterations. Meanwhile, after 100 iterations, the DE algorithm did not converge and the current optimal value of the cost function was just 1.0523E-05 (84.86% of the true maximum 1.2351E-05). At that moment, the current optimal initial perturbation calculated by the DE algorithm and the CNOP had similar spatial distributions except for differences in amplitudes and opposite phases, between which the correlation coefficient was -0.8660. In the subsequent iteration process, the current optimal perturbation obtained by the DE algorithm was converging to the CNOP slowly. In other words, many iterations (meaning more computation time) compensate for the lack of derivative information about the cost function during the search process of the DE algorithm. Nevertheless, with the rapid development of computers and the employment of parallel computing, the computation cost of the DE algorithm will be reduced greatly.

In this study, we only applied the DE algorithm to solve NOPs with dimensions of a magnitude of 100. Practically, the dimensions of NOPs in atmospheric and oceanic sciences always possess a magnitude of 10³. In future, we will try to explore the application of the DE algorithm in these aspects.

Acknowledgments. The authors thank the anonymous reviewers for their suggestions to improve the manuscript. Funding was provided by grants from the LASG State Key Laboratory Special Fund and the National Natural Science Foundation of China (Grant Nos. 40905050, 40830955, and 41375111).