# Optimal Combination of Aircraft Maintenance Tasks by a Novel Simplex Optimization Method

Authors: Huaiyuan Li; Hongfu Zuo; Dan Lei; Kun Liang; Tingting Lu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Combining maintenance tasks into work packages is not only necessary for arranging maintenance activities, but also critical for the reduction of maintenance cost. In order to optimize the combination of maintenance tasks by fuzzy C-means clustering algorithm, an improved fuzzy C-means clustering model is introduced in this paper. In order to reduce the dimension, variables representing clustering centers are eliminated in the improved cluster model. So the improved clustering model can be directly solved by the optimization method. To optimize the clustering model, a novel nonlinear simplex optimization method is also proposed in this paper. The novel method searches along all rays emitting from the center to each vertex, and those search directions are rightlyn+1 positive basis. The algorithm has both theoretical convergence and good experimental effect. Taking the optimal combination of some maintenance tasks of a certain aircraft as an instance, the novel simplex optimization method and the clustering model both exhibit excellent performance.


---

## Body

## 1. Introduction

Combining maintenance tasks into work packages is not only necessary for arranging maintenance activities, but also critical for the reduction of maintenance cost. In the field of civil aircraft maintenance, in order to reduce the maintenance cost, the number of character check is constantly increased. For example, the number ofA check is increased from 64 to 128. Recently, character checks have been completely cancelled, instead, the maintenance tasks executed in character check are combined into several small maintenance work packages and are equally allocated to line maintenance, so as to balance the maintenance and reduce the maintenance cost. However, it is still inevitable to optimally combine maintenance tasks into work packages no matter whether character check is utilized or not. As combinatorial optimizations usually belong to NP-hard problems that are difficult to be exactly solved, only optimal approximation solutions are obtained on the basis of fuzzy clustering of maintenance tasks in the paper. First of all, a modification is made on the fuzzy C-means clustering model to reduce the dimension. Then, a novel nonlinear simplex algorithm is put forward. Lastly, the optimal combination of maintenance tasks is transformed into a clustering problem which can be solved by a novel simplex algorithm.In the traditional fuzzy C-means (FCM) clustering method, not only the degree of fuzzy membership is taken as an independent variable but also the center of every class is unknown. For this reason, the FCM model is a high-dimension optimization problem that is difficult to be directly solved. Usually, only its local optimal solution can be found out by the alternative optimization algorithm. Therefore, the traditional FCM may present high error rate and poor adaptability. In this paper, class centers are regarded as a function of fuzzy membership degree and samples; thus, an improved FCM model is developed. The improved model does not explicitly include clustering centers, so the number of unknown variables is reduced, which not only results in the decrease in the computational cost but also makes it possible to directly solve the clustering model by the optimization method. According to the experimental results, if we directly solve the improved FCM model, the results of the improved FCM model are better than those of traditional FCM model in most cases.As a direct optimization method, the traditional nonlinear simplex method mainly has the following serious defects: poor efficiency in solving high dimension problems and rather slow calculating speed [1]. In high-dimension problems, a sufficient number of vertexes are needed to form a simplex, and the descending directions are more complex. However, the traditional simplex optimization method only searches in a single one direction and merely updates one vertex at each iteration, when it optimizes by reflection, expansion and contraction [2]. Even though the new vertex is much better than the former one, only a small number of vertexes are updated. Thus, the number of vertexes updated at each iteration is too small, which is not enough for the rapid revolution of the entire simplex. Namely, update of one vertex hardly contributes to the evolution of high-dimension simplex comprised with a large amount of vertexes. Hence, in the paper, center ray operation is introduced, which can search along all rays emitting from the center to each vertex. By the center ray operation, more search directions are available and most of vertexes in the simplex can be updated, which will greatly accelerate the optimization speed and improve capacity of solving high-dimension problems. Experiments show that the slow computing speed of the traditional simplex method is mainly caused by several reasons: (1) too many shrinking operations consume great computation cost [3]; (2) fixed search speed and step can not adapt to different conditions [4]; (3) the search direction is not necessarily accurate descent direction, which results in lower solving efficiencies [5]. In view of the above mentioned reasons, shrinkage operation is removed in the paper. In case of slow evolution of simplex, the novel simplex method regenerates a new small simplex according to the current estimated descent direction. This not only accelerates the optimization speed and reduces computational expense but also avoids the degeneration problem that easily occurs in the evolution course of traditional simplex method.The paper consists of three aspects. Firstly, efficiency of the traditional simplex method is low, especially when it is utilized to optimize high-dimension problems [6]. Accordingly, in the paper, a novel nonlinear simplex optimization method is primarily put forward, which can not only search along all rays emitting from the center to each vertex but also simultaneously update numerous vertexes. Secondly, concerning the problems of too many variables and large calculated amount in the traditional FCM, the FCM model is improved in the paper. Finally, optimal combination of maintenance tasks is concluded as a problem of constrained clustering, the novel simplex method is employed to solve the improved FCM model so as to achieve the optimal combination of maintenance tasks.


## 2. Improvement of Fuzzy C-Means Clustering

### 2.1. Introduction to Traditional FCM

Clustering means adaptively classifying the sample settj,j=0,…,n into c classes according to the similarity among samples. One of the most common methods of clustering is fuzzy C-means clustering [7]. If uij∈0,1 indicates fuzzy membership that the sample tj affiliates the ith class and ci is the center of the sample of the ith class, FCM takes dij2=tj-ci2 as quantitative index to measure the difference or dissimilarity between tj and ci. Assuming the fuzzy control parameter m>1, the FCM model is as (1)min⁡Jm=∑i=1c∑j=1nuijmdij2s.t.∑i=1cuij=1,∀j=1,…,nuij≥0,∀i=1,…,c,j=1,…,n.The key to the traditional FCM algorithm is to solve the class center matrixcic×d and the fuzzy membership matrix uijc×n. However, clustering is a NP problem, and (1) is usually a high-dimension optimization problem with numerous local minimum values [8]. Traditionally, Model (1) is alternatively solved by the steepest descent method. Firstly, assuming that the current clustering center ci is known, and the stationary point of Model (1) with respect to uij is taken as the next uij. Then, the new uij is supposed to be constant, and the optimum ci is achieved with the stationary point of Model (1) with respect to ci, as shown in (2)ci=∑j=1nuijm∑k=1nuikmxj.In the whole course of computation, the unknownci,uij,i=1,…,c, j=1,…,n are not simultaneously computed as a whole, and also, the correlation between ci and uij is ignored. Therefore, only the local minimum point can be found out. The traditional FCM model must repeatedly compute the distance between each sample and all centers of class at each iteration, which consumes a large amount of computation.


### 2.2. Improved Model of FCM

Actually, there is a close correlation betweenci and uij. According to (2), if the class center ci is regarded as the function of uij and xj, and then ci can be eliminated and the number of variables will be reduced. From the property of norm and inner product, dij2 can be denoted as dij2=tj,tj-2∑l=1nuilmtj,tl/∑k=1nuikm+∑p=1n∑q=1nuipmuiqmtp,tq/∑k=1nuikm2. As such, in the paper, (2) is substituted into Model (1) directly and the FCM model is obtained as (3), whose variables only contain uij not ci:(3)min⁡Jm=∑i=1c∑j=1nuijmdij2=∑i=1c∑j=1nuijmtj,tj-2∑l=1nuilmtj,tl∑k=1nuikm+∑p=1n∑q=1nuipmuiqmtp,tq∑k=1nuikm2s.t.∑i=1cuij=1,∀j=1,…,nuij≥0,∀i=1,…,c,j=1,…,n.Different from (1), (3) does not contain class center ci. Instead, it is composed of inner products among samples. If they can be computed and stored in advance, the computed amount of Model (3) can be reduced especially when sample size is large. Compared with (1), (3) contains much less unknown variables, which makes it possible to directly solve (3). As (3) has numerical stationary points with respect to uij, the optimization method which needs derivative is inappropriate for (3). Inspired by the traditional simplex algorithm, a novel direct optimization method that can be utilized to solve nonlinear programming, namely, a novel simplex optimization method based on center ray, is presented to solve (3) in the paper.


### 2.3. The Verification of Improved Model of FCM

To solve the constrained optimization Model (3), (3) is firstly transformed into the unconstrained Model (4) by the exterior penalty function method. Then, the novel optimization method proposed in the paper is employed to solve the unconstrained problem, namely, (4):(4)min⁡Jm=∑i=1c∑j=1nuijmtj,tj-2∑l=1nuilmtj,tl∑k=1nuikm+∑p=1n∑q=1nuipmuiqmtp,tq∑k=1nuikm2+δk∑i=1c∑j=1nmax2⁡0,-uij+∑j=1n∑i=1cuij-1.Ifci is the center of the ith class of simulation sample and rl is the random vector obeying certain distribution, then Ci=ci+rl,l=1,…,Ci can be regarded as the set of the ith class simulation sample. In the experiment, rl, respectively, follows normal distribution and uniform distribution, in other words, simulate, respectively, under rl~N0,(ci+1-ci)/22 and rl~U-((ci+1-ci)/2),(ci+1-ci)/2 circumstance. Three classes of eight-dimensional random vectors are randomly generated, which comprise sample set C1,C2,C3, where, respectively, set c1=1, c2=2, and c3=3. Before the sample is generated, the number of every class sample is randomly determined by Ci=κi+ϑri-0.5, where, ϑ=5, ri~U0,1 and κ1=15, κ2=20, κ3=25.Generally, the indexVpe=-(1/n)∑i=1c∑j=1nuijlog⁡uij, Vmpc=1-(c/(c-1))1-(1/n)∑i=1c∑j=1nuij2, Vxb=∑i=1c∑j=1nuijmtj-ci2/(nmini≠j⁡ci-cj2), Vk=(∑i=1c∑j=1nuijmtj-ci2+(1/c)∑i=1cci-c-2)/mini≠j⁡ci-cj2, and the error rate of classification is utilized to evaluate the effect of the classification and compare the performance of the clustering model. Model (4) is solved by Algorithm 5. In the experiment, numerical simulation is repeated 300 times under the condition of same experimental parameters and some average results are totaled, which are tabulated in Table 1. Jm in Table 1 is the practical optimum value of the clustering object function, which can indicate whether the clustering model has advantages and disadvantages. In the experiment, the main parameters of Algorithm 5 are as follows: λ⃛h=0.985, tmax=10, ρnei=0.1, hsmall=0.1, h0=0.5, αinitforward=2.2, αendforward=0.75, αinitreverse=-2.2, αendreverse=-1.2, μinner=0.5, μouter=1.3, αfe=αre=μe=1.35, αfc=μoc=0.7, μic=1.5, αs=μs=-0.5, λ˙h=3, λ¨h=0.95, αic=0.5, αfs=-0.4, αrs=0.4, and μs=0.35.Table 1
Comparison of clustering results.
Distribution
Model
V
p
e
V
m
p
c
V
x
b
V
k
Error rate
J
m
Uniform
Model (1)
0.050742
0.974330
0.079669
0.095152
0.007974
27.780529
Uniform
Model (3)
0.372381
0.729550
0.073406
0.088955
0.000000
25.252153
Norm
Model (1)
0.224796
0.836828
0.209994
0.225994
0.005346
72.118062
Norm
Model (3)
0.638111
0.467340
0.203394
0.220157
0.004615
61.697758After analyzing the results in Table1, the following conclusion is obtained. Jm in Model (3) is always smaller than that in Model (1), which implies the result of Model (3) is better than that of Model (1) from the view of mathematical optimization of object function. Vxb,Vk of Model (3) are also better than those of Model (1), which shows that, after clustering of Model (3), the interior of each class becomes more compact and the separability and dissimilarity between classes are more significant apparent. Besides, in terms of the classification error rate, Model (3) is also more accurate than Model (1). However, if Vpe and Vmpc are taken into consideration, Model (1) seems to be better, indicating that making-decisions of Model (1) are clearer and more arbitrary while the membership degree of Model (3) is more uniform than that of Model (1). Thus, compared with traditional FCM, Model (3) is stronger in adaptability and higher in recognition accuracy.


## 2.1. Introduction to Traditional FCM

Clustering means adaptively classifying the sample settj,j=0,…,n into c classes according to the similarity among samples. One of the most common methods of clustering is fuzzy C-means clustering [7]. If uij∈0,1 indicates fuzzy membership that the sample tj affiliates the ith class and ci is the center of the sample of the ith class, FCM takes dij2=tj-ci2 as quantitative index to measure the difference or dissimilarity between tj and ci. Assuming the fuzzy control parameter m>1, the FCM model is as (1)min⁡Jm=∑i=1c∑j=1nuijmdij2s.t.∑i=1cuij=1,∀j=1,…,nuij≥0,∀i=1,…,c,j=1,…,n.The key to the traditional FCM algorithm is to solve the class center matrixcic×d and the fuzzy membership matrix uijc×n. However, clustering is a NP problem, and (1) is usually a high-dimension optimization problem with numerous local minimum values [8]. Traditionally, Model (1) is alternatively solved by the steepest descent method. Firstly, assuming that the current clustering center ci is known, and the stationary point of Model (1) with respect to uij is taken as the next uij. Then, the new uij is supposed to be constant, and the optimum ci is achieved with the stationary point of Model (1) with respect to ci, as shown in (2)ci=∑j=1nuijm∑k=1nuikmxj.In the whole course of computation, the unknownci,uij,i=1,…,c, j=1,…,n are not simultaneously computed as a whole, and also, the correlation between ci and uij is ignored. Therefore, only the local minimum point can be found out. The traditional FCM model must repeatedly compute the distance between each sample and all centers of class at each iteration, which consumes a large amount of computation.


## 2.2. Improved Model of FCM

Actually, there is a close correlation betweenci and uij. According to (2), if the class center ci is regarded as the function of uij and xj, and then ci can be eliminated and the number of variables will be reduced. From the property of norm and inner product, dij2 can be denoted as dij2=tj,tj-2∑l=1nuilmtj,tl/∑k=1nuikm+∑p=1n∑q=1nuipmuiqmtp,tq/∑k=1nuikm2. As such, in the paper, (2) is substituted into Model (1) directly and the FCM model is obtained as (3), whose variables only contain uij not ci:(3)min⁡Jm=∑i=1c∑j=1nuijmdij2=∑i=1c∑j=1nuijmtj,tj-2∑l=1nuilmtj,tl∑k=1nuikm+∑p=1n∑q=1nuipmuiqmtp,tq∑k=1nuikm2s.t.∑i=1cuij=1,∀j=1,…,nuij≥0,∀i=1,…,c,j=1,…,n.Different from (1), (3) does not contain class center ci. Instead, it is composed of inner products among samples. If they can be computed and stored in advance, the computed amount of Model (3) can be reduced especially when sample size is large. Compared with (1), (3) contains much less unknown variables, which makes it possible to directly solve (3). As (3) has numerical stationary points with respect to uij, the optimization method which needs derivative is inappropriate for (3). Inspired by the traditional simplex algorithm, a novel direct optimization method that can be utilized to solve nonlinear programming, namely, a novel simplex optimization method based on center ray, is presented to solve (3) in the paper.


## 2.3. The Verification of Improved Model of FCM

To solve the constrained optimization Model (3), (3) is firstly transformed into the unconstrained Model (4) by the exterior penalty function method. Then, the novel optimization method proposed in the paper is employed to solve the unconstrained problem, namely, (4):(4)min⁡Jm=∑i=1c∑j=1nuijmtj,tj-2∑l=1nuilmtj,tl∑k=1nuikm+∑p=1n∑q=1nuipmuiqmtp,tq∑k=1nuikm2+δk∑i=1c∑j=1nmax2⁡0,-uij+∑j=1n∑i=1cuij-1.Ifci is the center of the ith class of simulation sample and rl is the random vector obeying certain distribution, then Ci=ci+rl,l=1,…,Ci can be regarded as the set of the ith class simulation sample. In the experiment, rl, respectively, follows normal distribution and uniform distribution, in other words, simulate, respectively, under rl~N0,(ci+1-ci)/22 and rl~U-((ci+1-ci)/2),(ci+1-ci)/2 circumstance. Three classes of eight-dimensional random vectors are randomly generated, which comprise sample set C1,C2,C3, where, respectively, set c1=1, c2=2, and c3=3. Before the sample is generated, the number of every class sample is randomly determined by Ci=κi+ϑri-0.5, where, ϑ=5, ri~U0,1 and κ1=15, κ2=20, κ3=25.Generally, the indexVpe=-(1/n)∑i=1c∑j=1nuijlog⁡uij, Vmpc=1-(c/(c-1))1-(1/n)∑i=1c∑j=1nuij2, Vxb=∑i=1c∑j=1nuijmtj-ci2/(nmini≠j⁡ci-cj2), Vk=(∑i=1c∑j=1nuijmtj-ci2+(1/c)∑i=1cci-c-2)/mini≠j⁡ci-cj2, and the error rate of classification is utilized to evaluate the effect of the classification and compare the performance of the clustering model. Model (4) is solved by Algorithm 5. In the experiment, numerical simulation is repeated 300 times under the condition of same experimental parameters and some average results are totaled, which are tabulated in Table 1. Jm in Table 1 is the practical optimum value of the clustering object function, which can indicate whether the clustering model has advantages and disadvantages. In the experiment, the main parameters of Algorithm 5 are as follows: λ⃛h=0.985, tmax=10, ρnei=0.1, hsmall=0.1, h0=0.5, αinitforward=2.2, αendforward=0.75, αinitreverse=-2.2, αendreverse=-1.2, μinner=0.5, μouter=1.3, αfe=αre=μe=1.35, αfc=μoc=0.7, μic=1.5, αs=μs=-0.5, λ˙h=3, λ¨h=0.95, αic=0.5, αfs=-0.4, αrs=0.4, and μs=0.35.Table 1
Comparison of clustering results.
Distribution
Model
V
p
e
V
m
p
c
V
x
b
V
k
Error rate
J
m
Uniform
Model (1)
0.050742
0.974330
0.079669
0.095152
0.007974
27.780529
Uniform
Model (3)
0.372381
0.729550
0.073406
0.088955
0.000000
25.252153
Norm
Model (1)
0.224796
0.836828
0.209994
0.225994
0.005346
72.118062
Norm
Model (3)
0.638111
0.467340
0.203394
0.220157
0.004615
61.697758After analyzing the results in Table1, the following conclusion is obtained. Jm in Model (3) is always smaller than that in Model (1), which implies the result of Model (3) is better than that of Model (1) from the view of mathematical optimization of object function. Vxb,Vk of Model (3) are also better than those of Model (1), which shows that, after clustering of Model (3), the interior of each class becomes more compact and the separability and dissimilarity between classes are more significant apparent. Besides, in terms of the classification error rate, Model (3) is also more accurate than Model (1). However, if Vpe and Vmpc are taken into consideration, Model (1) seems to be better, indicating that making-decisions of Model (1) are clearer and more arbitrary while the membership degree of Model (3) is more uniform than that of Model (1). Thus, compared with traditional FCM, Model (3) is stronger in adaptability and higher in recognition accuracy.


