# Optimization on Black Box Function Optimization Problem

Authors: Jin-ke Xiao; Wei-min Li; Wei Li; Xin-rong Xiao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

There are a large number of engineering optimization problems in real world, whose input-output relationships are vague and indistinct. Here, they are called black box function optimization problem (BBFOP). Then, inspired by the mechanism of neuroendocrine system regulating immune system, BP neural network modified immune optimization algorithm (NN-MIA) is proposed. NN-MIA consists of two phases: the first phase is training BP neural network with expected precision to confirm input-output relationship and the other phase is immune optimization phase, whose aim is to search global optima. BP neural network fitting without expected fitting precision could be replaced with polynomial fitting or other fitting methods within expected fitting precision. Experimental simulation confirms global optimization capability of MIA and the practical application of BBFOP optimization method.


---

## Body

## 1. Introduction

Black box function optimization problem (BBFOP), inspired from black box concept which originated from cybernetics, is one kind of mathematic optimization functions. Interior function and characteristic of BBFOP could only be understood through exterior observation and experimentation for some condition restriction reasons. Usually, the exterior impacts on BBFOP are called inputs and the function feedback from inputs is called outputs.A large number of complex engineering BBFOPs present ambiguity relationships between inputs and outputs, much optimization computational complexity issue in real manufacturing industry, in aerospace industry, and in computing industry. BBFOP is hard to be solved for the reason that input-output relationships are vague and indistinct. Therefore, finding the global optima or the local optima nearest to the global optima (if we cannot find the global optima) would provide the best solutions to the decision-makers.When we solve BBFOP, two difficulties must be given top priority: one is input-output relationship and the other one is the optimization method. There are two solutions to input-output relationship problem: one is giving BBFOP expression directly through studying interior structure and exploring interior controlling mechanism, which is almost impossible, and the other one is using fitting function as an indirect description of input-output relationship. The latter method enjoys wide application in the engineering practice for fitting function could be born from plentiful input-output sample data. The representative fitting method includes BP neural network fitting and polynomial fitting.For BBFOP optimization method, population-based optimization algorithms enjoy high popularity and have extensive applications in recent years in engineering optimization problems, including genetic algorithms (GAs) [1, 2], ant colony optimization (ACO) [3, 4], particle swarm optimization (PSO) [5, 6], artificial bee colony (ABC) [7, 8], bacterial foraging optimization (BFO) [9, 10], big bang-big crunch algorithm [11, 12], and teaching-learning-based optimization (TLBO) [13–15] and immune optimization algorithm (IA) [16–19]. These population-based optimization algorithms have one or more populations consisting of a certain number of individuals, which presents a solution of the problem to be solved, respectively. Fitness of individuals in population moves towards better one by applying some operators in the search space. Immune optimization algorithm (IA), a vital research and promising direction of Artificial Immune System (AIS), is a novel population-based intelligent algorithm based on the theoretical immune principles [20] inspired by adaptive immune system. Due to the advantages of simplicity and ease of implementation and the remarkable optimization capability, IA and its variants have attracted the attention of researchers and have been widely exploited to solve a great variety of related engineering optimization problems [19, 21–25].Body neural system, immune system, and endocrine system (Immune-Neuro-Endocrine, INE) have super information processing mechanism with characteristic of distribution and robustness. The interaction mechanism among neural system, endocrine system, and immune system gives us some inspirations for intelligent application on optimization [26, 27] while the mechanism where neural system and endocrine system enhance or suppress immune system function may be another useful source of inspirations for population-based algorithm for BBFOP.Here, the single object optimization BBFOP is under research for the reason that the single object optimization problem is the base for multiobject optimization problem. BBFOP optimization method is inspired from the mechanism of how immune system is regulated by neural system and endocrine system. BP neural network or other fitting methods are used to fit input-output sample data with ideal precision while IA and its variants are used to optimize input-output relationship.The remaining text in this paper is organized as follows. Section2 gives a detailed biological mechanism description and its inspiration of neuroendocrine system regulating immune system mechanism; Section 3 develops BBFOP algorithm statement and design; in Section 4, three experiments will be done to attest method validity; in Section 5, conclusions on solving BBFOP and modified optimization method directions are given.


## 2. Neuroendocrine System Regulating Immune System Mechanism

INE system is stereonetwork regulation system with high precision and high complexity through cytokine, hormone, and other chemistry substances [28, 29]. Physical and chemical biological characteristic of medium and hormone secreted by neuroendocrine system and lymph gene, monocyte secreted by immune system, and existing receptor in INE system cell surface are gradually discovered with the molecule-biology development.Immune system could eliminate the invading antigen through the immune response at first once the body is attacked by exterior antigen. If one antibody matches the antigen with higher affinity, it reproduces and eliminates antigen; else, other antibody populations coevolve to produce antibody matching the antigen with high affinity. Meanwhile, the real-time status of immune system is reflected to neuroendocrine system; then neuroendocrine system coregulates immune system to eliminate antigen faster through neurotransmitter and multihormone. Visual concept view is in Figure1.Figure 1
INE system regulation concept view.Neuroendocrine system strengthens immune effect through neurotransmitter and multihormone. This mechanism gives a source inspiration for solving BBFOP, which inspires us to propose the optimization method in Figure2.Figure 2
BBFOP optimization flow.In Figure2, black box function processing unit consists of input-output relationship which is calculated by neural network; optimization unit is the mechanism simulation where neuroendocrine system regulates immune system eliminating invading antigen and searches ideal solution based on input-output relationship.


