# A Modified Hybrid Genetic Algorithm for Solving Nonlinear Optimal Control Problems

Authors: Saeed Nezhadhosein; Aghileh Heydari; Reza Ghanbari
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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 et al. [12] developed an ecologically inspired optimization technique, called invasive weed optimization (IWO), for solving optimal control problems. The other well-known metaheuristic algorithms used for solving NOCPs are genetic programming (GP) [13], particle swarm optimization (PSO) [14, 15], ant colony optimization (ACO) [16], and DE [17, 18].To increase the quality of solutions and decrease the running time, hybrid methods were introduced, which used a local search in the implementation of a population-based metaheuristics, [19]. Modares and Naghibi-Sistani [20] proposed a hybrid algorithm by integrating an improved PSO with successive quadratic programming (SQP) for solving NOCPs. Recently, Sun et al. [21] proposed a hybrid improved GA, which used simplex method (SM) to perform a local search, for solving NOCPs and applied it for chemical processes.Based on the success of the hybrid methods for solving NOCPs, mentioned above, we here use a modified hybrid genetic algorithm (MHGA), which combines GA with SQP, see [22], as a local search. SQP is an iterative algorithm for solving nonlinear programming (NLP) problems, which uses gradient information. It can moreover be used for solving NOCPs, see [23–25]. For decreasing the running time in the early generations (iterations) of MHGA, a less number of iterations for SQP was used and then, when the promising region of search space was found, we increase the number of iterations of SQP, gradually.To perform MHGA for solving an NOCP, the time interval is uniformly divided by using a constant number of time nodes. Next, in each of these time nodes, the control variable is approximated by a scaler vector of control input values. Thus, an infinite dimensional NOCP is changed to a finite dimensional NLP. Now, we encounter two conflict situations: the quality of the global solution and the needed computational time. In other words, when the number of time nodes is increased, then we expect that the quality of the global solution is also increased, but we know that in this situation the computational time is increased dramatically. In other situation, if we consider less number of time nodes, then the computational time is decreased but we may find a poor local solution. To conquer these problems, MHGA performs in two phases. In the first phase (exploration phase), to decrease the computational time and to find a promising region of search space, MHGA uses a less number of time nodes. After phase 1, to increase the quality of solutions obtained from phase 1, the number of time nodes is increased. Using the population obtained in phase 1, the values of the new control inputs are estimated by Linear or Spline interpolations. Next, in the second phase (exploitation phase), MHGA uses the solutions constructed by the above procedure, as an initial population.The paper is organized as follows: in Section2, the formulation of the problem is introduced. In Section 3, the proposed MHGA is presented. In Section 4, we introduce our algorithm for solving NOCP. In Section 5, we provide 20 numerical benchmark examples, to compare the proposed algorithm with the other recently proposed algorithms. In Section 6, we consider two statistical approaches for the comparison of the LI and SI methods and investigation of sensitivity analysis of the algorithm parameters. The impact of SQP, as local search, in the proposed algorithm is surveyed in Section 7. We conclude in Section 8.


## 2. Formulation of Problem

The bounded continuous-time NOCP is considered as finding the control inputu
(
t
)
∈
R
m, over the planning horizon [
t
0
,
t
f
], which minimizes the cost functional:(1)
J
=
ϕ
x
t
f
,
t
f
+
∫
t
0
t
f
‍
g
x
t
,
u
t
,
t
d
tsubject to(2)
x
˙
t
=
f
x
t
,
u
t
,
t
,
(3)
c
x
,
u
,
t
=
0
,
(4)
d
x
,
u
,
t
≤
0
,
(5)
ψ
x
t
f
,
t
f
=
0
,
(6)
x
t
0
=
x
0
,where x
(
t
)
∈
R
n denotes the state vector for the system and x
0
∈
R
n is the initial state. The functions f
:
R
n
×
R
m
×
R
→
R
n, g
:
R
n
×
R
m
×
R
→
R, c
:
R
n
×
R
m
×
R
→
R
n
c, d
:
R
n
×
R
m
×
R
→
R
n
d, ψ
:
R
n
×
R
→
R
n
ψ, and ϕ
:
R
n
×
R
→
R are assumed to be sufficiently smooth on appropriate open sets. 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 x and u. The minimum time problems, tracking problem, terminal control problem and minimum energy are another special case of NOCPs.


## 3. 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 N
t
-
1 subintervals using time nodes t
0
,
…
,
t
N
t
-
1 and then the control input values are computed (or selected randomly). This can be done by the following stages.(1)
Lett
j
=
t
0
+
j
h, where h
=
t
f
-
t
0
/
N
t
-
1, j
=
0,1
,
…
,
N
t
-
1, be time nodes, where t
0 and t
f are the initial and final times, respectively.
(2)
The corresponding control input value at each time nodet
j is an m
×
1 vector, u
j, which can be calculated randomly, with the following components.(7)
u
i
j
=
u
j
i
=
u
left
i
+
u
right
i
-
u
left
i
·
r
i
j
,
iiiiiiiiiiiiiii
i
=
1,2
,
…
,
m
,
j
=
0,1
,
…
,
N
t
-
1
,where r
i
j is a random number in [
0,1
] with a uniform distribution and u
left
,
u
right
∈
R
m 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 m
×
N
t matrix as U
=
(
u
i
j
)
m
×
N
t
=
[
u
j
]
j
=
0
N
t
-
1. Next, we let U
(
k
)
=
(
u
i
j
)
m
×
N
t, k
=
1,2
,
…
,
N
p as kth individual of the population, which N
p is the size of the population.Evaluation. For each control input matrix, U
(
k
), k
=
1,2
,
…
,
N
p, the corresponding state variable is an n
×
N
t matrix, X
(
k
), 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, J
(
U
(
k
)
), 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 I
(
U
(
k
)
) to U
(
k
) as the fitness value as follows:(8)
I
U
k
=
J
~
+
∑
l
=
1
n
d
‍
∑
j
=
0
N
t
-
1
‍
M
1
l
max
⁡
⁡
0
,
d
l
x
j
,
u
j
,
t
j
+
∑
h
=
1
n
c
‍
∑
j
=
0
N
t
-
1
‍
M
2
h
c
h
2
x
j
,
u
j
,
t
j
+
∑
i
=
p
n
ψ
‍
M
3
p
ψ
p
2
x
N
t
-
1
,
t
N
t
-
1
,where M
1
=
[
M
11
,
…
,
M
1
n
d
]
T, M
2
=
[
M
21
,
…
,
M
2
n
c
]
T and M
3
=
[
M
31
,
…
,
M
3
n
ψ
]
T are big numbers, as the penalty parameters, c
h
(
·
,
·
), h
=
1,2
,
…
,
n
c, d
l
(
·
,
·
), l
=
1,2
,
…
,
n
d, and ψ
p
(
·
,
·
), p
=
1,2
,
…
,
n
ψ are defined in (3), (4), and (5), respectively.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 N
tour. The tournament algorithm is given in Algorithm 1.Algorithm 1:Tournament algorithm.
{Initialization} Input the the size of tournament, N
tour and select N
tour
individuals from population, randomly. Letk
i
,
i
=
1,2
,
…
,
N
tour be the indexes of them.
repeat
if
N
tour
=
2
then
if
I
(
U
k
1
)
<
I
(
U
k
2
)
then
Letk
=
k
1.
else
Letk
=
k
2.
end if
end if
Perform Tournament algorithm with sizeN
tour
/
2 with output k
1.
Perform Tournament algorithm with sizeN
tour
/
2 with output k
2.
until
k
1
≠
k
2
Return the kth individual of population.Crossover. When two parents U
(
1
) and U
(
2
) are selected, we use the following stages to construct an offspring.(1)
Select the following numbers(9)
λ
1
∈
0,1
,
λ
2
∈
[
-
λ
max
⁡
,
0
]
,
λ
3
∈
1,1
+
λ
max
⁡randomly, where λ
max
⁡ is a random number in [
0,1
].
(2)
Let(10)
of
k
=
λ
k
U
(
1
)
+
1
-
λ
k
U
(
2
)
,
k
=
1,2
,
3
,where λ
k, k
=
1,2
,
3 is defined in (9). For i
=
1
,
…
,
m and j
=
1
,
…
,
N
t, if (
of
k
)
i
j
>
u
right
i, then let (
of
k
)
i
j
=
u
right
i and if (
of
k
)
i
j
<
u
left
i, then let (
of
k
)
i
j
=
u
left
i.
(3)
Letof
=
of
*, where of
* is the best of
k, k
=
1,2
,
3 constructed by (10).Mutation. We apply a perturbation on each component of the offspring as follows:(11)
of
i
j
=
of
i
j
+
r
i
j
·
α
,
i
=
1,2
,
…
,
m
,
j
=
1,2
,
…
,
N
t
,where r
i
j is selected randomly in {
-
1,1
} and α is a random number, with a uniform distribution in [
0,1
]. If (
of
)
i
j
>
u
right
i, then let (
of
)
i
j
=
u
right
i and if (
of
)
i
j
<
u
left
i, then let (
of
)
i
j
=
u
left
i.Replacement. Here, in the underling GA, we use a traditional replacement strategy. The replacement is done, if the new offspring has two properties: first, it is better than the worst person in the population, I
(
of
)
<
max
⁡
1
≤
i
≤
N
p
I
(
U
(
i
)
), second, it is not very similar to a person in the population; that is, for each i
=
1
,
…
,
N
p, at least one of the following conditions is satisfied:(12)
I
of
-
I
U
i
>
ɛ
´
,
of
-
U
(
i
)
>
ɛ
´
,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,N
g, is reached.
(2)
Over a specified number of generations,N
i, 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,(13)
φ
f
=
ψ
2
<
ɛ
,where ψ
=
[
ψ
1
,
ψ
2
,
…
,
ψ
n
ψ
]
T is the vector of final state constraints, defined in (5).