## 3. A Novel Simplex Optimization Algorithm

### 3.1. Description of the Novel Simplex Method

The convex hull, namely, the hyper polyhedron composed ofn+1 vertexes in Rn is called a simplex. Assuming vertexes x0,x1,…,xn in the simplex are sorted in ascending order by function values, the simplex of the kth iteration is denoted as Sk=x0,…,xn. Obviously, n side vectors xi-xj, j=0,…,n∧j≠i, connecting with the ith vertex xi are linearly independent, from which the degeneration of the simplex can be judged. The whole iteration of the traditional nonlinear simplex algorithm mainly consists of four steps: reflection, expansion, contraction, and shrinkage. In contrast, the novel algorithm in the paper primarily solves problems by generating a simplex, raying outward from the simplex center and reflecting operations. Taking the unconstrained optimization problem min⁡fx as an example, the novel simplex algorithm can be described as below.


#### 3.1.1. Generation of a Simplex

Due to the tendency of deformation and degeneration in the process of the simplex’s iteration, it is usually better that the initial simplex is constructed as an orthogonal one. When one side of the orthogonal simplex gets close to the negative gradient direction-ρk, the simplex will search much rapidly. If Ai=τ1,…,τi and z is any of solution of the linear equations AiTZ=0, then z will be orthogonal to τ1,…,τi. Firstly, we set τ1=-ρk and A1=τ1. If z, that is, any of the solutions of AiTZ=0, is taken as the last column of Ai to generate a new matrix Ai+1, then the orthogonal matrix An=τ1,…,τn can be obtained after repeatedly solving the linear equations n-1 times, where, τ1,…,τn represent n orthogonal vectors. Subsequently, we choose the point whose function value is better from zk+hmτi,zk-hmτi as the vertex of the simplex. In this paper, the flow of simplex generation is shown as Algorithm 1.Algorithm 1.
Step  1. Let τ1=-ρk/ρk, i=1, Sk=zk, A1=τ1.
Step  2. Let(5)Sk=Sk∪zk+hmτifzk+hmτi<fzk-hmτiSk∪zk-hmτiothers.
Step  3. Solve the linear equations AiTZ=0 and choose any of the solutions as z and let (6)τi+1=zz,Ai+1=τ1,…,τi,τi+1.
Step  4. When i<n-1, then i=i+1 and go to Step  2; otherwise exit.Following the procedure, the simplexSk is generated at zk. Edges connecting with zk are orthogonal vectors whose lengths are all hm. Besides, the first edge is parallel to the negative gradient direction -ρk.


#### 3.1.2. Center Ray Operation

Theorem 2.
If the center of the simplex is denoted asx-=(1/(n+1))∑i=0nxi and the direction from the center to the ith vertex is denoted as ui=xi-x-, i=0,…,n, then U=ui∣i=0,…,n is a group of positive base.Proof.
Obviously,ui=xi-x-=xi-(1/(n+1))∑j=0nxj=-(1/(n+1))∑j=0,j≠inxj-xi. Assuming there exist constants l1,…,li,…,ln who are not all zero such that ∑i=1nliui=0 is tenable, then ∑i=1nliui=(-1/(n+1))∑i=1nli∑j=0,j≠inxj-xi=(-1/(n+1))∑i=1nnli-∑j=1,j≠inljx0-xi=0. As the side vectors x0-xi, i=1,…,n, are linearly independent, we get nli-∑j=1,j≠inlj=0, i=1,…,n, from which li=0, i=1,…,n is obtained. Therefore, the assumption is contradictory and untenable. Thus, ui,i=1,…,n is linearly independent, and rank⁡U=n is obtained. Therefore, U is a spanning set of Rn. On the other hand, ∑i=0nui=∑i=0nxi-x-=∑i=0nxi-n+1x-=0. In addition, the number of vectors is exactly right n+1. According to literatures [9–11], ui, i=0,…,n are positive bases.
The center ray operation searches along the direction ofui, i=0,…,n by comparing fx- and fxi, i=0,…,n. Firstly, as shown in (7), the first search is along the direction ui, i=0,…,n, where, μi must satisfy (8), and it is advised to calculate μi by (9) or (10). To be specific, when fx->fxi, we let μi=μouter>1, which means search outside the simplex. If fx-≤fxi, it is reasonable we assume better points locate between x- and xi, so we let μi=μinner∈0,1 in order to search inside the simplex:(7)x^ie=x-+μixi-x-,i=0,…,n,(8)μi=μinner∈0,1fx-≤fxiμouter∈1,∞others,(9)μouter=μ˙outer,μinner=μ˙inner,(10)μouter=μouterinit+μouterend-μouterinittmaxt,μinner=μinnerinit+μinnerend-μinnerinittmaxt.
Next, in order to look for a better point, the second detection is attempted according tofx^ie, indicated in (11). The second detection point is located by parameter μq whose scope is determined by comparing function values of x^ie,xi,x-, as indicated in (12):(11)x~ie=x-+μqμixi-x-,i=0,…,n,(12)μq=μe∈1,∞fx^ie<fxi<fx-μoc∈μi-1,1fxi≤fx^ie≤fx-μic∈1,μi-1fx-≤fxiμs∈0,μi-1fxi<fx-≤fx^ie.
Ifx^ie, which is found out by the first search, is the minimum point, then the direction is supposed as the descent direction. Naturally, a better point may be acquired after extending forward along the same direction. So let μq>1 in the case. If expansion point x~ie is better than x^ie, x~ie will be accepted; or else, x^ie will be accepted. When fxi≤fx^e≤fx-∨fx-<fxi, better points are supposed to locate between xi and x^e. Thus, set μi-1<μq<1, and the algorithm goes back for a certain distance to look for a better point x~ie. Then, xi is replaced with the better point from x~ie and x^ie, as a result, a new simplex will be formed. When fxi<fx-≤fx^e, assumption that better points locate between xi and x- is reasonable; thus, let μi-1<μq<0. If x~ie is better than xi, xi is replaced with x~ie. Otherwise, vertex xi still retains in the new simplex. Through the comparison between fxi and fx-, center ray operation can search both inside and outside the simplex, so the algorithm can get close to the extreme point no matter it is inside or outside the simplex. Some illustrations of the center ray operation under 2-dimensional circumstance are shown in Figure 1.
In the process of the center ray, the simplex can become larger or smaller. The adaptive change can reduce calculated amount and accelerate the convergence. If search direction is appropriate, the second search can expand outwards, which can enlarge the simplex and speed up the search. If the extreme point is inside the simplex, the simplex can be automatically downsized by the comparison of function values among the center, the vertex and the detected point. Consequently, the convergence will be accelerated. The search of traditional simplex method relies mainly on reflection operations; unfortunately, its search is only along one direction and only one vertex is updated. Computation of the traditional method is slow and its effect is poor in solving high-dimension problems. In comparison, the center ray operation searches along directions ofn+1 positive bases, and one operation may update several vertexes. As a result, not only the computation speed is improved but also the omission of solution is avoided.Figure 1
Illustration of center ray under 2-dimensional circumstance.


#### 3.1.3. Reflection Operation

In the simplex, the center of all vertexes except the worst vertexxn is denoted as x-′, as shown (13). Firstly, like the traditional simplex method, xn is reflected about x-′. As indicated in (14), x^r, namely, the reflection point of xn about x-′, is computed. Different from the traditional simplex method, the algorithm in the paper can search along two inverse reflection directions by comparison between fxn and fx-′, as shown in (15). Besides, like center ray operation, there are two attempts in reflection operation, as shown in (14) and (18), respectively:(13)x-′=∑i=0n-1xi,(14)x^r=x-′+αx-′-xn,(15)α=aforward>0fxn≥fx-′areverse<-1fxn<fx-′.The reflection coefficientα can be computed by either (16) or (17). As shown in (16), the reflection coefficient α is usually constant. When α decreases with iterations as shown in (17), the computation speed can be accelerated. In (17), tmax indicates the preset maximum iterations, αinitforward, αinitreverse are the reflection coefficients when t=0; and αendforward and αendreverse are the reflection coefficients when t=tmax:(16)αforward=α˙forward,αreverse=α˙reverse,(17)αforward=αinitforward-αinitforward-αendforwardtmaxt;0≤αendforward≤αinitforward,αreverse=αinitreverse-αinitreverse-αendreversetmaxt;0≤αinitreverse≤αendreverse.After that, by comparison amongfx^r, fx-′, fxn, x~r is computed by (18) to the second search. When fxn≥fx-′, coefficient αq must satisfy (19); when fxn<fx-′, αq must satisfy (20). The final reflection point is denoted as xr, as shown in (21). If fxr<fxn, xn is replaced with xr to form a new simplex. Otherwise, the reflection will fail. Figure 2 illustrates some of the process of reflection operations under 2-dimensional circumstance:(18)x~r=x-′+αqαx-′-xn,(19)αq=αfe∈1,∞fx^r<fx-′≤fxnαfc∈0,1fx-′≤fx^r<fxnαfs∈-α-1,0fx-′≤fxn<fx^r,(20)αq=αre∈1,∞fx^r<fxn<fx-′αrc∈-α-1,1fxn≤fx^r<fx-′αrs∈0,-α-1fxn<fx-′≤fx^r,(21)xr=x^rfx^r<fx~rx~rothers.Figure 2
Illustration of the reflection operation under 2-dimensional circumstance.Other than the traditional simplex method, both the center ray and the reflection in the paper twice searches along a same direction, likewise, object function is also computed twice. If results of the first detection are better, the algorithm outward expands along the same direction with a longer step to start the second search. Thus, computation speed will be accelerated. Otherwise, the algorithm will try to make certain contraction on the basis of the first search or search along inverse direction. This will improve the accuracy, and what is more important, this also can cause adaptive changes in both the size and the shape of the simplex in the search process. It means to automatically and continuously change the direction and step of the search, which enhance adaptability of the algorithm. One operation in the improved algorithm can achieve effects of several operations of the traditional algorithm, including reflection, expansion, contraction, and shrinkage. As long as the appropriate parameters are chosen, the computation efficiency can be greatly improved.


#### 3.1.4. Estimation of Gradient with Simplex

As derivative information is not utilized in the traditional simplex algorithm, the search is slow. Function values of at leastn+1 points are known in the simplex of n-dimensional problems. Therefore, the gradient even the Hessian matrix at each iteration can be estimated by these n+1 points. When one side of the simplex is in accordance with the negative gradient, the simplex may search along the descent direction.If theith side vectors xi-x0i=1,…,n connecting with the optimum point x0 is denoted as vi, then vi=viv→i, where vi=xi-x0 is its side length and v→i is its direction of unit vector, as shown in (22). Since the directional derivative of fx in the direction of v→i at x0 is (∂fx/∂v→i)x=x0=limt→0+⁡((fx0+tv→i-fx0)/t), di in (23) is approximate to (∂fx/∂v→i)x=x0 when the simplex is relatively small:(22)v→i=xi-x0xi-x0,(23)di=fxi-fx0xi-x0≈∂fx∂v→ix=x0.Suppose the gradient offx at x0 is ∇fx0; then, v→iT∇fx0=(∂fx/∂v→i)x=x0≈di. Let the matrix Vk=v→1,v→2,…,v→nT and the vector Dk=d1,d2,…,dnT; then, Vk∇fx0≈Dk, especially, when the simplex is relatively small. As v→1,v→2,…,v→n are linearly independent, (24) can be regarded as the estimated value of the gradient:(24)ρk=Vk-1Dk.In (24), the gradient can be accurately estimated without additionally calculating the function value. Especially, when the simplex is small or approximate to the extreme point, estimation is very accurate. Furthermore, Hessian matrix, Newton direction as well as the conjugate gradient at each iteration can be approximately estimated according to recent two successive gradients in the process of iterations.


#### 3.1.5. Flow of the Novel Simplex Method

