Conditional nonlinear optimal perturbations based on the particle swarm optimization and their applications to the predictability problems

In predictability problem research, the conditional nonlinear optimal perturbation (CNOP) describes the initial perturbation that satisfies a certain constraint condition and causes the largest prediction error at the prediction time. The CNOP has been successfully applied in estimation of the lower bound of maximum predictable time (LBMPT). Generally, CNOPs are calculated by a gradient descent algorithm based on the adjoint model, which is called ADJ-CNOP. This study, through the two-dimensional Ikeda model, investigates the impacts of the nonlinearity on ADJ-CNOP and the corresponding precision problems when using ADJ-CNOP to estimate the LBMPT. Our conclusions are that (1) when the initial perturbation is large or the prediction time is long, the strong nonlinearity of the dynamical model in the prediction variable will lead to failure of the ADJ-CNOP method, and (2) when the objective function has multiple extreme values, ADJ-CNOP has a large probability of producing local CNOPs, hence making a false estimation of the LBMPT. Furthermore, the particle swarm optimization (PSO) algorithm, one kind of intelligent algorithm, is introduced to solve this problem. The method using PSO to compute CNOP is called PSO-CNOP. The results of numerical experiments show that even with a large initial perturbation and long prediction time, or when the objective function has multiple extreme values, PSO-CNOP can always obtain the global CNOP. Since the PSO algorithm is a heuristic search algorithm based on the population, it can overcome the impact of nonlinearity and the disturbance from multiple extremes of the objective function. In addition, to check the estimation accuracy of the LBMPT presented by PSO-CNOP and ADJ-CNOP, we partition the constraint domain of initial perturbations into sufficiently fine grid meshes and take the LBMPT obtained by the filtering method as a benchmark. The result shows that the estimation presented by PSO-CNOP is closer to the true value than the one by ADJ-CNOP with the forecast time increasing.