## 3. Algorithm Statement and Design

### 3.1. Algorithm Formulation

As the minimal optimization problem could be converted into corresponding maximal one, the following nonlinear optimization problem of form in (1) is considered: (1)
max
⁡
f
X
=
f
x
1
,
x
2
,
…
,
x
n
,where f
X is unknown single function optimization BBFOP and also is a time continuous real function with n-dimension argument X
=
x
1
,
x
2
,
…
,
x
n; x
i (i
=
1,2
,
…
,
n) is in the interval l
i
,
h
i (i
=
1,2
,
…
,
n), respectively.Here, a flowchart of NN-MIA is given in Figure3. Figure 3 consists of BP neural network training and fitting and modified immune optimization algorithm (MIA), where MIA has a loop deciding the procedure ends.Figure 3
Algorithm flow.In Box1, Steps  1 and 2 are BP neural network training and fitting and a loop of optimization, respectively. Particularly, in Step  2.5, mutation operator is to be implemented on antibody population C
c
,
P
r instead of antibody population P
r, which could make full use of antibody information before combination; in Step  2.6, selecting the best antibody with the highest affinity from antibody population C
c
,
P
r
,
P
m
,
P
d and adding the best antibody into memory antibody population could make full use of antibody information before combination and mutation.Box 1:Pseudoprocedure flow.
Step 1. BP neural network training and fitting
Step 1.1. draw out mathematic optimization model from engineering problem and produce plenty of sample data
Step 1.2. construct appropriate BP neural network based on the sample data and the model complexity
Step 1.3. initialize BP neural network, train and conserve neural network
Step 1.4. test BP neural network through some sample data to forecast the output, that is affinity in the following MIA
Step 2. modified immune optimization algorithm
Step 2.1. population initialization: initialize antibody population scale P, memory population scale M, clone population scale C
and antibody populationP
P
Step 2.2. affinity calculating: calculate antibody-antigen affinity based on the trained BP neural network
Step 2.3. clone operator: clone antibody population P
P to population scale C and produce antibody population C
c
Step 2.4. combination operator: combine antibody population C
c to produce antibody population P
r
Step 2.5.  mutation operator: mutate antibody population C
c
,
P
r to produce antibody population P
m
Step 2.6. antibody population updating: produce antibody population P
d randomly with population scale d
Step 2.7. add in memory antibody population: select M antibody with the highest affinity from C
c
,
P
r
,
P
m
,
P
d
and add in memory antibody population
Step 2.8. if constrained condition satisfaction is true, skip to Step 8; else skip to Step 2
Step 2.9. select the antibody with highest affinity from the antibody population as the satisfying output result


### 3.2. Algorithm Illustration

#### 3.2.1. BP Neural Network Construction Phase

Input-output layer and hidden network layer of BP neural network should be constructed according to experimental conditions and “2
N
+
1” principle proposed by Hecht-Nielsen [30].


#### 3.2.2. Immune Optimization Phase

