function [ Mse_max ] = mse_max_computation( b, dmodel )
%UNTITLED3 此处显示有关此函数的摘要
%   此处显示详细说明

LB=[b(1)-3 b(2)-3]; %lower bounds of variables
UB=[b(1)+3 b(2)+3]; %upper bounds of variables
% pso parameters values
m=1; % number of variables，搜索空间维数
n=100; % population size，代表100个粒子
wmax=0.9; % inertia weight，惯性权重最大值
wmin=0.4; % inertia weight，惯性权重最小值
c1=2; % acceleration factor，局部加速因子
c2=2; % acceleration factor，全局加速因子
% pso main program----------------------------------------------------start
maxite=1000; % set maximum number of iteration，最大迭代次数

maxrun=1; % set maximum number of runs need to be，最大运行次数
% 添加目标函数---------start
% load('data_VS_4')
% S=A(:,1:2);
% Y=A(:,3);
% theta=[2 2]; lob = [1e-1 1e-1]; upb = [10 10];
% [dmodel, perf] = dacefit(S, Y, @regpoly0, @corrgauss, theta, lob, upb)

for run=1:maxrun
    % pso initialization----------------------------------------------start
    for i=1:n% i=1，代表第一个粒子，生成初始种群的坐标阵
        for j=1:m% j=1，代表第一个维度
            x0(i,j)=LB(j)+rand()*(UB(j)-LB(j));% 产生第i个粒子的第j维度的坐标值
        end
    end
    x=x0; % initial population，初始种群
    v=0.1*x0; % initial velocity，初始速度
    for i=1:n% 共n个初始粒子
        [y, dy, mse] = predictor(x0(i,:), dmodel);
        f0(i,1) = -mse;
%         f0(i,1) = -predictor(x0(i,:), dmodel);% 需输入dmodel，即目标函数
%         f0(i,1)=sum_motion_error(x0(i,:));% n*1,设计变量初值代入目标函数得到目标初值，将第i个粒子的坐标阵代入计算，返回适应度值
    end
    [fmin0,index0]=min(f0);% 从一系列的适应度值中选取最小的适应度值和相应的粒子
    pbest=x0; % initial pbest，初始局部最优解，即为初始粒子的坐标阵
    gbest=x0(index0,:); % initial gbest，初始全局最优解，即在初始粒子群中适应度函数值最小所对应的粒子坐标阵
    % pso initialization------------------------------------------------end
    % pso algorithm---------------------------------------------------start
    ite=1;% 迭代次数初值
    tolerance=1;% 误差容忍度初值
    while ite<=maxite && tolerance>10^-12% 如果迭代次数未达到上限且没达到误差容忍度
        w=wmax-(wmax-wmin)*ite/maxite; % update inertial weight，更新惯性权重
        % pso velocity updates
        for i=1:n% 外层循环为第i个粒子
            for j=1:m% 内层循环为第i个粒子对应的第j个坐标值，共有m个坐标值
                v(i,j)=w*v(i,j)+c1*rand()*(pbest(i,j)-x(i,j))...
                    +c2*rand()*(gbest(1,j)-x(i,j));% 对第i个粒子的第j维度速度进行更新，循环m次后，完成第i个粒子的所有维度的速度更新
            end
        end
        % pso position update
        for i=1:n
            for j=1:m
                x(i,j)=x(i,j)+v(i,j);% 对第i个粒子的第j维度位置(即位置分量)进行更新
            end
        end
        % handling boundary violations
        for i=1:n
            for j=1:m
                if x(i,j)<LB(j)
                    x(i,j)=LB(j);% 如果粒子坐标阵超过下边界，则将下边界值赋给粒子
                elseif x(i,j)>UB(j)
                    x(i,j)=UB(j);% 如果粒子坐标阵超过上边界，则将上边界值赋给粒子
                end
            end
        end
        % evaluating fitness
        for i=1:n
            [y, dy, mse] = predictor(x0(i,:), dmodel);
            f(i,1) = -mse;
%             f(i,1)=-predictor(x0(i,:), dmodel);
%             f(i,1)=sum_motion_error(x(i,:));% 计算新的适应度函数值
        end
        
        % updating pbest and fitness
        for i=1:n
            if f(i,1)<f0(i,1)
                pbest(i,:)=x(i,:);% 更新局部最优粒子
                f0(i,1)=f(i,1);% 更新适应度函数值
            end
        end
        [fmin,index]=min(f0); % finding out the best particle，返回新粒子群中最小的适应度函数值和对应的粒子索引值
        ffmin(ite,run)=fmin; % storing best fitness，存储迭代最优的适应度函数值到第ite行run列
        ffite(run)=ite; % storing iteration count，存储迭代次数
        hold on;
        % updating gbest and best fitness
        if fmin<fmin0
            gbest=pbest(index,:);% 从局部最优粒子群中，设定适应度函数值最小者为全局最优粒子
            fmin0=fmin;% 将新的最小适应度函数值赋给fmin0
        end
        % calculating tolerance
        if ite>100
            tolerance=abs(ffmin(ite-100,run)-fmin0);% 
        end
        % displaying iterative results
%         if ite==1
%             fprintf('Iteration Best particle Objective fun\n');
%         end
%         fprintf('%8g %8g %10.5f\n',ite,index,fmin0);
        ite=ite+1;
    end
    % pso algorithm-----------------------------------------------------end
    gbest;
    [y, dy, mse] = predictor(gbest, dmodel);
    fvalue= -mse;
%     fvalue=-predictor(gbest,dmodel);
    fff(run)=fvalue;
    rgbest(run,:)=gbest;
   
    
end
% pso main program------------------------------------------------------end
[bestfun,bestrun]=min(fff);
best_variables=rgbest(bestrun,:);
Mse_max =best_variables;

end