For each simplex, the algorithm looks for a better point along rays emitting from the centerx- to the vertex xi to replace the current vertex xi. By continuously replacing the former simplex with a new one, iterations will be repeated until the evolution of the simplex stagnates. Once no better points are worked out or the simplex collapses, a new smaller simplex will be generated at the current optimal point in line with the current estimated descent direction and resume iterations again. The flow of the algorithm is shown in Figure 3.Figure 3
The flow of Algorithm3.Algorithm 3.
Step  1. Initialize the parameters l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  4.
Step  4. Set k=k+1, t=t+1.
Step  5. Estimate the gradient ρk of the simplex Sk.
Step  5.1. Compute each side vector v→i by (22) and construct matrix Vk.
Step  5.2. Approximately compute the directional derivative in the side vector by (23) and construct the matrix D=d1,d2,…,dnT.
Step  5.3. Compute ρk with (24).
Step  6. Generate simplex T after center ray operations of simplex Sk.
Step  6.1. Compute x-,fx- and set i=1.
Step  6.2. Compute μi by (9) or (10), and compute x^ie and fx^ie with (7) and (8).
Step  6.3. Determine μq that satisfies (12), and then compute x~ie,fx~ie with (11) to the second detect.
Step  6.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of simplex T.
Step  6.5. If i<n, then i=i+1 and go to Step  6.2; otherwise, go to Step  6.6.
Step  6.6. Sort vertexes of simplex T in ascending order by their function values.
Step  7. Whenfx0<fk-εm, then let fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T and go to Step  3. Otherwise, let fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  8.
Step  8. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, then quit the algorithm; otherwise, go to Step  2 and repeat iterations.In Algorithm3, the mth execution of the inner loop from Step  3 to Step  7 generates the subsequence zk,fk, km≤k<km+1, which satisfies fk≤fk-1-εm<fk-1. As the outer loop from Step  2 to Step  8 repeatedly performs the inner loop, the algorithm generates the infinite sequence zk,fk,ym,gm, which satisfies fk≤fk-1, gm≤gm-1.The center ray only searches alongn+1 fixed directions. If the direction is not descending, the center ray may fail, and the search will stop for the simplex tends to beinfinitely small. For some ill-conditioned problems, even though search directions are descent, it is theoretically feasible to find optimum value; however, this can be realized only when the step is extremely small. Since such extremely small simplex is likely to collapse or the computation is too slow, it is difficult to realize in practice. Although Algorithm 3 is theoretically convergent, sometimes, it actually tends to converge to saddle point and nonoptimum point in case of complex functions. The performance only can be improved through continuous changes of searching directions. A new simplex will be forced to generate when the number of successive execution of ray operations exceeds certain threshold, which is equivalent to changes of the search direction to some degree. Whereby Algorithm 4 is formed with improvements of Algorithm 3, as shown in Figure 4.Figure 4
The flow of Algorithm4.Algorithm 4.
Step  1. Initialize parameters l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. If t<tmax, go to Step  4; otherwise, go to Step  9.
Step  4. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  5.
Step  5. Set k=k+1, t=t+1.
Step  6. Estimate the gradient ρk of the simplex Sk.
Step  6.1. Compute each side vector v→i by (22) and create the matrix Vk.
Step  6.2. Approximately compute the directional derivative in the side vector by (23) to construct the matrix D=d1,d2,…,dnT.
Step  6.3. Compute ρk by (24).
Step  7. Generate simplex T after the application of center ray operations of simplex Sk.
Step  7.1. Compute x-,fx- and let i=1.
Step  7.2. Compute μi by (9) or (10) and compute x^ie and fx^ie with (7) and (8).
Step  7.3. Determine μq which satisfies (12), compute x~ie,fx~ie with (11) to the second search.
Step  7.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of T.
Step  7.5. If i<n, then set i=i+1 and go to Step  7.2; otherwise, go to Step  7.6.
Step  7.6. Sort vertexes of simplex T in ascending order by their function values.
Step  8. When fx0<fk-εm, let fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T and go to Step  3. Otherwise, let fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  9.
Step  9. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, quit the algorithm; otherwise, go to Step  2 and repeat computation.Different from Algorithm3, Algorithm 4 makes restrictions on the times of the repetitive execution of ray operations, denoted as t, which is not allowed to surpass the maximum value tmax. If t≥tmax, even though ray operations are successful, a new simplex at the current optimal point zk is generated according to the current ρk so as to change the search direction.Considering the better performance of reflection operations in the searching process of traditional simplex, in order to further improve the search abilities, a reflection operation will be performed in case of failure of ray operations. Consequently, search directions are more abundant and Algorithm5 is formed, as shown in Figure 5.Figure 5
The flow of Algorithm5.Algorithm 5.
Step  1. Initialize parameters as l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. If t<tmax, go to Step  4; otherwise, go to Step  11.
Step  4. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  5.
Step  5. Set k=k+1, t=t+1.
Step  6. Estimate the gradient ρk of the simplex Sk.
Step  6.1. Compute each side vector v→i by (22) to construct matrix Vk.
Step  6.2. Approximately compute the directional derivative in the side vector by (23), and construct matrix D=d1,d2,…,dnT.
Step  6.3. Compute ρk by (24).
Step  7. Generate simplex T1 after the application of center ray operations of simplex Sk.
Step  7.1. Compute x-,fx- and set i=1.
Step  7.2. Determine μi by (9) or (10) and compute x^ie and fx^ie by (7) and (8).
Step  7.3. Determine μq which satisfies (12) and compute x~ie,fx~ie by (11) to the second detection.
Step  7.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T1 to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of simplex T1.
Step  7.5. If i<n, then set i=i+1 and go to Step  7.2; otherwise, go to Step  7.6.
Step  7.6. Sort vertexes of simplex T1 in ascending order by their function value.
Step  8. Whenfx0<fk-εm, set fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T1 and go to Step  3. Otherwise, go to Step  9.
Step  9. Form novel simplex T2 after the reflection of simplex T1.
Step  9.1. Compute x-′ by (13).
Step  9.2. The reflection coefficient α is computed by (16) or (17) and x^r and fx^r are computed by (14).
Step  9.3. Determine αq which satisfies (19) or (20), and then solve x~r by (18) and compute fx~r.
Step  9.4. Determine reflection point xr by (21). If fxr<fxn, then xn is replaced with xr to form new simplex T2 and sort vertexes of T2 in ascending order by function value.
Step  10. Whenfx0<fk-εm, set fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T2 and go to Step  3. Otherwise, set fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  11.
Step  11. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, quit the algorithm; otherwise, go to Step  2 to repeatedly do loop computations.Similarly, as for Algorithm5, Step  3 to Step  10 construct inner loop while Step  2 to Step  11 form the outer loop.Algorithms3, 4, and 5 all have two layers of loops. The sequence zk is composed of optimal points of simplex Sk. The mth execution of the inner loop generates the finite sequence zk,km≤k<km+1, which satisfies fzk≤fzk-1-εm<fzk-1. The outer loop will be repeatedly executed, which may generate the infinite sequence zk satisfying fzk+1≤fzk. That is to say, the sequence fzk generated by the algorithm must be nonincreasing. The sequence ym is composed of best vertexes of the simplex when both ray operations and reflection operations fail, which satisfies fym+1≤fym. Upon the failure of ray and reflection operations, a new smaller orthogonal simplex will be generated. The algorithm can ensure the current minimum point to be one of the vertexes of the new simplex. Besides, it also ensures that the current descent direction is parallel to an edge of the new simplex. The generation of a new smaller simplex achieves better effect than shrinkage operations, which can reduce the possibility of the simplex degeneration. lk indicates the number of simplexes generation till the kth iteration. km means that the mth simplex generation occurs at the kth iteration. hm means the initial edge length of the mth generated simplex. For the sake of enhancing the efficiency, hm can increase in the process of iterations. However, it is certain that there are decreasing subsequences in hm. εm is the control error, which is a monotonically decreasing sequence in the process of iterations. In order to converge, hm,εm must satisfy limm→∞⁡hm=0 and limm→∞⁡(εm/hm)=0. In the paper, hm,εm are solved by (25) or (26):(25)hm=hm-1λh,εm=εm-1λε,(26)λh=λ˙ht≥tmaxλ¨hρm≤ρnei∨t≤tdec∧hm-1>hsmallλ⃛hothers,λε=λ¨ερm≤ρneiλ⃛ερm>ρnei.The variation ofhm has important influences on the performance of the algorithm. The rate of change of hm,εm must satisfy 1<λ˙h, 0<λ¨h<λ⃛h<1, 0<λ¨ε<λ⃛ε<1 in the paper.Theorem 6.
As for Algorithm3, assuming the sequence zk is bounded, the objective function fx is continuously differentiable and f′x is Lipschitz in any bounded region of Rn; then, ym converges at the stationary point of fx.Proof.
It can be learned from Theorem2 that the search direction ui=xi-x-, i=0,…,n, is positive bases. From the flow of the algorithms, the finite sequence fzk, km≤k<km+1, generated by the inner loop is monotonically decreasing as fzk>fzk+1, km≤k<km+1. Upon the exiting from the inner loop, it is certain fzk+1+εm≥fzk will hold. So fzk≥fzk+1 will be obtained by the algorithms. Similar to the frame of Algorithms  3.1 and 3.2 described in Coope and Price (2000), zk constitutes a quasiframe [11]. Comparing those two algorithms in Coope and Price (2000) with Algorithms 3, 4, and 5 in the paper, the algorithm proposed in the paper can be regarded as a concrete instance of Algorithms  3.1 and 3.2 in Coope and Price (2000). The side vectors xi-x0, i=1,…,n, and x--x0 of the simplex are regarded as positive bases in Coope and Price (2000) [11]. However, xi-x-, i=0,…,n, are taken as positive bases in the paper. Algorithms in Coope and Price (2000) have same convergence properties as algorithms in the paper, whose convergence is proved in Coope and Price (2000) in details [11]. Compared with Algorithm 3, some acceleration methods are adopted in Algorithms 4 and 5, so their convergence is similar to Algorithm 3.


### 3.2. Experimental Verification of the Novel Simplex Algorithm

Schwefel and Rosenbrock functions are taken as benchmarks to test the performance of the improved simplex algorithm [12]. In the experiments, the initial point at each computation is all x0=50,…,50T. Once the number of iterations exceeds the predetermined maximum value or hm approaches zero or the solution is obtained, the algorithm will immediately be terminated. μouter and μinner are computed by (9) and, respectively, set μouter=1.2, μinner=0.9, μe=1.25, μoc=0.92, μic=1.06, and μs=0.4. Besides, reflection coefficient α is computed by (16) and set αforward=1.3, αreverse=-2.3, αfe=1.25, αfc=0.7, αfs=-0.4, αre=1.25, αrc=0.72, and αrs=0.15. In the experiments, the other parameters are h0=5,ε0=0.2,ρnei=0.1,λ¨ε=0.5,λ⃛ε=0.8,λ˙h=1.5. Assuming that actual number of iterations is Ω1, the number of function computation is Ω2 and the optimal solution is fx^, then Ω3=x^-x0/Ω1, Ω4=fx^-fx0/Ω1, Ω5=x^-x0/Ω2, Ω6=fx^-fx0/Ω2, Ω7=x^-x∗ can be utilized to evaluate the performance of the algorithm. Smaller Ω1, Ω2, Ω7 and larger Ω3, Ω4, Ω5, Ω6 indicate the performance of the algorithm is better.Price and Coope (2002) has improved NM Simplex method and introduced ghost simplex [13]. The algorithm revised by PRICE is tested in the paper. It is concluded that either Price’s algorithm or NM Simplex method can hardly solve when the dimension of Rosenbrock functions is more than 20. In other words, in solving high-dimension problems, N M Simplex method has presented very poor effect and the algorithm proposed by PRICE also makes little improvement. Taking 5-dimensional Schwefel function as the benchmark, the algorithm is terminated once Ω7=4e-5. Then, the results of this paper are compared with those of Price’s method. Some results are listed in Table 2. In order to verify the performance of the algorithm proposed in the paper in high-dimension optimization problems, 90-dimensional Schwefel function and 45-dimensional Rosenbrock function are solved by Algorithms 4, 5, and 7 (Figure 6). Some results are shown in Tables 3~8. In order to compare effects of all algorithms, average values are totaled in Tables 3~8, as shown in Figure 5. The experiment proves that the effect and performance of improved simplex algorithm proposed in the paper is much better than that of traditional simplex methods.Table 2
Comparison between the algorithm proposed in the paper and Price.
Algorithm
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Price
1.361376
E
-
09
2457
3569
0.045504
55.96256
0.031326
38.5262
Algorithm3
2.161297
E
-
09
602
11226
0.18572
228.4053
0.009959
12.24835
Algorithm4
1.703474
E
-
09
526
11305
0.212554
261.4068
0.00989
12.16276
Algorithm5
1.314375
E
-
09
448
9220
0.249561
306.9196
0.012126
14.91323Table 3
Parts of results of solving 90-dimensional Schwefel function by Algorithm3.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.99
6.96
E
-
08
60683
13153005
0.007817
10178.51
3.61
E
-
05
46.9598
0.000516
0.9
0.98
7.37
E
-
08
53731
11575442
0.008828
11495.46
4.1
E
-
05
53.35973
0.000465
0.9
0.975
8.98
E
-
08
48455
10472012
0.009789
12747.14
4.53
E
-
05
58.98222
0.000593
0.85
0.99
6.5
E
-
08
59400
12882921
0.007986
10398.36
3.68
E
-
05
47.94429
0.000504
0.85
0.98
7.38
E
-
08
49668
10766210
0.00955
12435.82
4.41
E
-
05
57.37047
0.000429
0.85
0.975
8.98
E
-
08
48455
10472012
0.009789
12747.14
4.53
E
-
05
58.98222
0.000593Table 4
Parts of results of solving 45-dimensional Rosenbrock function by Algorithm3.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.995
3.32
E
-
09
299760
30953026
0.001097
88107.5
1.06
E
-
05
853.2641
0.000115
0.9
0.985
3.32
E
-
09
446429
46260909
0.000736
59160.82
7.11
E
-
06
570.9163
0.000115
0.85
0.995
3.32
E
-
09
316215
32404535
0.001039
83522.62
1.01
E
-
05
815.0435
0.000115
0.85
0.985
3.32
E
-
09
454935
47117032
0.000723
58054.68
6.98
E
-
06
560.5426
0.000115
0.8
0.995
3.32
E
-
09
303212
31274881
0.001084
87104.42
1.05
E
-
05
844.483
0.000115
0.8
0.985
3.32
E
-
09
460166
47576852
0.000714
57394.74
6.91
E
-
06
555.1251
0.000115Table 5
Parts of results of solving 90-dimensional Schwefel function by Algorithm4.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.99
8.34
E
-
08
67054
14486323
0.007074
9211.419
3.27
E
-
05
42.63763
0.000567
0.9
0.98
8.78
E
-
08
130898
27575059
0.003624
4718.655
1.72
E
-
05
22.39932
0.000589
0.9
0.975
7.08
E
-
08
194194
40244951
0.002443
3180.647
1.18
E
-
05
15.34758
0.00052
0.85
0.99
8.75
E
-
08
73219
15845655
0.006478
8435.823
2.99
E
-
05
38.97993
0.000584
0.85
0.98
8.52
E
-
08
138163
29067816
0.003433
4470.535
1.63
E
-
05
21.24902
0.00058
0.85
0.975
4.68
E
-
08
195527
40501017
0.002426
3158.963
3.23
E
-
05
15.25054
0.000505Table 6
Parts of results of solving 45-dimensional Rosenbrock function by Algorithm4.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.995
3.32
E
-
09
285070
29702517
0.001153
92647.79
1.11
E
-
05
889.1875
1.72
E
-
05
0.9
0.985
3.32
E
-
09
287837
29884227
0.001142
91757.16
1.10
E
-
05
883.7808
1.72
E
-
05
0.85
0.995
3.32
E
-
09
283558
29541881
0.001159
93141.81
1.11
E
-
05
894.0225
1.72
E
-
05
0.85
0.985
3.32
E
-
09
285806
29672230
0.00115
92409.21
1.11
E
-
05
890.0951
1.72
E
-
05
0.8
0.995
3.32
E
-
09
287742
29952651
0.001142
91787.45
1.10
E
-
05
881.7619
1.72
E
-
05
0.8
0.985
3.32
E
-
09
288202
29917444
0.001141
91640.95
1.10
E
-
05
882.7995
1.72
E
-
05Table 7
Parts of results of solving 90-dimensional Schwefel function by Algorithm5.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.99
4.06
E
-
08
74783
21242674
0.006343
8259.397
2.23
E
-
05
29.07649
0.0004
0.9
0.98
7.08
E
-
08
75482
16388224
0.006284
8182.911
2.89
E
-
05
37.68941
0.000516
0.9
0.975
8.86
E
-
08
81503
17609179
0.00582
7578.402
2.69
E
-
05
35.07617
0.000581
0.85
0.99
7.57
E
-
08
59474
13368291
0.007976
10385.42
3.55
E
-
05
46.20355
0.000542
0.85
0.98
7.50
E
-
08
72768
15772677
0.006519
8488.106
3.01
E
-
05
39.16028
0.000532
0.85
0.975
2.38
E
-
08
84097
18157608
0.00564
7344.644
1.42
E
-
05
34.01673
0.000611Table 8
Parts of results of solving 45-dimensional Rosenbrock function by Algorithm5.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.995
3.32
E
-
09
391063
52265101
0.000841
67536.7
6.29
E
-
06
505.3297
1.72
E
-
05
0.9
0.985
3.32
E
-
09
323520
40051422
0.001016
81636.7
8.21
E
-
06
659.4299
1.72
E
-
05
0.85
0.995
3.32
E
-
09
394208
52604873
0.000834
66997.89
6.25
E
-
06
502.0658
1.72
E
-
05
0.85
0.985
3.33
E
-
09
323368
39892530
0.001016
81675.07
8.24
E
-
06
662.0564
1.72
E
-
05
0.8
0.995
3.34
E
-
09
395380
52745863
0.000831
66799.3
6.23
E
-
06
500.7237
1.72
E
-
05
0.8
0.985
3.33
E
-
09
325750
40212144
0.001009
81077.84
8.17
E
-
06
656.7943
1.72
E
-
05Figure 6
Comparison of different effects of Algorithms3, 4, and 5.In terms of optimization of low-dimensional functions, though the algorithm proposed by the paper and by PRICE can solve such problems, the improved simplex method in the paper demonstrates higher efficiency. It can be known from Table2 that the algorithm proposed in the paper computes much rapidly although the times of object function computation by Algorithms 3, 4, and 5 is more than that of Price’s algorithm. The object function is calculated twice every detection in iteration, so not only Ω2 of the algorithm proposed by the paper may be great than that of other simplex method but also Ω5, Ω6 of algorithm in my paper may be less than that of other simplex method under low dimensional circumstance. In addition, as Price’s algorithm and the traditional simplex algorithm are quite inefficient in high-dimension optimization problems, they can hardly solve any problem as long as the number of dimensions is over 20. In comparison, the algorithm proposed in the paper still behaves more effectively in solving high-dimension problems. It indicates that the improved simplex method, which can update several vertexes at each iteration by center ray operations, is effective and efficient.With respect to the performance, generally, Algorithm5 is the best, while Algorithm 4 is better than Algorithm 3. For both Rosenbrock function and Schwefel function, Algorithm 5 can obtain the smallest Ω1, Ω2 and the largest Ω3, Ω4, Ω5, Ω6; this implies that increasing of search directions can greatly improve the efficiency. Ω1, Ω2 of Algorithm 4 are smaller than those of Algorithm 3, which indicates that the former costs less in computation. In general, Ω3, Ω4, Ω5, Ω6 of Algorithm 3 is smaller than that of Algorithm 4, which reflects that efficiency and speed of Algorithm 4 is higher than that of Algorithm 3.The parameters of the algorithm, especially the parameters controlling changes ofhm have great influence on the performance of the algorithm. According to the experiment, it is usually possible for the algorithm to get good effects when main parameters are within the following ranges: (27)μinner∈0.5,1,μouter∈1.2,2.2,μe∈1.1,1.6,μoc∈0.5,0.8,μic∈1.1,1.8,μs∈0.1,0.5,λ˙h∈1.5,3,λ¨h∈0.5,0.95,λ⃛h∈0.95,0.999,tmax∈6,16,ρnei∈0.1,0.2,hsmall∈0.5,0.05,h0∈0.5,5,λ¨ε∈0.5,0.95,λ⃛ε∈0.8,0.9,ε0∈0.2,0.5,αforward∈0.5,2.5,αreverse∈-1.1,-3,αinitforward∈2,3,αendforward∈0.3,0.6,αinitreverse∈-2.5,-3.5,αendreverse∈-1.1,-1.5,αfe∈1.2,2.5,αfc∈0.2,0.8,αfs∈-2,-0.4,αre∈1.2,2.5,αrc∈0.3,0.9,αrs∈0,0.4.


## 3.1. Description of the Novel Simplex Method

The convex hull, namely, the hyper polyhedron composed ofn+1 vertexes in Rn is called a simplex. Assuming vertexes x0,x1,…,xn in the simplex are sorted in ascending order by function values, the simplex of the kth iteration is denoted as Sk=x0,…,xn. Obviously, n side vectors xi-xj, j=0,…,n∧j≠i, connecting with the ith vertex xi are linearly independent, from which the degeneration of the simplex can be judged. The whole iteration of the traditional nonlinear simplex algorithm mainly consists of four steps: reflection, expansion, contraction, and shrinkage. In contrast, the novel algorithm in the paper primarily solves problems by generating a simplex, raying outward from the simplex center and reflecting operations. Taking the unconstrained optimization problem min⁡fx as an example, the novel simplex algorithm can be described as below.


### 3.1.1. Generation of a Simplex