(1) Initialization Operator. The algorithm generates populations corresponding to the solutions in the search space. The solutions are produced randomly within the range of the boundaries of the variables in(2)
x
i
,
j
=
x
j
m
i
n
+
rand
⁡
·
x
j
m
a
x
-
x
j
m
i
n
,
i
=
1,2
,
…
,
P
;
j
=
1,2
,
…
,
n
.P is the individual number in population, n is the dimension of optimization problem, representing the number of parameters to be optimized, and x
j
m
i
n and x
j
m
a
x are lower and upper bounds of the jth parameter.(2) Clone Operator. Antibody clone scale is adjusted adaptively and dynamically according to the affinity between antibodies, which is inspired from the inhibitory characteristic mechanism of immune response.Here is the definition of antibody affinityΦ
i between antibody i and the other antibody j in(3)
Φ
i
=
min
⁡
exp
⁡
X
i
-
X
j
=
min
⁡
exp
⁡
∑
i
i
=
1
n
x
i
i
i
-
y
i
i
j
2
,
i
≠
j
,
i
,
j
=
1,2
,
…
,
P
,where X
i
-
X
j is the Euclidean distance between antibody i and antibody j, which is in the interval 0,1 after normalization processing. The higher the affinity between antibody and the other antibody is, the higher their comparability value is and the smaller the value of Φ
i is, which means the stronger inhibitory effect. C
i stands for the clone scale of antibody i and is calculated in(4)
C
i
=
ceil
C
∗
Φ
i
∗
f
X
i
∑
i
=
1
n
f
X
i
,where C is the clone scale initialized before and function c
e
i
l
x stands for the elements of x to the nearest integers which can be seen in MATLAB 7.11.0.(3) Recombination Operator. Suppose X
1 and X
2 are two vectors in n-dimension space. Then, select a dot between two vector regions randomly to produce offspring antibody X
n
e
w, where X
n
e
w
=
a
·
X
1
+
1
-
a
·
X
2, a
=
r
a
n
d.Antibody combination operator executing flow is in Box2.Box 2:Antibody combination operator executing flow.
Step 1. initialization j
:
=
1
Step 2.
Step 2.1. produce random number a
∈
0,1
Step 2.2. produce X
n
e
w
=
a
⋅
X
1
+
1
-
a
⋅
X
2 after combination
Step 3. j
:
=
j
+
1, if termination condition reach, iterative terminates;
else, skip to Step 2.(4) Mutation Operator. Hypermutation, altering on some genes with light probability, leads to affinity maturation. Usually, uniform mutation is used, a mutation combining the result of mutation operator with iterative times. Mutation range is larger relatively at the early stage while the latter is; mutation range is smaller and smaller, which will have fine-tuning impact on evolvement system. Uniform mutation is described as follows.Suppose father antibodyX
=
x
1
,
x
2
,
…
,
x
n
x
k
∈
a
k
,
b
k is selected to mutate, where a
k and b
k are lower and upper bounds of x
k; offspring after mutation is in(5)
X
′
=
x
1
,
…
,
x
k
-
1
,
x
k
′
,
…
,
x
n
,
x
k
′
=
x
k
+
Δ
i
t
,
b
k
-
x
k
,
r
n
d
2
=
0
,
x
k
-
Δ
i
t
,
b
k
-
x
k
,
r
n
d
2
=
1
,where r
n
d
2 is the modulus after division 2, i
t is the current iteration, and range of function Δ
i
t
,
y is 0
,
y. The bigger i
t is, the larger probability that Δ
i
t
,
y is close to 0, which enables the algorithm search across a large scale in the early evolvement stage and search across a local scale in the later evolvement stage. r
∈
0,1 is random number, T is the maximum iteration, and λ
∈
2,5 is the parameter that decides uniform mutation degree, which adjusts the local search region.Main steps of uniform mutation operator are in Box3, where b
u
i
,
k and b
d
i
,
k are lower and upper bounds of kth component in ith antibody, respectively.Box 3:Main steps of uniform mutation operator.
Step 1. initialization i
:
=
1
Step 2. k
=
r
o
u
n
d
r
a
n
d
⋅
n
-
1
+
1, where n is the dimension of decision-making vector
Step 2.1. select kth component in ith antibody A
i
t
i
,
k from antibody population A
i
t
Step 2.2. select positive integer a randomly and calculate its modulus after division 2
Step 2.3. if r
n
d
2
=
=
1
A
i
t
i
,
k
=
A
i
t
i
,
k
+
b
u
i
,
k
-
A
i
t
i
,
k
1
-
r
1
-
i
t
/
T
2
else
A
i
t
i
,
k
=
A
i
t
i
,
k
-
A
i
t
i
,
k
-
b
d
i
,
k
1
-
r
1
-
i
t
/
T
2
Step 3. i
:
=
i
+
1, if termination condition reach, iterative terminates;
else, skip to Step 2.


## 3.1. Algorithm Formulation

As the minimal optimization problem could be converted into corresponding maximal one, the following nonlinear optimization problem of form in (1) is considered: (1)
max
⁡
f
X
=
f
x
1
,
x
2
,
…
,
x
n
,where f
X is unknown single function optimization BBFOP and also is a time continuous real function with n-dimension argument X
=
x
1
,
x
2
,
…
,
x
n; x
i (i
=
1,2
,
…
,
n) is in the interval l
i
,
h
i (i
=
1,2
,
…
,
n), respectively.Here, a flowchart of NN-MIA is given in Figure3. Figure 3 consists of BP neural network training and fitting and modified immune optimization algorithm (MIA), where MIA has a loop deciding the procedure ends.Figure 3
Algorithm flow.In Box1, Steps  1 and 2 are BP neural network training and fitting and a loop of optimization, respectively. Particularly, in Step  2.5, mutation operator is to be implemented on antibody population C
c
,
P
r instead of antibody population P
r, which could make full use of antibody information before combination; in Step  2.6, selecting the best antibody with the highest affinity from antibody population C
c
,
P
r
,
P
m
,
P
d and adding the best antibody into memory antibody population could make full use of antibody information before combination and mutation.Box 1:Pseudoprocedure flow.
Step 1. BP neural network training and fitting
Step 1.1. draw out mathematic optimization model from engineering problem and produce plenty of sample data
Step 1.2. construct appropriate BP neural network based on the sample data and the model complexity
Step 1.3. initialize BP neural network, train and conserve neural network
Step 1.4. test BP neural network through some sample data to forecast the output, that is affinity in the following MIA
Step 2. modified immune optimization algorithm
Step 2.1. population initialization: initialize antibody population scale P, memory population scale M, clone population scale C
and antibody populationP
P
Step 2.2. affinity calculating: calculate antibody-antigen affinity based on the trained BP neural network
Step 2.3. clone operator: clone antibody population P
P to population scale C and produce antibody population C
c
Step 2.4. combination operator: combine antibody population C
c to produce antibody population P
r
Step 2.5.  mutation operator: mutate antibody population C
c
,
P
r to produce antibody population P
m
Step 2.6. antibody population updating: produce antibody population P
d randomly with population scale d
Step 2.7. add in memory antibody population: select M antibody with the highest affinity from C
c
,
P
r
,
P
m
,
P
d
and add in memory antibody population
Step 2.8. if constrained condition satisfaction is true, skip to Step 8; else skip to Step 2
Step 2.9. select the antibody with highest affinity from the antibody population as the satisfying output result


## 3.2. Algorithm Illustration

### 3.2.1. BP Neural Network Construction Phase