### 3.2. SQP

The fitness value in (8) can be viewed as a nonlinear objective function with the decision variable as u
=
[
u
0
,
u
1
,
…
,
u
N
t
-
1
]. This cost function with upper and lower bounds of input signals construct a finite dimensional NLP problem as following(14)

min
⁡
I
u
=
I
u
0
,
u
1
,
…
,
u
N
t
-
1
(15)
s.t.
iiiiiii
u
left
≤
u
j
≤
u
right
,
j
=
0,1
,
…
,
N
t
-
1SQP algorithm [22, 28] is performed on the NLP (14)-(15), using u
¯
0
=
of and constructed in (11), as the initial solution when the maximum number of iteration is s
q
p
m
a
x
i
t
e
r.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,u
¯
k, k
=
0,1
,
…, a sequence of solutions as u
¯
k
+
1
=
u
¯
k
+
d
k is constructed, which d
k is the optimal solution of the constructed quadratic programming (QP) that approximates NLP in the iteration k based on u
¯
k, 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 L
(
u
,
λ
)
=
I
(
u
)
+
λ
T
h
(
u
), where the vector λ is Lagrangian multiplier and h
(
u
) return the vector of, inequality constraints evaluated at u. The QP is obtained by linearizing the nonlinear functions as follows:(16)
min
⁡
1
2
d
T
H
u
¯
k
d
+
∇
I
u
¯
k
T
d
(17)
∇
h
u
¯
k
T
d
+
h
u
¯
k
≤
0
.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:(18)
∂
I
∂
u
j
=
I
(
⋯
u
j
+
δ
⋯
)
-
I
(
u
)
δ
,
j
=
0,1
,
…
,
N
t
-
1
,where δ is the double precision of machine. So, the gradient vector is ∇
I
=
[
∂
I
/
∂
u
0
,
…
,
∂
I
/
∂
u
N
t
-
1
]
T. Also, at each major iteration a positive definite quasi-Newton approximation of the Hessian of the Lagrangian function, H, is calculated using the BFGS method [28], where λ
i, i
=
1
,
…
,
m, is an estimate of the Lagrange multipliers. The general procedure for SQP, for NLP (14)-(15), is as follows.(1)
Given an initial solutionu
¯
0. Let k
=
0.
(2)
Construct the QP subproblem (16)-(17), based on u
¯
0, using the approximations of the gradient and the Hessian of the the Lagrangian function.
(3)
Compute the new point asu
¯
k
+
1
=
u
¯
k
+
d
k, where d
k is the optimal solution of the current QP.
(4)
Letk
=
k
+
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.


### 3.3. 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 Figure1.Figure 1
Flowchart of the MHGA algorithm.


## 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 N
t
-
1 subintervals using time nodes t
0
,
…
,
t
N
t
-
1 and then the control input values are computed (or selected randomly). This can be done by the following stages.(1)
Lett
j
=
t
0
+
j
h, where h
=
t
f
-
t
0
/
N
t
-
1, j
=
0,1
,
…
,
N
t
-
1, be time nodes, where t
0 and t
f are the initial and final times, respectively.
(2)
The corresponding control input value at each time nodet
j is an m
×
1 vector, u
j, which can be calculated randomly, with the following components.(7)
u
i
j
=
u
j
i
=
u
left
i
+
u
right
i
-
u
left
i
·
r
i
j
,
iiiiiiiiiiiiiii
i
=
1,2
,
…
,
m
,
j
=
0,1
,
…
,
N
t
-
1
,where r
i
j is a random number in [
0,1
] with a uniform distribution and u
left
,
u
right
∈
R
m 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 m
×
N
t matrix as U
=
(
u
i
j
)
m
×
N
t
=
[
u
j
]
j
=
0
N
t
-
1. Next, we let U
(
k
)
=
(
u
i
j
)
m
×
N
t, k
=
1,2
,
…
,
N
p as kth individual of the population, which N
p is the size of the population.Evaluation. For each control input matrix, U
(
k
), k
=
1,2
,
…
,
N
p, the corresponding state variable is an n
×
N
t matrix, X
(
k
), 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, J
(
U
(
k
)
), 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 I
(
U
(
k
)
) to U
(
k
) as the fitness value as follows:(8)
I
U
k
=
J
~
+
∑
l
=
1
n
d
‍
∑
j
=
0
N
t
-
1
‍
M
1
l
max
⁡
⁡
0
,
d
l
x
j
,
u
j
,
t
j
+
∑
h
=
1
n
c
‍
∑
j
=
0
N
t
-
1
‍
M
2
h
c
h
2
x
j
,
u
j
,
t
j
+
∑
i
=
p
n
ψ
‍
M
3
p
ψ
p
2
x
N
t
-
1
,
t
N
t
-
1
,where M
1
=
[
M
11
,
…
,
M
1
n
d
]
T, M
2
=
[
M
21
,
…
,
M
2
n
c
]
T and M
3
=
[
M
31
,
…
,
M
3
n
ψ
]
T are big numbers, as the penalty parameters, c
h
(
·
,
·
), h
=
1,2
,
…
,
n
c, d
l
(
·
,
·
), l
=
1,2
,
…
,
n
d, and ψ
p
(
·
,
·
), p
=
1,2
,
…
,
n
ψ are defined in (3), (4), and (5), respectively.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 N
tour. The tournament algorithm is given in Algorithm 1.Algorithm 1:Tournament algorithm.
{Initialization} Input the the size of tournament, N
tour and select N
tour
individuals from population, randomly. Letk
i
,
i
=
1,2
,
…
,
N
tour be the indexes of them.
repeat
if
N
tour
=
2
then
if
I
(
U
k
1
)
<
I
(
U
k
2
)
then
Letk
=
k
1.
else
Letk
=
k
2.
end if
end if
Perform Tournament algorithm with sizeN
tour
/
2 with output k
1.
Perform Tournament algorithm with sizeN
tour
/
2 with output k
2.
until
k
1
≠
k
2
Return the kth individual of population.Crossover. When two parents U
(
1
) and U
(
2
) are selected, we use the following stages to construct an offspring.(1)
Select the following numbers(9)
λ
1
∈
0,1
,
λ
2
∈
[
-
λ
max
⁡
,
0
]
,
λ
3
∈
1,1
+
λ
max
⁡randomly, where λ
max
⁡ is a random number in [
0,1
].
(2)
Let(10)
of
k
=
λ
k
U
(
1
)
+
1
-
λ
k
U
(
2
)
,
k
=
1,2
,
3
,where λ
k, k
=
1,2
,
3 is defined in (9). For i
=
1
,
…
,
m and j
=
1
,
…
,
N
t, if (
of
k
)
i
j
>
u
right
i, then let (
of
k
)
i
j
=
u
right
i and if (
of
k
)
i
j
<
u
left
i, then let (
of
k
)
i
j
=
u
left
i.
(3)
Letof
=
of
*, where of
* is the best of
k, k
=
1,2
,
3 constructed by (10).Mutation. We apply a perturbation on each component of the offspring as follows:(11)
of
i
j
=
of
i
j
+
r
i
j
·
α
,
i
=
1,2
,
…
,
m
,
j
=
1,2
,
…
,
N
t
,where r
i
j is selected randomly in {
-
1,1
} and α is a random number, with a uniform distribution in [
0,1
]. If (
of
)
i
j
>
u
right
i, then let (
of
)
i
j
=
u
right
i and if (
of
)
i
j
<
u
left
i, then let (
of
)
i
j
=
u
left
i.Replacement. Here, in the underling GA, we use a traditional replacement strategy. The replacement is done, if the new offspring has two properties: first, it is better than the worst person in the population, I
(
of
)
<
max
⁡
1
≤
i
≤
N
p
I
(
U
(
i
)
), second, it is not very similar to a person in the population; that is, for each i
=
1
,
…
,
N
p, at least one of the following conditions is satisfied:(12)
I
of
-
I
U
i
>
ɛ
´
,
of
-
U
(
i
)
>
ɛ
´
,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,N
g, is reached.
(2)
Over a specified number of generations,N
i, 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,(13)
φ
f
=
ψ
2
<
ɛ
,where ψ
=
[
ψ
1
,
ψ
2
,
…
,
ψ
n
ψ
]
T is the vector of final state constraints, defined in (5).