Due to the tendency of deformation and degeneration in the process of the simplex’s iteration, it is usually better that the initial simplex is constructed as an orthogonal one. When one side of the orthogonal simplex gets close to the negative gradient direction-ρk, the simplex will search much rapidly. If Ai=τ1,…,τi and z is any of solution of the linear equations AiTZ=0, then z will be orthogonal to τ1,…,τi. Firstly, we set τ1=-ρk and A1=τ1. If z, that is, any of the solutions of AiTZ=0, is taken as the last column of Ai to generate a new matrix Ai+1, then the orthogonal matrix An=τ1,…,τn can be obtained after repeatedly solving the linear equations n-1 times, where, τ1,…,τn represent n orthogonal vectors. Subsequently, we choose the point whose function value is better from zk+hmτi,zk-hmτi as the vertex of the simplex. In this paper, the flow of simplex generation is shown as Algorithm 1.Algorithm 1.
Step  1. Let τ1=-ρk/ρk, i=1, Sk=zk, A1=τ1.
Step  2. Let(5)Sk=Sk∪zk+hmτifzk+hmτi<fzk-hmτiSk∪zk-hmτiothers.
Step  3. Solve the linear equations AiTZ=0 and choose any of the solutions as z and let (6)τi+1=zz,Ai+1=τ1,…,τi,τi+1.
Step  4. When i<n-1, then i=i+1 and go to Step  2; otherwise exit.Following the procedure, the simplexSk is generated at zk. Edges connecting with zk are orthogonal vectors whose lengths are all hm. Besides, the first edge is parallel to the negative gradient direction -ρk.


### 3.1.2. Center Ray Operation

Theorem 2.
If the center of the simplex is denoted asx-=(1/(n+1))∑i=0nxi and the direction from the center to the ith vertex is denoted as ui=xi-x-, i=0,…,n, then U=ui∣i=0,…,n is a group of positive base.Proof.
Obviously,ui=xi-x-=xi-(1/(n+1))∑j=0nxj=-(1/(n+1))∑j=0,j≠inxj-xi. Assuming there exist constants l1,…,li,…,ln who are not all zero such that ∑i=1nliui=0 is tenable, then ∑i=1nliui=(-1/(n+1))∑i=1nli∑j=0,j≠inxj-xi=(-1/(n+1))∑i=1nnli-∑j=1,j≠inljx0-xi=0. As the side vectors x0-xi, i=1,…,n, are linearly independent, we get nli-∑j=1,j≠inlj=0, i=1,…,n, from which li=0, i=1,…,n is obtained. Therefore, the assumption is contradictory and untenable. Thus, ui,i=1,…,n is linearly independent, and rank⁡U=n is obtained. Therefore, U is a spanning set of Rn. On the other hand, ∑i=0nui=∑i=0nxi-x-=∑i=0nxi-n+1x-=0. In addition, the number of vectors is exactly right n+1. According to literatures [9–11], ui, i=0,…,n are positive bases.
The center ray operation searches along the direction ofui, i=0,…,n by comparing fx- and fxi, i=0,…,n. Firstly, as shown in (7), the first search is along the direction ui, i=0,…,n, where, μi must satisfy (8), and it is advised to calculate μi by (9) or (10). To be specific, when fx->fxi, we let μi=μouter>1, which means search outside the simplex. If fx-≤fxi, it is reasonable we assume better points locate between x- and xi, so we let μi=μinner∈0,1 in order to search inside the simplex:(7)x^ie=x-+μixi-x-,i=0,…,n,(8)μi=μinner∈0,1fx-≤fxiμouter∈1,∞others,(9)μouter=μ˙outer,μinner=μ˙inner,(10)μouter=μouterinit+μouterend-μouterinittmaxt,μinner=μinnerinit+μinnerend-μinnerinittmaxt.
Next, in order to look for a better point, the second detection is attempted according tofx^ie, indicated in (11). The second detection point is located by parameter μq whose scope is determined by comparing function values of x^ie,xi,x-, as indicated in (12):(11)x~ie=x-+μqμixi-x-,i=0,…,n,(12)μq=μe∈1,∞fx^ie<fxi<fx-μoc∈μi-1,1fxi≤fx^ie≤fx-μic∈1,μi-1fx-≤fxiμs∈0,μi-1fxi<fx-≤fx^ie.
Ifx^ie, which is found out by the first search, is the minimum point, then the direction is supposed as the descent direction. Naturally, a better point may be acquired after extending forward along the same direction. So let μq>1 in the case. If expansion point x~ie is better than x^ie, x~ie will be accepted; or else, x^ie will be accepted. When fxi≤fx^e≤fx-∨fx-<fxi, better points are supposed to locate between xi and x^e. Thus, set μi-1<μq<1, and the algorithm goes back for a certain distance to look for a better point x~ie. Then, xi is replaced with the better point from x~ie and x^ie, as a result, a new simplex will be formed. When fxi<fx-≤fx^e, assumption that better points locate between xi and x- is reasonable; thus, let μi-1<μq<0. If x~ie is better than xi, xi is replaced with x~ie. Otherwise, vertex xi still retains in the new simplex. Through the comparison between fxi and fx-, center ray operation can search both inside and outside the simplex, so the algorithm can get close to the extreme point no matter it is inside or outside the simplex. Some illustrations of the center ray operation under 2-dimensional circumstance are shown in Figure 1.
In the process of the center ray, the simplex can become larger or smaller. The adaptive change can reduce calculated amount and accelerate the convergence. If search direction is appropriate, the second search can expand outwards, which can enlarge the simplex and speed up the search. If the extreme point is inside the simplex, the simplex can be automatically downsized by the comparison of function values among the center, the vertex and the detected point. Consequently, the convergence will be accelerated. The search of traditional simplex method relies mainly on reflection operations; unfortunately, its search is only along one direction and only one vertex is updated. Computation of the traditional method is slow and its effect is poor in solving high-dimension problems. In comparison, the center ray operation searches along directions ofn+1 positive bases, and one operation may update several vertexes. As a result, not only the computation speed is improved but also the omission of solution is avoided.Figure 1
Illustration of center ray under 2-dimensional circumstance.


### 3.1.3. Reflection Operation

In the simplex, the center of all vertexes except the worst vertexxn is denoted as x-′, as shown (13). Firstly, like the traditional simplex method, xn is reflected about x-′. As indicated in (14), x^r, namely, the reflection point of xn about x-′, is computed. Different from the traditional simplex method, the algorithm in the paper can search along two inverse reflection directions by comparison between fxn and fx-′, as shown in (15). Besides, like center ray operation, there are two attempts in reflection operation, as shown in (14) and (18), respectively:(13)x-′=∑i=0n-1xi,(14)x^r=x-′+αx-′-xn,(15)α=aforward>0fxn≥fx-′areverse<-1fxn<fx-′.The reflection coefficientα can be computed by either (16) or (17). As shown in (16), the reflection coefficient α is usually constant. When α decreases with iterations as shown in (17), the computation speed can be accelerated. In (17), tmax indicates the preset maximum iterations, αinitforward, αinitreverse are the reflection coefficients when t=0; and αendforward and αendreverse are the reflection coefficients when t=tmax:(16)αforward=α˙forward,αreverse=α˙reverse,(17)αforward=αinitforward-αinitforward-αendforwardtmaxt;0≤αendforward≤αinitforward,αreverse=αinitreverse-αinitreverse-αendreversetmaxt;0≤αinitreverse≤αendreverse.After that, by comparison amongfx^r, fx-′, fxn, x~r is computed by (18) to the second search. When fxn≥fx-′, coefficient αq must satisfy (19); when fxn<fx-′, αq must satisfy (20). The final reflection point is denoted as xr, as shown in (21). If fxr<fxn, xn is replaced with xr to form a new simplex. Otherwise, the reflection will fail. Figure 2 illustrates some of the process of reflection operations under 2-dimensional circumstance:(18)x~r=x-′+αqαx-′-xn,(19)αq=αfe∈1,∞fx^r<fx-′≤fxnαfc∈0,1fx-′≤fx^r<fxnαfs∈-α-1,0fx-′≤fxn<fx^r,(20)αq=αre∈1,∞fx^r<fxn<fx-′αrc∈-α-1,1fxn≤fx^r<fx-′αrs∈0,-α-1fxn<fx-′≤fx^r,(21)xr=x^rfx^r<fx~rx~rothers.Figure 2
Illustration of the reflection operation under 2-dimensional circumstance.Other than the traditional simplex method, both the center ray and the reflection in the paper twice searches along a same direction, likewise, object function is also computed twice. If results of the first detection are better, the algorithm outward expands along the same direction with a longer step to start the second search. Thus, computation speed will be accelerated. Otherwise, the algorithm will try to make certain contraction on the basis of the first search or search along inverse direction. This will improve the accuracy, and what is more important, this also can cause adaptive changes in both the size and the shape of the simplex in the search process. It means to automatically and continuously change the direction and step of the search, which enhance adaptability of the algorithm. One operation in the improved algorithm can achieve effects of several operations of the traditional algorithm, including reflection, expansion, contraction, and shrinkage. As long as the appropriate parameters are chosen, the computation efficiency can be greatly improved.


### 3.1.4. Estimation of Gradient with Simplex

As derivative information is not utilized in the traditional simplex algorithm, the search is slow. Function values of at leastn+1 points are known in the simplex of n-dimensional problems. Therefore, the gradient even the Hessian matrix at each iteration can be estimated by these n+1 points. When one side of the simplex is in accordance with the negative gradient, the simplex may search along the descent direction.If theith side vectors xi-x0i=1,…,n connecting with the optimum point x0 is denoted as vi, then vi=viv→i, where vi=xi-x0 is its side length and v→i is its direction of unit vector, as shown in (22). Since the directional derivative of fx in the direction of v→i at x0 is (∂fx/∂v→i)x=x0=limt→0+⁡((fx0+tv→i-fx0)/t), di in (23) is approximate to (∂fx/∂v→i)x=x0 when the simplex is relatively small:(22)v→i=xi-x0xi-x0,(23)di=fxi-fx0xi-x0≈∂fx∂v→ix=x0.Suppose the gradient offx at x0 is ∇fx0; then, v→iT∇fx0=(∂fx/∂v→i)x=x0≈di. Let the matrix Vk=v→1,v→2,…,v→nT and the vector Dk=d1,d2,…,dnT; then, Vk∇fx0≈Dk, especially, when the simplex is relatively small. As v→1,v→2,…,v→n are linearly independent, (24) can be regarded as the estimated value of the gradient:(24)ρk=Vk-1Dk.In (24), the gradient can be accurately estimated without additionally calculating the function value. Especially, when the simplex is small or approximate to the extreme point, estimation is very accurate. Furthermore, Hessian matrix, Newton direction as well as the conjugate gradient at each iteration can be approximately estimated according to recent two successive gradients in the process of iterations.


### 3.1.5. Flow of the Novel Simplex Method

For each simplex, the algorithm looks for a better point along rays emitting from the centerx- to the vertex xi to replace the current vertex xi. By continuously replacing the former simplex with a new one, iterations will be repeated until the evolution of the simplex stagnates. Once no better points are worked out or the simplex collapses, a new smaller simplex will be generated at the current optimal point in line with the current estimated descent direction and resume iterations again. The flow of the algorithm is shown in Figure 3.Figure 3
The flow of Algorithm3.Algorithm 3.
Step  1. Initialize the parameters l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  4.
Step  4. Set k=k+1, t=t+1.
Step  5. Estimate the gradient ρk of the simplex Sk.
Step  5.1. Compute each side vector v→i by (22) and construct matrix Vk.
Step  5.2. Approximately compute the directional derivative in the side vector by (23) and construct the matrix D=d1,d2,…,dnT.
Step  5.3. Compute ρk with (24).
Step  6. Generate simplex T after center ray operations of simplex Sk.
Step  6.1. Compute x-,fx- and set i=1.
Step  6.2. Compute μi by (9) or (10), and compute x^ie and fx^ie with (7) and (8).
Step  6.3. Determine μq that satisfies (12), and then compute x~ie,fx~ie with (11) to the second detect.
Step  6.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of simplex T.
Step  6.5. If i<n, then i=i+1 and go to Step  6.2; otherwise, go to Step  6.6.
Step  6.6. Sort vertexes of simplex T in ascending order by their function values.
Step  7. Whenfx0<fk-εm, then let fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T and go to Step  3. Otherwise, let fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  8.
Step  8. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, then quit the algorithm; otherwise, go to Step  2 and repeat iterations.In Algorithm3, the mth execution of the inner loop from Step  3 to Step  7 generates the subsequence zk,fk, km≤k<km+1, which satisfies fk≤fk-1-εm<fk-1. As the outer loop from Step  2 to Step  8 repeatedly performs the inner loop, the algorithm generates the infinite sequence zk,fk,ym,gm, which satisfies fk≤fk-1, gm≤gm-1.The center ray only searches alongn+1 fixed directions. If the direction is not descending, the center ray may fail, and the search will stop for the simplex tends to beinfinitely small. For some ill-conditioned problems, even though search directions are descent, it is theoretically feasible to find optimum value; however, this can be realized only when the step is extremely small. Since such extremely small simplex is likely to collapse or the computation is too slow, it is difficult to realize in practice. Although Algorithm 3 is theoretically convergent, sometimes, it actually tends to converge to saddle point and nonoptimum point in case of complex functions. The performance only can be improved through continuous changes of searching directions. A new simplex will be forced to generate when the number of successive execution of ray operations exceeds certain threshold, which is equivalent to changes of the search direction to some degree. Whereby Algorithm 4 is formed with improvements of Algorithm 3, as shown in Figure 4.Figure 4
The flow of Algorithm4.Algorithm 4.
Step  1. Initialize parameters l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. If t<tmax, go to Step  4; otherwise, go to Step  9.
Step  4. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  5.
Step  5. Set k=k+1, t=t+1.
Step  6. Estimate the gradient ρk of the simplex Sk.
Step  6.1. Compute each side vector v→i by (22) and create the matrix Vk.
Step  6.2. Approximately compute the directional derivative in the side vector by (23) to construct the matrix D=d1,d2,…,dnT.
Step  6.3. Compute ρk by (24).
Step  7. Generate simplex T after the application of center ray operations of simplex Sk.
Step  7.1. Compute x-,fx- and let i=1.
Step  7.2. Compute μi by (9) or (10) and compute x^ie and fx^ie with (7) and (8).
Step  7.3. Determine μq which satisfies (12), compute x~ie,fx~ie with (11) to the second search.
Step  7.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of T.
Step  7.5. If i<n, then set i=i+1 and go to Step  7.2; otherwise, go to Step  7.6.
Step  7.6. Sort vertexes of simplex T in ascending order by their function values.
Step  8. When fx0<fk-εm, let fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T and go to Step  3. Otherwise, let fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  9.
Step  9. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, quit the algorithm; otherwise, go to Step  2 and repeat computation.Different from Algorithm3, Algorithm 4 makes restrictions on the times of the repetitive execution of ray operations, denoted as t, which is not allowed to surpass the maximum value tmax. If t≥tmax, even though ray operations are successful, a new simplex at the current optimal point zk is generated according to the current ρk so as to change the search direction.Considering the better performance of reflection operations in the searching process of traditional simplex, in order to further improve the search abilities, a reflection operation will be performed in case of failure of ray operations. Consequently, search directions are more abundant and Algorithm5 is formed, as shown in Figure 5.Figure 5
The flow of Algorithm5.Algorithm 5.
Step  1. Initialize parameters as l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. If t<tmax, go to Step  4; otherwise, go to Step  11.
Step  4. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  5.
Step  5. Set k=k+1, t=t+1.
Step  6. Estimate the gradient ρk of the simplex Sk.
Step  6.1. Compute each side vector v→i by (22) to construct matrix Vk.
Step  6.2. Approximately compute the directional derivative in the side vector by (23), and construct matrix D=d1,d2,…,dnT.
Step  6.3. Compute ρk by (24).
Step  7. Generate simplex T1 after the application of center ray operations of simplex Sk.
Step  7.1. Compute x-,fx- and set i=1.
Step  7.2. Determine μi by (9) or (10) and compute x^ie and fx^ie by (7) and (8).
Step  7.3. Determine μq which satisfies (12) and compute x~ie,fx~ie by (11) to the second detection.
Step  7.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T1 to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of simplex T1.
Step  7.5. If i<n, then set i=i+1 and go to Step  7.2; otherwise, go to Step  7.6.
Step  7.6. Sort vertexes of simplex T1 in ascending order by their function value.
Step  8. Whenfx0<fk-εm, set fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T1 and go to Step  3. Otherwise, go to Step  9.
Step  9. Form novel simplex T2 after the reflection of simplex T1.
Step  9.1. Compute x-′ by (13).
Step  9.2. The reflection coefficient α is computed by (16) or (17) and x^r and fx^r are computed by (14).
Step  9.3. Determine αq which satisfies (19) or (20), and then solve x~r by (18) and compute fx~r.
Step  9.4. Determine reflection point xr by (21). If fxr<fxn, then xn is replaced with xr to form new simplex T2 and sort vertexes of T2 in ascending order by function value.
Step  10. Whenfx0<fk-εm, set fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T2 and go to Step  3. Otherwise, set fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  11.
Step  11. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, quit the algorithm; otherwise, go to Step  2 to repeatedly do loop computations.Similarly, as for Algorithm5, Step  3 to Step  10 construct inner loop while Step  2 to Step  11 form the outer loop.Algorithms3, 4, and 5 all have two layers of loops. The sequence zk is composed of optimal points of simplex Sk. The mth execution of the inner loop generates the finite sequence zk,km≤k<km+1, which satisfies fzk≤fzk-1-εm<fzk-1. The outer loop will be repeatedly executed, which may generate the infinite sequence zk satisfying fzk+1≤fzk. That is to say, the sequence fzk generated by the algorithm must be nonincreasing. The sequence ym is composed of best vertexes of the simplex when both ray operations and reflection operations fail, which satisfies fym+1≤fym. Upon the failure of ray and reflection operations, a new smaller orthogonal simplex will be generated. The algorithm can ensure the current minimum point to be one of the vertexes of the new simplex. Besides, it also ensures that the current descent direction is parallel to an edge of the new simplex. The generation of a new smaller simplex achieves better effect than shrinkage operations, which can reduce the possibility of the simplex degeneration. lk indicates the number of simplexes generation till the kth iteration. km means that the mth simplex generation occurs at the kth iteration. hm means the initial edge length of the mth generated simplex. For the sake of enhancing the efficiency, hm can increase in the process of iterations. However, it is certain that there are decreasing subsequences in hm. εm is the control error, which is a monotonically decreasing sequence in the process of iterations. In order to converge, hm,εm must satisfy limm→∞⁡hm=0 and limm→∞⁡(εm/hm)=0. In the paper, hm,εm are solved by (25) or (26):(25)hm=hm-1λh,εm=εm-1λε,(26)λh=λ˙ht≥tmaxλ¨hρm≤ρnei∨t≤tdec∧hm-1>hsmallλ⃛hothers,λε=λ¨ερm≤ρneiλ⃛ερm>ρnei.The variation ofhm has important influences on the performance of the algorithm. The rate of change of hm,εm must satisfy 1<λ˙h, 0<λ¨h<λ⃛h<1, 0<λ¨ε<λ⃛ε<1 in the paper.Theorem 6.
As for Algorithm3, assuming the sequence zk is bounded, the objective function fx is continuously differentiable and f′x is Lipschitz in any bounded region of Rn; then, ym converges at the stationary point of fx.Proof.
It can be learned from Theorem2 that the search direction ui=xi-x-, i=0,…,n, is positive bases. From the flow of the algorithms, the finite sequence fzk, km≤k<km+1, generated by the inner loop is monotonically decreasing as fzk>fzk+1, km≤k<km+1. Upon the exiting from the inner loop, it is certain fzk+1+εm≥fzk will hold. So fzk≥fzk+1 will be obtained by the algorithms. Similar to the frame of Algorithms  3.1 and 3.2 described in Coope and Price (2000), zk constitutes a quasiframe [11]. Comparing those two algorithms in Coope and Price (2000) with Algorithms 3, 4, and 5 in the paper, the algorithm proposed in the paper can be regarded as a concrete instance of Algorithms  3.1 and 3.2 in Coope and Price (2000). The side vectors xi-x0, i=1,…,n, and x--x0 of the simplex are regarded as positive bases in Coope and Price (2000) [11]. However, xi-x-, i=0,…,n, are taken as positive bases in the paper. Algorithms in Coope and Price (2000) have same convergence properties as algorithms in the paper, whose convergence is proved in Coope and Price (2000) in details [11]. Compared with Algorithm 3, some acceleration methods are adopted in Algorithms 4 and 5, so their convergence is similar to Algorithm 3.


