A Modified Hybrid Genetic Algorithm for Solving Nonlinear Optimal Control Problems

Here, a two-phase algorithm is proposed for solving bounded continuous-time nonlinear optimal control problems (NOCP). In each phase of the algorithm, a modified hybrid genetic algorithm (MHGA) is applied, which performs a local search on offsprings. In first phase, a random initial population of control input values in time nodes is constructed. Next, MHGA starts with this population. After phase 1, to achieve more accurate solutions, the number of time nodes is increased. The values of the associated new control inputs are estimated by Linear interpolation (LI) or Spline interpolation (SI), using the curves obtained from the phase 1. In addition, to maintain the diversity in the population, some additional individuals are added randomly. Next, in the second phase, MHGA restarts with the new population constructed by above procedure and tries to improve the obtained solutions at the end of phase 1. We implement our proposed algorithm on 20 well-known benchmark and real world problems; then the results are compared with some recently proposed algorithms. Moreover, two statistical approaches are considered for the comparison of the LI and SI methods and investigation of sensitivity analysis for the MHGA parameters.


Introduction
NOCPs are dynamic optimization problems with many applications in industrial processes such as airplane, robotic arm, bio-process system, biomedicine, electric power systems, and plasma physics [1].
High-quality solutions and the less required computational time are main issues for solving NOCPs.The numerical methods, direct [2] or indirect [3], usually have two main deficiencies, less accuracy and convergence to a local solution.In direct methods, the quality of solution depends on discretization resolution.Since these methods, using control parametrization, convert the continuous problem to discrete problem, they have less accuracy.However, the adaptive strategies [4,5] can overcome these defects, but they may be trapped by a local optimal, yet.In indirect approaches, the problem, through the use of the Pontryagins minimum principle (PMP), is converted into a two-boundary value problem (TBVP) that can be solved by numerical methods such as shooting method [6].These methods need the good initial guesses that lie within the domain of convergence.Therefore, the numerical methods, usually are not suitable for solving NOCPs, especially for large-scale and multimodal models.
Metaheuristics as the global optimization methods can overcome these problems, but they usually need more computational time, though they do not really need good initial guesses and deterministic rules.Several researchers used metaheuristics to solve optimal control problems.For instance, Michalewicz et al. [7] applied floating-point Genetic algorithms (GA) to solve discrete time optimal control problems; Yamashita and Shima [8] used the classical GAs to solve the free final time optimal control problems with terminal constraints.Abo-Hammour et al. [9] used continuous GA for solving NOCPs.Moreover, the other usages of GA for optimal control problems can be found in [6,10].Lopez Cruz et al. [11] applied differential evolution (DE) algorithms for solving the multimodal optimal control problems.Recently, Ghosh
Cost function (1) must be minimized subject to dynamic (2), control and state equality constraints (3) and control and state inequality constraints (4), the final state constraints (5), and the initial condition (6).A special case of the NOCPs is the linear quadratic regulator (LQR) problem where the dynamic equations are linear and the objective function is a quadratic function of  and .The minimum time problems, tracking problem, terminal control problem and minimum energy are another special case of NOCPs.

Modified Hybrid Genetic Algorithm
In this section, first MHGA, as a subprocedure for the main algorithm, is introduced.To perform MHGA, the control variables are discretized.Next, NOCP is changed into a finite dimensional NLP; see [21,26].Now, we can imply a GA to find the global solution of the corresponding NLP.In the following, we introduce GA operators.
3.1.Underlying GA.GAs, introduced by Holland in 1975, are heuristics and probabilistic methods [27].These algorithms start with an initial population of solutions.This population is evaluated by using genetic operators that include selection, crossover, and mutation.Here, in MHGA, the underlying GA has the following steps.
Initialization.The time interval is divided into   −1 subintervals using time nodes  0 , . . .,    −1 and then the control input values are computed (or selected randomly).This can be done by the following stages.
(2) The corresponding control input value at each time node   is an  × 1 vector,   , which can be calculated randomly, with the following components.