## 3.2. SQP

The fitness value in (8) can be viewed as a nonlinear objective function with the decision variable as u
=
[
u
0
,
u
1
,
…
,
u
N
t
-
1
]. This cost function with upper and lower bounds of input signals construct a finite dimensional NLP problem as following(14)

min
⁡
I
u
=
I
u
0
,
u
1
,
…
,
u
N
t
-
1
(15)
s.t.
iiiiiii
u
left
≤
u
j
≤
u
right
,
j
=
0,1
,
…
,
N
t
-
1SQP algorithm [22, 28] is performed on the NLP (14)-(15), using u
¯
0
=
of and constructed in (11), as the initial solution when the maximum number of iteration is s
q
p
m
a
x
i
t
e
r.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,u
¯
k, k
=
0,1
,
…, a sequence of solutions as u
¯
k
+
1
=
u
¯
k
+
d
k is constructed, which d
k is the optimal solution of the constructed quadratic programming (QP) that approximates NLP in the iteration k based on u
¯
k, 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 L
(
u
,
λ
)
=
I
(
u
)
+
λ
T
h
(
u
), where the vector λ is Lagrangian multiplier and h
(
u
) return the vector of, inequality constraints evaluated at u. The QP is obtained by linearizing the nonlinear functions as follows:(16)
min
⁡
1
2
d
T
H
u
¯
k
d
+
∇
I
u
¯
k
T
d
(17)
∇
h
u
¯
k
T
d
+
h
u
¯
k
≤
0
.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:(18)
∂
I
∂
u
j
=
I
(
⋯
u
j
+
δ
⋯
)
-
I
(
u
)
δ
,
j
=
0,1
,
…
,
N
t
-
1
,where δ is the double precision of machine. So, the gradient vector is ∇
I
=
[
∂
I
/
∂
u
0
,
…
,
∂
I
/
∂
u
N
t
-
1
]
T. Also, at each major iteration a positive definite quasi-Newton approximation of the Hessian of the Lagrangian function, H, is calculated using the BFGS method [28], where λ
i, i
=
1
,
…
,
m, is an estimate of the Lagrange multipliers. The general procedure for SQP, for NLP (14)-(15), is as follows.(1)
Given an initial solutionu
¯
0. Let k
=
0.
(2)
Construct the QP subproblem (16)-(17), based on u
¯
0, using the approximations of the gradient and the Hessian of the the Lagrangian function.
(3)
Compute the new point asu
¯
k
+
1
=
u
¯
k
+
d
k, where d
k is the optimal solution of the current QP.
(4)
Letk
=
k
+
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.


## 3.3. 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 Figure1.Figure 1
Flowchart of the MHGA algorithm.


## 4. 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 fromN
t
1 in phase 1 to N
t
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.Algorithm 2:The proposed algorithm.
{Initialization} Input the desired precision, ε, in (13), the penalty parameters,
M
i
,
i
=
1,2
,
3, in (8), the bounds of control input values in (7), u
left and u
right.
{Phase 1} Perform MHGA with a random population and N
t
1
,
N
p
1
,
N
i
1
,
N
g
1
,
P
m
1 and s
q
p
m
a
x
i
t
e
r
1.
{Construction of the initial population of the phase 2} Increase time nodes uniformly to N
t
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.
CreateN
p
2
-
N
p
1 new different individuals with N
t
2 time nodes, randomly.
{Phase 2} Perform MHGA with the constructed population and
N
t
2
,
N
p
2
,
N
i
2
,
N
g
2, P
m
2 and s
q
p
m
a
x
i
t
e
r
2.


## 5. 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. LetJ be the obtained performance index by an algorithm, φ
f, defined in (13), be the error of final state constraints, and J
* be the best obtained solution among all implementations, or the exact solution (when exists). Now the relative error of J, E
J, of the algorithm can be defined as(19)
E
J
=
J
-
J
*
J
*
.To more accurate study, we now define a new criterion, called factor, to compare the algorithms as follows:(20)
K
ψ
=
E
J
+
φ
f
.Note that K
ψ shows the summation of two important errors. Thus, based on K
ψ 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; includingN
t, N
p, N
i, N
g, P
m and s
q
p
m
a
x
i
t
e
r, in both phases in Algorithm 2, and the problem’s parameters including ɛ, in (13), M
i, i
=
1,2
,
3, in (8), u
left and u
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 s
q
p
m
a
x
i
t
e
r
i
=
4 and P
m
i
=
0.8, for i
=
1
,
2. Also, we consider N
g
1
=
N
g
2, N
i
1
=
N
i
2, N
p
1
=
9 and N
p
2
=
12. 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.Table 1
The best numerical results in 12 different runs of the SI and LI methods for NOCPs in the Appendix.
Problem
Algorithm
J
E
J
φ
f
K
ψ
Time
VDP
CGA
1.7404
0.0912
2.67
×
1
0
-
11
0.0912
501.28
LI
1.5949
0
1.08
×
1
0
-
13
1.08
×
1
0
-
13
282.82
SI
1.5950
6.27
×
1
0
-
5
9.99
×
1
0
-
14
6.27
×
1
0
-
5
239.14
CRP
CGA
1.63
×
1
0
-
2
0.2835
7.57
×
1
0
-
10
0.2835
501.28
LI
1.27
×
1
0
-
2
0
1.15
×
1
0
-
9
1.15
×
1
0
-
9
171.94
SI
1.27
×
1
0
-
2
0
5.99
×
1
0
-
9
5.99
×
1
0
-
9
124.05
FFRP
CGA
83.63
1.3256
4.65
×
1
0
-
3
1.3302
1413
LI
35.96
0
1.22
×
1
0
-
5
1.22
×
1
0
-
5
1268.58
SI
35.99
8.34
×
1
0
-
4
9.01
×
1
0
-
6
8.43
×
1
0
-
4
1343.40
MSNIC
IPSO
0.1727
0.0165
—
—
126.26
LI
0.1699
0
—
—
138.82
SI
0.1700
5.89
×
1
0
-
4
—
—
138.68
CSTCR
IPSO
0.1355
0.2098
—
—
31.34
LI
0.1120
0
—
—
103.58
SI
0.1120
0
—
—
106.73
No. 6
Bézier
-
5.3898
0.0251
—
—
NRa
LI
-
5.4309
0.0177
—
—
138.38
SI
-
5.4309
0.0177
—
—
107.34
Number 7
HPM
0.2353
0.1677
4.20
×
1
0
-
6
0.1677
NR
LI
0.2015
0
2.35
×
1
0
-
9
2.35
×
1
0
-
9
213.92
SI
0.2015
0
2.82
×
1
0
-
10
2.82
×
1
0
-
10
206.96
Number 8
SQP
6.36
×
1
0
-
6
2.20
×
1
0
9
—
—
NR
SUMT
5.15
×
1
0
-
6
1.78
×
1
0
9
—
—
NR
LI
2.89
×
1
0
-
15
0
—
—
72.82
SI
3.64
×
1
0
-
15
0.2595
—
—
68.68
Number 9
SQP
1.7950
0.0873
—
—
NR
SUMT
1.7980
0.0891
—
—
NR
LI
1.6509
0
—
—
639.10
SI
1.6509
0
—
—
595.0
Number 10
SQP
0.2163
0.3964
—
—
NR
SUMT
0.1703
0.0994
—
—
NR
LI
0.1549
0
—
—
625.70
SI
0.1549
0
—
—
676.62
Number 11
SQP
3.25
0.1648
0
0.1648
NR
SUMT
3.25
0.1648
0
0.1648
NR
LI
2.7901
0
1.62
×
1
0
-
9
1.62
×
1
0
-
9
147.78
SI
2.7901
0
5.86
×
1
0
-
10
5.86
×
1
0
-
10
219.46
Number 12
SQP
-
0.2490
4.0
×
1
0
-
3
0
4.0
×
1
0
-
3
NR
SUMT
-
0.2490
4.0
×
1
0
-
3
0
4.0
×
1
0
-
3
NR
LI
-
0.2500
0
2.85
×
1
0
-
8
2.85
×
1
0
-
8
547.71
SI
-
0.2500
0
3.90
×
1
0
-
10
3.90
×
1
0
-
10
545.48
Number 13
SQP
1.68
×
1
0
-
2
0.1748
0
0.1748
NR
SUMT
1.67
×
1
0
-
2
0.1678
0
0.1678
NR
LI
1.43
×
1
0
-
2
0
1.18
×
1
0
-
9
1.18
×
1
0
-
9
425.99
SI
1.44
×
1
0
-
2
7.0
×
1
0
-
3
6.32
×
1
0
-
9
7.0
×
1
0
-
3
490.27
Number 14
SQP
3.7220
0.0961
0
0.0961
NR
SUMT
3.7700
0.1103
0
0.1103
NR
LI
3.3956
0
7.86
×
1
0
-
7
7.86
×
1
0
-
7
1041.24
SI
3.3965
2.65
×
1
0
-
4
2.76
×
1
0
-
6
2.67
×
1
0
-
4
1151.44
Number 15
SQP
1.03
×
1
0
-
3
4.3368
0
4.3368
NR
SUMT
9.29
×
1
0
-
4
3.8135
0
3.8135
NR
LI
1.93
×
1
0
-
4
0
8.10
×
1
0
-
10
8.10
×
1
0
-
10
333.54
SI
1.94
×
1
0
-
4
5.20
×
1
0
-
3
1.64
×
1
0
-
9
5.20
×
1
0
-
3
340.14
Number 16
SQP
2.2120
0.0753
0
0.0753
NR
SUMT
2.2080
0.0734
0
0.0734
NR
LI
2.0571
0
3.94
×
1
0
-
12
3.94
×
1
0
-
12
469.82
SI
2.0571
0
5.60
×
1
0
-
13
5.60
×
1
0
-
13
704.23
Number 17
SQP
-
8.8690
2.25
×
1
0
-
5
0
2.25
×
1
0
-
5
NR
SUMT
-
8.8690
2.25
×
1
0
-
5
0
2.25
×
1
0
-
5
NR
LI
-
8.8692
0
1.45
×
1
0
-
10
1.45
×
1
0
-
10
321.83
SI
-
8.8692
0
2.79
×
1
0
-
10
2.79
×
1
0
-
10
309.59
Number 18
SQP
0.0368
0.1288
—
—
NR
SUMT
0.0386
0.1840
—
—
NR
LI
0.0326
0
—
—
551.63
SI
0.0326
0
—
—
606.71
Number 19
SQP
0.3439
4.9293
0
4.9293
NR
SUMT
0.3428
4.9103
0
4.9103
NR
LI
0.0689
0.1879
4.10
×
1
0
-
3
0.1919
1645.09
SI
0.0580
0
4.05
×
1
0
-
4
4.05
×
1
0
-
4
1787.83
Number 20
SQP
77.52
1.2554
0
1.2554
NR
SUMT
76.83
1.2353
0
1.2353
NR
LI
34.3716
2.32
×
1
0
-
5
2.25
×
1
0
-
4
2.48
×
1
0
-
4
815.19
SI
34.3708
0
1.57
×
1
0
-
4
1.57
×
1
0
-
4
804.20
aNot reported.Table 2
The problem parameters for NOCPs in the Appendix.
Parameters
Problem number
1
2
3
4
5
6
7
8
9
10
u
left
-
0.5
-
1.5
-
15
-
20
0
-
2
0
-
2
-
1
-
20
u
right
2
2
10
20
5
2
1
2
1
20
M
i
1
0
3
1
0
2
70
1
—
1
—
1
0
2
—
1
0
2
ε
1
0
-
12
1
0
-
10
1
0
-
3
—
—
—
—
—
—
1
0
-
9
Parameters
Problem no.
11
12
13
14
15
16
17
18
19
20
u
left
-
5
-
1
-
2
-
π
-
1
-
3
-
30
-
1
-
15
-
15
u
right
5
1
2
π
1
3
30
1
10
10
M
i
10
1
0
2
—
1
0
2
1
0
2
1
0
2
—
10
70
1
0
-
3
ε
1
0
-
9
1
0
-
11
1
0
-
11
—
1
0
-
10
1
0
-
10
1
0
-
10
—
1
0
-
3
1
0
-
4The reported numerical results of the proposed algorithm, for each NOCP, include the value of performance index,J, the relative error of J, E
J, defined in (19), the required computational time, Time, the norm of final state constraints, φ
f, defined in (13) and the factor, K
ψ, 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.
(2)
SI: spline interpolation.For comparing the numerical results of the proposed algorithm two subsections are considered, comparison with some metaheuristic algorithms, in Section5.1, and comparison with some numerical methods, in Section 5.2. We give more details of these comparisons in the following subsections.