## 3.1.1. Generation of a Simplex

Due to the tendency of deformation and degeneration in the process of the simplex’s iteration, it is usually better that the initial simplex is constructed as an orthogonal one. When one side of the orthogonal simplex gets close to the negative gradient direction-ρk, the simplex will search much rapidly. If Ai=τ1,…,τi and z is any of solution of the linear equations AiTZ=0, then z will be orthogonal to τ1,…,τi. Firstly, we set τ1=-ρk and A1=τ1. If z, that is, any of the solutions of AiTZ=0, is taken as the last column of Ai to generate a new matrix Ai+1, then the orthogonal matrix An=τ1,…,τn can be obtained after repeatedly solving the linear equations n-1 times, where, τ1,…,τn represent n orthogonal vectors. Subsequently, we choose the point whose function value is better from zk+hmτi,zk-hmτi as the vertex of the simplex. In this paper, the flow of simplex generation is shown as Algorithm 1.Algorithm 1.
Step  1. Let τ1=-ρk/ρk, i=1, Sk=zk, A1=τ1.
Step  2. Let(5)Sk=Sk∪zk+hmτifzk+hmτi<fzk-hmτiSk∪zk-hmτiothers.
Step  3. Solve the linear equations AiTZ=0 and choose any of the solutions as z and let (6)τi+1=zz,Ai+1=τ1,…,τi,τi+1.
Step  4. When i<n-1, then i=i+1 and go to Step  2; otherwise exit.Following the procedure, the simplexSk is generated at zk. Edges connecting with zk are orthogonal vectors whose lengths are all hm. Besides, the first edge is parallel to the negative gradient direction -ρk.


## 3.1.2. Center Ray Operation

Theorem 2.
If the center of the simplex is denoted asx-=(1/(n+1))∑i=0nxi and the direction from the center to the ith vertex is denoted as ui=xi-x-, i=0,…,n, then U=ui∣i=0,…,n is a group of positive base.Proof.
Obviously,ui=xi-x-=xi-(1/(n+1))∑j=0nxj=-(1/(n+1))∑j=0,j≠inxj-xi. Assuming there exist constants l1,…,li,…,ln who are not all zero such that ∑i=1nliui=0 is tenable, then ∑i=1nliui=(-1/(n+1))∑i=1nli∑j=0,j≠inxj-xi=(-1/(n+1))∑i=1nnli-∑j=1,j≠inljx0-xi=0. As the side vectors x0-xi, i=1,…,n, are linearly independent, we get nli-∑j=1,j≠inlj=0, i=1,…,n, from which li=0, i=1,…,n is obtained. Therefore, the assumption is contradictory and untenable. Thus, ui,i=1,…,n is linearly independent, and rank⁡U=n is obtained. Therefore, U is a spanning set of Rn. On the other hand, ∑i=0nui=∑i=0nxi-x-=∑i=0nxi-n+1x-=0. In addition, the number of vectors is exactly right n+1. According to literatures [9–11], ui, i=0,…,n are positive bases.
The center ray operation searches along the direction ofui, i=0,…,n by comparing fx- and fxi, i=0,…,n. Firstly, as shown in (7), the first search is along the direction ui, i=0,…,n, where, μi must satisfy (8), and it is advised to calculate μi by (9) or (10). To be specific, when fx->fxi, we let μi=μouter>1, which means search outside the simplex. If fx-≤fxi, it is reasonable we assume better points locate between x- and xi, so we let μi=μinner∈0,1 in order to search inside the simplex:(7)x^ie=x-+μixi-x-,i=0,…,n,(8)μi=μinner∈0,1fx-≤fxiμouter∈1,∞others,(9)μouter=μ˙outer,μinner=μ˙inner,(10)μouter=μouterinit+μouterend-μouterinittmaxt,μinner=μinnerinit+μinnerend-μinnerinittmaxt.
Next, in order to look for a better point, the second detection is attempted according tofx^ie, indicated in (11). The second detection point is located by parameter μq whose scope is determined by comparing function values of x^ie,xi,x-, as indicated in (12):(11)x~ie=x-+μqμixi-x-,i=0,…,n,(12)μq=μe∈1,∞fx^ie<fxi<fx-μoc∈μi-1,1fxi≤fx^ie≤fx-μic∈1,μi-1fx-≤fxiμs∈0,μi-1fxi<fx-≤fx^ie.
Ifx^ie, which is found out by the first search, is the minimum point, then the direction is supposed as the descent direction. Naturally, a better point may be acquired after extending forward along the same direction. So let μq>1 in the case. If expansion point x~ie is better than x^ie, x~ie will be accepted; or else, x^ie will be accepted. When fxi≤fx^e≤fx-∨fx-<fxi, better points are supposed to locate between xi and x^e. Thus, set μi-1<μq<1, and the algorithm goes back for a certain distance to look for a better point x~ie. Then, xi is replaced with the better point from x~ie and x^ie, as a result, a new simplex will be formed. When fxi<fx-≤fx^e, assumption that better points locate between xi and x- is reasonable; thus, let μi-1<μq<0. If x~ie is better than xi, xi is replaced with x~ie. Otherwise, vertex xi still retains in the new simplex. Through the comparison between fxi and fx-, center ray operation can search both inside and outside the simplex, so the algorithm can get close to the extreme point no matter it is inside or outside the simplex. Some illustrations of the center ray operation under 2-dimensional circumstance are shown in Figure 1.
In the process of the center ray, the simplex can become larger or smaller. The adaptive change can reduce calculated amount and accelerate the convergence. If search direction is appropriate, the second search can expand outwards, which can enlarge the simplex and speed up the search. If the extreme point is inside the simplex, the simplex can be automatically downsized by the comparison of function values among the center, the vertex and the detected point. Consequently, the convergence will be accelerated. The search of traditional simplex method relies mainly on reflection operations; unfortunately, its search is only along one direction and only one vertex is updated. Computation of the traditional method is slow and its effect is poor in solving high-dimension problems. In comparison, the center ray operation searches along directions ofn+1 positive bases, and one operation may update several vertexes. As a result, not only the computation speed is improved but also the omission of solution is avoided.Figure 1
Illustration of center ray under 2-dimensional circumstance.


## 3.1.3. Reflection Operation

In the simplex, the center of all vertexes except the worst vertexxn is denoted as x-′, as shown (13). Firstly, like the traditional simplex method, xn is reflected about x-′. As indicated in (14), x^r, namely, the reflection point of xn about x-′, is computed. Different from the traditional simplex method, the algorithm in the paper can search along two inverse reflection directions by comparison between fxn and fx-′, as shown in (15). Besides, like center ray operation, there are two attempts in reflection operation, as shown in (14) and (18), respectively:(13)x-′=∑i=0n-1xi,(14)x^r=x-′+αx-′-xn,(15)α=aforward>0fxn≥fx-′areverse<-1fxn<fx-′.The reflection coefficientα can be computed by either (16) or (17). As shown in (16), the reflection coefficient α is usually constant. When α decreases with iterations as shown in (17), the computation speed can be accelerated. In (17), tmax indicates the preset maximum iterations, αinitforward, αinitreverse are the reflection coefficients when t=0; and αendforward and αendreverse are the reflection coefficients when t=tmax:(16)αforward=α˙forward,αreverse=α˙reverse,(17)αforward=αinitforward-αinitforward-αendforwardtmaxt;0≤αendforward≤αinitforward,αreverse=αinitreverse-αinitreverse-αendreversetmaxt;0≤αinitreverse≤αendreverse.After that, by comparison amongfx^r, fx-′, fxn, x~r is computed by (18) to the second search. When fxn≥fx-′, coefficient αq must satisfy (19); when fxn<fx-′, αq must satisfy (20). The final reflection point is denoted as xr, as shown in (21). If fxr<fxn, xn is replaced with xr to form a new simplex. Otherwise, the reflection will fail. Figure 2 illustrates some of the process of reflection operations under 2-dimensional circumstance:(18)x~r=x-′+αqαx-′-xn,(19)αq=αfe∈1,∞fx^r<fx-′≤fxnαfc∈0,1fx-′≤fx^r<fxnαfs∈-α-1,0fx-′≤fxn<fx^r,(20)αq=αre∈1,∞fx^r<fxn<fx-′αrc∈-α-1,1fxn≤fx^r<fx-′αrs∈0,-α-1fxn<fx-′≤fx^r,(21)xr=x^rfx^r<fx~rx~rothers.Figure 2
Illustration of the reflection operation under 2-dimensional circumstance.Other than the traditional simplex method, both the center ray and the reflection in the paper twice searches along a same direction, likewise, object function is also computed twice. If results of the first detection are better, the algorithm outward expands along the same direction with a longer step to start the second search. Thus, computation speed will be accelerated. Otherwise, the algorithm will try to make certain contraction on the basis of the first search or search along inverse direction. This will improve the accuracy, and what is more important, this also can cause adaptive changes in both the size and the shape of the simplex in the search process. It means to automatically and continuously change the direction and step of the search, which enhance adaptability of the algorithm. One operation in the improved algorithm can achieve effects of several operations of the traditional algorithm, including reflection, expansion, contraction, and shrinkage. As long as the appropriate parameters are chosen, the computation efficiency can be greatly improved.


## 3.1.4. Estimation of Gradient with Simplex

As derivative information is not utilized in the traditional simplex algorithm, the search is slow. Function values of at leastn+1 points are known in the simplex of n-dimensional problems. Therefore, the gradient even the Hessian matrix at each iteration can be estimated by these n+1 points. When one side of the simplex is in accordance with the negative gradient, the simplex may search along the descent direction.If theith side vectors xi-x0i=1,…,n connecting with the optimum point x0 is denoted as vi, then vi=viv→i, where vi=xi-x0 is its side length and v→i is its direction of unit vector, as shown in (22). Since the directional derivative of fx in the direction of v→i at x0 is (∂fx/∂v→i)x=x0=limt→0+⁡((fx0+tv→i-fx0)/t), di in (23) is approximate to (∂fx/∂v→i)x=x0 when the simplex is relatively small:(22)v→i=xi-x0xi-x0,(23)di=fxi-fx0xi-x0≈∂fx∂v→ix=x0.Suppose the gradient offx at x0 is ∇fx0; then, v→iT∇fx0=(∂fx/∂v→i)x=x0≈di. Let the matrix Vk=v→1,v→2,…,v→nT and the vector Dk=d1,d2,…,dnT; then, Vk∇fx0≈Dk, especially, when the simplex is relatively small. As v→1,v→2,…,v→n are linearly independent, (24) can be regarded as the estimated value of the gradient:(24)ρk=Vk-1Dk.In (24), the gradient can be accurately estimated without additionally calculating the function value. Especially, when the simplex is small or approximate to the extreme point, estimation is very accurate. Furthermore, Hessian matrix, Newton direction as well as the conjugate gradient at each iteration can be approximately estimated according to recent two successive gradients in the process of iterations.


## 3.1.5. Flow of the Novel Simplex Method