𝑢
where   is a random number in [0, 1] with a uniform distribution and  left ,  right ∈ R  are the lower and the upper bound vectors of control input values, which can be given by the problem's definition or the user (e.g., see the NOCPs numbers (6) and (5) in the Appendix, resp.).So, each individual of the population is an  ×   matrix as  = (  ) . Next, we let  () = (  ) ×  ,  = 1, 2, . . .,   as th individual of the population, which   is the size of the population.
Evaluation.For each control input matrix,  () ,  = 1, 2, . . .,   , the corresponding state variable is an  ×   matrix,  () , and it can be computed by the forth Runge-Kutta method on dynamic system (2) with the initial condition (6), approximately.Then, the performance index, ( () ), is approximated by a numerical method (denoted by J).If NOCP includes equality or inequality constraints (3) or ( 4), then we add some penalty terms to the corresponding fitness value of the solution.Finally, we assign ( () ) to  () as the fitness value as follows: where Selection.To select two parents, we use a tournament selection [27].It can be applied for parallel GA.The tournament operator applies competition among the same individuals, and the best of them is selected for next generation.At first we select a specified number of individuals from population, randomly.This number is tournament selection parameter, which is denoted by  tour .The tournament algorithm is given in Algorithm 1.
Crossover.When two parents  (1) and  (2) are selected, we use the following stages to construct an offspring.
Mutation.We apply a perturbation on each component of the offspring as follows: where where έ is the machine epsilon.
Termination Conditions.Underlying GA is terminated when at least one of the following conditions is occurred.
(1) The maximum number of generations,   , is reached.
(2) Over a specified number of generations,   , we do not have any improvement (the best individual is not changed) or the two-norm, or error, of final state constraints will reach a small number as the desired precision, ; that is, where  = [ 1 ,  2 , . . .,    ]  is the vector of final state constraints, defined in (5).

SQP.
SQP algorithm [22,28] is performed on the NLP ( 14)-( 15), using  0 = of and constructed in (11), as the initial solution when the maximum number of iteration is .
SQP, is an effective and iterative algorithm for the numerical solution of the constrained NLP problem.This technique is based on finding a solution to the system of nonlinear equations that arise from the first-order necessary conditions for an extremum of the NLP problem.Using an initial solution of NLP,   ,  = 0, 1, . .., a sequence of solutions as  +1 =   +   is constructed, which   is the optimal solution of the constructed quadratic programming (QP) that approximates NLP in the iteration  based on   , as the search direction in the line search procedure.For the NLP (15), the principal idea is the formulation of a QP subproblem based on a quadratic approximation of the Lagrangian function as (, ) = () +   ℎ(), where the vector  is Lagrangian multiplier and ℎ() return the vector of, inequality constraints evaluated at .The QP is obtained by linearizing the nonlinear functions as follows: Similar to [26], here a finite difference approximation is applied to compute the gradient of the cost function and the constraints, with the following components: where  is the double precision of machine.So, the gradient vector is ∇ = [/ 0 , . . ., /   −1 ]  .Also, at each major iteration a positive definite quasi-Newton approximation of the Hessian of the Lagrangian function, , is calculated using the BFGS method [28], where   ,  = 1, . . ., , is an estimate of the Lagrange multipliers.The general procedure for SQP, for NLP ( 14)- (15), is as follows.
(2) Construct the QP subproblem ( 16)- (17), based on  0 , using the approximations of the gradient and the Hessian of the the Lagrangian function.
(3) Compute the new point as  +1 =   +   , where   is the optimal solution of the current QP.(4) Let  =  + 1 and go to step (2).
Here, in MHGA, SQP is used as the local search, and we use the maximum number of iterations as the main criterion for stopping SQP.In other words, we terminate SQP when it converges either to local solution or the maximum number of SQP's iterations is reached.

MHGA.
In MHGA, GA uses a local search method to improve solutions.Here, we use SQP as a local search.Using SQP as a local search in the hybrid metaheuristic is common, for example, see [20].
MHGA can be seen as a multi start local search where initial solutions are constructed by GA.From another perspective, MHGA can be seen as a GA that the quality of its population is intensified by SQP.In the beginning of MHGA, a less number of iterations for SQP was used.Then, when the promising regions of search space were found by GA operators, we increase the number of iterations of SQP gradually.Using this approach, we may decrease the needed running time (in [19] the philosophy of this approach is discussed).
Finally, we give our modified MHGA, to find the global solution, by the flowchart in Figure 1.