Input-output layer and hidden network layer of BP neural network should be constructed according to experimental conditions and “2
N
+
1” principle proposed by Hecht-Nielsen [30].


### 3.2.2. Immune Optimization Phase

(1) Initialization Operator. The algorithm generates populations corresponding to the solutions in the search space. The solutions are produced randomly within the range of the boundaries of the variables in(2)
x
i
,
j
=
x
j
m
i
n
+
rand
⁡
·
x
j
m
a
x
-
x
j
m
i
n
,
i
=
1,2
,
…
,
P
;
j
=
1,2
,
…
,
n
.P is the individual number in population, n is the dimension of optimization problem, representing the number of parameters to be optimized, and x
j
m
i
n and x
j
m
a
x are lower and upper bounds of the jth parameter.(2) Clone Operator. Antibody clone scale is adjusted adaptively and dynamically according to the affinity between antibodies, which is inspired from the inhibitory characteristic mechanism of immune response.Here is the definition of antibody affinityΦ
i between antibody i and the other antibody j in(3)
Φ
i
=
min
⁡
exp
⁡
X
i
-
X
j
=
min
⁡
exp
⁡
∑
i
i
=
1
n
x
i
i
i
-
y
i
i
j
2
,
i
≠
j
,
i
,
j
=
1,2
,
…
,
P
,where X
i
-
X
j is the Euclidean distance between antibody i and antibody j, which is in the interval 0,1 after normalization processing. The higher the affinity between antibody and the other antibody is, the higher their comparability value is and the smaller the value of Φ
i is, which means the stronger inhibitory effect. C
i stands for the clone scale of antibody i and is calculated in(4)
C
i
=
ceil
C
∗
Φ
i
∗
f
X
i
∑
i
=
1
n
f
X
i
,where C is the clone scale initialized before and function c
e
i
l
x stands for the elements of x to the nearest integers which can be seen in MATLAB 7.11.0.(3) Recombination Operator. Suppose X
1 and X
2 are two vectors in n-dimension space. Then, select a dot between two vector regions randomly to produce offspring antibody X
n
e
w, where X
n
e
w
=
a
·
X
1
+
1
-
a
·
X
2, a
=
r
a
n
d.Antibody combination operator executing flow is in Box2.Box 2:Antibody combination operator executing flow.
Step 1. initialization j
:
=
1
Step 2.
Step 2.1. produce random number a
∈
0,1
Step 2.2. produce X
n
e
w
=
a
⋅
X
1
+
1
-
a
⋅
X
2 after combination
Step 3. j
:
=
j
+
1, if termination condition reach, iterative terminates;
else, skip to Step 2.(4) Mutation Operator. Hypermutation, altering on some genes with light probability, leads to affinity maturation. Usually, uniform mutation is used, a mutation combining the result of mutation operator with iterative times. Mutation range is larger relatively at the early stage while the latter is; mutation range is smaller and smaller, which will have fine-tuning impact on evolvement system. Uniform mutation is described as follows.Suppose father antibodyX
=
x
1
,
x
2
,
…
,
x
n
x
k
∈
a
k
,
b
k is selected to mutate, where a
k and b
k are lower and upper bounds of x
k; offspring after mutation is in(5)
X
′
=
x
1
,
…
,
x
k
-
1
,
x
k
′
,
…
,
x
n
,
x
k
′
=
x
k
+
Δ
i
t
,
b
k
-
x
k
,
r
n
d
2
=
0
,
x
k
-
Δ
i
t
,
b
k
-
x
k
,
r
n
d
2
=
1
,where r
n
d
2 is the modulus after division 2, i
t is the current iteration, and range of function Δ
i
t
,
y is 0
,
y. The bigger i
t is, the larger probability that Δ
i
t
,
y is close to 0, which enables the algorithm search across a large scale in the early evolvement stage and search across a local scale in the later evolvement stage. r
∈
0,1 is random number, T is the maximum iteration, and λ
∈
2,5 is the parameter that decides uniform mutation degree, which adjusts the local search region.Main steps of uniform mutation operator are in Box3, where b
u
i
,
k and b
d
i
,
k are lower and upper bounds of kth component in ith antibody, respectively.Box 3:Main steps of uniform mutation operator.
Step 1. initialization i
:
=
1
Step 2. k
=
r
o
u
n
d
r
a
n
d
⋅
n
-
1
+
1, where n is the dimension of decision-making vector
Step 2.1. select kth component in ith antibody A
i
t
i
,
k from antibody population A
i
t
Step 2.2. select positive integer a randomly and calculate its modulus after division 2
Step 2.3. if r
n
d
2
=
=
1
A
i
t
i
,
k
=
A
i
t
i
,
k
+
b
u
i
,
k
-
A
i
t
i
,
k
1
-
r
1
-
i
t
/
T
2
else
A
i
t
i
,
k
=
A
i
t
i
,
k
-
A
i
t
i
,
k
-
b
d
i
,
k
1
-
r
1
-
i
t
/
T
2
Step 3. i
:
=
i
+
1, if termination condition reach, iterative terminates;
else, skip to Step 2.


## 3.2.1. BP Neural Network Construction Phase

Input-output layer and hidden network layer of BP neural network should be constructed according to experimental conditions and “2
N
+
1” principle proposed by Hecht-Nielsen [30].


## 3.2.2. Immune Optimization Phase