For each simplex, the algorithm looks for a better point along rays emitting from the centerx- to the vertex xi to replace the current vertex xi. By continuously replacing the former simplex with a new one, iterations will be repeated until the evolution of the simplex stagnates. Once no better points are worked out or the simplex collapses, a new smaller simplex will be generated at the current optimal point in line with the current estimated descent direction and resume iterations again. The flow of the algorithm is shown in Figure 3.Figure 3
The flow of Algorithm3.Algorithm 3.
Step  1. Initialize the parameters l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  4.
Step  4. Set k=k+1, t=t+1.
Step  5. Estimate the gradient ρk of the simplex Sk.
Step  5.1. Compute each side vector v→i by (22) and construct matrix Vk.
Step  5.2. Approximately compute the directional derivative in the side vector by (23) and construct the matrix D=d1,d2,…,dnT.
Step  5.3. Compute ρk with (24).
Step  6. Generate simplex T after center ray operations of simplex Sk.
Step  6.1. Compute x-,fx- and set i=1.
Step  6.2. Compute μi by (9) or (10), and compute x^ie and fx^ie with (7) and (8).
Step  6.3. Determine μq that satisfies (12), and then compute x~ie,fx~ie with (11) to the second detect.
Step  6.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of simplex T.
Step  6.5. If i<n, then i=i+1 and go to Step  6.2; otherwise, go to Step  6.6.
Step  6.6. Sort vertexes of simplex T in ascending order by their function values.
Step  7. Whenfx0<fk-εm, then let fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T and go to Step  3. Otherwise, let fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  8.
Step  8. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, then quit the algorithm; otherwise, go to Step  2 and repeat iterations.In Algorithm3, the mth execution of the inner loop from Step  3 to Step  7 generates the subsequence zk,fk, km≤k<km+1, which satisfies fk≤fk-1-εm<fk-1. As the outer loop from Step  2 to Step  8 repeatedly performs the inner loop, the algorithm generates the infinite sequence zk,fk,ym,gm, which satisfies fk≤fk-1, gm≤gm-1.The center ray only searches alongn+1 fixed directions. If the direction is not descending, the center ray may fail, and the search will stop for the simplex tends to beinfinitely small. For some ill-conditioned problems, even though search directions are descent, it is theoretically feasible to find optimum value; however, this can be realized only when the step is extremely small. Since such extremely small simplex is likely to collapse or the computation is too slow, it is difficult to realize in practice. Although Algorithm 3 is theoretically convergent, sometimes, it actually tends to converge to saddle point and nonoptimum point in case of complex functions. The performance only can be improved through continuous changes of searching directions. A new simplex will be forced to generate when the number of successive execution of ray operations exceeds certain threshold, which is equivalent to changes of the search direction to some degree. Whereby Algorithm 4 is formed with improvements of Algorithm 3, as shown in Figure 4.Figure 4
The flow of Algorithm4.Algorithm 4.
Step  1. Initialize parameters l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. If t<tmax, go to Step  4; otherwise, go to Step  9.
Step  4. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  5.
Step  5. Set k=k+1, t=t+1.
Step  6. Estimate the gradient ρk of the simplex Sk.
Step  6.1. Compute each side vector v→i by (22) and create the matrix Vk.
Step  6.2. Approximately compute the directional derivative in the side vector by (23) to construct the matrix D=d1,d2,…,dnT.
Step  6.3. Compute ρk by (24).
Step  7. Generate simplex T after the application of center ray operations of simplex Sk.
Step  7.1. Compute x-,fx- and let i=1.
Step  7.2. Compute μi by (9) or (10) and compute x^ie and fx^ie with (7) and (8).
Step  7.3. Determine μq which satisfies (12), compute x~ie,fx~ie with (11) to the second search.
Step  7.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of T.
Step  7.5. If i<n, then set i=i+1 and go to Step  7.2; otherwise, go to Step  7.6.
Step  7.6. Sort vertexes of simplex T in ascending order by their function values.
Step  8. When fx0<fk-εm, let fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T and go to Step  3. Otherwise, let fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  9.
Step  9. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, quit the algorithm; otherwise, go to Step  2 and repeat computation.Different from Algorithm3, Algorithm 4 makes restrictions on the times of the repetitive execution of ray operations, denoted as t, which is not allowed to surpass the maximum value tmax. If t≥tmax, even though ray operations are successful, a new simplex at the current optimal point zk is generated according to the current ρk so as to change the search direction.Considering the better performance of reflection operations in the searching process of traditional simplex, in order to further improve the search abilities, a reflection operation will be performed in case of failure of ray operations. Consequently, search directions are more abundant and Algorithm5 is formed, as shown in Figure 5.Figure 5
The flow of Algorithm5.Algorithm 5.
Step  1. Initialize parameters as l1=1, k1=1, m=1, k=1, z1=y1, f1=f(y1), g1=f1,ρ1,ε1.
Step  2. Generate simplex Sk according to zk,hm,-ρk.
Step  2.1. Generate the orthogonal simplex with the side length hm at zk in accordance with Algorithm 1.
Step  2.2. Sort vertexes xi of simplex Sk in ascending order by f(xi) to form Sk=x0,x1,…,xi,…,xn.
Step  2.3. Set t=0.
Step  3. If t<tmax, go to Step  4; otherwise, go to Step  11.
Step  4. Determine whether Sk has degenerated. In case of degeneration, go to Step  2; or else, go to Step  5.
Step  5. Set k=k+1, t=t+1.
Step  6. Estimate the gradient ρk of the simplex Sk.
Step  6.1. Compute each side vector v→i by (22) to construct matrix Vk.
Step  6.2. Approximately compute the directional derivative in the side vector by (23), and construct matrix D=d1,d2,…,dnT.
Step  6.3. Compute ρk by (24).
Step  7. Generate simplex T1 after the application of center ray operations of simplex Sk.
Step  7.1. Compute x-,fx- and set i=1.
Step  7.2. Determine μi by (9) or (10) and compute x^ie and fx^ie by (7) and (8).
Step  7.3. Determine μq which satisfies (12) and compute x~ie,fx~ie by (11) to the second detection.
Step  7.4. If fx->fx^e, choose the best point from x~ie,x^ie as a new vertex of simplex T1 to replace xi; otherwise, choose the best point from xi,x~ie as a new vertex of simplex T1.
Step  7.5. If i<n, then set i=i+1 and go to Step  7.2; otherwise, go to Step  7.6.
Step  7.6. Sort vertexes of simplex T1 in ascending order by their function value.
Step  8. Whenfx0<fk-εm, set fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T1 and go to Step  3. Otherwise, go to Step  9.
Step  9. Form novel simplex T2 after the reflection of simplex T1.
Step  9.1. Compute x-′ by (13).
Step  9.2. The reflection coefficient α is computed by (16) or (17) and x^r and fx^r are computed by (14).
Step  9.3. Determine αq which satisfies (19) or (20), and then solve x~r by (18) and compute fx~r.
Step  9.4. Determine reflection point xr by (21). If fxr<fxn, then xn is replaced with xr to form new simplex T2 and sort vertexes of T2 in ascending order by function value.
Step  10. Whenfx0<fk-εm, set fk+1=fx0, zk+1=x0, lk+1=lk, Sk+1=T2 and go to Step  3. Otherwise, set fk+1=fk, zk+1=zk, gm=fk, ym=zk, lk+1=lk+1, km=k, m=m+1 and go to Step  11.
Step  11. Compute hm+1, εm+1 by (25) or (26). If the result satisfies termination conditions, quit the algorithm; otherwise, go to Step  2 to repeatedly do loop computations.Similarly, as for Algorithm5, Step  3 to Step  10 construct inner loop while Step  2 to Step  11 form the outer loop.Algorithms3, 4, and 5 all have two layers of loops. The sequence zk is composed of optimal points of simplex Sk. The mth execution of the inner loop generates the finite sequence zk,km≤k<km+1, which satisfies fzk≤fzk-1-εm<fzk-1. The outer loop will be repeatedly executed, which may generate the infinite sequence zk satisfying fzk+1≤fzk. That is to say, the sequence fzk generated by the algorithm must be nonincreasing. The sequence ym is composed of best vertexes of the simplex when both ray operations and reflection operations fail, which satisfies fym+1≤fym. Upon the failure of ray and reflection operations, a new smaller orthogonal simplex will be generated. The algorithm can ensure the current minimum point to be one of the vertexes of the new simplex. Besides, it also ensures that the current descent direction is parallel to an edge of the new simplex. The generation of a new smaller simplex achieves better effect than shrinkage operations, which can reduce the possibility of the simplex degeneration. lk indicates the number of simplexes generation till the kth iteration. km means that the mth simplex generation occurs at the kth iteration. hm means the initial edge length of the mth generated simplex. For the sake of enhancing the efficiency, hm can increase in the process of iterations. However, it is certain that there are decreasing subsequences in hm. εm is the control error, which is a monotonically decreasing sequence in the process of iterations. In order to converge, hm,εm must satisfy limm→∞⁡hm=0 and limm→∞⁡(εm/hm)=0. In the paper, hm,εm are solved by (25) or (26):(25)hm=hm-1λh,εm=εm-1λε,(26)λh=λ˙ht≥tmaxλ¨hρm≤ρnei∨t≤tdec∧hm-1>hsmallλ⃛hothers,λε=λ¨ερm≤ρneiλ⃛ερm>ρnei.The variation ofhm has important influences on the performance of the algorithm. The rate of change of hm,εm must satisfy 1<λ˙h, 0<λ¨h<λ⃛h<1, 0<λ¨ε<λ⃛ε<1 in the paper.Theorem 6.
As for Algorithm3, assuming the sequence zk is bounded, the objective function fx is continuously differentiable and f′x is Lipschitz in any bounded region of Rn; then, ym converges at the stationary point of fx.Proof.
It can be learned from Theorem2 that the search direction ui=xi-x-, i=0,…,n, is positive bases. From the flow of the algorithms, the finite sequence fzk, km≤k<km+1, generated by the inner loop is monotonically decreasing as fzk>fzk+1, km≤k<km+1. Upon the exiting from the inner loop, it is certain fzk+1+εm≥fzk will hold. So fzk≥fzk+1 will be obtained by the algorithms. Similar to the frame of Algorithms  3.1 and 3.2 described in Coope and Price (2000), zk constitutes a quasiframe [11]. Comparing those two algorithms in Coope and Price (2000) with Algorithms 3, 4, and 5 in the paper, the algorithm proposed in the paper can be regarded as a concrete instance of Algorithms  3.1 and 3.2 in Coope and Price (2000). The side vectors xi-x0, i=1,…,n, and x--x0 of the simplex are regarded as positive bases in Coope and Price (2000) [11]. However, xi-x-, i=0,…,n, are taken as positive bases in the paper. Algorithms in Coope and Price (2000) have same convergence properties as algorithms in the paper, whose convergence is proved in Coope and Price (2000) in details [11]. Compared with Algorithm 3, some acceleration methods are adopted in Algorithms 4 and 5, so their convergence is similar to Algorithm 3.


## 3.2. Experimental Verification of the Novel Simplex Algorithm

Schwefel and Rosenbrock functions are taken as benchmarks to test the performance of the improved simplex algorithm [12]. In the experiments, the initial point at each computation is all x0=50,…,50T. Once the number of iterations exceeds the predetermined maximum value or hm approaches zero or the solution is obtained, the algorithm will immediately be terminated. μouter and μinner are computed by (9) and, respectively, set μouter=1.2, μinner=0.9, μe=1.25, μoc=0.92, μic=1.06, and μs=0.4. Besides, reflection coefficient α is computed by (16) and set αforward=1.3, αreverse=-2.3, αfe=1.25, αfc=0.7, αfs=-0.4, αre=1.25, αrc=0.72, and αrs=0.15. In the experiments, the other parameters are h0=5,ε0=0.2,ρnei=0.1,λ¨ε=0.5,λ⃛ε=0.8,λ˙h=1.5. Assuming that actual number of iterations is Ω1, the number of function computation is Ω2 and the optimal solution is fx^, then Ω3=x^-x0/Ω1, Ω4=fx^-fx0/Ω1, Ω5=x^-x0/Ω2, Ω6=fx^-fx0/Ω2, Ω7=x^-x∗ can be utilized to evaluate the performance of the algorithm. Smaller Ω1, Ω2, Ω7 and larger Ω3, Ω4, Ω5, Ω6 indicate the performance of the algorithm is better.Price and Coope (2002) has improved NM Simplex method and introduced ghost simplex [13]. The algorithm revised by PRICE is tested in the paper. It is concluded that either Price’s algorithm or NM Simplex method can hardly solve when the dimension of Rosenbrock functions is more than 20. In other words, in solving high-dimension problems, N M Simplex method has presented very poor effect and the algorithm proposed by PRICE also makes little improvement. Taking 5-dimensional Schwefel function as the benchmark, the algorithm is terminated once Ω7=4e-5. Then, the results of this paper are compared with those of Price’s method. Some results are listed in Table 2. In order to verify the performance of the algorithm proposed in the paper in high-dimension optimization problems, 90-dimensional Schwefel function and 45-dimensional Rosenbrock function are solved by Algorithms 4, 5, and 7 (Figure 6). Some results are shown in Tables 3~8. In order to compare effects of all algorithms, average values are totaled in Tables 3~8, as shown in Figure 5. The experiment proves that the effect and performance of improved simplex algorithm proposed in the paper is much better than that of traditional simplex methods.Table 2
Comparison between the algorithm proposed in the paper and Price.
Algorithm
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Price
1.361376
E
-
09
2457
3569
0.045504
55.96256
0.031326
38.5262
Algorithm3
2.161297
E
-
09
602
11226
0.18572
228.4053
0.009959
12.24835
Algorithm4
1.703474
E
-
09
526
11305
0.212554
261.4068
0.00989
12.16276
Algorithm5
1.314375
E
-
09
448
9220
0.249561
306.9196
0.012126
14.91323Table 3
Parts of results of solving 90-dimensional Schwefel function by Algorithm3.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.99
6.96
E
-
08
60683
13153005
0.007817
10178.51
3.61
E
-
05
46.9598
0.000516
0.9
0.98
7.37
E
-
08
53731
11575442
0.008828
11495.46
4.1
E
-
05
53.35973
0.000465
0.9
0.975
8.98
E
-
08
48455
10472012
0.009789
12747.14
4.53
E
-
05
58.98222
0.000593
0.85
0.99
6.5
E
-
08
59400
12882921
0.007986
10398.36
3.68
E
-
05
47.94429
0.000504
0.85
0.98
7.38
E
-
08
49668
10766210
0.00955
12435.82
4.41
E
-
05
57.37047
0.000429
0.85
0.975
8.98
E
-
08
48455
10472012
0.009789
12747.14
4.53
E
-
05
58.98222
0.000593Table 4
Parts of results of solving 45-dimensional Rosenbrock function by Algorithm3.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.995
3.32
E
-
09
299760
30953026
0.001097
88107.5
1.06
E
-
05
853.2641
0.000115
0.9
0.985
3.32
E
-
09
446429
46260909
0.000736
59160.82
7.11
E
-
06
570.9163
0.000115
0.85
0.995
3.32
E
-
09
316215
32404535
0.001039
83522.62
1.01
E
-
05
815.0435
0.000115
0.85
0.985
3.32
E
-
09
454935
47117032
0.000723
58054.68
6.98
E
-
06
560.5426
0.000115
0.8
0.995
3.32
E
-
09
303212
31274881
0.001084
87104.42
1.05
E
-
05
844.483
0.000115
0.8
0.985
3.32
E
-
09
460166
47576852
0.000714
57394.74
6.91
E
-
06
555.1251
0.000115Table 5
Parts of results of solving 90-dimensional Schwefel function by Algorithm4.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.99
8.34
E
-
08
67054
14486323
0.007074
9211.419
3.27
E
-
05
42.63763
0.000567
0.9
0.98
8.78
E
-
08
130898
27575059
0.003624
4718.655
1.72
E
-
05
22.39932
0.000589
0.9
0.975
7.08
E
-
08
194194
40244951
0.002443
3180.647
1.18
E
-
05
15.34758
0.00052
0.85
0.99
8.75
E
-
08
73219
15845655
0.006478
8435.823
2.99
E
-
05
38.97993
0.000584
0.85
0.98
8.52
E
-
08
138163
29067816
0.003433
4470.535
1.63
E
-
05
21.24902
0.00058
0.85
0.975
4.68
E
-
08
195527
40501017
0.002426
3158.963
3.23
E
-
05
15.25054
0.000505Table 6
Parts of results of solving 45-dimensional Rosenbrock function by Algorithm4.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.995
3.32
E
-
09
285070
29702517
0.001153
92647.79
1.11
E
-
05
889.1875
1.72
E
-
05
0.9
0.985
3.32
E
-
09
287837
29884227
0.001142
91757.16
1.10
E
-
05
883.7808
1.72
E
-
05
0.85
0.995
3.32
E
-
09
283558
29541881
0.001159
93141.81
1.11
E
-
05
894.0225
1.72
E
-
05
0.85
0.985
3.32
E
-
09
285806
29672230
0.00115
92409.21
1.11
E
-
05
890.0951
1.72
E
-
05
0.8
0.995
3.32
E
-
09
287742
29952651
0.001142
91787.45
1.10
E
-
05
881.7619
1.72
E
-
05
0.8
0.985
3.32
E
-
09
288202
29917444
0.001141
91640.95
1.10
E
-
05
882.7995
1.72
E
-
05Table 7
Parts of results of solving 90-dimensional Schwefel function by Algorithm5.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.99
4.06
E
-
08
74783
21242674
0.006343
8259.397
2.23
E
-
05
29.07649
0.0004
0.9
0.98
7.08
E
-
08
75482
16388224
0.006284
8182.911
2.89
E
-
05
37.68941
0.000516
0.9
0.975
8.86
E
-
08
81503
17609179
0.00582
7578.402
2.69
E
-
05
35.07617
0.000581
0.85
0.99
7.57
E
-
08
59474
13368291
0.007976
10385.42
3.55
E
-
05
46.20355
0.000542
0.85
0.98
7.50
E
-
08
72768
15772677
0.006519
8488.106
3.01
E
-
05
39.16028
0.000532
0.85
0.975
2.38
E
-
08
84097
18157608
0.00564
7344.644
1.42
E
-
05
34.01673
0.000611Table 8
Parts of results of solving 45-dimensional Rosenbrock function by Algorithm5.
λ
¨
h
λ
⃛
h
f
x
^
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
0.9
0.995
3.32
E
-
09
391063
52265101
0.000841
67536.7
6.29
E
-
06
505.3297
1.72
E
-
05
0.9
0.985
3.32
E
-
09
323520
40051422
0.001016
81636.7
8.21
E
-
06
659.4299
1.72
E
-
05
0.85
0.995
3.32
E
-
09
394208
52604873
0.000834
66997.89
6.25
E
-
06
502.0658
1.72
E
-
05
0.85
0.985
3.33
E
-
09
323368
39892530
0.001016
81675.07
8.24
E
-
06
662.0564
1.72
E
-
05
0.8
0.995
3.34
E
-
09
395380
52745863
0.000831
66799.3
6.23
E
-
06
500.7237
1.72
E
-
05
0.8
0.985
3.33
E
-
09
325750
40212144
0.001009
81077.84
8.17
E
-
06
656.7943
1.72
E
-
05Figure 6
Comparison of different effects of Algorithms3, 4, and 5.In terms of optimization of low-dimensional functions, though the algorithm proposed by the paper and by PRICE can solve such problems, the improved simplex method in the paper demonstrates higher efficiency. It can be known from Table2 that the algorithm proposed in the paper computes much rapidly although the times of object function computation by Algorithms 3, 4, and 5 is more than that of Price’s algorithm. The object function is calculated twice every detection in iteration, so not only Ω2 of the algorithm proposed by the paper may be great than that of other simplex method but also Ω5, Ω6 of algorithm in my paper may be less than that of other simplex method under low dimensional circumstance. In addition, as Price’s algorithm and the traditional simplex algorithm are quite inefficient in high-dimension optimization problems, they can hardly solve any problem as long as the number of dimensions is over 20. In comparison, the algorithm proposed in the paper still behaves more effectively in solving high-dimension problems. It indicates that the improved simplex method, which can update several vertexes at each iteration by center ray operations, is effective and efficient.With respect to the performance, generally, Algorithm5 is the best, while Algorithm 4 is better than Algorithm 3. For both Rosenbrock function and Schwefel function, Algorithm 5 can obtain the smallest Ω1, Ω2 and the largest Ω3, Ω4, Ω5, Ω6; this implies that increasing of search directions can greatly improve the efficiency. Ω1, Ω2 of Algorithm 4 are smaller than those of Algorithm 3, which indicates that the former costs less in computation. In general, Ω3, Ω4, Ω5, Ω6 of Algorithm 3 is smaller than that of Algorithm 4, which reflects that efficiency and speed of Algorithm 4 is higher than that of Algorithm 3.The parameters of the algorithm, especially the parameters controlling changes ofhm have great influence on the performance of the algorithm. According to the experiment, it is usually possible for the algorithm to get good effects when main parameters are within the following ranges: (27)μinner∈0.5,1,μouter∈1.2,2.2,μe∈1.1,1.6,μoc∈0.5,0.8,μic∈1.1,1.8,μs∈0.1,0.5,λ˙h∈1.5,3,λ¨h∈0.5,0.95,λ⃛h∈0.95,0.999,tmax∈6,16,ρnei∈0.1,0.2,hsmall∈0.5,0.05,h0∈0.5,5,λ¨ε∈0.5,0.95,λ⃛ε∈0.8,0.9,ε0∈0.2,0.5,αforward∈0.5,2.5,αreverse∈-1.1,-3,αinitforward∈2,3,αendforward∈0.3,0.6,αinitreverse∈-2.5,-3.5,αendreverse∈-1.1,-1.5,αfe∈1.2,2.5,αfc∈0.2,0.8,αfs∈-2,-0.4,αre∈1.2,2.5,αrc∈0.3,0.9,αrs∈0,0.4.