Proposed Algorithm
Here, we give a new algorithm, which is a direct approach, based on MHGA, for solving NOCPs.The proposed algorithm has two main phases.In the first phase, we perform MHGA with a completely random initial population constructed by (7).In the first phase, to find the promising regions of the search space, in a less running time, we use a few numbers of time nodes.In addition, to have a faster MHGA, the size of the population in the first phase is usually less than the size of the population in the second phase.
After phase 1, to maintain the property of individuals in the last population of phase 1 and to increase the accurately of solutions, we add some additional time nodes.When the number of time nodes is increased, it is estimated that the quality of solution obtained by numerical methods (e.g., Runge-Kutta and Simpson) is increased.Thus, we increase time nodes from   1 in phase 1 to   2 in phase 2. The corresponding control input values of the new time nodes are added to individuals.To use the information of the obtained solutions from phase 1 in the construction of the initial population of phase 2, we use either Linear or Spline interpolation to estimate the value of the control inputs in the new time nodes in each individual of the last population of phase 1.Moreover, to maintain the diversity in the initial population of phase 2, we add new random individuals to the population using (7).In the second phase, MHGA starts with this population and new value of parameters.Finally, the proposed algorithm is given in Algorithm 2.

Numerical Experiments
In this section, to investigate the efficiency of the proposed algorithm, 20 well-known and real world NOCPs, as benchmark problems, are considered which are presented in terms of ( 1)-( 6) in the Appendix.These NOCPs are selected with single control signal and multi control signals, which will be demonstrated in a general manner.
The numerical behaviour of algorithms can be studied from two view of points, the relative error of the performance index and the status of the final state constraints.Let  be the obtained performance index by an algorithm,   , defined in {Initialization} Input the desired precision, , in (13), the penalty parameters,   ,  = 1, 2, 3, in (8), the bounds of control input values in (7) Algorithm 2: The proposed algorithm.(13), be the error of final state constraints, and  * be the best obtained solution among all implementations, or the exact solution (when exists).Now the relative error of ,   , of the algorithm can be defined as To more accurate study, we now define a new criterion, called factor, to compare the algorithms as follows: Note that   shows the summation of two important errors.Thus, based on   we can study the behaviour of the algorithms on the quality and feasibility of given solutions, simultaneously.
To solve any NOCP described in the Appendix, we must know the algorithm's parameters including: MHGA's parameters; including   ,   ,   ,   ,   and , in both phases in Algorithm 2, and the problem's parameters including , in (13),   ,  = 1, 2, 3, in (8),  left and  right , in (7).To estimate the best value of the algorithm's parameters, we ran the proposed algorithm with different values of parameters and then select the best.However, the sensitivity of MHGA parameters are studied in the next section.In both of phases of Algorithm 2, in MHGA, we let   = 4 and    = 0.8, for  = 1, 2. Also, we consider The other MHGA parameters are given in the associated subsection and the problem's parameters in Table 2.For each NOCP, 12 different runs were done and the best results are reported in Table 1, which the best value of each column is seen in the bold.
The reported numerical results of the proposed algorithm, for each NOCP, include the value of performance index, , the relative error of ,   , defined in (19), the required computational time, Time, the norm of final state constraints,   , defined in (13) and the factor,   , defined in (20).
The algorithm was implemented in Matlab R2011a environment on a Notebook with Windows 7 Ultimate, CPU 2.53 GHz and 4.00 GB RAM.Also, to implement SQP in our proposed algorithm, we used "fmincon" in Matlab when the "Algorithm" was set to "SQP".Moreover, we use composite Simpson's method [29] to approximate integrations.
Remark 1.We use the following abbreviations to show the used interpolation method in our proposed algorithm: (1) LI: linear interpolation.
For comparing the numerical results of the proposed algorithm two subsections are considered, comparison with some metaheuristic algorithms, in Section 5.1, and comparison with some numerical methods, in Section 5.2.We give more details of these comparisons in the following subsections.