### 5.1. 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].


#### 5.1.1. VDP Problem [9]

The first NOCP in the Appendix is Van Der Pol Problem, VDP, which has two state variables and one control variable. VDP problem has a final state constraint, which isψ
=
x
1
(
t
f
)
-
x
2
(
t
f
)
+
1
=
0. The results of the proposed algorithm with the MHGA’s parameters as N
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
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 K
ψ.


#### 5.1.2. CRP Problem [9]

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 asN
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
200, are shown in the second row of Table 1. CRP problem has two final state constraints, ψ
=
[
x
1
,
x
2
]
T. Although, from Table 1, the norm of final state constraints, φ
f, for the CGA, equals φ
f
*
=
7.57
×
1
0
-
10, is less than φ
f’s of LI and SI methods, which equals 1.15
×
1
0
-
9 and 5.99
×
1
0
-
9, respectively, but the performance index, the relative error of J and the factor of the proposed algorithm is better, and so the proposed algorithm is more robust than CGA.


#### 5.1.3. FFRP Problem [9]

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 asN
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
200. This problem has six final state constraints, ψ
=
[
x
1
-
4
,
x
2
,
x
3
-
4
,
x
4
,
x
5
,
x
6
]
T. The numerical results are shown in Table 1. The values of J, E
J, φ
f and K
ψ 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.


#### 5.1.4. MSNIC Problem [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,d
(
x
,
t
)
=
x
2
(
t
)
+
0.5
-
8
(
t
-
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 N
t
1
=
31, N
t
2
=
91, N
g
=
100 and N
i
=
60. From the forth row of Table 1, the absolute error of J, E
J, for LI and SI methods equal 0 and 5.89
×
1
0
-
4, respectively, which are less than IPSO’s, 0.0165.Subplots (a) and (b) in Figure2, show the graphs of the convergence rate for the performance index and the inequality constraint, respect to the number of iteration, respectively.Figure 2
Graphical results of MSNIC problem using SI method: (a) convergence rate of the performance index and (b) the inequality constraint,d
(
x
,
t
)
≤
0, in respect to the number of iterations.
(a)
(b)


#### 5.1.5. CSTCR Problem [20]

The fifth NOCP in the Appendix is a model of a nonlinear Continuous Stirred-tank Chemical Reactor, CSTCR. It has two state variablesx
1
(
t
) and x
2
(
t
), as the deviation from the steady-state temperature and concentration, and one control variable u
(
t
), 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 N
t
1
=
31, N
t
2
=
51, N
g
=
100, and N
i
=
50, are compared with IPSO. From Table 1, the performance index, J, for LI and SI methods is equal to J
*
=
0.1120 which is less than IPSOs 0.1355.


### 5.2. 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 N
t
1
=
31, N
t
2
=
51, N
g
=
100, N
i
=
50, with the problem parameters in Table 2.


#### 5.2.1. Compared with Bézier [34]

The NOCP number(
6
), in the Appendix, has exact solution, which has an inequality constraint as d
(
x
1
,
t
)
=
-
6
-
x
1
(
t
)
≤
0. The exact value of performance index equals J
*
=
-
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.Figure3, 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.Figure 3
Graphical results for NOCP number(
6
) using SI method: (a) convergence rate for the performance index and (b) the inequality constraint, d
(
x
1
,
t
)
≤
0, in respect to the number of iterations.
(a)
(b)


#### 5.2.2. Compared with HPM [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 ψ
=
x
-
0.5
=
0. From [36], the norm of final state constraint for HPM equals 4.2
×
1
0
-
6, however, this criterion for the LI and SI methods equals 2.35
×
1
0
-
9 and φ
f
*
=
2.82
×
1
0
-
10, respectively. From Table 1, it is obvious that the obtained values of the performance index, the norm of final state constraint and K
ψ for the SI method are more accurate than LI than HPM methods.


#### 5.2.3. Compared with SQP and SUMT [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 φ
f
=
0 to calculate the factor K
ψ.Figures4–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 Figures 7–13 and 15-16. 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.Figure 4
Convergence rate for the performance index of SI method for NOCP number(
8
).Figure 5
(a) Convergence rate for the performance index and (b) the inequality constraint of SI method, for NOCP number(
9
).
(a)
(b)Figure 6
(a) Convergence rate of the performance index and (b) the inequality constraint for SI method, for NOCP number(
10
).
(a)
(b)Figure 7
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
11
).
(a)
(b)Figure 8
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
12
).
(a)
(b)Figure 9
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
13
).
(a)
(b)Figure 10
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
14
).
(a)
(b)Figure 11
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
15
).
(a)
(b)Figure 12
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
16
).
(a)
(b)Figure 13
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
17
).
(a)
(b)Figure 14
Convergence rate for the performance index of SI method for NOCP numbers(
18
).Figure 15
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
19
).
(a)
(b)Figure 16
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
20
).
(a)
(b)Table1 shows that the proposed algorithm, LI and SI methods, was 100 percent successful in point of views the performance index, J, and the factor, K
ψ, numerically. So, the proposed algorithm provides robust solutions with respect to the other mentioned, numerical or metaheuristic, methods. To compare J 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 φ
f 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.To compare with CGA, the mean of relative error ofJ, E
J, for CGA, LI and SI methods, in NOCPs (
1
)
–
(
3
), equal 0.5668, 0 and 2.98
×
1
0
-
4, respectively. Also the mean for the error of the final state constraints, φ
f, are 1.6
×
1
0
-
3, 4.07
×
1
0
-
6 and 3.01
×
1
0
-
6, respectively, and for K
ψ, these values equal 0.5883, 4.07
×
1
0
-
6 and 3.02
×
1
0
-
4. Thus, we can say that the feasibility of the solutions given by the proposed algorithm and CGA are competitive.From Table1, to compare with numerical methods, SQP and SUMT, in NOCPs (
8
)
–
(
20
), the mean of E
J for LI, SI, SQP and SUMT equals 0.0145, 0.0209, 1.69
×
1
0
8 and 1.36
×
1
0
8, respectively. Also, the mean for the error of final state constraints, for these NOCPs, equal 3.34
×
1
0
-
4, 4.41
×
1
0
-
5, 0 and 0, respectively. For K
ψ, these values are 0.0213, 0.0014, 1.2263 and 1.1644. Therefore, the performance index, J, and the factor, K
ψ, 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 Table1.