Introduction
Weather and climate predictability problems are attractive and significant in atmospheric and oceanic sciences.The goal of studying predictability problems is to investigate and understand the reasons for and mechanisms of the prediction uncertainty.Early in 1975, Lorenz divided the predictability problem into two types based on the consideration that prediction uncertainties were mainly caused by the initial and model errors.The first type of predictability problem is concerned with the uncertainty in the forecast results induced by initial errors, and the second type deals with the uncertainty caused by model errors.He further introduced the singular vector (SV) into the predictability problem study.The leading singular vector is the initial perturbation with the greatest linear growth, and, accordingly, can be used to estimate the evolution of initial errors in the tangent linear regime during the course of a forecast (Buizza and Palmer, 1995;Lacarra and Talagrand, 1988;Farrell, 1990;Borges and Hartmann, 1992;Kalnay, 2003).Since the atmospheric and oceanic movement is a nonlinear physical process, an SV based on linear theory cannot effectively express the procedure.This Q. Zheng et al.: Conditional nonlinear optimal perturbations limits its applications to the predictability problems (Mu and Duan, 2003).To this end, Mu et al. (2003) presented the conditional nonlinear optimal perturbation (CNOP) method.CNOP refers to the initial perturbation that satisfies a certain constraint condition and has the maximum nonlinear evolution at the prediction time.Therefore, it stands for the initial uncertainty which leads to the largest prediction error.CNOP has been widely applied to weather and climate predictability studies (Duan et al., 2004(Duan et al., , 2009;;Mu et al., 2007aMu et al., , b, 2009;;Terwisscha and Dijkstra, 2008;Yu et al., 2012;Wang et al., 2012Wang et al., , 2013)).Riviere et al. (2010) used an extension of the CNOP approach, i.e., the nonlinear singular vectors, to estimate the predictability of atmospheric moisture processes so as to reveal the effects of nonlinear processes.
Based on actual demands, Mu et al. (2002) separated the predictability problem into three sub-problems, i.e., the problem of the LBMPT, the problem of the upper bound of maximum prediction error, and the problem of the lower bound of maximum allowable initial error and parameter error.Duan and Luo (2010) formulated these three problems into three constrained nonlinear optimization problems.Meanwhile, they used the CNOP method to find the solutions of the three sub-problems.
Capturing CNOPs is a kind of constraint optimization problem, and optimization algorithms commonly used in solving CNOPs are based on the gradient descent method, including the spectral projected gradient 2 (SPG2; Brigin et al., 2000), sequential quadratic programming (SQP; Powell, 1982), and limited memory BFGS (L-BFGS; Liu and Nocedal, 1989).Among these algorithms, the gradient information is always provided by the backward integral of the corresponding adjoint model of the prediction model (Duan et al., 2004(Duan et al., , 2008;;Mu and Zhang, 2006;Mu et al., 2009;Jiang and Wang, 2010;Yu et al., 2012;Wang et al., 2012Wang et al., , 2013)).But the optimal algorithms based on gradient information involve the forward integral of the tangent model and backward integration of the adjoint model.It might cause the following two problems: (1) for the numerical prediction model with complex physical processes, the validity of the tangent linear approximation cannot be guaranteed when the forecast period is long; (2) for the actual prediction model, it is quite difficult and time consuming to develop the adjoint model.Recently, Zheng et al. (2012Zheng et al. ( , 2014) ) attempted to apply genetic algorithms (GAs) to capture CNOPs of the dynamical model containing discontinuous "on-off" switches.They concluded that GAs, with proper genetic operator configuration, can overcome the non-smooth influences and obtain the global CNOP with high probability.Thus, in non-smooth cases, using GAs to solve predictability problems is more effective than using the conventional optimization algorithm.
The particle swarm optimization (PSO) algorithm is an intelligent algorithm proposed by Kennedy and Eberhart (1995) which imitated the process of bird foraging.In the PSO algorithm, each particle, as a vector in solution space, represents a potential solution of the optimal prob-lem.Compared to the gradient descent algorithm based on the adjoint model, the PSO algorithm has better effectiveness in searching the global optimal solution for nonlinear and non-smooth optimal problems; PSO does not require a gradient of the objective function.PSO is similar to the GA in the sense that they are both population-based search approaches and that they both depend on information sharing among their population members to enhance their search processes using a combination of deterministic and probabilistic rules (Hassan et al., 2005).PSO has the same effectiveness as the GA but with significantly better computational efficiency; it has memory and constructive cooperation between particles (Fang and Zheng, 2009).This paper tries to use PSO to calculate the CNOP of the forecast model, and obtains a highly precise estimation of the lower bound of maximum predictable time.
The Ikeda model was originally proposed by Ikeda (1979) as a model describing light going across a nonlinear optical resonator.The Ikeda model has strong nonlinearity, and the two-dimensional difference scheme is its most common form.Li et al. (2016) investigated the stability of solutions of the Ikeda model and tested the dependence of the solutions on the model parameter.In addition, they provided the solution description of various shapes corresponding to parameter values of different regions.This paper takes the twodimensional Ikeda model as the prediction model to reveal how the nonlinearity impacts the precision when estimating the LBMPT based on the ADJ-CNOP method.A new method, PSO-CNOP, is presented to solve this problem.
This paper is organized as follows: Sect. 2 is devoted to describing the three predictability sub-problems, the definition of CNOP and the two-dimensional Ikeda model.Section 3 provides the knowledge about the particle swarm optimization (PSO) algorithm.In Sect.4, the performances of ADJ-CNOP and PSO-CNOP are compared when solving CNOPs through numerical experiments.The impacts on the estimation precision of the lower bound of maximum predictable time are also demonstrated in this section.The conclusion and discussion are presented in Sect. 5.

Related conceptions and the forecast model
The three predictability sub-problems, the definition of CNOP and the two-dimensional Ikeda model are briefly described as follows, and more detailed introductions can be found in Mu et al. (2002Mu et al. ( , 2003) ) and Li and Zheng (2016).