Comparison with Metaheuristic
Algorithms.The numerical results for the NOCPs numbers (1)-(3), in the Appendix, are compared with a continuous GA, CGA, as a metaheuristic, proposed in [9], which gave better solutions than shooting method and gradient algorithm, from the indirect methods category [2,30], and SUMT from the direct methods category [26].For NOCPs numbers (4) and (5), the results are compared with another metaheuristic, which is a hybrid improved PSO, called IPSO, proposed in [20].[9].The first NOCP in the Appendix is Van  the relative error of  and the factor of the proposed algorithm is better, and so the proposed algorithm is more robust than CGA.[9].  1.The values of ,   ,   and   for LI and SI methods, separately, are less than CGA.Therefore, the proposed algorithm can achieve much better quality solutions than the CGA, with reasonable computational time.[20].For the forth NOCP in the Appendix, which is a Mathematical System with Nonlinear Inequality Constraint, NSNIC, the numerical results are compared with IPSO.MSNIC contains an inequality constraint, (, ) =  2 () + 0.5 − 8( − 0.5) 2 ≤ 0. The problem solved by several numerical methods as [24,31].From [20], IPSO method could achieved more accurate results than mentioned numerical methods.Also, MSNIC can be solved by the proposed algorithm, with the MHGA's parameters as   1 = 31,   2 = 91,   = 100 and   = 60.From the forth row of Table 1, the absolute error of ,   , for LI and SI methods equal 0 and 5.89 × 10 −4 , respectively, which are less than IPSO's, 0.0165.Subplots (a) and (b) in Figure 2, show the graphs of the convergence rate for the performance index and the inequality constraint, respect to the number of iteration, respectively.[20].The fifth NOCP in the Appendix is a model of a nonlinear Continuous Stirred-tank Chemical Reactor, CSTCR.It has two state variables  1 () and  2 (), as the deviation from the steady-state temperature and concentration, and one control variable (), which represent the effect of the flow rate of cooling fluid on chemical reactor.The objective is to maintain the temperature and concentration close to steady-state values without expending large amount of control effort.Also, this is a benchmark problem in the handbook of test problems in local and global optimization [32], which is a multimodal optimal control problem [33].It involves two different local minima.The values of the performance indexes, for these solutions, equal 0.244 and 0.133.Similarly to the MSNIC, the numerical results of the proposed algorithm, with the MHGA's parameters as   1 = 31,   2 = 51,   = 100, and   = 50, are compared with IPSO.From Table 1, the performance index, , for LI and SI methods is equal to  * = 0.1120 which is less than IPSOs 0.1355.

Comparison with Numerical Methods.
For NOCPs numbers (6)-( 20), the comparison are done with some numerical methods.Unfortunately, for these methods, usually, the final state constraints and the required computational time are not reported, which are shown with NR in Table 1, but these values are reported for both LI and SI methods in Table 1.For all NOCPs, in this section, the MHGA's parameters are considered as   1 = 31,   2 = 51,   = 100,   = 50, with the problem parameters in Table 2. [34].The NOCP number (6), in the Appendix, has exact solution, which has an inequality constraint as ( 1 , ) = −6 −  1 () ≤ 0. The exact value of performance index equals  * = −5.5285[35].This problem has been solved by a numerical method proposed in [34], called Bézier.From sixth row of Table 1, the absolute error of the LI and SI methods equal 0.0177, which is less than Bézier's, 0.0251.

Compared with Bézier
Figure 3, shows the graphs of the convergence rate of the performance index, subplot (a), and inequality constraint, subplot (b), respect to the number of iteration, using SI method.[36].For NOCP number (7) in the Appendix, which is a constraint nonlinear model, the numerical results of the proposed algorithm are compared with HPM, proposed in [36].This NOCP has a final state constraint as  =  − 0.5 = 0. From [36], the norm of final state constraint for HPM equals 4.2 × 10 −6 , however, this criterion for the LI and SI methods equals 2.35 × 10 −9 and  *  = 2.82 × 10 −10 , respectively.From Table 1, it is obvious that the obtained values of the performance index, the norm of final state constraint and   for the SI method are more accurate than LI than HPM methods.[26].For the NOCPs numbers (8)- (20), in the Appendix, the numerical results of LI and SI methods are compared with two numerical methods contain: SQP and SUMT, proposed in [26].Among these NOCPs, only two problems numbers (8) and (18) are unconstrained, and the others have at least one constraint, final state constraint or inequality constraint.All of these NOCPs solved by the proposed algorithm, with the problem parameters in Table 2, and their results are summarized in Table 1.Because the final state constraints, in these methods, are not reported, we let   = 0 to calculate the factor   .Figures 4-16 show the graphical results of NOCPs numbers (8)-( 20) in the Appendix, using SI method.For unconstrained NOCPs, numbers (8) and (18), only the graph of convergent rate of performance index is shown; see Figures 4 and 14.For constraint NOCPs and NOCPs numbers (11)-( 17) and numbers ( 19)- (20), with final state constraint, the graphs of convergent rates of performance index and the error of final state constraint are shown; see  For the constraint NOCPs with inequality constraints, NOCPs numbers (9) and (10), the graphs of convergent rate of performance index and inequality constraint are shown; see Figures 5 and 6.