(1) Initialization Operator. The algorithm generates populations corresponding to the solutions in the search space. The solutions are produced randomly within the range of the boundaries of the variables in(2)
x
i
,
j
=
x
j
m
i
n
+
rand
⁡
·
x
j
m
a
x
-
x
j
m
i
n
,
i
=
1,2
,
…
,
P
;
j
=
1,2
,
…
,
n
.P is the individual number in population, n is the dimension of optimization problem, representing the number of parameters to be optimized, and x
j
m
i
n and x
j
m
a
x are lower and upper bounds of the jth parameter.(2) Clone Operator. Antibody clone scale is adjusted adaptively and dynamically according to the affinity between antibodies, which is inspired from the inhibitory characteristic mechanism of immune response.Here is the definition of antibody affinityΦ
i between antibody i and the other antibody j in(3)
Φ
i
=
min
⁡
exp
⁡
X
i
-
X
j
=
min
⁡
exp
⁡
∑
i
i
=
1
n
x
i
i
i
-
y
i
i
j
2
,
i
≠
j
,
i
,
j
=
1,2
,
…
,
P
,where X
i
-
X
j is the Euclidean distance between antibody i and antibody j, which is in the interval 0,1 after normalization processing. The higher the affinity between antibody and the other antibody is, the higher their comparability value is and the smaller the value of Φ
i is, which means the stronger inhibitory effect. C
i stands for the clone scale of antibody i and is calculated in(4)
C
i
=
ceil
C
∗
Φ
i
∗
f
X
i
∑
i
=
1
n
f
X
i
,where C is the clone scale initialized before and function c
e
i
l
x stands for the elements of x to the nearest integers which can be seen in MATLAB 7.11.0.(3) Recombination Operator. Suppose X
1 and X
2 are two vectors in n-dimension space. Then, select a dot between two vector regions randomly to produce offspring antibody X
n
e
w, where X
n
e
w
=
a
·
X
1
+
1
-
a
·
X
2, a
=
r
a
n
d.Antibody combination operator executing flow is in Box2.Box 2:Antibody combination operator executing flow.
Step 1. initialization j
:
=
1
Step 2.
Step 2.1. produce random number a
∈
0,1
Step 2.2. produce X
n
e
w
=
a
⋅
X
1
+
1
-
a
⋅
X
2 after combination
Step 3. j
:
=
j
+
1, if termination condition reach, iterative terminates;
else, skip to Step 2.(4) Mutation Operator. Hypermutation, altering on some genes with light probability, leads to affinity maturation. Usually, uniform mutation is used, a mutation combining the result of mutation operator with iterative times. Mutation range is larger relatively at the early stage while the latter is; mutation range is smaller and smaller, which will have fine-tuning impact on evolvement system. Uniform mutation is described as follows.Suppose father antibodyX
=
x
1
,
x
2
,
…
,
x
n
x
k
∈
a
k
,
b
k is selected to mutate, where a
k and b
k are lower and upper bounds of x
k; offspring after mutation is in(5)
X
′
=
x
1
,
…
,
x
k
-
1
,
x
k
′
,
…
,
x
n
,
x
k
′
=
x
k
+
Δ
i
t
,
b
k
-
x
k
,
r
n
d
2
=
0
,
x
k
-
Δ
i
t
,
b
k
-
x
k
,
r
n
d
2
=
1
,where r
n
d
2 is the modulus after division 2, i
t is the current iteration, and range of function Δ
i
t
,
y is 0
,
y. The bigger i
t is, the larger probability that Δ
i
t
,
y is close to 0, which enables the algorithm search across a large scale in the early evolvement stage and search across a local scale in the later evolvement stage. r
∈
0,1 is random number, T is the maximum iteration, and λ
∈
2,5 is the parameter that decides uniform mutation degree, which adjusts the local search region.Main steps of uniform mutation operator are in Box3, where b
u
i
,
k and b
d
i
,
k are lower and upper bounds of kth component in ith antibody, respectively.Box 3:Main steps of uniform mutation operator.
Step 1. initialization i
:
=
1
Step 2. k
=
r
o
u
n
d
r
a
n
d
⋅
n
-
1
+
1, where n is the dimension of decision-making vector
Step 2.1. select kth component in ith antibody A
i
t
i
,
k from antibody population A
i
t
Step 2.2. select positive integer a randomly and calculate its modulus after division 2
Step 2.3. if r
n
d
2
=
=
1
A
i
t
i
,
k
=
A
i
t
i
,
k
+
b
u
i
,
k
-
A
i
t
i
,
k
1
-
r
1
-
i
t
/
T
2
else
A
i
t
i
,
k
=
A
i
t
i
,
k
-
A
i
t
i
,
k
-
b
d
i
,
k
1
-
r
1
-
i
t
/
T
2
Step 3. i
:
=
i
+
1, if termination condition reach, iterative terminates;
else, skip to Step 2.


## 4. Experiments