Three sub-problems of the predictability problem
The three predictability sub-problems associated with the first kind of predictability problem are introduced, and the forecast model is supposed to be perfect in the following.
Problem 1: the lower bound of maximum predictable time (LBMPT) Assuming that there is an error in the initial condition (IC) of the forecast model, it will lead to a prediction error when integrating forward the model from the IC to predict the atmospheric or oceanic states in the future.
Let u 0 be the IC, u t T the true state at time T , and M T the nonlinear propagator of the numerical forecast model from 0 to time T ; then, under the assumption of the perfect model, u t T = M T (u t 0 ), where u t 0 is the true state at the initial time.With a given prediction precision ε > 0, the maximum predictable time T ε is defined as follows: Since the true value u t T cannot be obtained exactly, it is impossible to get T ε by solving the nonlinear optimization problem (1).Inspired by the fact that the IC u 0 is often provided by an analysis field, and the associated analysis error can generally be controlled in a specified range, Mu et al. (2002) reduced the maximum predictable time problem to the following LBMPT problem.
If we have an estimation of the uncertainty in the IC as follows, then the LBMPT T l is defined as where σ > 0 denotes the accuracy of the IC in terms of the norm || • ||, and δu 0 is an initial perturbation superposed on the IC.According to Eq. ( 2), the true initial state is within the constraint region; we have Problem 2: the upper bound of maximum prediction error When a forecast is produced from an incorrect initial IC u 0 , the prediction error at the prediction time T is Similarly to problem 1, since the true value u t T cannot be obtained precisely, Mu et al. (2002) instead introduced the upper bound of the maximum prediction error within the given initial error limitation as follows: (5) Note that u t 0 satisfies Eq. ( 2) and u t T = M T (u t 0 ) under the assumption of perfect model; we have

Problem 3: the lower bound of maximum allowable initial error
Given the prediction time T > 0 and prediction precision ε > 0, the maximum allowable initial error is Similar to problems 1 and 2, the above problem was reduced by Mu et al. (2002) to the following lower bound of maximum allowable initial error:

Conditional nonlinear optimal perturbation (CNOP)
In consideration of the nonlinearity impacts, Mu et al. (2003) introduced CNOPs into the study of predictability problems.Suppose the atmospheric or oceanic motions can be described by the following dynamic system: T is the basic state, which is an n-dimensional vector; the superscript T represents the transpose, U 0 is the initial basic state, and ∈ ⊂ R m and t are the spatial and temporal variables, respectively; t = 0 is the initial time; and F is a nonlinear partial differential operator.
Suppose M τ is the nonlinear transmission propagator from the initial time t = 0 to the forecast time t = τ ; thus, the state of model ( 8) at time τ is If u 0 stands for the initial perturbation of the basic state U(t) and u I (τ ) is the development of u 0 at time τ , that is, then the initial perturbation u * 0 is called the conditional nonlinear optimal perturbation (CNOP) if and only if u * 0 is the solution of the following optimization problem: where B σ = {u 0 u 0 ∈ R n , u 0 ≤ σ } is the constraint domain on the initial perturbation.In terms of the L 2 norm, B σ is a ball with the center at the origin and the radius σ .In addition, J is called the objective function in the context of optimal control theory.

Estimation of the LBMPT
Duan and Luo (2010) designed a numerical method to calculate the LBMPT in their research of predictability (Fig. 1).
It should be noted that CNOPs stand for the initial uncertainty in the given constrained domain which leads to the largest prediction error.Therefore, the maximum prediction error can be estimated by solving CNOPs.
In detail, for a given first guess value T G of the prediction time, one can use a constraint nonlinear optimization algorithm to capture the CNOP so as to estimate the maximum prediction error E T G at time T G caused by the initial error in a constraint domain B σ .
If E T G > E m (E m stands for the allowable prediction error), we try to reduce the integral time step T G = T G − T , where T is a certain constant, and calculate the maximum prediction error at the reduced time.If E T G < E m , we will increase the integral time step T G = T G + T and calculate the maximum prediction error at the new T G .
If T G satisfies the conditions E T G + T > E m and E T G − T ≤ E m , then T G is considered the lower bound of maximum predictable time which satisfies the prediction precision E m under the given initial error.The operation flow chart is shown as Fig. 1.

