Complex algorithm and matlab code
1. Complex algorithm
Complex algorithm is an optimization algorithm, which does not necessarily require the derivability of objective function or constraint function when optimizing non-linear problems. This method is widely used in low-dimensional engineering optimization problems.
The basic idea is to use the statistical information of the complex vertices (the worst point and the best point as well as the shape center), and to gradually replace the worst solution in the complex and form a new complex by iteration, until the complex shrinks to the specified size before terminating the algorithm.
The main steps are:
Step 1: Initialize the complex, then calculate the objective function value of each vertex of the complex, and determine the best point Xl and the worst point Xh in the complex.
Step 2: Find out the centroid Xc corresponding to all points except the worst point in the current complex.
Step 3: Judge the feasibility of Xc centroid. If the centroid is not in the feasible region, the boundary of the design variable is redefined according to the best point and centroid of the current complex and turned to step 1. If the centroid is in the feasible region, step 4 is continued.
Step 4: Map the worst point Xh in the compound type according to the formula: Xr = Xc + alpha*(Xc - Xh); alpha is generally 1.3.
Step 5: Then determine whether the mapping point Xr is in the feasible region or, if not, move to step 8. Otherwise proceed to step 6.
Step 6: Compare the objective function values of the mapping point Xr and the worst point Xh (when the mapping point Xr is in the feasible region), and if the objective function values of the mapping point are better, replace the worst point with the mapping point. Then turn to step 1.
Step 7: If the objective function value of the mapping point is not better than the function value of the worst point, first determine whether alpha is less than a specified good number (delta), if it means that the search fails, then redefine the boundary of the design variable according to the best point and center of the current complex and turn to step 1. If alpha is not less than delta, turn to step 8.
Step 8: Generate T random points on the annular hypersphere, and then determine whether there are feasible points: if there are, replace the worst point in the current composite type, then turn to step 1, if not, reduce alpha proportionally: alpha= alpha * 0.5, and then turn to step 7.
2. matlab code
function [Dopti eopti Popti] = ComplexMethod()
%Generate 6 vertices to form a complex
D = 400 + unifrnd (0,1,1,6).*(1600 - 400);
e = 1.3 + unifrnd (0,1,1,6).*(2.0 - 1.3);
P = 3 + unifrnd (0,1,1,6).*(10 - 3);
Fitness = zeros(1,6);
plan = zeros(6,4);%Save all 6 vertices and corresponding objective function values: D e P Fitness
for i = 1:6
point = [D(i) e(i) P(i)];
Fitness(i) = ComputeScore(point);
plan(i,:) = [record Fitness(i)];
end
[worstPoint,bestPoint,xingxin] =updataInfo(plan);
alpha = 1.3;
while (~JudgeIsEnd(plan))
if JudgeInRange(xingxin)%If the centroid is within the feasible range
Xr = ComputeReflectPoint(xingxin,plan,alpha); %Mapping the worst point to get the mapping point
if JudgeInRange(Xr) %If the reflection point Xr Within the feasible region
while(1)
if Xr(:,4) < worstPoint(:,4) % Reflection point is better than worst point in complex
[~,dim] = max(plan(:,4));
plan(dim,:) = Xr(:,4);%Renewal Complex
[worstPoint,bestPoint,xingxin] = updataInfo(plan);
goon = true;
break;
else % Reflection points are worse than the worst points in a complex.
alpha = 0.5*alpha;
if alpha < 10^(-5) %alpha To a certain extent
plan = regenerate(xingxin,bestPoint);
[worstPoint,bestPoint,xingxin] = updataInfo(plan);
goon = true;
break;
else %Calculating Spherical Random Points
Xr = ComputeReflectPoint(xingxin,plan,alpha);
end
end
end
if goon == true
goon = flase;
continue;%Re-evaluation of the new complex
end
else %The mapping point is not in the feasible region, and the centroid has been ensured in the feasible region. (S8)
%while
Xr = ObtainReflectPoint(xingxin,alpha,Xr);%Newly generated Xr It must be a point in the feasible domain.
if Xr(:,4) < worstPoint(:,4) % Reflection point is better than worst point in complex
[~,dim] = max(plan(:,4));
plan(dim,:) = Xr(:,4);%Renewal Complex
[worstPoint,bestPoint,xingxin] = updataInfo(plan);
goon = true;
continue;
else % Reflection points are worse than the worst points in a complex.
while (1)
alpha = alpha*0.5;
if alpha < 10^(-5) %alpha To a certain extent
plan = regenerate(xingxin,bestPoint);
[worstPoint,bestPoint,xingxin] = updataInfo(plan);
goon = true;
break;
else %Spherical randomization is used to find feasible points and replace the worst points of the current composite type.
Vr = alpha.* (xingxin - worstPoint);
Xrand = findGlobalPoint(Vr,xingxin,plan);
if JudgeInRange(Xrand)
temp = [Xrand ComputeScore(Xrand)];
[~,dim] = max(plan(:,4));
plan(dim,:) = temp;
[worstPoint,bestPoint,xingxin] = updataInfo(plan);
goon = true;
break;
else
continue;
end
end
end
if goon == true
goon = flase;
continue;
end
end
end
else %When the centroid is not in the feasible region, the complex is regenerated
plan = regenerate(xingxin,bestPoint);
[worstPoint,bestPoint,xingxin] = updataInfo(plan);
continue;
end
end
[~,dim] = min(plan(:,4));
optimal = plan(dim,:);
Dopti = optimal(1);
eopti = optimal(2);
Popti = optimal(3);
function [worstPoint,bestPoint,xingxin] = updataInfo(plan)
xingxin = ComputeXingxin(plan);
[~,dim] = max(plan(:,4));
worstPoint = plan(dim,:);
[~,dim1] = min(plan(:,4));
bestPoint = plan(dim1,:);
function Xrand = findGlobalPoint(Xr,xingxin,alpha)
D = Xr(1);
e = Xr(2);
P = Xr(3);
ro = rand;
beta = 10/180*pi;
Dv = 400 + unifrnd (0,1,1,1).*(1600 - 400);
ev = 1.3 + unifrnd (0,1,1,1).*(2.0 - 1.3);
Pv = (Dv*D + ev*e)/(P);
Dn = Dv/sqrt(Dv^2 + ev^2 + Pv^2);
en = ev/sqrt(Dv^2 + ev^2 + Pv^2);
Pn = Pv/sqrt(Dv^2 + ev^2 + Pv^2);
Vn = [Dn en Pn];
Xrand = xingxin(1,1:3) + alpha*norm(Xr)*cos(ro*beta).* Vn + alpha*tan(ro.beta).*(Xr(1,1:3));
function Xr = ObtainReflectPoint(xingxin,alpha,Xr)
while(~JudgeInRange(Xr))
alpha = alpha/2;
Xr = xingxin + alpha.*(xingxin - Xworst);
end
function plan = regenerate(xingxin,bestPoint)
%Using New Boundary to Generate New Complex
plan = zeros(6,4);1.3
for i = 1:6
if xingxin(1) > bestPoint(1)
D = xingxin(1) + unifrnd (0,1,1,1).*(xingxin(1) - bestPoint(1));
else
D = bestPoint(1) + unifrnd (0,1,1,1).*(bestPoint(1) - xingxin(1));
end
if D > 1600 || D< 400
D = 400 + unifrnd (0,1,1,1).*(1600 - 400);
end
if xingxin(2) > bestPoint(2)
e = bestPoint(2) + unifrnd (0,1,1,1).*(xingxin(2) - bestPoint(2));
else
e = xingxin(2) + unifrnd (0,1,1,1).*(bestPoint(2) - xingxin(2));
end
if e > 2 || e< 1.3
e = 1.3 + unifrnd (0,1,1,1).*(2 - 1.3);
end
if xingxin(3) > bestPoint(3)
P = bestPoint(3) + unifrnd (0,1,1,1).*(xingxin(3) - bestPoint(3));
else
P = xingxin(3) + unifrnd (0,1,1,1).*(bestPoint(3) - xingxin(3));
end
if P > 10 || P< 3
P = 3 + unifrnd (0,1,1,1).*(10 - 3);
end
record = [D e P];
Fitness = ComputeScore(record);
plan(i,:) = [D e P Fitness];
end
function ret = JudgeIsEnd(plan)
plan = plan(:,1:3);
Xmean = sum(plan)/6;
%ret =false;
for i = 1:6
temp = Xmean - plan(i,:);
if norm(temp) >= 10^(-3)
ret = false;
return;
end
end
ret = true;
function Xr = ComputeReflectPoint(xingxin,plan,alpha)
[~,dim] = max(plan(:,4));
Xworst = plan(dim,:);
[~,index] = sort(plan(:,4),'descend');
Xworst2 = plan(index(2),:);
if alpha > 10^(-5)
Xr = xingxin + alpha.*(xingxin - Xworst);
else
Xr = xingxin + alpha.*(xingxin - Xworst2);
end
function ret = JudgeInRange(xingxin)
D = xingxin(1);
e = xingxin(2);
P = xingxin(3);
ret = false;
if (D >= 400) && (D <= 1600) && (e >= 1.3) && (e <= 2.0) && (P >= 3 && P <= 10)
ret = true;
end
function xingxin = ComputeXingxin(plan)
[~,dim] = max(plan(:,4));%Worst point
xingxin = (sum(plan) - plan(dim,:))./5;
xingxin(1,4) = ComputeScore(xingxin);