1 简介
2 部分代码
%% MOICA Main Function
close all
clc; clear
tic;
global myGlob;
myGlob = 0;
%% Problem Statement
ProblemParams.CostFuncName = 'BenchmarkFunction'; % cost function
ProblemParams.CostFuncExtraParams = 30; % change function number according to test problem number in BenchmarkFunction.m file
ProblemParams.NPar = 30; % Number of optimization variables of objective function
ProblemParams.VarMin = 0; % lower bound
ProblemParams.VarMax = 1; % upper bound
ProblemParams.M = 3; % number of objective functions
switch ProblemParams.CostFuncExtraParams
case 13
ProblemParams.VarMin=zeros(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
case 14
ProblemParams.VarMin=zeros(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
case 15
ProblemParams.VarMin=zeros(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
case 16
ProblemParams.VarMin=zeros(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
ProblemParams.VarMin(2:end)=ProblemParams.VarMin(2:end)-5;
ProblemParams.VarMax(2:end)=ProblemParams.VarMax(2:end)*5;
case 17
ProblemParams.VarMin=zeros(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
case 18
ProblemParams.NPar = 1;
ProblemParams.VarMin=-1000;
ProblemParams.VarMax=1000;
case 19
ProblemParams.NPar = 3;
ProblemParams.VarMin = repmat(-5,1,ProblemParams.NPar);
ProblemParams.VarMax = repmat(5,1,ProblemParams.NPar);
case 20
ProblemParams.NPar = 3;
ProblemParams.VarMin = repmat(-4,1,ProblemParams.NPar);
ProblemParams.VarMax = repmat(4,1,ProblemParams.NPar);
case 21
ProblemParams.VarMin=-1*ones(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1)=0;
case 22
ProblemParams.VarMin=-1*ones(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1)=0;
case 23
ProblemParams.VarMin=zeros(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
case 24
ProblemParams.VarMin=-2*ones(1,ProblemParams.NPar);
ProblemParams.VarMax=2*ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1)=0;
ProblemParams.VarMax(1)=1;
case 25
ProblemParams.VarMin=-1*ones(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1)=0;
case 26
ProblemParams.VarMin=-1*ones(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1)=0;
case 27
ProblemParams.VarMin=-1*ones(1,ProblemParams.NPar);
ProblemParams.VarMax=ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1)=0;
case 28
ProblemParams.VarMin = -2*ones(1,ProblemParams.NPar);
ProblemParams.VarMax = 2*ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1) = 0;
ProblemParams.VarMin(2) = 0;
ProblemParams.VarMax(1) = 1;
ProblemParams.VarMax(2) = 1;
ProblemParams.M = 3;
case 29
ProblemParams.VarMin = -2*ones(1,ProblemParams.NPar);
ProblemParams.VarMax = 2*ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1) = 0;
ProblemParams.VarMin(2) = 0;
ProblemParams.VarMax(1) = 1;
ProblemParams.VarMax(2) = 1;
ProblemParams.M = 3;
case 30
ProblemParams.VarMin = -2*ones(1,ProblemParams.NPar);
ProblemParams.VarMax = 2*ones(1,ProblemParams.NPar);
ProblemParams.VarMin(1) = 0;
ProblemParams.VarMin(2) = 0;
ProblemParams.VarMax(1) = 1;
ProblemParams.VarMax(2) = 1;
ProblemParams.M = 3;
otherwise
error('Invalid function number is used.');
end
ProblemParams.SearchSpaceSize = ProblemParams.VarMax - ProblemParams.VarMin;
%% Algorithmic Parameter Setting
AlgorithmParams.NumOfCountries = 100; % Number of initial countries.
AlgorithmParams.NumOfInitialImperialists = 8; % Number of Initial Imperialists.
AlgorithmParams.NumOfAllColonies = AlgorithmParams.NumOfCountries - AlgorithmParams.NumOfInitialImperialists;
AlgorithmParams.NumOfDecades = 2e+9;
AlgorithmParams.NumOfFunEvals = 25000;
AlgorithmParams.RevolutionRate = 0.3; % Revolution is the process in which the socio-political characteristics of a country change suddenly.
AlgorithmParams.UnitingThreshold = 0.02; % The percent of Search Space Size, which enables the uniting process of two Empires.
AlgorithmParams.ImperialistPercentage = 0.3; % The maximum percent of Imperialists that can be contained in the Empire.
%% Creation of Initial Empires
InitialCountries = GenerateNewCountry(AlgorithmParams.NumOfCountries , ProblemParams);
% Calculates the cost of each country. The less the cost is, the more is the power.
InitialCost = feval(ProblemParams.CostFuncName,InitialCountries,ProblemParams.CostFuncExtraParams);
[InitialCost, SortInd] = NonDominationSort(InitialCost,ProblemParams.M); % apply non-domination sorting...
InitialCost = InitialCost(:,1:ProblemParams.M);
InitialCountries = InitialCountries(SortInd,:); % Sort the population with respect to their cost.
Empires = CreateInitialEmpires(InitialCountries,InitialCost,AlgorithmParams, ProblemParams);
Solutions = [];
SolutionPositions = [];
ArchInd = 1;
for Decade = 1:AlgorithmParams.NumOfDecades
for ii = 1:numel(Empires)
Empires(ii) = AssimilateColonies(Empires(ii),ProblemParams, Empires);
Empires(ii) = RevolveColonies(Empires(ii),AlgorithmParams,ProblemParams);
%% New Cost Evaluation
Empires(ii).ColoniesCost = feval(ProblemParams.CostFuncName,Empires(ii).ColoniesPosition,ProblemParams.CostFuncExtraParams);
%% Empire Possession (****** Power Possession, Empire Possession) % exchange Imperialist & Best Colony positions if any...
Empires(ii) = PossesEmpire(Empires(ii), ProblemParams, AlgorithmParams);
%% Computation of Total Cost for Empires
empirePop = [Empires(ii).ImperialistCost
Empires(ii).ColoniesCost];
[TotalCosts, ~] = NonDominationSort(empirePop,ProblemParams.M);
FirstFront = find(TotalCosts(:,ProblemParams.M+1) == 1);
Empires(ii).TotalCost = size(FirstFront, 1);
end
%uncomment below line to see the decrease of number of empires in each
%generation...
%numberOfEmpires = numel(Empires)
%% Uniting Similiar Empires
Empires = UniteSimilarEmpires(Empires,AlgorithmParams,ProblemParams);
%% Imperialistic Competition
Empires = ImperialisticCompetition(Empires);
%display number of function evaluations performed...
FunEvals = myGlob
%% Get Non-dominated Solutions
AllCosts = [];zeros(AlgorithmParams.NumOfCountries, ProblemParams.M);
AllPositions = [];zeros(AlgorithmParams.NumOfCountries, ProblemParams.NPar);
AllCostsIndex = 1;
for u=1:numel(Empires)
AllCosts(AllCostsIndex:AllCostsIndex + size(Empires(u).ImperialistCost,1) - 1,:) = Empires(u).ImperialistCost;
AllPositions(AllCostsIndex:AllCostsIndex + size(Empires(u).ImperialistPosition,1) - 1,:) = Empires(u).ImperialistPosition;
AllCostsIndex = AllCostsIndex + size(Empires(u).ImperialistCost,1);
end
tempArr = [];
tempArrPositions = [];
tempInd = 1;
for y=1:size(AllCosts,1)
if isempty(find(sum(ismember(tempArr,AllCosts(y,:)),2) == ProblemParams.M))
tempArr(tempInd,:) = AllCosts(y,:);
tempArrPositions(tempInd, :) = AllPositions(y,:);
tempInd = tempInd + 1;
end
end
[tempArr, NewIn] = NonDominationSort(tempArr,ProblemParams.M); % apply non-domination sorting...
tempArrPositions = tempArrPositions(NewIn,:);
Solutions(ArchInd:ArchInd+size(tempArr(find(tempArr(:,ProblemParams.M+1) == 1),1:ProblemParams.M),1)-1,:) = tempArr(find(tempArr(:,ProblemParams.M+1) == 1),1:ProblemParams.M);
SolutionPositions(ArchInd:ArchInd+size(tempArr(find(tempArr(:,ProblemParams.M+1) == 1),1:ProblemParams.M),1)-1,:) = tempArrPositions(1:size(tempArr(find(tempArr(:,ProblemParams.M+1) == 1),1:ProblemParams.M),1),:);
ArchInd = ArchInd + size(tempArr(find(tempArr(:,ProblemParams.M+1) == 1),1:ProblemParams.M),1);
% Terminate if Max Fun Evals reached...
if myGlob > AlgorithmParams.NumOfFunEvals
break;
end
end % End of Algorithm
Solutions = [SolutionPositions Solutions];
[Solutions, I] = unique(Solutions, 'rows');
Solutions = Get_Non_Dominated_Solutions(Solutions,ProblemParams.M);
% [Solutions, ~] = NonDominationSort(Solutions,ProblemParams.M);
% Solutions = Solutions(find(Solutions(:,end-1) == 1),ProblemParams.NPar+1:ProblemParams.NPar + ProblemParams.M);
dlmwrite('Solutions',Solutions,'delimiter','\t','precision', '%.12f');
f = Solutions;
if ProblemParams.M == 2
uf = load('ZDT4.pf');
plot(uf(:,1),uf(:,2),'r*');
hold on;
plot(Solutions(:,1),Solutions(:,2),'*');
title('MOICA Two Functions');
xlabel('f(x_1)');
ylabel('f(x_2)');
else
uf = load('UF10.pf');
plot3(uf(:,1),uf(:,2),uf(:,3),'r*');
hold on;
plot3(Solutions(:,1),Solutions(:,2),Solutions(:,3),'*');
title('MOICA Three Functions');
xlabel('f(x_1)');
ylabel('f(x_2)');
zlabel('f(x_3)');
end
3 仿真结果
4 参考文献
[1]曲倩雯. 基于ICA和GA混合算法的装配序列规划研究. Diss. 山东大学, 2016.
5 MATLAB代码与数据下载地址
见博客主页