The two-dimensional Ikeda model
The following two-dimensional Ikeda model is adopted as the prediction model: where 0 ≤ µ ≤ 1, a = 0.4 and b = 6.
From the expression of the model we find that the trigonometric functions appear in Eq. ( 12), and Eq. ( 13) is a frac- tion whose denominator includes two quadratic components.Thus, the two-dimensional Ikeda model has fairly strong nonlinearity (Li and Zheng, 2016).The solutions of the model present different behaviors with the change in the model parameter µ.When the parameter varies from 0 to 1, the numerical solutions change from a point attractor to periodic solutions, then to chaos, and end up with a limit cycle (Li and Zheng, 2016).The predictability problems are always launching under chaos.According to the conclusions given by Li and Zheng (2016), the model solution appeared chaotic when µ ∈ [0.700, 0.902].Figure 2 shows the numerical solution of the last 5000 steps in 10 000 integral steps, while the initial value is set to (x 0 , y 0 ) = (0.25, −0.325) and the model parameter is µ = 0.83.

The particle swarm optimization (PSO) algorithm
The PSO algorithm has recently got more and more attention (Eberhart and Shi, 2001;Banks et al., 2007Banks et al., , 2008;;Poli et al., 2007).The PSO algorithm was originally proposed by social psychologists Kennedy and Eberhart (1995).It simulates the collective behavior of birds foraging.Each particle represents a potential solution in the PSO algorithm, flies with specific velocity, respectively, and adjusts its trajectories according to the flying experiences of its own and the companion's, finally finding the optimal location in the solution space.To avoid the particles rapidly flying out of the solving region and to improve the ability of PSO to search for a global optimal solution, Shi and Eberhart (1998) introduced the maximum velocity and the inertia weight into PSOs, so as to restrain the particle's behaviors in the searching process.Clerc and Kennedy (2002) proposed a limiting factor for the two acceleration coefficients after they analyzed theoretically the convergence of PSO.
The basic PSO algorithm consists of three processes, namely, generating particles' positions and velocities, assessing particles, and updating particles' positions and velocities.
The mathematical description of the classic PSO algorithm is as follows: in an n-dimensional search space, each particle of PSO represents a potential solution of the optimization problem.We denote M the swarm size, the position and the velocity of the ith particle at the kth generation, respectively, P i (k) = (p i 1 (k), p i 2 (k), • • •, p i n (k)) the personal historical best position of the ith particle found so far, and ) the best position that the whole swarm attained so far; then the particle i's velocity and position in the next k + 1 generation can be updated according to the following formula: 1 and c 2 are the acceleration coefficients, which make particles having the ability to self-summarize and learn from excellent particles among the group to approach their own and group historical optimal points.r 1 and r 2 are two random numbers that are subject to uniform distribution on the interval [0, 1].ω is the inertia weight.It can be set as a fixed constant or a linear reduction function with the increase in the evolutional generations.The flow chart of the PSO algorithm is shown as Fig. 3.When using a PSO to search CNOPs for the estimation of the LBMPT, the prediction error at the specified forecast time is the associated objective function J .The initial perturbation δu, which is a two-dimensional vector in the search space in our situation, is the optimization variable.
4 Numerical experiments and their results analyses