## 5.1. 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].


### 5.1.1. VDP Problem [9]

The first NOCP in the Appendix is Van Der Pol Problem, VDP, which has two state variables and one control variable. VDP problem has a final state constraint, which isψ
=
x
1
(
t
f
)
-
x
2
(
t
f
)
+
1
=
0. The results of the proposed algorithm with the MHGA’s parameters as N
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
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 K
ψ.


### 5.1.2. CRP Problem [9]

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 asN
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
200, are shown in the second row of Table 1. CRP problem has two final state constraints, ψ
=
[
x
1
,
x
2
]
T. Although, from Table 1, the norm of final state constraints, φ
f, for the CGA, equals φ
f
*
=
7.57
×
1
0
-
10, is less than φ
f’s of LI and SI methods, which equals 1.15
×
1
0
-
9 and 5.99
×
1
0
-
9, respectively, but the performance index, the relative error of J and the factor of the proposed algorithm is better, and so the proposed algorithm is more robust than CGA.


### 5.1.3. FFRP Problem [9]

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 asN
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
200. This problem has six final state constraints, ψ
=
[
x
1
-
4
,
x
2
,
x
3
-
4
,
x
4
,
x
5
,
x
6
]
T. The numerical results are shown in Table 1. The values of J, E
J, φ
f and K
ψ 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.


### 5.1.4. MSNIC Problem [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,d
(
x
,
t
)
=
x
2
(
t
)
+
0.5
-
8
(
t
-
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 N
t
1
=
31, N
t
2
=
91, N
g
=
100 and N
i
=
60. From the forth row of Table 1, the absolute error of J, E
J, for LI and SI methods equal 0 and 5.89
×
1
0
-
4, respectively, which are less than IPSO’s, 0.0165.Subplots (a) and (b) in Figure2, show the graphs of the convergence rate for the performance index and the inequality constraint, respect to the number of iteration, respectively.Figure 2
Graphical results of MSNIC problem using SI method: (a) convergence rate of the performance index and (b) the inequality constraint,d
(
x
,
t
)
≤
0, in respect to the number of iterations.
(a)
(b)


### 5.1.5. CSTCR Problem [20]

The fifth NOCP in the Appendix is a model of a nonlinear Continuous Stirred-tank Chemical Reactor, CSTCR. It has two state variablesx
1
(
t
) and x
2
(
t
), as the deviation from the steady-state temperature and concentration, and one control variable u
(
t
), 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 N
t
1
=
31, N
t
2
=
51, N
g
=
100, and N
i
=
50, are compared with IPSO. From Table 1, the performance index, J, for LI and SI methods is equal to J
*
=
0.1120 which is less than IPSOs 0.1355.


## 5.1.1. VDP Problem [9]

The first NOCP in the Appendix is Van Der Pol Problem, VDP, which has two state variables and one control variable. VDP problem has a final state constraint, which isψ
=
x
1
(
t
f
)
-
x
2
(
t
f
)
+
1
=
0. The results of the proposed algorithm with the MHGA’s parameters as N
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
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 K
ψ.


## 5.1.2. CRP Problem [9]

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 asN
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
200, are shown in the second row of Table 1. CRP problem has two final state constraints, ψ
=
[
x
1
,
x
2
]
T. Although, from Table 1, the norm of final state constraints, φ
f, for the CGA, equals φ
f
*
=
7.57
×
1
0
-
10, is less than φ
f’s of LI and SI methods, which equals 1.15
×
1
0
-
9 and 5.99
×
1
0
-
9, respectively, but the performance index, the relative error of J and the factor of the proposed algorithm is better, and so the proposed algorithm is more robust than CGA.


## 5.1.3. FFRP Problem [9]

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 asN
t
1
=
31, N
t
2
=
71, N
g
=
300 and N
i
=
200. This problem has six final state constraints, ψ
=
[
x
1
-
4
,
x
2
,
x
3
-
4
,
x
4
,
x
5
,
x
6
]
T. The numerical results are shown in Table 1. The values of J, E
J, φ
f and K
ψ 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.


## 5.1.4. MSNIC Problem [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,d
(
x
,
t
)
=
x
2
(
t
)
+
0.5
-
8
(
t
-
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 N
t
1
=
31, N
t
2
=
91, N
g
=
100 and N
i
=
60. From the forth row of Table 1, the absolute error of J, E
J, for LI and SI methods equal 0 and 5.89
×
1
0
-
4, respectively, which are less than IPSO’s, 0.0165.Subplots (a) and (b) in Figure2, show the graphs of the convergence rate for the performance index and the inequality constraint, respect to the number of iteration, respectively.Figure 2
Graphical results of MSNIC problem using SI method: (a) convergence rate of the performance index and (b) the inequality constraint,d
(
x
,
t
)
≤
0, in respect to the number of iterations.
(a)
(b)


## 5.1.5. CSTCR Problem [20]

The fifth NOCP in the Appendix is a model of a nonlinear Continuous Stirred-tank Chemical Reactor, CSTCR. It has two state variablesx
1
(
t
) and x
2
(
t
), as the deviation from the steady-state temperature and concentration, and one control variable u
(
t
), 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 N
t
1
=
31, N
t
2
=
51, N
g
=
100, and N
i
=
50, are compared with IPSO. From Table 1, the performance index, J, for LI and SI methods is equal to J
*
=
0.1120 which is less than IPSOs 0.1355.


## 5.2. 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 N
t
1
=
31, N
t
2
=
51, N
g
=
100, N
i
=
50, with the problem parameters in Table 2.


### 5.2.1. Compared with Bézier [34]

The NOCP number(
6
), in the Appendix, has exact solution, which has an inequality constraint as d
(
x
1
,
t
)
=
-
6
-
x
1
(
t
)
≤
0. The exact value of performance index equals J
*
=
-
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.Figure3, 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.Figure 3
Graphical results for NOCP number(
6
) using SI method: (a) convergence rate for the performance index and (b) the inequality constraint, d
(
x
1
,
t
)
≤
0, in respect to the number of iterations.
(a)
(b)


### 5.2.2. Compared with HPM [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 ψ
=
x
-
0.5
=
0. From [36], the norm of final state constraint for HPM equals 4.2
×
1
0
-
6, however, this criterion for the LI and SI methods equals 2.35
×
1
0
-
9 and φ
f
*
=
2.82
×
1
0
-
10, respectively. From Table 1, it is obvious that the obtained values of the performance index, the norm of final state constraint and K
ψ for the SI method are more accurate than LI than HPM methods.


### 5.2.3. Compared with SQP and SUMT [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 φ
f
=
0 to calculate the factor K
ψ.Figures4–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 Figures 7–13 and 15-16. 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.Figure 4
Convergence rate for the performance index of SI method for NOCP number(
8
).Figure 5
(a) Convergence rate for the performance index and (b) the inequality constraint of SI method, for NOCP number(
9
).
(a)
(b)Figure 6
(a) Convergence rate of the performance index and (b) the inequality constraint for SI method, for NOCP number(
10
).
(a)
(b)Figure 7
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
11
).
(a)
(b)Figure 8
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
12
).
(a)
(b)Figure 9
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
13
).
(a)
(b)Figure 10
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
14
).
(a)
(b)Figure 11
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
15
).
(a)
(b)Figure 12
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
16
).
(a)
(b)Figure 13
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
17
).
(a)
(b)Figure 14
Convergence rate for the performance index of SI method for NOCP numbers(
18
).Figure 15
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
19
).
(a)
(b)Figure 16
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
20
).
(a)
(b)Table1 shows that the proposed algorithm, LI and SI methods, was 100 percent successful in point of views the performance index, J, and the factor, K
ψ, numerically. So, the proposed algorithm provides robust solutions with respect to the other mentioned, numerical or metaheuristic, methods. To compare J 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 φ
f 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.To compare with CGA, the mean of relative error ofJ, E
J, for CGA, LI and SI methods, in NOCPs (
1
)
–
(
3
), equal 0.5668, 0 and 2.98
×
1
0
-
4, respectively. Also the mean for the error of the final state constraints, φ
f, are 1.6
×
1
0
-
3, 4.07
×
1
0
-
6 and 3.01
×
1
0
-
6, respectively, and for K
ψ, these values equal 0.5883, 4.07
×
1
0
-
6 and 3.02
×
1
0
-
4. Thus, we can say that the feasibility of the solutions given by the proposed algorithm and CGA are competitive.From Table1, to compare with numerical methods, SQP and SUMT, in NOCPs (
8
)
–
(
20
), the mean of E
J for LI, SI, SQP and SUMT equals 0.0145, 0.0209, 1.69
×
1
0
8 and 1.36
×
1
0
8, respectively. Also, the mean for the error of final state constraints, for these NOCPs, equal 3.34
×
1
0
-
4, 4.41
×
1
0
-
5, 0 and 0, respectively. For K
ψ, these values are 0.0213, 0.0014, 1.2263 and 1.1644. Therefore, the performance index, J, and the factor, K
ψ, 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 Table1.