Experimental procedures and results will be presented and discussed in this section. All algorithms were executed in MATLAB 7.11.0 using computer with Intel Pentium CPU G640, 2.80 GHZ, 2 GB RAM. The operating system of the computer is Windows 7.The global optimization performance of MIA is tested and compared with other optimization algorithms on several well-known benchmark functions including (A) Rosenbrock function, (B) Step function, (C) Quadric function, (D) Schwefel’s function, (E) Rastrigin function, (F) Ackley’s function, (G) Griewank function, and (H) Rotate hyperellipsoid function in the Appendix.In addition, it should be mentioned that, as these population optimization algorithms have stochastic search, which would lead to some changing results in different iterations to some degree, we observe the box plot of optimization function with 30 dimensions (in Section4.1) obtained by each algorithm based on 60 independent runs.Section4.1 reports a comparison between MIA experimental results and three known immune algorithms: SAMCCTLBO, CLONALG, and HPGA; Section 4.2 reports BP neural network-MIA solving BBFO; Section 4.3 reports polynomial fitting-MIA solving BBFO.


### 4.1. Experiment  1: Performance Comparisons of MIA with SAMCCTLBO, CLONALG, and HPGA

In order to examine the proposed MIA, three known algorithms (a multiclass cooperative teaching-learning-based optimization algorithm with simulated annealing (SAMCCTLBO) [31], Clone Selection Algorithm (CLONALG) [32], and Hierarchical Parallel Genetic Algorithms (HPGA) [33]) are selected for direct comparisons.The number of function evaluations (FEs) is becoming a more and more applausive and reasonable measurement termination criterion [34–36] for measuring different algorithms as computational complexity taken in each iteration may exhibit difference in different algorithms. Therefore, combined with search result precision, the termination criterion of the algorithms is one of the objectives, f
b
e
s
t
≤
ε, ε
=
1
0
-
5, where  f
b
e
s
t is the function value obtained by the optimization algorithms or a maximum number of function evaluations (F
E
s) that have been reached, F
e
s
=
10000; here ε
=
1
0
-
5 is used for all functions. If f
b
e
s
t
≤
ε, this algorithm is considered the best result that has been obtained.Box plots [37] are used to illustrate the distribution of these samples obtained from 60 independent runs. The upper and lower ends of the box are the upper and lower quartiles while a line within the box denotes the median and thin appendages summarize the spread and shape of the distribution. The notches represent a robust estimate of the uncertainty about the medians for box-to-box comparison. Symbol “+” denotes outliers. Here, box plots are used to illustrate the distribution of the function values obtained from 60 independent runs.Detailed parameter settings of SAMCCTLBO, CLONALG, and HPGA can be seen in the respective literature and in MIA. The higher the parameters are, the more excellent the results are, the more complex the algorithm’s space and time complexity are, and the longer the running time is. Therefore, parameters are set toP
=
100, M
=
50, C
=
200, r
=
3, λ
=
3, d
=
10, and T
=
500 according to several experiments’ testing. Once they have been all experimentally determined, they are kept the same through all experiments except that MIA parameters are set again. Then, for a summarized performance comparison for different algorithms solving unimodal function, unrotated multimodal function problems are presented in Figure 4.Figure 4
Performance comparison among optimization algorithms for the test functions.
(a)
(b)
(c)
(d)Then, NN-MIA formulation pseudoprocedure is in Box1 correspondingly.From the results in Figure4 box plot, the upper quartiles, lower quartiles, and medians obtained by MIA are much less than those obtained by CLONALG and HPGA and, in addition, are much less or not much more than those obtained by SAMCCTLBO. Furthermore, for (A) Rosenbrock function, (B) Step function, (E) Rastrigin function, (F) Ackley’s function, and (G) Griewank function, MIA and SAMCCTLBO, CLONALG, and HPGA all can obtain global optimal 0 in all the 60 runs with accepted standard deviation, but MIA generates the best results while SAMCCTLBO performs better than CLONALG and HPGA. For (C) Quadric function and (D) Schwefel’s function, MIA and SAMCCTLBO can obtain global optimal 0 in all the 60 runs, and CLONALG and HPGA could obtain global results near to global optima and beyond the results precision, while MIA obtains global optima with smaller standard deviation than that of SAMCCTLBO. For (H) Rotate hyperellipsoid function, only MIA can obtain the best optima near the results precision. Here, it can be easily deduced from the box plots that MIA has a much higher capability of finding the global optima.


### 4.2. Experiment  2: BP Neural Network-MIA Solving BBFO

Here, NN-MIA is used to optimize the following nonline function in(6)
max
⁡
f
x
=
1
1
+
∑
l
=
1
n
x
l
2
,
-
5
≤
x
l
≤
5
,
∀
l
=
1,2
,
…
,
n
,where f
x
∗
=
1 with x
1
⋯
x
n
=
0
⋯
0.The function graphic is in Figure5 when n
=
2.Figure 5
Function graphic whenn
=
2.It is easy to conclude that the function global maximum is 1 and its corresponding coordinate is0,0 from the above function expression in (6) and its function graphic. Function extremum and its corresponding coordinate are easy to be concluded in a known instead of an unknown function expression case.Supposen
=
8 in (6) which means the corresponding black box function has 8 input parameters and 1 output parameter. Then, BP neural network is constructed: suppose BP neural network has single hidden network, the neuron number is set to be 17 according to “2
N
+
1” method, and transferring function is Sigmoid function in every neuron. Finally, topology structure 8-17-1 BP neural network with 8 inputs, 17 hidden layers, and 1 output is demonstrated in Figure 6.Figure 6
BP neural network designed to model in experimentation 1.Here, 8000 inputs and outputs are produced according to (6) and 7000 pieces of data are selected randomly to train BP neural network; the other 1000 pieces of data are used to test the fitting accuracy of BP neural network. Then, comparison of forecasting outputs and expected outputs is demonstrated in Figure 7. The accumulative error between forecasting outputs and expected outputs is limited within 1.1495, which is within the error requirement. Therefore, the input-output relationship could be solved by topology structure 8-17-1 BP neural network. Then, the following is optimizing results.Figure 7
Comparison of forecasting output and expected output.The efficiency of real number code is much higher than that of binary code usually; therefore, real coding is applied in MIA and antibody fitness is the BP neural network output; the larger the fitness is, the super the antibody is. Then, the best fitness is0.999794 and its best antibody individual is [
0.0002
,
-
0.00009
,
-
0.0003,0.0001,0.0003, -
0.0004
,
-
0.0003
,
-
0.0009
], which is very close to real maximum 1 and its corresponding coordinate 0,0. The method is confirmed to be valid.


