基于粒子群算法的非合作博弈的matlab程序

clc;
clear;
%%
%设定均衡点初值                                                                                                                                                                                                       
x0=0;
y0=0;
%%
%各博弈参与者或联盟依次进行独立优化决策
%此处选用粒子群算法作为优化算法
%此处以双目标优化问题为例，从博弈的角度来求解最小值，fx=（x-1）^2+（x-y）^2fy=（y-3）^2+（x-y）^2

%参数初始化
c1=3;                %个体学习因子
c2=3;                %社会学习因子
w_max=0.9;           %设置最大惯性权重为0.9
w_min=0.4;           %设置最小惯性权重为0.4

nums_zibianliang=1;   %每个函数有1个自变量
max_die_dai=50;      %迭代次数设置为50
size_zhong_qun=50;   %粒子种群规模设置为50

x=20.*rands（size_zhong_qunnums_zibianliang）;    %初始化粒子x的位置
y=20.*rands（size_zhong_qunnums_zibianliang）;    %初始化粒子y的位置

vx=3.*rands（size_zhong_qunnums_zibianliang）;    %初始化粒子x的飞翔速度
vy=3.*rands（size_zhong_qunnums_zibianliang）;    %初始化粒子y的飞翔速度

%定义适应度函数
fitness_1=inline（‘（x-1）^2+（x-y）^2‘‘x‘‘y‘）;    %适应度函数fitness_1
for i=1:size_zhong_qun
	fx（i:）=fitness_1（x（i:）y0）;	
end

%这个地方逻辑上存在点问题，需要改一下
fitness_2=inline（‘（y-3）^2+（x-y）^2‘‘x‘‘y‘）;    %适应度函数fitness_2
for i=1:size_zhong_qun
	fy（i:）=fitness_2（x0y（i:））;	
end

%最大速度值过大，容易导致越过最优区域，值过小，容易陷入局部最优
V_max=3;                                        %速度最大值为2
V_min=-3;                                       %速度最小值为-2

x_max=20;                                       %种群x中个体位置的最大值
x_min=-20;                                      %种群x中个体位置的最小值
y_max=20;                                       %种群y中个体位置的最大值
y_min=-20;                                      %种群y中个体位置的最小值

%群体中个体的最佳位置随着迭代的进行而不断发生变化
personal_best_local_x=x;                               %种群x中每个个体的最佳位置，随着迭代的进行，数值将会发生改变
personal_best_value_x=fx;                              %种群x中每个个体的最佳适应度值
[global_best_value_xi]=min（personal_best_value_x）;    %计算出当前的全局最佳值和最佳值所对应的位置
global_best_local_x=personal_best_local_x（i:）;        %当前的全局最佳位置

personal_best_local_y=y;                               %种群y中每个个体的最佳位置
personal_best_value_y=fy;                              %种群y中每个个体的最佳适应度值
[global_best_value_yi]=min（personal_best_value_y）;    %计算出当前的全局最佳值和最佳值所对应的位置
global_best_local_y=personal_best_local_y（i:）;        %当前的全局最佳位置
%%
%迭代寻优
for i = 1:max_die_dai
    %更新惯性权重因子
    w=w_max-（（w_max-w_min）/max_die_dai）*i;
      %粒子群x的更新
      for j = 1:size_zhong_qun
        % 速度更新
        vx（j:） = w*vx（j:） + c1*rand*（personal_best_local_x（j:）- x（j:）） + c2*rand*（global_best_local_x - x（j:））;
        if（vx（j:）> V_max）
            vx（j:） = V_max;
        else
            if（vx（j:）< V_min）
               vx（j:） = V_min;
            end
        end
        
        % 种群更新
        x（j:） = x（j:） + vx（j:）;
        if（ x（j:）> x_max）
            x（j:） = x_max;
        else
            if（x（j:）< x_min）
              x（j:）= x_min;
            end
        end        
        
        % 适应度值更新
        fx_update（