Compared with SQP and SUMT
Table 1 shows that the proposed algorithm, LI and SI methods, was 100 percent successful in point of views the performance index, , and the factor,   , numerically.So, the proposed algorithm provides robust solutions with respect to the other mentioned, numerical or metaheuristic, methods.To compare  in LI and SI methods, in 35 percent of NOCPs, LI is more accurate than SI, in 10 percent SI is more accurate than LI and in 55 percent are same.In point of view the final conditions, 65 percent of NOCPs in the Appendix have final state constraints.In all of them   in the proposed algorithm is improved, except CRP problem, however the factor of the proposed algorithm is the best, yet.In 54 percent of these NOCPs, LI is better performance and in 46 percent SI is better.Also, in point of the running time, Time, LI is same as SI, that is, in 50 percent of NOCPs LI and in 50 percent SI is better than another.So, the proposed algorithm could provide very suitable solutions, in a reasonable computational time.Also, for more accurate comparison of LI and SI methods a statistical approach will done in next section.
From Table 1, to compare with numerical methods, SQP and SUMT, in NOCPs (8)- (20), the mean of   for LI, SI, SQP and SUMT equals 0.0145, 0.0209, 1.69 × 10 8 and 1.36 × 10 8 , respectively.Also, the mean for the error of final state constraints, for these NOCPs, equal 3.34 × 10 −4 , 4.41 × 10 −5 , 0 and 0, respectively.For   , these values are 0.0213, 0.0014, 1.2263 and 1.1644.Therefore, the performance index, , and the factor,   , for the LI and SI methods are more accurate than SQP and SUMT.So the proposed algorithm gave more better solution in comparison with the numerical methods.Therefore, based on this numerical study, we can conclude that the proposed algorithm outperform well-known numerical method.Since, the algorithms were not implemented on the same PC, the computational times of them are not competitive.Therefore, we did not give the computational times in bold in Table 1.

Sensitivity Analysis and Comparing LI and SI
In this section, two statistical analysis, based on the oneway analysis of variance (ANOVA), used for investigating the sensitivity of MHGA parameters, and Mann-Whitney,  applied comparing LI and SI methods, are done, by the statistical software IBM SPSS version 21.

Sensitivity Analysis.
In order to survey the sensitivity of the MHGA parameters from the proposed algorithm, the VDP problem is selected, for example.The influence of these parameters are investigated for this NOCP on the dependent outputs consist of the performance index, , the relative error of ,   , the required computational time, Time, the error of final state constraints,   and the factor,   .The independent parameters are consist of the number of time nodes in both two phases,   1 and   2 , the size of population in both two phases,   1 and   2 , the maximum number of generations without improvement,   ,  and the maximum number of generations,   .Because the mutation implementation probability,   , has a less influence in numerical results, it is not considered.Also  is changed in each iteration of the proposed algorithm, so it is not considered too.
At first, we selected at least four constant value for each of parameters and then in each pose, 35 different runs were made, independently.The statistical analysis is done based on ANOVA.The descriptive statistics, which contains the number of different runs (), the mean of each   output in  different runs, (Mean), standard deviation, (S.d), maximum, (Max) and minimum, (Min), of each output, could be achieved.Among of the MHGA parameters, we present the descriptive statistics of the ANOVA, only for the   1 parameter, which is reported in Table 3. From above cases, it is obvious that all parameters can be effect on the required computational time, except   .Moreover, intuitions shows the norm of final state constraints,   is independent output with respect to all parameters, that is, any of the parameters could not effect on this output and it is not sensitive with respect to any of the MHGA parameters.