## 5.2.1. Compared with Bézier [34]

The NOCP number(
6
), in the Appendix, has exact solution, which has an inequality constraint as d
(
x
1
,
t
)
=
-
6
-
x
1
(
t
)
≤
0. The exact value of performance index equals J
*
=
-
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.Figure3, 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.Figure 3
Graphical results for NOCP number(
6
) using SI method: (a) convergence rate for the performance index and (b) the inequality constraint, d
(
x
1
,
t
)
≤
0, in respect to the number of iterations.
(a)
(b)


## 5.2.2. Compared with HPM [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 ψ
=
x
-
0.5
=
0. From [36], the norm of final state constraint for HPM equals 4.2
×
1
0
-
6, however, this criterion for the LI and SI methods equals 2.35
×
1
0
-
9 and φ
f
*
=
2.82
×
1
0
-
10, respectively. From Table 1, it is obvious that the obtained values of the performance index, the norm of final state constraint and K
ψ for the SI method are more accurate than LI than HPM methods.


## 5.2.3. Compared with SQP and SUMT [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 φ
f
=
0 to calculate the factor K
ψ.Figures4–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 Figures 7–13 and 15-16. 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.Figure 4
Convergence rate for the performance index of SI method for NOCP number(
8
).Figure 5
(a) Convergence rate for the performance index and (b) the inequality constraint of SI method, for NOCP number(
9
).
(a)
(b)Figure 6
(a) Convergence rate of the performance index and (b) the inequality constraint for SI method, for NOCP number(
10
).
(a)
(b)Figure 7
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
11
).
(a)
(b)Figure 8
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
12
).
(a)
(b)Figure 9
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
13
).
(a)
(b)Figure 10
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
14
).
(a)
(b)Figure 11
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
15
).
(a)
(b)Figure 12
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
16
).
(a)
(b)Figure 13
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
17
).
(a)
(b)Figure 14
Convergence rate for the performance index of SI method for NOCP numbers(
18
).Figure 15
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
19
).
(a)
(b)Figure 16
(a) Convergence rates of the performance index and (b) the error of final state constraint for SI method, for NOCP number(
20
).
(a)
(b)Table1 shows that the proposed algorithm, LI and SI methods, was 100 percent successful in point of views the performance index, J, and the factor, K
ψ, numerically. So, the proposed algorithm provides robust solutions with respect to the other mentioned, numerical or metaheuristic, methods. To compare J 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 φ
f 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.To compare with CGA, the mean of relative error ofJ, E
J, for CGA, LI and SI methods, in NOCPs (
1
)
–
(
3
), equal 0.5668, 0 and 2.98
×
1
0
-
4, respectively. Also the mean for the error of the final state constraints, φ
f, are 1.6
×
1
0
-
3, 4.07
×
1
0
-
6 and 3.01
×
1
0
-
6, respectively, and for K
ψ, these values equal 0.5883, 4.07
×
1
0
-
6 and 3.02
×
1
0
-
4. Thus, we can say that the feasibility of the solutions given by the proposed algorithm and CGA are competitive.From Table1, to compare with numerical methods, SQP and SUMT, in NOCPs (
8
)
–
(
20
), the mean of E
J for LI, SI, SQP and SUMT equals 0.0145, 0.0209, 1.69
×
1
0
8 and 1.36
×
1
0
8, respectively. Also, the mean for the error of final state constraints, for these NOCPs, equal 3.34
×
1
0
-
4, 4.41
×
1
0
-
5, 0 and 0, respectively. For K
ψ, these values are 0.0213, 0.0014, 1.2263 and 1.1644. Therefore, the performance index, J, and the factor, K
ψ, 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 Table1.


## 6. Sensitivity Analysis and Comparing LI and SI

In this section, two statistical analysis, based on the one-way 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.


### 6.1. 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,J, the relative error of J, E
J, the required computational time, Time, the error of final state constraints, φ
f and the factor, K
ψ. The independent parameters are consist of the number of time nodes in both two phases, N
t
1 and N
t
2, the size of population in both two phases, N
p
1 and N
p
2, the maximum number of generations without improvement, N
i, and the maximum number of generations, N
g. Because the mutation implementation probability, P
m, has a less influence in numerical results, it is not considered. Also s
q
p
m
a
x
i
t
e
r 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(
N
), the mean of each output in N 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 N
t
1 parameter, which is reported in Table 3.Table 3
Descriptive statistics for the sensitivity analysis ofN
t
1, for VDP problem.
Output
N
t
1
N
Mean
S.d
Min
Max
J
11
35
1.5333
0.0069
1.5286
1.5598
21
35
1.5328
0.0097
1.5285
1.5849
31
35
1.5329
0.0076
1.5285
1.5699
41
35
1.5455
0.0607
1.5285
1.8449
Total
140
1.5357
0.0294
1.5285
1.8449
φ
f
11
35
2.81
×
1
0
-
8
9.52
×
1
0
-
8
1
×
1
0
-
12
4.52
×
1
0
-
7
21
35
4.67
×
1
0
-
8
2.35
×
1
0
-
7
1
×
1
0
-
12
1.41
×
1
0
-
6
31
35
1.63
×
1
0
-
8
4.91
×
1
0
-
8
1
×
1
0
-
12
2.21
×
1
0
-
7
41
35
8.17
×
1
0
-
8
2.85
×
1
0
-
7
1
×
1
0
-
12
1.26
×
1
0
-
6
Total
140
4.26
×
1
0
-
8
1.91
×
1
0
-
7
1
×
1
0
-
12
1.41
×
1
0
-
6
Time
11
35
89.7838
27.2081
45.6614
155.4394
21
35
109.8147
46.0264
60.4192
284.0310
31
35
116.3591
30.0227
77.2673
179.5103
41
35
144.8290
39.3091
81.0425
284.9514
Total
140
114.2131
41.0891
45.6615
284.9514
E
J
11
35
0.0031
0.0045
0
0.0204
21
35
0.0028
0.0064
0
0.0369
31
35
0.0029
0.0050
0
0.0271
41
35
0.0111
0.0397
0
0.2070
Total
140
0.0047
0.0192
0
0.2070
K
ψ
11
35
0.0031
0.0045
4.75
×
1
0
-
9
0.0204
21
35
0.0027
0.0064
5.50
×
1
0
-
10
0.0369
31
35
0.0029
0.0050
1.54
×
1
0
-
13
0.0271
41
35
0.0111
0.0397
1.93
×
1
0
-
10
0.2070
Total
140
0.0047
0.0192
1.54
×
1
0
-
13
0.2070Table4, summarized the statistical data, contain the test statistics (F) and P-values, of ANOVA tests. Sensitivity analysis for each parameter, separately, are done based on this table, as follows.N
t
1
:
From Table4, the significant level, or P-value, for K
ψ is equal to 0.279, which is greater than 0.05. So, from ANOVA, N
t
1 parameter has no significant effect on the factor K
ψ. With similar analysis, this parameter has no effect on the other outputs, except required computational time, Time, because its P-value is equal to 0, which is less than 0.05, that is, this output sensitive to the parameter N
t
1.
N
t
2
:
From the second row of Table4, the outputs J, E
J, K
ψ and Time are sensitive to the parameter N
t
2, and the only output φ
f is independent.
N
p
1
:
From the third row of Table4, only the computational required time, Time, is sensitive, because the its P-value is equal to 0.006, which is less than 0.05, and the other outputs, contain J, E
J, φ
f, K
ψ, are independent.
N
p
2
:
From the forth row of Table4, the output φ
f is independent respect to the parameter N
p
2, but other outputs, contain J, E
J, Time, K
ϕ are sensitive to this parameter.
N
g
:
From the fifth row of Table4, the outputs J, E
J and K
ψ are independent to the parameter N
g, and other outputs, contain φ
f and Time are sensitive respect to this parameter.
N
i
:
The sensitivity analysis is similar toN
p
1. From above cases, it is obvious that all parameters can be effect on the required computational time, except N
g. Moreover, intuitions shows the norm of final state constraints, φ
f 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.Table 4
Summary statistical data of ANOVA test for the parameters,N
t
1
,
N
t
2
,
N
p
1
,
N
p
2
,
N
g
,
N
i.
Parameters
J
E
J
φ
f
Time
K
ψ
Test statistic (F)
N
t
1
1.294
1.296
0.624
10.79
1.296
N
t
2
1105.72
3.93
0.868
59.39
3.93
N
p
1
1.945
1.835
1.271
4.317
1.835
N
p
2
2.740
2.478
1.011
23.02
2.478
N
g
3.541
3.591
0.491
0.890
3.591
N
i
0.495
0.309
0.301
5.849
1.046
P value
N
t
1
0.280
0.279
0.601
0
0.279
N
t
2
0
0.005
0.489
0
0.005
N
p
1
0.125
0.143
0.286
0.006
0.143
N
p
2
0.021
0.034
0.413
0
0.034
N
g
0.009
0.008
0.743
0.472
0.008
N
i
0.739
0.819
0.877
0
0.386