The numerical experiments solving CNOPs by different optimization algorithms
In order to compare the performances of the ADJ-CNOP and PSO-CNOP in solving CNOPs, the CNOPs yielded by the filtering method are taken as the benchmark after finedividing the constraint domain of initial perturbations.The filtering method is implemented as follows.The corresponding circumscribed square of a constraint region of the CNOP is considered; four square meshes of a certain size are used to discretize the circumscribed square.For any mesh point outside the region, it is connected with the center of the region; the intersection point of this line with the boundary of the region is obtained.Integrating the Ikeda model from the initial basic state superimposed each of these intersection The particle swarm initialization scheme in PSO-CNOP is as follows: X i = (x i,1 (0), x i,2 (0)) are random vectors obeying a uniform distribution on B σ .
The first guess of the perturbation δu = (x 0 , y 0 ) for ADJ-CNOP is randomly picked from B σ .The constrained optimization algorithm used in ADJ-CNOP is the SPG2.
With the prediction time increasing, there would appear simultaneously many CNOPs for the two-dimensional Ikeda model because of the impact of the strong nonlinearity.Hence, different forecast times are adopted to test the ability of ADJ-CNOP and PSO-CNOP to attack the nonlinearity obstacles.For each forecast time, the numerical experiment using ADJ-CNOP or PSO-CNOP to obtain CNOPs is conducted 40 times, respectively; 40 CNOPs are clustered by the fuzzy c-means clustering (FCM) method and the accuracy of the CNOPs is statistically analyzed.The numerical experiment results show that when the integration time of the prediction model is short, the corresponding objective function (Eq.11) presents good behavior with the change in the initial perturbation.It has only two extreme values in the constraint ball of the initial perturbation.One is the global maximum, and the extreme point corresponds to the global CNOP.The other is the local maximum, and the extreme point is a local CNOP.Figure 4 shows the distribution of objective function values (OFVs) when the prediction times are on 6th unit time steps, i.e., 6 t (left) and 13 t (right), respectively, in which the global maximum point is located at point a and point b is the position of the local maximum.
Tables 1 and 2 demonstrate the analysis results of the CNOPs produced by ADJ-CNOP and PSO-CNOP when capturing CNOPs 40 times at the forecast times 6 t and 13 t, and the related results generated by the filtering method.Through the FCM method, we find that CNOPs obtained by the ADJ-CNOP method are divided into two categories: one is related to the global CNOP that accounted for 47.5 % (70 %) for the forecast time 6 t (13 t) of the total; the other is the local CNOP that makes up 52.5 % (30 %) of the total.However, 40 CNOPs captured by PSO-CNOP are completely the same, and they are coincident with the CNOP yielded by the filtering method.
From Tables 1 and 2, we can see that although the prediction time is short, the ADJ-CNOP method still has a large probability of capturing local CNOPs, while PSO-CNOP can always catch the global CNOP.Actually, we can draw the same conclusion with the prediction time being increased to 13 t.
When the forecast time increases to 14 t, the 40 CNOPs yielded by ADJ-CNOP and PSO-CNOP are demonstrated in the following Fig. 5.
Figure 6 indicates all CNOPs generated by the filtering method with fine division of the constraint domain of initial perturbations (the circumscribed squares of the constraint ball of CNOPs are separated into 1001×1001 small quadrate patches with very small side length 1.6402 × 10 −5 ).
Figure 6 shows that when the prediction time reaches 14 t, there exist many global CNOPs, and all of them are located in a line.Based on this, we find that ADJ-CNOP not only gets global or local CNOPs, but also captures false CNOPs at the prediction time 14 t, since no matter how small an area we take around one of these "CNOPs", there always exists one point whose objective function value is larger than the objective function value of the "CNOP".According to the definition of CNOPs, these "CNOPs" are not true CNOPs.Hence, we call them false CNOPs.Additionally, comparing the right panel of Fig. 5 with Fig. 6, it is easy to know that although many CNOPs are produced by PSO-CNOP in 40 repeated numerical experiments, all of the CNOP points are located on the same line as the one presented by the filtering method.Therefore, the CNOPs yielded by PSO-CNOP are all global.
When we keep extending the prediction time, the behavior of the objective function will get much worse, and more extreme points will appear.In order to verify the performance of the PSO-CONP method in solving CNOPs in the strong nonlinear case, the mean value and variance of the OFVs of the 40 CNOPs at different forecast times are calculated and compared with the maximal OFV (MOFV) obtained by the filtering method.
According to Table 3, the OFVs of 40 CNOPs calculated by the PSO-CNOP method at each forecast time are almost consistent with the maximum of the objective function gotten by the filtering method at the same forecast time.Therefore, PSO-CNOP is still capable of solving global CNOPs of the two-dimensional Ikeda model for long forecast times.Figure 7 demonstrates the distributions of OFVs at the prediction times 15 t, 18 t and 22 t, as well as the locations of all CNOPs generated by PSO-CNOP.
From the upper two panels of Fig. 7, we can see clearly that the CNOPs are all located in the maximal OFV region, which indicates that all of the CNOPs captured by the PSO-CNOP method are global.Because of the complexity of the  distributions of OFVs at the prediction time 22 t, it cannot be confirmed directly from the lower panel of Fig. 7 whether or not the CNOPs are global.Therefore, we select one CNOP randomly from each cluster of the 40 CNOPs and zoom into the graph nearby the CNOP point to look at the OFV distri-bution. Figure 8 gives one of the results, from which we can see that the CNOP is still located in the maximal OFV area.With further increasing of the prediction time, the strong nonlinearity deteriorates the behavior of the objective function seriously.In this situation, a predictability study based on the CNOP method becomes no longer meaningful because the CNOPs are too dispersive.