## 4. Combinatorial Optimization of Aircraft Maintenance Tasks

### 4.1. Analysis of the Combinatorial Optimization of Aircraft Maintenance Tasks

Taking the actual requirements of the combination of aircraft maintenance tasks into consideration, a constraint optimization model of the combination of maintenance tasks is created on the basis of model (3). Algorithms 3, 4, and 5 in the paper can be applied to solve the optimization model, and subsequently the approximate optimal maintenance work packages are obtained [14].


#### 4.1.1. Combinatorial Optimization Model of Aircraft Maintenance Tasks

The basic properties of thejth preventive maintenance task of aircraft includes ATA code, maintenance interval, zone, check and, man-hours, respectively, denoted as ATAj, Periodj, Zonej, Checkj, and Hourj. ATA code is a six-digit code indicating the system, subsystem and components of the aircraft. The last two digits of ATA code are not uniformly standardized, so only the preceding four codes are adopted in the paper. Being interval-valued, maintenance interval refers to the floating range of the maintenance cycle and is denoted as Periodj=P_j,P-j. In terms of ATA standards, the structure and the space of a civil aircraft are usually divided into several zones. Every task is executed in some zones. There are great similarities among maintenance tasks in adjacent zones. In line with MSG-3, there are five types of checks, namely, LU/SV, OP/VC, IN/FC, RS, and DS. If we define set I=LU/SV,OP/VC,IN/FC,RS,DS, then Checkj∈I. Man-hours represents how long time execution of the maintenance task need, which can be utilized to measure the size of the work package. ATAj, Periodj, Zonej, and Checkj are utilized to measure similarities between maintenance tasks. However, roles and influences of the four properties on similarity evaluation are very different. Weights wperiod, wzone, wcheck, and wATA are, respectively, quantified those differences.Maintenance interval is interval-valued. It is necessary to redefine its distance and inner product to adapt to model (3) [15, 16]. Supposing that set V is composed of interval-valued a,b,a,b∈Rn, and s1=a1,b1∈V, s2=a2,b2∈V, k∈R, if addition is defined as s1+s2=a1+a2,b1+b2 and the scalar multiplication is defined as ks1=ka1,b1=ka1,kb1, it can be proved that V is the linear space on R and the zero element is 0,0. In addition, if the inner product is defined as s1,s2=a1,b1,a2,b2=a1Ta2+b1Tb2, V is proved to be the inner product space. Hence, the distance between s1 and s2 can be defined as ds1s22=a1,b1-a2,b2,a1,b1-a2,b2 = a1-a2Ta1-a2+b1-b2Tb1-b2 which is in consistent with the definition of the city block distance [17, 18]. Obviously, it also meets the requirements of metric distance. Consequently, the distance between Periodi and Periodj is defined as dPeriodi,Periodj2=P_i-P_j2+P-i-P-j2.Check merely involves several discrete values. So the measurement of the similarity or distance is usually specified subjectively from expertise. After detailed investigations and comprehensive researches, similarities of those checks are summarized in Table9. Table 9 develops the similarity matrix Iij5×5 of Checki and Checkj; thus, the distance matrix can be obtained as dChecki,Checkj2=1-Iij [19].Table 9
Similarities of checks.
LU/SV
OP/VC
IN/FC
RS
DS
LU/SV
1
0.5
0.5
0.2
0.2
OP/VC
0.5
1
0.6
0.25
0.25
IN/FC
0.5
0.6
1
0.3
0.3
RS
0.2
0.25
0.3
1
0.6
DS
0.2
0.25
0.3
0.6
1Zones are symbolized by digital symbols from 100 to 800. These zone codes can be regarded as digits to measure similarities although they are not numeric values in the real sense. Estimating the distance between two zones by the common numerical computation can approximately meet the requirements of engineering. Suppose that theith task associates with Zonei zones, denoted as Zonei=CZi,1,…,CZi,Zonei, and the jth task involves Zonej zones, expressed as Zonej=CZj,1,…,CZj,Zonej; then, the distance between Zonei and Zonej is defined as dZonei,Zonej2=(1/ZoneiZonej)∑m=1Zonei∑n=1ZonejCZi,m-CZj,n2.Similar to zone codes, it is acceptable in practice to evaluate similarities of ATA codes by numerical computations if they are treated as numerical values. Thus, the distance betweenATAi and ATAj is defined as dATAi,ATAj2=ATAi-ATAj,ATAi-ATAj.Supposing theith maintenance work package is Packagei, man-hours of the package can be computed by (28) if maintenance tasks are executed in sequence. Its man-hours shall be computed by (29) if maintenance tasks are performed concurrently. To realize the successful execution of maintenance tasks, the workload of the maintenance work package must meet the requirement of actual maintenance capability, namely, PHi<PHmax, where PHmax indicates the upper limit of every work package man-hours, which is required by maintenance conditions:(28)PHi=∑j∈PackageiHourj,(29)PHi=maxj∈PackageiHourj.To satisfy requirements of the reliability and economy, when those work packages are executed, maintenance intervals of work packages should be acceptable to all maintenance tasks in the package. Namely, there shall be a public interval among the intersection of all tasks in the package. If the maintenance intervalPPi of the ith maintenance work package Packagei is defined as (30), then P~i-P~i≥0, ∀i is required:(30)PPi=P~i,P~i=∩j∈PackageiP_j,P-j.Correspondingly, the combinatorial optimization model of the maintenance task will be(31)JmU=wperiod∑i=1c∑j=1nuijmP_j2+P-j2-2∑l=1nuilmP_jP_l+P-jP-l∑k=1nuikm+∑p=1n∑q=1nuipmuiqmP_pP_q+P-pP-q∑k=1nuikm2+wzone∑i=1c∑j=1nuijmdZonei,Zonej2+wATA∑i=1c∑j=1nuijmATAj2-2∑l=1nuilmATAjATAl∑k=1nuikm+∑p=1n∑q=1nuipmuiqmATApATAq∑k=1nuikm2+wcheck∑i=1c∑j=1nuijmdChecki,Checkj2s.t.∑i=1cuij=1,∀j=1,…,nuij≥0,∀i=1,…,c,j=1,…,nPHi<PHmax,∀i=1,…,cP~i>P~,∀i=1,…,c.


#### 4.1.2. The Flow of Optimal Combination of Maintenance Tasks

Based on the clustering model, Model (31) is supplemented with two additional constraints. Besides, the constraints are difficult to be explicitly represented by uij,ci. It is not suitable to directly solve the maintenance tasks combination problem as a normal clustering problem. However, maintenance tasks can be clustered first, and then the result of clustering should be examined or checked whether clustering results satisfy constraints. If clustering results do not meet the requirements, we increase the number of clustering and recomputation until the results satisfy the constraint [20]. Its detailed steps are described in Algorithm 7.Algorithm 7.
Step  1. Initialize the parameter and set c=2.
Step  2. Solve clustering model (31).
Step  2.1. Set k=1, κ=1.5, σ1=500.
Step  2.2. By the penalty function, the clustering model is transformed into a unconstrained optimization problem as yU=JmU+σkPU, where PU=∑i=1c∑j=1nmax20,-uij+∑j=1n∑i=1cuij-12, JmU is (31).
Step  2.3. Optimize the unconstrained optimization problem min⁡yU by Algorithms 3 or 4 or 5.
Step  2.4. When σkPU<ε, go to Step  3. Otherwise, set σk+1=κσk, k=k+1 and go to Step  2.3.
Step  3. Tasks in the same class are taken as same maintenance work package. The man-hours PHi of the ith maintenance work package Packagei is totaled by (28) or (29). Summarize the maintenance interval PPi=P~i,P~i by (30).
Step  4. Determine whether PHi<PHmax, P~i>P~, ∀i=1,…,c are satisfied. When PHi<PHmax, P~i>P~, ∀i=1,…,c, it indicates that the current clustering result is a reasonable maintenance task combination. Otherwise, set c=c+1, go to Step  2 and compute again.
Step  5. Analyze and compare maintenance cost between before and after the combination.Step  2 to Step  4 form an inner loop, in which the preset number of classes is increased constantly until a combination of maintenance task satisfies the constraints.


### 4.2. Instances of Combination Optimization of Aircraft Maintenance Tasks

Some scheduled maintenance tasks of certain aircraft are investigated from an airline company, shown in Table10, where the maintenance interval is measured by flight hours. Judging from engineering experience, we set wperiod=0.8, wzone=0.05, wcheck=0.1, and wATA=0.05. PHi is computed by (28), and set PHmax=2. Algorithm 7 is employed to find out optimal combination of scheduled maintenance tasks in Table 10 in order to form maintenance work package. When penalty function is utilized to solve unconstraint clustering problems, set κ=1.5 and σ1=500. Algorithm 5 is applied to solve unconstraint optimization problem, where μouter, μinner are computed by (9) and α is determined with (16). Let other parameters as the following μinner=0.5, μouter=1.4, μe=1.2, μoc=0.7, μic=1.5, μs=0.3, λ˙h=2, λ¨h=0.9, λ⃛h=0.985, tmax=10, ρnei=0.1, hsmall=0.1, αforward=2, αreverse=-2, αfe=2, αfc=0.65, αfs=-1, and αre=2, αrc=0.6. The result of Algorithm 7 is tabulated in Table 11, from which it can be known that seventeen maintenance tasks are combined into thirteen maintenance work packages.Table 10
Some scheduled maintenance tasks of an airline.
Number
ATA
Maintenance interval
Check
Zone
Man-hours
Lower boundary
Upper boundary
1
2410
338
394
IN/FC
210
0.1
2
3201
240
278
RS
734; 744
0.05
3
3201
221
258
RS
734; 744
0.05
4
5381
1543
1779
IN/FC
121; 122
0.3
5
5780
2746
3176
IN/FC
510; 520; 530; 540; 550; 560; 570
0.3
6
3243
1231
1468
IN/FC
310
0.05
7
3301
1228
1417
OP/VC
200
0.02
8
5280
2359
2760
IN/FC
820; 830; 840
0.4
9
2304
895
1014
OP/VC
211
0.1
10
2702
12331
14161
LU/SV
133; 570
0.5
11
2702
7920
9242
LU/SV
134; 672
0.4
12
2703
16116
18200
IN/FC
210; 572; 672
0.4
13
2704
23413
27600
IN/FC
133; 134; 325
0.4
14
2709
9788
11258
IN/FC
333; 334; 343; 344
1
15
2704
28497
32543
IN/FC
131; 134; 325
0.7
16
2704
24885
29716
IN/FC
131; 134; 325
0.4
17
2713
11671
13050
DS
133; 134; 211; 541; 542; 543; 544; 553
0.5Table 11
Optimal work package after the combination of maintenance tasks.
Number
Maintenance interval
Man-hours of package
Tasks
Lower boundary
Upper boundary
1
338
394
0.1
1
2
240
258
0.1
2, 3
3
1543
1779
0.3
4
4
2746
2760
0.7
5, 8
5
1231
1417
0.07
6, 7
6
895
1014
0.1
9
7
12331
14161
0.5
10
8
7920
9242
0.4
11
9
16116
18200
0.4
12
10
23413
27600
0.4
13
11
9788
11258
1
14
12
28497
29716
1.1
15, 16
13
11671
13050
0.5
17Suppose costs of every maintenance includes four parts: labor costs, device costs, halt costs, and fixed costs, which are, respectively, represented asCostl, Costd, Costh, and Costf. Specifically, the fixed cost typically refers to the depreciation of maintenance tools and fixed assets and the management fees. Costh and Costf are always the same for any repair. Costd, which is determined by maintenance items, is the total device costs of all maintenance items repaired and is irrelative with the way of combination. Thereby, maintenance cost before the combination is Costb=Flabour∑j=117Hourj+∑j=117Devicej+17Costh+17Costf. Flabour refers to the labor costs per unit time. Devicej indicates the device cost of the jth maintenance item. After the combination, the maintenance cost will be Costa=Flabour∑j=113PHj+∑j=117Devicej+13Costh+13Costf. The maintenance cost is compared between before and after combination is tabulated in Table 12. As Costb-Costa=4Costh+4Costf, it can be concluded that the combination helps saving cost and improving efficiency. Halt cost will cause great cost, which accounts for a large proportion in the total maintenance cost. Therefore, optimal combination of maintenance tasks can greatly increase profits of the airlines.Table 12
Optimal work package after the combination of maintenance tasks.
Before combination
After combination
Cost saved
Labor cost
F
labour
∑
j
=
1
17
Hou
r
j
F
labour
∑
j
=
1
13
P
H
j
0
Device cost
∑
j
=
1
17
Devic
e
j
∑
j
=
1
17
Devic
e
j
0
Halt cost
17
Cost
h
13
Cost
h
4
Cost
h
Fixed cost
17
Cost
f
13
Cost
f
4
Cost
f


## 4.1. Analysis of the Combinatorial Optimization of Aircraft Maintenance Tasks

Taking the actual requirements of the combination of aircraft maintenance tasks into consideration, a constraint optimization model of the combination of maintenance tasks is created on the basis of model (3). Algorithms 3, 4, and 5 in the paper can be applied to solve the optimization model, and subsequently the approximate optimal maintenance work packages are obtained [14].


### 4.1.1. Combinatorial Optimization Model of Aircraft Maintenance Tasks

The basic properties of thejth preventive maintenance task of aircraft includes ATA code, maintenance interval, zone, check and, man-hours, respectively, denoted as ATAj, Periodj, Zonej, Checkj, and Hourj. ATA code is a six-digit code indicating the system, subsystem and components of the aircraft. The last two digits of ATA code are not uniformly standardized, so only the preceding four codes are adopted in the paper. Being interval-valued, maintenance interval refers to the floating range of the maintenance cycle and is denoted as Periodj=P_j,P-j. In terms of ATA standards, the structure and the space of a civil aircraft are usually divided into several zones. Every task is executed in some zones. There are great similarities among maintenance tasks in adjacent zones. In line with MSG-3, there are five types of checks, namely, LU/SV, OP/VC, IN/FC, RS, and DS. If we define set I=LU/SV,OP/VC,IN/FC,RS,DS, then Checkj∈I. Man-hours represents how long time execution of the maintenance task need, which can be utilized to measure the size of the work package. ATAj, Periodj, Zonej, and Checkj are utilized to measure similarities between maintenance tasks. However, roles and influences of the four properties on similarity evaluation are very different. Weights wperiod, wzone, wcheck, and wATA are, respectively, quantified those differences.Maintenance interval is interval-valued. It is necessary to redefine its distance and inner product to adapt to model (3) [15, 16]. Supposing that set V is composed of interval-valued a,b,a,b∈Rn, and s1=a1,b1∈V, s2=a2,b2∈V, k∈R, if addition is defined as s1+s2=a1+a2,b1+b2 and the scalar multiplication is defined as ks1=ka1,b1=ka1,kb1, it can be proved that V is the linear space on R and the zero element is 0,0. In addition, if the inner product is defined as s1,s2=a1,b1,a2,b2=a1Ta2+b1Tb2, V is proved to be the inner product space. Hence, the distance between s1 and s2 can be defined as ds1s22=a1,b1-a2,b2,a1,b1-a2,b2 = a1-a2Ta1-a2+b1-b2Tb1-b2 which is in consistent with the definition of the city block distance [17, 18]. Obviously, it also meets the requirements of metric distance. Consequently, the distance between Periodi and Periodj is defined as dPeriodi,Periodj2=P_i-P_j2+P-i-P-j2.Check merely involves several discrete values. So the measurement of the similarity or distance is usually specified subjectively from expertise. After detailed investigations and comprehensive researches, similarities of those checks are summarized in Table9. Table 9 develops the similarity matrix Iij5×5 of Checki and Checkj; thus, the distance matrix can be obtained as dChecki,Checkj2=1-Iij [19].Table 9
Similarities of checks.
LU/SV
OP/VC
IN/FC
RS
DS
LU/SV
1
0.5
0.5
0.2
0.2
OP/VC
0.5
1
0.6
0.25
0.25
IN/FC
0.5
0.6
1
0.3
0.3
RS
0.2
0.25
0.3
1
0.6
DS
0.2
0.25
0.3
0.6
1Zones are symbolized by digital symbols from 100 to 800. These zone codes can be regarded as digits to measure similarities although they are not numeric values in the real sense. Estimating the distance between two zones by the common numerical computation can approximately meet the requirements of engineering. Suppose that theith task associates with Zonei zones, denoted as Zonei=CZi,1,…,CZi,Zonei, and the jth task involves Zonej zones, expressed as Zonej=CZj,1,…,CZj,Zonej; then, the distance between Zonei and Zonej is defined as dZonei,Zonej2=(1/ZoneiZonej)∑m=1Zonei∑n=1ZonejCZi,m-CZj,n2.Similar to zone codes, it is acceptable in practice to evaluate similarities of ATA codes by numerical computations if they are treated as numerical values. Thus, the distance betweenATAi and ATAj is defined as dATAi,ATAj2=ATAi-ATAj,ATAi-ATAj.Supposing theith maintenance work package is Packagei, man-hours of the package can be computed by (28) if maintenance tasks are executed in sequence. Its man-hours shall be computed by (29) if maintenance tasks are performed concurrently. To realize the successful execution of maintenance tasks, the workload of the maintenance work package must meet the requirement of actual maintenance capability, namely, PHi<PHmax, where PHmax indicates the upper limit of every work package man-hours, which is required by maintenance conditions:(28)PHi=∑j∈PackageiHourj,(29)PHi=maxj∈PackageiHourj.To satisfy requirements of the reliability and economy, when those work packages are executed, maintenance intervals of work packages should be acceptable to all maintenance tasks in the package. Namely, there shall be a public interval among the intersection of all tasks in the package. If the maintenance intervalPPi of the ith maintenance work package Packagei is defined as (30), then P~i-P~i≥0, ∀i is required:(30)PPi=P~i,P~i=∩j∈PackageiP_j,P-j.Correspondingly, the combinatorial optimization model of the maintenance task will be(31)JmU=wperiod∑i=1c∑j=1nuijmP_j2+P-j2-2∑l=1nuilmP_jP_l+P-jP-l∑k=1nuikm+∑p=1n∑q=1nuipmuiqmP_pP_q+P-pP-q∑k=1nuikm2+wzone∑i=1c∑j=1nuijmdZonei,Zonej2+wATA∑i=1c∑j=1nuijmATAj2-2∑l=1nuilmATAjATAl∑k=1nuikm+∑p=1n∑q=1nuipmuiqmATApATAq∑k=1nuikm2+wcheck∑i=1c∑j=1nuijmdChecki,Checkj2s.t.∑i=1cuij=1,∀j=1,…,nuij≥0,∀i=1,…,c,j=1,…,nPHi<PHmax,∀i=1,…,cP~i>P~,∀i=1,…,c.