### 6.2. 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. Table 5, shows the mean of numerical results, separately. Here, we apply the fix MHGA parameters as N
p
1
=
9, N
p
2
=
12, N
t
1
=
31, N
t
2
=
91, N
g
=
100 and N
i
=
50 with the previous problem parameters in Table 2. Comparison criterions contain: J, φ
f, K
ψ and Time. The results of Mann-Whitney test are shown in Tables 6–9. From Table 6, since Z
=
0
>
Z
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 J. Similarly, the results of Tables 7–9 indicate that LI and SI methods, from the perspective of errors for φ
f, Time, and K
ψ, have same behaviour.Table 5
The mean of numerical results of 15 different runs for NOCPs in the Appendix, using LI and SI methods.
Problem
LI
SI
J
φ
f
K
ψ
Time
J
φ
f
K
ψ
Time
VDP
1.5977
1.57
×
1
0
-
9
1.70
×
1
0
-
3
177.20
1.6002
6.67
×
1
0
-
8
3.26
×
1
0
-
3
157.57
CRP
0.0126
3.07
×
1
0
-
8
3.07
×
1
0
-
8
197.14
0.0119
1.61
×
1
0
-
8
1.61
×
1
0
-
8
195.32
FFPP
47.74
8.05
×
1
0
-
5
0.3277
1332.80
53.51
5.92
×
1
0
-
5
0.4868
1379.89
MSNIC
0.1702
—
—
136.78
0.1703
—
—
140.92
CSTCR
0.1120
—
—
173.85
0.1120
—
—
176.06
Number 6
-
5.4307
—
—
105.41
-
5.4308
—
—
103.72
Number 7
0.2015
3.26
×
1
0
-
9
6.95
×
1
0
-
6
208.07
0.2020
3.13
×
1
0
-
9
6.62
×
1
0
-
6
210.19
Number 8
6.30
×
1
0
-
15
—
—
63.76
1.24
×
1
0
-
14
—
—
58.37
Number 9
1.6512
—
—
592.39
1.6512
—
—
581.25
Number 10
0.1550
—
—
631.53
0.1550
—
—
655.97
Number 11
2.9701
2.70
×
1
0
-
9
2.72
×
1
0
-
7
228.60
2.9701
3.37
×
1
0
-
9
8.11
×
1
0
-
7
232.76
Number 12
-
0.2499
2.01
×
1
0
-
8
4.64
×
1
0
-
5
547.18
-
0.2499
1.15
×
1
0
-
9
8.27
×
1
0
-
5
527.65
Number 13
0.0147
1.15
×
1
0
-
8
0.0252
577.97
0.0151
1.39
×
1
0
-
8
0.0472
660.98
Number 14
3.4761
6.72
×
1
0
-
6
0.0237
1045.85
3.4290
7.57
×
1
0
-
6
9.50
×
1
0
-
3
1087.67
Number 15
1.95
×
1
0
-
4
1.11
×
1
0
-
8
5.69
×
1
0
-
3
363
1.94
×
1
0
-
4
7.04
×
1
0
-
9
4.86
×
1
0
-
3
365
Number 16
2.0571
1.66
×
1
0
-
10
1.66
×
1
0
-
10
543.91
2.0571
1.15
×
1
0
-
10
1.15
×
1
0
-
10
544.31
Number 17
-
8.8692
8.03
×
1
0
-
8
8.03
×
1
0
-
8
227.44
-
8.8692
6.61
×
1
0
-
8
6.61
×
1
0
-
8
232.54
Number 18
0.0326
—
—
659.57
0.0326
—
—
708.97
Number 19
0.11588
6.55
×
1
0
-
4
0.6828
1618.13
0.1091
8.84
×
1
0
-
4
0.8808
1599.90
Number 20
52.42
0.0176
0.5429
841.65
44.64
4.53
×
1
0
-
4
0.3765
868.91Table 6
Results of Mann-Whitney test on relative errors of the pair (LI, SI) forJ.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
20.50
410.0
Mann-WhitneyU
200
SI
20.50
410.0
WilcoxonW
410
—
—
—
Z
0Table 7
Results of Mann-Whitney test on relative errors of the pair (LI, SI) forφ
f.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
13.62
177
Mann-WhitneyU
83
SI
13.38
174
WilcoxonW
174
—
—
—
Z
-
0.077Table 8
Results of Mann-Whitney test on relative errors of the pair (LI, SI) for Time.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
20.25
405.0
Mann-WhitneyU
195
SI
20.75
415.0
WilcoxonW
405
—
—
—
Z
-
0.135Table 9
Results of Mann-Whitney test on relative errors of the pair (LI, SI) forK
ψ.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
13.54
176.0
Mann-WhitneyU
84
SI
13.46
175.0
WilcoxonW
175
—
—
—
Z
-
0.026