Comparison between PSO-CNOP and GA-CNOP
In the following, the GA is adopted to capture CNOPs of the two-dimensional Ikeda model, and the results are compared with the ones obtained by the PSO-CNOP.The method using the GA to compute CNOP is called GA-CNOP.The relevant operational flow chart of the GA is shown in Fig. 9.
The configuration of the genetic operators and the relevant parameter are the same as in Zheng et al. (2014).A more detailed description of the GA and numerical experiment scheme of GA-CNOP can be found in Zheng et al. (2014).
The performances of PSO-CNOP and GA-CNOP in solving the CNOPs are tested for different population sizes; the  results are statistically analyzed and presented in Tables 4  and 5.
It can clearly be seen from Tables 4 and 5 that the optimal population size of PSO-CNOP is about 15, but the optimal population size of GA-CNOP is about 40.Furthermore, the computational times and CFV calculated, respectively, by PSO-CNOP with a population size of 15 and GA-CNOP with a population size of 40, are compared with that of ADJ-CNOP.6 illustrates the mean CFV and the average computation time obtained by different methods in their 40 numerical experiments.From the numerical results we can see that the CFV of GA-CNOP is almost the same as PSO-CNOP; the GA is an effective optimal algorithm to obtain the optimal solution.However, the GA is more time consuming than PSO.At the same time, the computational time of PSO-CNOP and GA-CNOP is much greater than ADJ-CNOP, which is also the difference between stochastic searching algorithms and deterministic searching algorithms.Fortunately, in intelligent optimization algorithms, the parallel computation can be easily realized.The operators of different individuals in one generation are independent, and can be done in different CPUs, which thereby can take full advantage of fast developed parallel computation technology (Fang et al., 2009).

Estimation of the LBMPT
The CNOP method can be adopted to do a similar study on other two predictability sub-problems.Here we only focus on the estimation of the LBMPT to discuss the influence of nonlinearity.To demonstrate the effectiveness of PSO-CNOP in solving this problem, the filtering method, ADJ-CNOP and PSO-CNOP are used in the numerical experiments, respectively.In order to compare the estimation accuracy of the LBMPT generated by ADJ-CNOP and PSO-CNOP, the LBMPT computed by the filtering method with a fine division is taken as the benchmark.Different allowable prediction errors E m , 0.5, 0.8, 1.1, 1.4, and 1.7, and a different constrained radius δ of initial perturbations, 0.01, 0.02, 0.03, 0.04, and 0.05, are employed to verify the performance of the ADJ-CNOP and PSO-CNOP methods.The lattice spacing of the filtering method is 0.001.Tables 7, 8 and 9 illustrate the LBMPTs computed by the filtering method, PSO-CNOP and ADJ-CNOP, respectively.
Comparing Table 7 with Table 8, we find that the LBMPTs estimated by the PSO-CNOP approach are completely the same as the ones computed by the filtering method.It is www.nonlin-processes-geophys.net/24/101/2017/ Nonlin.Processes Geophys., 24, 101-112, 2017 The bold numbers in Table 9 are the LBMPTs that are different from the ones computed by the filtering method.The LBMPTs yielded by the ADJ-CNOP method are generally larger.The reason is that the CNOP given by the ADJ-CNOP method is often local, even false.Therefore the estimation of the maximum prediction error based on the CNOP is usually questionable and untrusted.
To investigate the probability that the ADJ-CNOP method generates incorrect LBMPTs, we operate the numerical experiment shown in Table 9 40 times with various first guesses of the initial perturbations.The statistical analysis results are given in Table 10.
From Table 10, we can see that the ratio of incorrect LBMPTs based on ADJ-CNOP is high.The highest one even reaches 95 % when δ = 0.02 and E m = 1.7.This problem is serious for real weather forecasts since it can mislead forecasts with a large probability.Furthermore, precision problems of using ADJ-CNOP to estimate the LBMPT are investigated.Results show that when the objective function has multiple extreme values, ADJ-CNOP has a large probability of producing the local CNOP, hence inducing false estimation of the LBMPT.As PSO-CNOP can always yield global CNOPs, therefore, the estimation of the LBMPT presented by PSO-CNOP is precise.It is consistent with the one yielded by the filtering method with fine division.
As we know, our numerical experiments focus on twodimensional prediction models only.When considering highdimensional and more complex models, whether or not the classic PSO algorithm used in this study can overcome the influence of high dimensions and the computation time meet the real requirement is still unknown.The problems of the curse of dimensionality and multimodal function are a big challenge for almost all intelligent optimization algorithms, also PSO.Whether it can be effective in higher-dimensional and more complicated models deserves further research.In short, the PSO-CNOP approach is an alternative method to study predictability problems in the case of strong nonlinearity.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Flow chart of solving the LBMPT.

