PSO-constrained optimization

    PSO.m
    clear all
    close all
    clc
    
    %define the details of this problem
    nVar = 2;
    ub = [10 10];
    lb = [-10 -10];
    fobj = @ObjectiveFunction;
    
    %define the PSO's parameters
    noP = 30;
    maxIter = 500;
    wMax = 0.9;
    wMin = 0.2;
    c1 = 2;
    c2 = 2;
    vMax = (ub-lb) .* 0.2;
    vMin = -vMax;
    
    %initialise the particles
    for k = 1 : noP
    Swarm.Particles(k).X = (ub-lb) .* rand(1,nVar) + lb;
    Swarm.Particles(k).V = zeros(1,nVar);
    Swarm.Particles(k).PBEST.X = zeros(1,nVar);
    Swarm.Particles(k).PBEST.O = inf;
    
    Swarm.GBEST.X = zeros(1,nVar);
    Swarm.GBEST.O = inf;
    end
    
    %main loop
    for t = 1 : maxIter
    %calculate the objective value
    for k= 1 : noP
        currentX = Swarm.Particles(k).X;
        Swarm.Particles(k).O = fobj(currentX);
        if Swarm.Particles(k).O < Swarm.Particles(k).PBEST.O   %update the PBEST
            Swarm.Particles(k).PBEST.X = currentX;
            Swarm.Particles(k).PBEST.O = Swarm.Particles(k).O; 
        end
        if Swarm.Particles(k).O < Swarm.GBEST.O                %update the GBEST
            Swarm.GBEST.X = currentX;
            Swarm.GBEST.O = Swarm.Particles(k).O;
        end
    end
    
    %update the X and V vectors
    w = wMax - t .* ((wMax-wMin)/maxIter);
    
    for k = 1 : noP
        Swarm.Particles(k).V = w .* Swarm.Particles(k).V + c1 .* rand(1,nVar) .* (Swarm.Particles(k).PBEST.X-Swarm.Particles(k).X)...
                                                         + c2 .* rand(1,nVar) .* (Swarm.GBEST.X-Swarm.Particles(k).X);
        index1 = find(Swarm.Particles(k).V > vMax);   %check velocities
        index2 = find(Swarm.Particles(k).V < vMin);
        Swarm.Particles(k).V(index1) = vMax(index1);
        Swarm.Particles(k).V(index2) = vMin(index2);
        Swarm.Particles(k).X = Swarm.Particles(k).X + Swarm.Particles(k).V;
        index1 = find(Swarm.Particles(k).X > ub);     %check positions
        index2 = find(Swarm.Particles(k).V < lb);
        Swarm.Particles(k).X(index1) = ub(index1);
        Swarm.Particles(k).X(index2) = lb(index2);
    end
    
    outmsg = ['Iteration#' , num2str(t) , 'Swarm.GBEST.O = ' , num2str(Swarm.GBEST.O)];
    disp(outmsg);
    cgCurve(t) = Swarm.GBEST.O;
    end
    
    semilogy(cgCurve);

    ObjectiveFunction.m
    function[o] = ObjectiveFunction(x)
    
    %List of your constraints
    const1 = x(2) <= 3.2 || x(2) >= 6.4;
    const2 = (x(1)^2 + x(2)^2) >= 1;
    const3 = x(1) ~= x(2);
    
    if const1 == 1 && const2 == 1 && const3 == 1
    o = sum(x.^2);
    else
    o = sum(x.^2) + 200;
    end
    
    end

    DrawLandscape.m
    %Visualise the landscape
    fobj = @ObjectiveFunction
    x = -10 : 0.05 : 10;
    y = -10 : 0.05 : 10;
    [x_new , y_new] = meshgrid(x,y);     %生成绘制3D图形所需的网格数据
    for i = 1 : size(x_new , 1)       
    for j =1 : size(x_new , 2)
        currentX = [ x_new(i,j), y_new(i,j) ];
        o(i,j) = fobj(currentX);
    end
    end
    
    surfc(x_new , y_new , o)
    shading interp
    camlight
    
    xlabel('x_1')
    ylabel('x_2')

    Iteration#491Swarm.GBEST.O = 1
    Iteration#492Swarm.GBEST.O = 1
    Iteration#493Swarm.GBEST.O = 1
    Iteration#494Swarm.GBEST.O = 1
    Iteration#495Swarm.GBEST.O = 1
    Iteration#496Swarm.GBEST.O = 1
    Iteration#497Swarm.GBEST.O = 1
    Iteration#498Swarm.GBEST.O = 1
    Iteration#499Swarm.GBEST.O = 1
    Iteration#500Swarm.GBEST.O = 1
    >> Swarm.GBEST.X
    
    ans =
    
    0.8900   -0.4559