## 6.1. 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,J, the relative error of J, E
J, the required computational time, Time, the error of final state constraints, φ
f and the factor, K
ψ. The independent parameters are consist of the number of time nodes in both two phases, N
t
1 and N
t
2, the size of population in both two phases, N
p
1 and N
p
2, the maximum number of generations without improvement, N
i, and the maximum number of generations, N
g. Because the mutation implementation probability, P
m, has a less influence in numerical results, it is not considered. Also s
q
p
m
a
x
i
t
e
r 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(
N
), the mean of each output in N 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 N
t
1 parameter, which is reported in Table 3.Table 3
Descriptive statistics for the sensitivity analysis ofN
t
1, for VDP problem.
Output
N
t
1
N
Mean
S.d
Min
Max
J
11
35
1.5333
0.0069
1.5286
1.5598
21
35
1.5328
0.0097
1.5285
1.5849
31
35
1.5329
0.0076
1.5285
1.5699
41
35
1.5455
0.0607
1.5285
1.8449
Total
140
1.5357
0.0294
1.5285
1.8449
φ
f
11
35
2.81
×
1
0
-
8
9.52
×
1
0
-
8
1
×
1
0
-
12
4.52
×
1
0
-
7
21
35
4.67
×
1
0
-
8
2.35
×
1
0
-
7
1
×
1
0
-
12
1.41
×
1
0
-
6
31
35
1.63
×
1
0
-
8
4.91
×
1
0
-
8
1
×
1
0
-
12
2.21
×
1
0
-
7
41
35
8.17
×
1
0
-
8
2.85
×
1
0
-
7
1
×
1
0
-
12
1.26
×
1
0
-
6
Total
140
4.26
×
1
0
-
8
1.91
×
1
0
-
7
1
×
1
0
-
12
1.41
×
1
0
-
6
Time
11
35
89.7838
27.2081
45.6614
155.4394
21
35
109.8147
46.0264
60.4192
284.0310
31
35
116.3591
30.0227
77.2673
179.5103
41
35
144.8290
39.3091
81.0425
284.9514
Total
140
114.2131
41.0891
45.6615
284.9514
E
J
11
35
0.0031
0.0045
0
0.0204
21
35
0.0028
0.0064
0
0.0369
31
35
0.0029
0.0050
0
0.0271
41
35
0.0111
0.0397
0
0.2070
Total
140
0.0047
0.0192
0
0.2070
K
ψ
11
35
0.0031
0.0045
4.75
×
1
0
-
9
0.0204
21
35
0.0027
0.0064
5.50
×
1
0
-
10
0.0369
31
35
0.0029
0.0050
1.54
×
1
0
-
13
0.0271
41
35
0.0111
0.0397
1.93
×
1
0
-
10
0.2070
Total
140
0.0047
0.0192
1.54
×
1
0
-
13
0.2070Table4, summarized the statistical data, contain the test statistics (F) and P-values, of ANOVA tests. Sensitivity analysis for each parameter, separately, are done based on this table, as follows.N
t
1
:
From Table4, the significant level, or P-value, for K
ψ is equal to 0.279, which is greater than 0.05. So, from ANOVA, N
t
1 parameter has no significant effect on the factor K
ψ. With similar analysis, this parameter has no effect on the other outputs, except required computational time, Time, because its P-value is equal to 0, which is less than 0.05, that is, this output sensitive to the parameter N
t
1.
N
t
2
:
From the second row of Table4, the outputs J, E
J, K
ψ and Time are sensitive to the parameter N
t
2, and the only output φ
f is independent.
N
p
1
:
From the third row of Table4, only the computational required time, Time, is sensitive, because the its P-value is equal to 0.006, which is less than 0.05, and the other outputs, contain J, E
J, φ
f, K
ψ, are independent.
N
p
2
:
From the forth row of Table4, the output φ
f is independent respect to the parameter N
p
2, but other outputs, contain J, E
J, Time, K
ϕ are sensitive to this parameter.
N
g
:
From the fifth row of Table4, the outputs J, E
J and K
ψ are independent to the parameter N
g, and other outputs, contain φ
f and Time are sensitive respect to this parameter.
N
i
:
The sensitivity analysis is similar toN
p
1. From above cases, it is obvious that all parameters can be effect on the required computational time, except N
g. Moreover, intuitions shows the norm of final state constraints, φ
f 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.Table 4
Summary statistical data of ANOVA test for the parameters,N
t
1
,
N
t
2
,
N
p
1
,
N
p
2
,
N
g
,
N
i.
Parameters
J
E
J
φ
f
Time
K
ψ
Test statistic (F)
N
t
1
1.294
1.296
0.624
10.79
1.296
N
t
2
1105.72
3.93
0.868
59.39
3.93
N
p
1
1.945
1.835
1.271
4.317
1.835
N
p
2
2.740
2.478
1.011
23.02
2.478
N
g
3.541
3.591
0.491
0.890
3.591
N
i
0.495
0.309
0.301
5.849
1.046
P value
N
t
1
0.280
0.279
0.601
0
0.279
N
t
2
0
0.005
0.489
0
0.005
N
p
1
0.125
0.143
0.286
0.006
0.143
N
p
2
0.021
0.034
0.413
0
0.034
N
g
0.009
0.008
0.743
0.472
0.008
N
i
0.739
0.819
0.877
0
0.386


## 6.2. 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. Table 5, shows the mean of numerical results, separately. Here, we apply the fix MHGA parameters as N
p
1
=
9, N
p
2
=
12, N
t
1
=
31, N
t
2
=
91, N
g
=
100 and N
i
=
50 with the previous problem parameters in Table 2. Comparison criterions contain: J, φ
f, K
ψ and Time. The results of Mann-Whitney test are shown in Tables 6–9. From Table 6, since Z
=
0
>
Z
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 J. Similarly, the results of Tables 7–9 indicate that LI and SI methods, from the perspective of errors for φ
f, Time, and K
ψ, have same behaviour.Table 5
The mean of numerical results of 15 different runs for NOCPs in the Appendix, using LI and SI methods.
Problem
LI
SI
J
φ
f
K
ψ
Time
J
φ
f
K
ψ
Time
VDP
1.5977
1.57
×
1
0
-
9
1.70
×
1
0
-
3
177.20
1.6002
6.67
×
1
0
-
8
3.26
×
1
0
-
3
157.57
CRP
0.0126
3.07
×
1
0
-
8
3.07
×
1
0
-
8
197.14
0.0119
1.61
×
1
0
-
8
1.61
×
1
0
-
8
195.32
FFPP
47.74
8.05
×
1
0
-
5
0.3277
1332.80
53.51
5.92
×
1
0
-
5
0.4868
1379.89
MSNIC
0.1702
—
—
136.78
0.1703
—
—
140.92
CSTCR
0.1120
—
—
173.85
0.1120
—
—
176.06
Number 6
-
5.4307
—
—
105.41
-
5.4308
—
—
103.72
Number 7
0.2015
3.26
×
1
0
-
9
6.95
×
1
0
-
6
208.07
0.2020
3.13
×
1
0
-
9
6.62
×
1
0
-
6
210.19
Number 8
6.30
×
1
0
-
15
—
—
63.76
1.24
×
1
0
-
14
—
—
58.37
Number 9
1.6512
—
—
592.39
1.6512
—
—
581.25
Number 10
0.1550
—
—
631.53
0.1550
—
—
655.97
Number 11
2.9701
2.70
×
1
0
-
9
2.72
×
1
0
-
7
228.60
2.9701
3.37
×
1
0
-
9
8.11
×
1
0
-
7
232.76
Number 12
-
0.2499
2.01
×
1
0
-
8
4.64
×
1
0
-
5
547.18
-
0.2499
1.15
×
1
0
-
9
8.27
×
1
0
-
5
527.65
Number 13
0.0147
1.15
×
1
0
-
8
0.0252
577.97
0.0151
1.39
×
1
0
-
8
0.0472
660.98
Number 14
3.4761
6.72
×
1
0
-
6
0.0237
1045.85
3.4290
7.57
×
1
0
-
6
9.50
×
1
0
-
3
1087.67
Number 15
1.95
×
1
0
-
4
1.11
×
1
0
-
8
5.69
×
1
0
-
3
363
1.94
×
1
0
-
4
7.04
×
1
0
-
9
4.86
×
1
0
-
3
365
Number 16
2.0571
1.66
×
1
0
-
10
1.66
×
1
0
-
10
543.91
2.0571
1.15
×
1
0
-
10
1.15
×
1
0
-
10
544.31
Number 17
-
8.8692
8.03
×
1
0
-
8
8.03
×
1
0
-
8
227.44
-
8.8692
6.61
×
1
0
-
8
6.61
×
1
0
-
8
232.54
Number 18
0.0326
—
—
659.57
0.0326
—
—
708.97
Number 19
0.11588
6.55
×
1
0
-
4
0.6828
1618.13
0.1091
8.84
×
1
0
-
4
0.8808
1599.90
Number 20
52.42
0.0176
0.5429
841.65
44.64
4.53
×
1
0
-
4
0.3765
868.91Table 6
Results of Mann-Whitney test on relative errors of the pair (LI, SI) forJ.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
20.50
410.0
Mann-WhitneyU
200
SI
20.50
410.0
WilcoxonW
410
—
—
—
Z
0Table 7
Results of Mann-Whitney test on relative errors of the pair (LI, SI) forφ
f.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
13.62
177
Mann-WhitneyU
83
SI
13.38
174
WilcoxonW
174
—
—
—
Z
-
0.077Table 8
Results of Mann-Whitney test on relative errors of the pair (LI, SI) for Time.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
20.25
405.0
Mann-WhitneyU
195
SI
20.75
415.0
WilcoxonW
405
—
—
—
Z
-
0.135Table 9
Results of Mann-Whitney test on relative errors of the pair (LI, SI) forK
ψ.
Method
Mean rank
Sum of ranks
Test statistics
Value
LI
13.54
176.0
Mann-WhitneyU
84
SI
13.46
175.0
WilcoxonW
175
—
—
—
Z
-
0.026


## 7. 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 Algorithm2, 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.Table 10
The numerical results of the without SQP algorithm, for NOCPs in the Appendix.
Problem no.
1
2
3
4
5
6
7
8
9
10
J
1.7296
0.0174
213.80
0.1702
0.1130
-
5.2489
0.2390
1.69
×
1
0
-
1
1.8606
0.2048
ϕ
f
3.30
×
1
0
-
9
1.35
×
1
0
-
5
1.1445
—
—
—
7.58
×
1
0
-
9
—
—
—
Problem no.
11
12
13
14
15
16
17
18
19
20
J
4.8816
-
0.1273
0.0153
4.5727
9.79
×
1
0
-
4
2.3182
-
9.0882
0.0336
25.13
202.74
ϕ
f
1.67
×
1
0
-
5
4.68
×
1
0
-
9
0.0161
4.33
×
1
0
-
4
5.94
×
1
0
-
3
1.56
×
1
0
-
5
1.64
—
17.06
2.5283


## 8. 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.



---

*Source: 139036-2015-03-05.xml*