Comparison of LI and SI.
To compare the efficiency of the LI and SI methods, for NOCPs in the Appendix, we used the Mann-Whitney nonparametric statistical test [37].Using the nonparametric statistical tests for comparing the performance of several algorithms presented in [38].For this purpose, At first, 15 different runs, for each of NOCP were made.6-9.From Table 6, since  = 0 >  0.01 = −2.34,then the average relative error of the LI and SI methods are same with the significant level 0.01.In other words, with probability 99%, we can say that LI and SI methods have same effectiveness form the perspective of error for .Similarly, the results of Tables 7-9 indicate that LI

The Impact of SQP
In this section, we investigation the impact of the SQP, on the proposed algorithm.For this purpose, we remove the SQP from Algorithm 2, as without SQP algorithm.For comparing the numerical results with the previous results, the required running time in each NOCP is considered fixed, which is the maximum of running time in 12 different runs, were done in Section 5. Also, the all of the parameters for each problem was set as same as Section 5.The numerical results of the without SQP algorithm is summarized in Table 10.By comparing the results of Tables 10 and 1, it is obvious that, for all NOCPs in the Appendix, the obtained values of the performance index and the norm of final state constraint for the proposed algorithm (Algorithm 2), are more accurate than the without SQP algorithm.

Conclusions
Here, a two-phase algorithm was proposed for solving bounded continuous-time NOCPs.In each phase of the algorithm, a MHGA was applied, which performed a local search on offsprings.In first phase, a random initial population of control input values in time nodes was constructed.
Next, MHGA started with this population.After phase 1, to achieve more accurate solutions, the number of time nodes was increased.The values of the associated new control inputs were estimated by Linear interpolation (LI) or Spline interpolation (SI), using the curves obtained from phase 1.In addition, to maintain the diversity in the population, some additional individuals were added randomly.Next, in the second phase, MHGA restarted with the new population constructed by above procedure and tried to improve the obtained solutions at the end of phase 1.We implemented our proposed algorithm on 20 well-known benchmark and real world problems; then the results were compared with some recently proposed algorithms.Moreover, two statistical approaches were considered for the comparison of the LI and SI methods and investigation of sensitivity analysis for the MHGA parameters.

Figure 2 :
Figure 2: Graphical results of MSNIC problem using SI method: (a) convergence rate of the performance index and (b) the inequality constraint, (, ) ≤ 0, in respect to the number of iterations.

Figure 3 :
Figure 3: Graphical results for NOCP number (6) using SI method: (a) convergence rate for the performance index and (b) the inequality constraint, ( 1 , ) ≤ 0, in respect to the number of iterations.

Figure 4 :
Figure 4: Convergence rate for the performance index of SI method for NOCP number (8).

Figure 5 :
Figure 5: (a) Convergence rate for the performance index and (b) the inequality constraint of SI method, for NOCP number (9).

Figure 6 :
Figure 6: (a) Convergence rate of the performance index and (b) the inequality constraint for SI method, for NOCP number (10).

Figure 7 :
Figure 7: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (11).

Figure 8 :
Figure 8: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (12).

Figure 9 :
Figure 9: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (13).

Figure 10 :
Figure 10: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (14).

Figure 11 :
Figure 11: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (15).

Figure 12 :
Figure 12: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (16).

Figure 13 :
Figure 13: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (17).

Figure 14 :
Figure 14: Convergence rate for the performance index of SI method for NOCP numbers (18).

Figure 15 :
Figure 15: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (19).

