Curve fitting using particle swarm optimization

Dear All I have a code for curve fitting using PSO but it's only for one variable how I can change it to multivariable inputs. clc; clear; close all;

%% Problem Definition

global NFE; NFE=0;

model.x=[0 1 1 4]; model.y=[1 0 2 2]; model.fhat=@(x,a) polyval(a,x); % model.fhat=@(x,a) a(1)+a(2)/(1+a(3)*x)+a(4)*x;

CostFunction=@(a) MyCost(a,model); % Cost Function

nVar=2; % Number of Decision Variables

VarSize=[1 nVar]; % Size of Decision Variables Matrix

VarMin=-10; % Lower Bound of Variables VarMax= 10; % Upper Bound of Variables

%% PSO Parameters

MaxIt=500; % Maximum Number of Iterations

nPop=100; % Population Size (Swarm Size)

w=1; % Inertia Weight wdamp=0.99; % Inertia Weight Damping Ratio c1=2; % Personal Learning Coefficient c2=2; % Global Learning Coefficient

% % Constriction Coefficients % phi1=2.05; % phi2=2.05; % phi=phi1+phi2; % chi=2/(phi-2+sqrt(phi^2-4*phi)); % w=chi; % Inertia Weight % wdamp=1; % Inertia Weight Damping Ratio % c1=chi*phi1; % Personal Learning Coefficient % c2=chi*phi2; % Global Learning Coefficient

% Velocity Limits VelMax=0.1*(VarMax-VarMin); VelMin=-VelMax;

mu=0.05;

%% Initialization

empty_particle.Position=[]; empty_particle.Cost=[]; empty_particle.Sol=[]; empty_particle.Velocity=[]; empty_particle.Best.Position=[]; empty_particle.Best.Cost=[]; empty_particle.Best.Sol=[];

particle=repmat(empty_particle,nPop,1);

GlobalBest.Cost=inf;

for i=1:nPop

    % Initialize Position
    particle(i).Position=unifrnd(VarMin,VarMax,VarSize);

    % Initialize Velocity
    particle(i).Velocity=zeros(VarSize);

    % Evaluation
    [particle(i).Cost particle(i).Sol]=CostFunction(particle(i).Position);

    % Update Personal Best
    particle(i).Best.Position=particle(i).Position;
    particle(i).Best.Cost=particle(i).Cost;
    particle(i).Best.Sol=particle(i).Sol;

    % Update Global Best
    if particle(i).Best.Cost<GlobalBest.Cost

        GlobalBest=particle(i).Best;

    end

end

BestCost=zeros(MaxIt,1);

nfe=zeros(MaxIt,1);

%% PSO Main Loop

for it=1:MaxIt

    for i=1:nPop

        % Update Velocity
        particle(i).Velocity = w*particle(i).Velocity ...
            +c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...
            +c2*rand(VarSize).*(GlobalBest.Position-particle(i).Position);

        % Apply Velocity Limits
        particle(i).Velocity = max(particle(i).Velocity,VelMin);
        particle(i).Velocity = min(particle(i).Velocity,VelMax);

        % Update Position
        particle(i).Position = particle(i).Position + particle(i).Velocity;

        % Velocity Mirror Effect
        IsOutside=(particle(i).Position<VarMin | particle(i).Position>VarMax);
        particle(i).Velocity(IsOutside)=-particle(i).Velocity(IsOutside);

        % Apply Position Limits
        particle(i).Position = max(particle(i).Position,VarMin);
        particle(i).Position = min(particle(i).Position,VarMax);

        % Evaluation
        [particle(i).Cost particle(i).Sol] = CostFunction(particle(i).Position);

        NewSol.Position=Mutate(particle(i).Position,mu,VarMin,VarMax);
        [NewSol.Cost NewSol.Sol]=CostFunction(NewSol.Position);
        if NewSol.Cost<=particle(i).Cost
            particle(i).Position=NewSol.Position;
            particle(i).Cost=NewSol.Cost;
            particle(i).Sol=NewSol.Sol;
        end

        % Update Personal Best
        if particle(i).Cost<particle(i).Best.Cost

            particle(i).Best.Position=particle(i).Position;
            particle(i).Best.Cost=particle(i).Cost;
            particle(i).Best.Sol=particle(i).Sol;

            % Update Global Best
            if particle(i).Best.Cost<GlobalBest.Cost

                GlobalBest=particle(i).Best;

            end

        end

    end

    NewSol.Position=Mutate(GlobalBest.Position,mu,VarMin,VarMax);
    [NewSol.Cost NewSol.Sol]=CostFunction(NewSol.Position);
    if NewSol.Cost<=GlobalBest.Cost
        GlobalBest=NewSol;
    end

    BestCost(it)=GlobalBest.Cost;

    nfe(it)=NFE;

    disp(['Iteration ' num2str(it) ': NFE = ' num2str(nfe(it)) ', Best Cost = ' num2str(BestCost(it))]);

    w=w*wdamp;

    figure(1);
    PlotSolution(GlobalBest.Position,model);

end

%% Results

figure; plot(nfe,BestCost,'LineWidth',2); xlabel('NFE'); ylabel('Best Cost');

function [z sol]=MyCost(a,model)

    global NFE;
    if isempty(NFE)
        NFE=0;
    end

    NFE=NFE+1;

    x=model.x;
    y=model.y;
    fhat=model.fhat;

    yhat=zeros(size(x));
    for i=1:numel(x)
        yhat(i)=fhat(x(i),a);
    end

    e=y-yhat;

    z=sum(e.^2);

    sol.yhat=yhat;
    sol.e=e;

end

function xnew=Mutate(x,mu,VarMin,VarMax)

    nVar=numel(x);

    nmu=ceil(mu*nVar);

    j=randsample(nVar,nmu);

    sigma=0.1*(VarMax-VarMin);

    xnew=x;
    xnew(j)=x(j)+sigma*randn(size(j));

    xnew=max(xnew,VarMin);
    xnew=min(xnew,VarMax);

end function PlotSolution(a,model)

    x=model.x;
    y=model.y;
    fhat=model.fhat;

    xmin=min(x);
    xmax=max(x);
    dx=xmax-xmin;
    xmin=xmin-0.1*dx;
    xmax=xmax+0.1*dx;

    xx=linspace(xmin,xmax,100);
    yy=zeros(size(xx));
    for i=1:numel(xx)
        yy(i)=fhat(xx(i),a);
    end

    plot(x,y,'ro','MarkerSize',10,'MarkerFaceColor','y');
    hold on;
    plot(xx,yy,'k','LineWidth',2);
    legend('Train Data','Fitted Curve');
    xlabel('x');
    ylabel('y');
    hold off;
    grid on;

end

Answer 1

[John D'Errico]

Answer: this is literally impossible to do, for us to show you how to change this mess of code to do something significantly different. DON'T DO IT. Don't write your own optimization tools. You will end up with buggy code that nobody can use, modify, or debug; including you.

Instead, use existing toolboxes to solve the problem. Find PSO tools on the File Exchange. Or use the optimization tools supplied by MathWorks.

[sadegh]

Thank you, can you add some link please.