### 4.3. Experiment  3: Polynomial Fitting-MIA Solving BBFO

However, in experiment  2, this method has to be modified to some extent. Here, modified Rosenbrock function is in (7) and the function graphic is in Figure 8 when n
=
2:(7)
max
⁡
f
x
=
1
1
+
∑
l
=
1
n
-
1
100
x
l
+
1
-
x
l
2
2
+
x
l
-
1
2
,
-
5
≤
x
l
≤
5
,
∀
l
=
1,2
,
…
,
n
,where f
x
∗
=
1 with x
1
⋯
x
n
=
1
⋯
1
.Figure 8
Function graphic whenn
=
2.The same BP neural network in experiment is used here and the comparison of forecasting output and expected output is demonstrated in Figure9. The accumulative error between forecasting output and expected output reaches 7.9212, which has gone beyond error requirement.Figure 9
Comparison of forecasting output and expected output.In this case, the fitness results of MIA will not be satisfied. Therefore, polynomial fitting method is used instead of BP neural network and the comparison of forecasting output and expected output is demonstrated in Figure10. The accumulative error between forecasting output and expected output is limited within 1.2032, which accords with error requirement.Figure 10
Comparison of forecasting output and expected output.Here, MIA is used for optimization and the best fitness0.9999981 is produced, which is very close to real maximum 1. The method of combining polynomial fitting with MIA is confirmed to be valid here.


## 4.1. Experiment  1: Performance Comparisons of MIA with SAMCCTLBO, CLONALG, and HPGA

In order to examine the proposed MIA, three known algorithms (a multiclass cooperative teaching-learning-based optimization algorithm with simulated annealing (SAMCCTLBO) [31], Clone Selection Algorithm (CLONALG) [32], and Hierarchical Parallel Genetic Algorithms (HPGA) [33]) are selected for direct comparisons.The number of function evaluations (FEs) is becoming a more and more applausive and reasonable measurement termination criterion [34–36] for measuring different algorithms as computational complexity taken in each iteration may exhibit difference in different algorithms. Therefore, combined with search result precision, the termination criterion of the algorithms is one of the objectives, f
b
e
s
t
≤
ε, ε
=
1
0
-
5, where  f
b
e
s
t is the function value obtained by the optimization algorithms or a maximum number of function evaluations (F
E
s) that have been reached, F
e
s
=
10000; here ε
=
1
0
-
5 is used for all functions. If f
b
e
s
t
≤
ε, this algorithm is considered the best result that has been obtained.Box plots [37] are used to illustrate the distribution of these samples obtained from 60 independent runs. The upper and lower ends of the box are the upper and lower quartiles while a line within the box denotes the median and thin appendages summarize the spread and shape of the distribution. The notches represent a robust estimate of the uncertainty about the medians for box-to-box comparison. Symbol “+” denotes outliers. Here, box plots are used to illustrate the distribution of the function values obtained from 60 independent runs.Detailed parameter settings of SAMCCTLBO, CLONALG, and HPGA can be seen in the respective literature and in MIA. The higher the parameters are, the more excellent the results are, the more complex the algorithm’s space and time complexity are, and the longer the running time is. Therefore, parameters are set toP
=
100, M
=
50, C
=
200, r
=
3, λ
=
3, d
=
10, and T
=
500 according to several experiments’ testing. Once they have been all experimentally determined, they are kept the same through all experiments except that MIA parameters are set again. Then, for a summarized performance comparison for different algorithms solving unimodal function, unrotated multimodal function problems are presented in Figure 4.Figure 4
Performance comparison among optimization algorithms for the test functions.
(a)
(b)
(c)
(d)Then, NN-MIA formulation pseudoprocedure is in Box1 correspondingly.From the results in Figure4 box plot, the upper quartiles, lower quartiles, and medians obtained by MIA are much less than those obtained by CLONALG and HPGA and, in addition, are much less or not much more than those obtained by SAMCCTLBO. Furthermore, for (A) Rosenbrock function, (B) Step function, (E) Rastrigin function, (F) Ackley’s function, and (G) Griewank function, MIA and SAMCCTLBO, CLONALG, and HPGA all can obtain global optimal 0 in all the 60 runs with accepted standard deviation, but MIA generates the best results while SAMCCTLBO performs better than CLONALG and HPGA. For (C) Quadric function and (D) Schwefel’s function, MIA and SAMCCTLBO can obtain global optimal 0 in all the 60 runs, and CLONALG and HPGA could obtain global results near to global optima and beyond the results precision, while MIA obtains global optima with smaller standard deviation than that of SAMCCTLBO. For (H) Rotate hyperellipsoid function, only MIA can obtain the best optima near the results precision. Here, it can be easily deduced from the box plots that MIA has a much higher capability of finding the global optima.