Figure 2 .
Figure 2. The distribution of solutions at the last 5000 steps when µ = 0.83.

Figure 3 .
Figure 3.The flow chart of the PSO algorithm.

Figure 4 .
Figure 4.The distributions of OFVs at the prediction times 6 t (left) and 13 t (right), in which dots a and b are the global and local maximum points.

Figure 5 .
Figure 5.The distributions of OFVs at the prediction time 14 t, where the dots in the left (right) panel are the CNOPs produced by ADJ-CNOP (PSO-CNOP).

Figure 6 .
Figure 6.The distributions of OFVs at the prediction time 14 t, where dots are the CNOPs produced by the filtering method.
Figure 7.The distribution of OFVs at the prediction times 15 t (the left of the upper panel), 18 t (the right of the upper panel) and 22 t (the lower panel), where dots denote CNOPs captured by PSO-CNOP.

Figure 8 .
Figure 8. Local distribution of the OFV at the prediction time 22 t nearby the CNOP point.

Figure 9 .
Figure 9.The flow chart of the GA.

Table 1 .
Statistical analysis of CNOPs produced by different methods at 6 t.

Table 2 .
Same as Table 1, except at 13 t.

Table 3 .
The precision analysis of the CNOPs produced by PSO-CNOP at different forecast times.

Table 4 .
Mean CFV of 40 CNOPs produced by PSO-CNOP for every given population size.

Table 5 .
Same as Table 4, except that 40 CNOPs are produced by GA-CNOP.

Table 6 .
Analysis of CFV and average computational time by different methods.

Table 7 .
The LBMPT computed by the filtering method.

Table 8 .
The LBMPT estimated by the PSO-CNOP method.

Table 10 .
Incorrect ratio in all 40 LBMPTs yielded by ADJ-CNOP.E m δ = 0.01 δ = 0.02 δ = 0.03 δ = 0.04 δ = 0.05 Since the two-dimensional Ikeda model has strong nonlinearity, when we utilize the ADJ-CNOP method to capture CNOPs, not only global or local CNOPs, but even false CNOPs, are obtained.The reason for this is that in the case of strong nonlinearity, the gradient provided by the adjoint model is incorrect.When the traditional optimization algorithm uses a wrong descent direction to search for extreme values of the objective function, false CNOPs are presented.PSO is a heuristic search algorithm based on population.It can overcome the nonlinear influences and produce global CNOPs with high probability.In addition, the operation of the PSO algorithm is simple.This study applies the PSO algorithm to capture the CNOP of the two-dimensional Ikeda model.Numerical experiment results with different forecast times demonstrate that although the objective function has awful behavior and multiple extreme values, PSO-CNOP can still capture global CNOPs.