Figure 16 :
Figure 16: (a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number (20).
{Initialization} Input the the size of tournament,  tour and select  tour individuals from population, randomly.Let   ,  = 1, 2, . . .,  tour be the indexes of them.repeat if  tour = 2 then if (  1 ) < (  2 ) then Let  =  1 .Perform Tournament algorithm with size  tour /2 with output  1 .Perform Tournament algorithm with size  tour /2 with output  2 .until  1 ̸ =  2 Return the th individual of population.
,  left and  right .{Phase 1} Perform MHGA with a random population and   1 ,   1 ,   1 ,   1 ,   1 and  1 .{Construction of the initial population of the phase 2} Increase time nodes uniformly to   2 and estimate the corresponding control input values of the new time nodes in each individual obtained from phase 1, using either Linear or Spline interpolation.Create   2 −   1 new different individuals with   2 time nodes, randomly.{Phase 2} Perform MHGA with the constructed population and   2 ,   2 ,   2 ,   2 ,   2 and  2 .

Table 1 :
Der Pol Problem, VDP, which has two state variables and one control variable.VDP problem has a final state constraint, which is  =  1 (  ) −  2 (  ) + 1 = 0.The results of the proposed algorithm with the MHGA's parameters as   1 = 31,   2 = 71,   = 300 and   = 200, are reported in Table 1.From Table 1, it is obvious that the numerical results of LI and SI methods are more accurate than CGA, with less amount of   .The second NOCP in the Appendix is Chemical Reactor Problem, CRP, which has two state variables and one control variable.The results of the proposed algorithm, with the MHGA's parameters as   1 = 31,   2 = 71,   = 300 and   = 200, are shown in the second row of Table 1.CRP problem has two final state constraints,  = [ 1 ,  2 ]  .Although, from Table 1, the norm of final state constraints,   , for the CGA, equals  *  = 7.57 × 10 −10 , is less than   's of LI and SI methods, which equals 1.15 × 10 −9 and 5.99 × 10 −9 , respectively, but the performance index, The best numerical results in 12 different runs of the SI and LI methods for NOCPs in the Appendix.
a Not reported.
The third NOCP in the Appendix is Free Floating Robot Problem, FFRP, which has six state variables and four control variables.FFRP has been solved by CGA and the proposed algorithm, with the MHGA's parameters as   1 = 31,   2 = 71,   = 300 and   = 200.This problem has six final state constraints,  = [ 1 − 4,  2 ,  3 − 4,  4 ,  5 ,  6 ]  .The numerical results are shown in Table

Table 2 :
The problem parameters for NOCPs in the Appendix.

Table 4 ,
summarized the statistical data, contain the test statistics () and -values, of ANOVA tests.Sensitivity analysis for each parameter, separately, are done based on this table, as follows.  1 : From Table 4, the significant level, or -value, for   is equal to 0.279, which is greater than 0.05.So, from ANOVA,   1 parameter has no significant effect on the factor   .With similar analysis, this parameter has no effect on the other outputs, except required computational time, Time, because its value is equal to 0, which is less than 0.05, that is, this output sensitive to the parameter   1 .  2 : From the second row of Table 4, the outputs ,   ,   and Time are sensitive to the parameter   2 , and the only output   is independent.  1 : From the third row of Table 4, only the computational required time, Time, is sensitive, because the its value is equal to 0.006, which is less than 0.05, and the other outputs, contain ,   ,   ,   , are independent.  2 : From the forth row of Table 4, the output   is independent respect to the parameter   2 , but other outputs, contain ,   , Time,   are sensitive to this parameter.  : From the fifth row of Table 4, the outputs ,   and   are independent to the parameter   , and other outputs, contain   and Time are sensitive respect to this parameter.  : The sensitivity analysis is similar to   1 .

Table 5 ,
shows the mean of numerical results, separately.Here, we apply the fix MHGA parameters as   1 = 9,   2 = 12,   1 = 31,   2 = 91,   = 100 and   = 50 with the previous problem parameters in Table 2. Comparison criterions contain: ,   ,   and Time.The results of Mann-Whitney test are shown in Tables

Table 3 :
Descriptive statistics for the sensitivity analysis of   1 , for VDP problem.

Table 5 :
The mean of numerical results of 15 different runs for NOCPs in the Appendix, using LI and SI methods.

Table 6 :
Results of Mann-Whitney test on relative errors of the pair (LI, SI) for .

Table 7 :
Results of Mann-Whitney test on relative errors of the pair (LI, SI) for   .

Table 8 :
Results of Mann-Whitney test on relative errors of the pair (LI, SI) for Time.

Table 9 :
Results of Mann-Whitney test on relative errors of the pair (LI, SI) for   .

Table 10 :
The numerical results of the without SQP algorithm, for NOCPs in the Appendix.