## 4.2. Experiment  2: BP Neural Network-MIA Solving BBFO

Here, NN-MIA is used to optimize the following nonline function in(6)
max
⁡
f
x
=
1
1
+
∑
l
=
1
n
x
l
2
,
-
5
≤
x
l
≤
5
,
∀
l
=
1,2
,
…
,
n
,where f
x
∗
=
1 with x
1
⋯
x
n
=
0
⋯
0.The function graphic is in Figure5 when n
=
2.Figure 5
Function graphic whenn
=
2.It is easy to conclude that the function global maximum is 1 and its corresponding coordinate is0,0 from the above function expression in (6) and its function graphic. Function extremum and its corresponding coordinate are easy to be concluded in a known instead of an unknown function expression case.Supposen
=
8 in (6) which means the corresponding black box function has 8 input parameters and 1 output parameter. Then, BP neural network is constructed: suppose BP neural network has single hidden network, the neuron number is set to be 17 according to “2
N
+
1” method, and transferring function is Sigmoid function in every neuron. Finally, topology structure 8-17-1 BP neural network with 8 inputs, 17 hidden layers, and 1 output is demonstrated in Figure 6.Figure 6
BP neural network designed to model in experimentation 1.Here, 8000 inputs and outputs are produced according to (6) and 7000 pieces of data are selected randomly to train BP neural network; the other 1000 pieces of data are used to test the fitting accuracy of BP neural network. Then, comparison of forecasting outputs and expected outputs is demonstrated in Figure 7. The accumulative error between forecasting outputs and expected outputs is limited within 1.1495, which is within the error requirement. Therefore, the input-output relationship could be solved by topology structure 8-17-1 BP neural network. Then, the following is optimizing results.Figure 7
Comparison of forecasting output and expected output.The efficiency of real number code is much higher than that of binary code usually; therefore, real coding is applied in MIA and antibody fitness is the BP neural network output; the larger the fitness is, the super the antibody is. Then, the best fitness is0.999794 and its best antibody individual is [
0.0002
,
-
0.00009
,
-
0.0003,0.0001,0.0003, -
0.0004
,
-
0.0003
,
-
0.0009
], which is very close to real maximum 1 and its corresponding coordinate 0,0. The method is confirmed to be valid.


## 4.3. Experiment  3: Polynomial Fitting-MIA Solving BBFO

However, in experiment  2, this method has to be modified to some extent. Here, modified Rosenbrock function is in (7) and the function graphic is in Figure 8 when n
=
2:(7)
max
⁡
f
x
=
1
1
+
∑
l
=
1
n
-
1
100
x
l
+
1
-
x
l
2
2
+
x
l
-
1
2
,
-
5
≤
x
l
≤
5
,
∀
l
=
1,2
,
…
,
n
,where f
x
∗
=
1 with x
1
⋯
x
n
=
1
⋯
1
.Figure 8
Function graphic whenn
=
2.The same BP neural network in experiment is used here and the comparison of forecasting output and expected output is demonstrated in Figure9. The accumulative error between forecasting output and expected output reaches 7.9212, which has gone beyond error requirement.Figure 9
Comparison of forecasting output and expected output.In this case, the fitness results of MIA will not be satisfied. Therefore, polynomial fitting method is used instead of BP neural network and the comparison of forecasting output and expected output is demonstrated in Figure10. The accumulative error between forecasting output and expected output is limited within 1.2032, which accords with error requirement.Figure 10
Comparison of forecasting output and expected output.Here, MIA is used for optimization and the best fitness0.9999981 is produced, which is very close to real maximum 1. The method of combining polynomial fitting with MIA is confirmed to be valid here.


## 5. Conclusion

BBFOP optimization algorithm is proposed in this paper, which is inspired by the mechanism of neuroendocrine system regulating immune system. In this algorithm, BP neural network is used to fit this input-output relationship based on sample data. If the fitness precision is not reached, the polynomial fitting and other fitting methods are adopted, and then MIA are to optimize the fitting function.(1) Engineering Value of Black Box Function and Its Optimization Method. BBFOP optimizing method enjoys important applications in some kinds of engineering problem. For example, one experiment aim is to conclude the experiment conditions according to the optimal experiment results; however, such experiment could only repeat numbered times and experiment conditions could not be reached just based on experiment results. Similar engineering problems could be considered as BBFOP.(2) Function Fitting Methods Analysis on BBFOP. The forecasting result precision of BP neural network determines final optimizing results to some extent. With higher forecasting result precision, the final optimizing results are closer to the real results, which requires more training sample data. However, not all black box function inputs and outputs could be expressed by BP neural network with expected precision. Therefore, other function fitting methods including polynomial fitting method could be applied to meet expected precision and the following steps are almost the same as the above.In order to promote the development on the theory and application of BBFOP and its optimization method, our future research will focus on basic theory of BBFOP: biology theory including neural system and immune system should be explored at length to enhance neural network fitting capability and MIA capability; therefore, deeper researches should be carried out by means of theoretical analysis and related experiments.



---

*Source: 647234-2015-10-29.xml*