### 4.1.2. The Flow of Optimal Combination of Maintenance Tasks

Based on the clustering model, Model (31) is supplemented with two additional constraints. Besides, the constraints are difficult to be explicitly represented by uij,ci. It is not suitable to directly solve the maintenance tasks combination problem as a normal clustering problem. However, maintenance tasks can be clustered first, and then the result of clustering should be examined or checked whether clustering results satisfy constraints. If clustering results do not meet the requirements, we increase the number of clustering and recomputation until the results satisfy the constraint [20]. Its detailed steps are described in Algorithm 7.Algorithm 7.
Step  1. Initialize the parameter and set c=2.
Step  2. Solve clustering model (31).
Step  2.1. Set k=1, κ=1.5, σ1=500.
Step  2.2. By the penalty function, the clustering model is transformed into a unconstrained optimization problem as yU=JmU+σkPU, where PU=∑i=1c∑j=1nmax20,-uij+∑j=1n∑i=1cuij-12, JmU is (31).
Step  2.3. Optimize the unconstrained optimization problem min⁡yU by Algorithms 3 or 4 or 5.
Step  2.4. When σkPU<ε, go to Step  3. Otherwise, set σk+1=κσk, k=k+1 and go to Step  2.3.
Step  3. Tasks in the same class are taken as same maintenance work package. The man-hours PHi of the ith maintenance work package Packagei is totaled by (28) or (29). Summarize the maintenance interval PPi=P~i,P~i by (30).
Step  4. Determine whether PHi<PHmax, P~i>P~, ∀i=1,…,c are satisfied. When PHi<PHmax, P~i>P~, ∀i=1,…,c, it indicates that the current clustering result is a reasonable maintenance task combination. Otherwise, set c=c+1, go to Step  2 and compute again.
Step  5. Analyze and compare maintenance cost between before and after the combination.Step  2 to Step  4 form an inner loop, in which the preset number of classes is increased constantly until a combination of maintenance task satisfies the constraints.


## 4.1.1. Combinatorial Optimization Model of Aircraft Maintenance Tasks

The basic properties of thejth preventive maintenance task of aircraft includes ATA code, maintenance interval, zone, check and, man-hours, respectively, denoted as ATAj, Periodj, Zonej, Checkj, and Hourj. ATA code is a six-digit code indicating the system, subsystem and components of the aircraft. The last two digits of ATA code are not uniformly standardized, so only the preceding four codes are adopted in the paper. Being interval-valued, maintenance interval refers to the floating range of the maintenance cycle and is denoted as Periodj=P_j,P-j. In terms of ATA standards, the structure and the space of a civil aircraft are usually divided into several zones. Every task is executed in some zones. There are great similarities among maintenance tasks in adjacent zones. In line with MSG-3, there are five types of checks, namely, LU/SV, OP/VC, IN/FC, RS, and DS. If we define set I=LU/SV,OP/VC,IN/FC,RS,DS, then Checkj∈I. Man-hours represents how long time execution of the maintenance task need, which can be utilized to measure the size of the work package. ATAj, Periodj, Zonej, and Checkj are utilized to measure similarities between maintenance tasks. However, roles and influences of the four properties on similarity evaluation are very different. Weights wperiod, wzone, wcheck, and wATA are, respectively, quantified those differences.Maintenance interval is interval-valued. It is necessary to redefine its distance and inner product to adapt to model (3) [15, 16]. Supposing that set V is composed of interval-valued a,b,a,b∈Rn, and s1=a1,b1∈V, s2=a2,b2∈V, k∈R, if addition is defined as s1+s2=a1+a2,b1+b2 and the scalar multiplication is defined as ks1=ka1,b1=ka1,kb1, it can be proved that V is the linear space on R and the zero element is 0,0. In addition, if the inner product is defined as s1,s2=a1,b1,a2,b2=a1Ta2+b1Tb2, V is proved to be the inner product space. Hence, the distance between s1 and s2 can be defined as ds1s22=a1,b1-a2,b2,a1,b1-a2,b2 = a1-a2Ta1-a2+b1-b2Tb1-b2 which is in consistent with the definition of the city block distance [17, 18]. Obviously, it also meets the requirements of metric distance. Consequently, the distance between Periodi and Periodj is defined as dPeriodi,Periodj2=P_i-P_j2+P-i-P-j2.Check merely involves several discrete values. So the measurement of the similarity or distance is usually specified subjectively from expertise. After detailed investigations and comprehensive researches, similarities of those checks are summarized in Table9. Table 9 develops the similarity matrix Iij5×5 of Checki and Checkj; thus, the distance matrix can be obtained as dChecki,Checkj2=1-Iij [19].Table 9
Similarities of checks.
LU/SV
OP/VC
IN/FC
RS
DS
LU/SV
1
0.5
0.5
0.2
0.2
OP/VC
0.5
1
0.6
0.25
0.25
IN/FC
0.5
0.6
1
0.3
0.3
RS
0.2
0.25
0.3
1
0.6
DS
0.2
0.25
0.3
0.6
1Zones are symbolized by digital symbols from 100 to 800. These zone codes can be regarded as digits to measure similarities although they are not numeric values in the real sense. Estimating the distance between two zones by the common numerical computation can approximately meet the requirements of engineering. Suppose that theith task associates with Zonei zones, denoted as Zonei=CZi,1,…,CZi,Zonei, and the jth task involves Zonej zones, expressed as Zonej=CZj,1,…,CZj,Zonej; then, the distance between Zonei and Zonej is defined as dZonei,Zonej2=(1/ZoneiZonej)∑m=1Zonei∑n=1ZonejCZi,m-CZj,n2.Similar to zone codes, it is acceptable in practice to evaluate similarities of ATA codes by numerical computations if they are treated as numerical values. Thus, the distance betweenATAi and ATAj is defined as dATAi,ATAj2=ATAi-ATAj,ATAi-ATAj.Supposing theith maintenance work package is Packagei, man-hours of the package can be computed by (28) if maintenance tasks are executed in sequence. Its man-hours shall be computed by (29) if maintenance tasks are performed concurrently. To realize the successful execution of maintenance tasks, the workload of the maintenance work package must meet the requirement of actual maintenance capability, namely, PHi<PHmax, where PHmax indicates the upper limit of every work package man-hours, which is required by maintenance conditions:(28)PHi=∑j∈PackageiHourj,(29)PHi=maxj∈PackageiHourj.To satisfy requirements of the reliability and economy, when those work packages are executed, maintenance intervals of work packages should be acceptable to all maintenance tasks in the package. Namely, there shall be a public interval among the intersection of all tasks in the package. If the maintenance intervalPPi of the ith maintenance work package Packagei is defined as (30), then P~i-P~i≥0, ∀i is required:(30)PPi=P~i,P~i=∩j∈PackageiP_j,P-j.Correspondingly, the combinatorial optimization model of the maintenance task will be(31)JmU=wperiod∑i=1c∑j=1nuijmP_j2+P-j2-2∑l=1nuilmP_jP_l+P-jP-l∑k=1nuikm+∑p=1n∑q=1nuipmuiqmP_pP_q+P-pP-q∑k=1nuikm2+wzone∑i=1c∑j=1nuijmdZonei,Zonej2+wATA∑i=1c∑j=1nuijmATAj2-2∑l=1nuilmATAjATAl∑k=1nuikm+∑p=1n∑q=1nuipmuiqmATApATAq∑k=1nuikm2+wcheck∑i=1c∑j=1nuijmdChecki,Checkj2s.t.∑i=1cuij=1,∀j=1,…,nuij≥0,∀i=1,…,c,j=1,…,nPHi<PHmax,∀i=1,…,cP~i>P~,∀i=1,…,c.


## 4.1.2. The Flow of Optimal Combination of Maintenance Tasks

Based on the clustering model, Model (31) is supplemented with two additional constraints. Besides, the constraints are difficult to be explicitly represented by uij,ci. It is not suitable to directly solve the maintenance tasks combination problem as a normal clustering problem. However, maintenance tasks can be clustered first, and then the result of clustering should be examined or checked whether clustering results satisfy constraints. If clustering results do not meet the requirements, we increase the number of clustering and recomputation until the results satisfy the constraint [20]. Its detailed steps are described in Algorithm 7.Algorithm 7.
Step  1. Initialize the parameter and set c=2.
Step  2. Solve clustering model (31).
Step  2.1. Set k=1, κ=1.5, σ1=500.
Step  2.2. By the penalty function, the clustering model is transformed into a unconstrained optimization problem as yU=JmU+σkPU, where PU=∑i=1c∑j=1nmax20,-uij+∑j=1n∑i=1cuij-12, JmU is (31).
Step  2.3. Optimize the unconstrained optimization problem min⁡yU by Algorithms 3 or 4 or 5.
Step  2.4. When σkPU<ε, go to Step  3. Otherwise, set σk+1=κσk, k=k+1 and go to Step  2.3.
Step  3. Tasks in the same class are taken as same maintenance work package. The man-hours PHi of the ith maintenance work package Packagei is totaled by (28) or (29). Summarize the maintenance interval PPi=P~i,P~i by (30).
Step  4. Determine whether PHi<PHmax, P~i>P~, ∀i=1,…,c are satisfied. When PHi<PHmax, P~i>P~, ∀i=1,…,c, it indicates that the current clustering result is a reasonable maintenance task combination. Otherwise, set c=c+1, go to Step  2 and compute again.
Step  5. Analyze and compare maintenance cost between before and after the combination.Step  2 to Step  4 form an inner loop, in which the preset number of classes is increased constantly until a combination of maintenance task satisfies the constraints.


## 4.2. Instances of Combination Optimization of Aircraft Maintenance Tasks

Some scheduled maintenance tasks of certain aircraft are investigated from an airline company, shown in Table10, where the maintenance interval is measured by flight hours. Judging from engineering experience, we set wperiod=0.8, wzone=0.05, wcheck=0.1, and wATA=0.05. PHi is computed by (28), and set PHmax=2. Algorithm 7 is employed to find out optimal combination of scheduled maintenance tasks in Table 10 in order to form maintenance work package. When penalty function is utilized to solve unconstraint clustering problems, set κ=1.5 and σ1=500. Algorithm 5 is applied to solve unconstraint optimization problem, where μouter, μinner are computed by (9) and α is determined with (16). Let other parameters as the following μinner=0.5, μouter=1.4, μe=1.2, μoc=0.7, μic=1.5, μs=0.3, λ˙h=2, λ¨h=0.9, λ⃛h=0.985, tmax=10, ρnei=0.1, hsmall=0.1, αforward=2, αreverse=-2, αfe=2, αfc=0.65, αfs=-1, and αre=2, αrc=0.6. The result of Algorithm 7 is tabulated in Table 11, from which it can be known that seventeen maintenance tasks are combined into thirteen maintenance work packages.Table 10
Some scheduled maintenance tasks of an airline.
Number
ATA
Maintenance interval
Check
Zone
Man-hours
Lower boundary
Upper boundary
1
2410
338
394
IN/FC
210
0.1
2
3201
240
278
RS
734; 744
0.05
3
3201
221
258
RS
734; 744
0.05
4
5381
1543
1779
IN/FC
121; 122
0.3
5
5780
2746
3176
IN/FC
510; 520; 530; 540; 550; 560; 570
0.3
6
3243
1231
1468
IN/FC
310
0.05
7
3301
1228
1417
OP/VC
200
0.02
8
5280
2359
2760
IN/FC
820; 830; 840
0.4
9
2304
895
1014
OP/VC
211
0.1
10
2702
12331
14161
LU/SV
133; 570
0.5
11
2702
7920
9242
LU/SV
134; 672
0.4
12
2703
16116
18200
IN/FC
210; 572; 672
0.4
13
2704
23413
27600
IN/FC
133; 134; 325
0.4
14
2709
9788
11258
IN/FC
333; 334; 343; 344
1
15
2704
28497
32543
IN/FC
131; 134; 325
0.7
16
2704
24885
29716
IN/FC
131; 134; 325
0.4
17
2713
11671
13050
DS
133; 134; 211; 541; 542; 543; 544; 553
0.5Table 11
Optimal work package after the combination of maintenance tasks.
Number
Maintenance interval
Man-hours of package
Tasks
Lower boundary
Upper boundary
1
338
394
0.1
1
2
240
258
0.1
2, 3
3
1543
1779
0.3
4
4
2746
2760
0.7
5, 8
5
1231
1417
0.07
6, 7
6
895
1014
0.1
9
7
12331
14161
0.5
10
8
7920
9242
0.4
11
9
16116
18200
0.4
12
10
23413
27600
0.4
13
11
9788
11258
1
14
12
28497
29716
1.1
15, 16
13
11671
13050
0.5
17Suppose costs of every maintenance includes four parts: labor costs, device costs, halt costs, and fixed costs, which are, respectively, represented asCostl, Costd, Costh, and Costf. Specifically, the fixed cost typically refers to the depreciation of maintenance tools and fixed assets and the management fees. Costh and Costf are always the same for any repair. Costd, which is determined by maintenance items, is the total device costs of all maintenance items repaired and is irrelative with the way of combination. Thereby, maintenance cost before the combination is Costb=Flabour∑j=117Hourj+∑j=117Devicej+17Costh+17Costf. Flabour refers to the labor costs per unit time. Devicej indicates the device cost of the jth maintenance item. After the combination, the maintenance cost will be Costa=Flabour∑j=113PHj+∑j=117Devicej+13Costh+13Costf. The maintenance cost is compared between before and after combination is tabulated in Table 12. As Costb-Costa=4Costh+4Costf, it can be concluded that the combination helps saving cost and improving efficiency. Halt cost will cause great cost, which accounts for a large proportion in the total maintenance cost. Therefore, optimal combination of maintenance tasks can greatly increase profits of the airlines.Table 12
Optimal work package after the combination of maintenance tasks.
Before combination
After combination
Cost saved
Labor cost
F
labour
∑
j
=
1
17
Hou
r
j
F
labour
∑
j
=
1
13
P
H
j
0
Device cost
∑
j
=
1
17
Devic
e
j
∑
j
=
1
17
Devic
e
j
0
Halt cost
17
Cost
h
13
Cost
h
4
Cost
h
Fixed cost
17
Cost
f
13
Cost
f
4
Cost
f


## 5. Conclusions and Discussions

Firstly, the improved fuzzy C-means clustering model is presented, which can be directly solved. Considering FCM model has many stationary points, the optimization method which needs derivatives is not applicable to solve it. Based on the positive bases theory, center ray operation is introduced and a novel nonlinear simplex optimization method is developed. Besides, its convergences are roughly analyzed. As a direct optimization method, it can be utilized to solve the clustering model effectively. Finally, the improved FCM model and the novel optimization method are employed to optimally combine maintenance tasks and exhibit good effects.However, there are still several problems to be further discussed. (1) The number of independent variables contained in the improved clustering model is still large. In case of a large number of samples, it is still inconvenient to solve FCM directly; while the clustering model, which is easier to be solved, needs to be further studied. (2) The simplex algorithm proposed in the paper is much more effective than the traditional algorithm. However, sometimes it is likely that the improved method fails to achieve good result because the searching direction can not be adjusted timely. From this view, the algorithm has to be improved. The most prominent point is that a certain rotation operation should be taken in the simplex algorithm to change the search directions. (3)hm has great influences on the performance of the improved algorithm proposed in the paper. Although it is proved to be feasible that hm is taken as a decreasing exponential function, the efficiency is still to be improved. Various factors can be associated to find out other way to control changes of hm so as to improve the performance. (4) The experimental results show that the algorithm proposed in the paper can rapidly approach to the nearby region of the extreme point from the far initial point. However, it takes a long time to converge to an accurate solution within the neighbourhood of the extreme point. For this reason, how to improve local convergence performance of the algorithm should be an important research issue. (5) For the optimal combination of maintenance tasks, maintenance tasks are simplified and only partial properties are selected to measure similarities. Although it is proved to be practical in most engineering problems, the result is acceptable but not optimal due to the ignorance of details of the execution of maintenance tasks such as maintenance tools and maintenance orders. A more economical combination will be achieved if complete and comprehensive factors can be taken into consideration to determine similarities of maintenance tasks.



---

*Source: 169310-2015-06-28.xml*