【优化求解】基于帝国主义竞争算法ICA求解多目标问题Matlab代码

秒速五厘米 2022-08-28 01:46 214阅读 0赞

## 1 简介 ##

![img][]

![img][img 1]

![img][img 2]

![img][img 3]

![img][img 4]

## 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 仿真结果 ##

![d7500bd96ebef34f6248cda4218dbad5.png][]

## 4 参考文献 ##

\[1\]曲倩雯. 基于ICA和GA混合算法的装配序列规划研究. Diss. 山东大学, 2016.

![图片][293db08c386dae29e57ec92e2a5ee0bc.png]

[img]: /images/20220828/da8749c3e208480380446dd04b96e1f0.png
[img 1]: /images/20220828/fae3a9e4430c494f92ee6775e66bdc33.png
[img 2]: /images/20220828/589c176d25914a6f99ac7eb4e399a805.png
[img 3]: /images/20220828/5235bc4bba1442349707279a4c9e0703.png
[img 4]: /images/20220828/1d2ed27745034c89bef81d7e26c61bf9.png
[d7500bd96ebef34f6248cda4218dbad5.png]: /images/20220828/3dab614a4ca143f4ab8149ec8a9bd69e.png
[293db08c386dae29e57ec92e2a5ee0bc.png]: /images/20220828/7886728478f94e669df54e463706bf66.png