博科园AI人工智能助手
图灵
闪电是因为天空中的大量正负电荷聚集在一起引发强烈放电才形成的。肉当雷雨云移到某处时，云的中下部是强大负电荷中心，云底相对的下垫面变成正电荷中心，在云底与地面间形成强大电场。在电荷越积越多，电场越来越强的情况下，云底首先出现大气被强烈电离的一段气柱，称梯级先导。这种电离气柱逐级向地面延伸，每级梯级先导是直径约5米、长50米、电流约100安培的暗淡光柱，它以平均约150000米/秒的高速度一级一级地伸向地面，在离地面5─50米左右时，地面便突然向上回击，回击的通道是从地面到云底，沿着上述梯级先导开辟出的电离通道回击以5万公里/秒的更高速度从地面驰向云底，发出光亮无比的光柱，历时40微秒，通过电流超过1万安培，这即第一次闪击。
相隔几秒之后，从云中一根暗淡光柱，携带巨大电流，沿第一次闪击的路径飞驰向地面，称直窜先导，当它离地面5─50米左右时，地面再向上回击，再形成光亮无比光柱，这即第二次闪击。接着又类似第二次那样产生第三、四次闪击。通常由3─4次闪击构成一次闪电过程。一次闪电过程历时约0.25秒，在此短时间内，窄狭的闪电通道上要释放巨大的电能，因而形成强烈的爆炸，产生冲击波，然后形成声波向四周传开，这就是雷声或说“打雷”。
![[s-70]](https://www.bokeyuan.net/pic/image/emoji/cas/70.png)
![请详细介绍这个算法：
clear
clc
close all
clear global
global para_D index_0 para_size
global para_D_max para_size
%% 数据读取
[temp_a,temp_b,temp_c]=xlsread(‘案例5-附件1-华东区域供应商分布及发货情况.xlsx’);
[temp_a_1,temp_b_1,temp_c_1]=xlsread(‘location（附件一）.xlsx’);
temp=[];
temp_name=[];
%1-序号 2-供应商序号 3-货量 4-经度 5-纬度
for i=2:size(temp_c,1)
temp(i-1,:)=[i-1 temp_c{i,1} temp_c{i,6} temp_c_1{i,2} temp_c_1{i,3}];
temp_name{i-1,1}=temp_c{i,2};
end
index_0=[];
place_name=[];
temp_x=unique(temp(:,2));
for i=1:length(temp_x)
temp_index=find(temp(:,2)==temp_x(i));
index_0(end+1,:)=temp(temp_index(1),:);
index_0(end,3)=sum(temp(temp_index,3));
place_name{end+1,1}=temp_name{temp_index(1)};
end
%% 数据输入
para_D_max=300;%辐射半径
para_size=[min(index_0(:,4)) max(index_0(:,4))
min(index_0(:,5)) max(index_0(:,5))];
index_0(:,3)=index_0(:,3)/sum(index_0(:,3));
%% 数据剔除
center_ori=[121.19329 31.31426];
temp_D=[];
for i=1:size(index_0,1)
temp_D(i,1)=Spherediff(index_0(i,4:5),center_ori)/1000;
end
error_index=find(temp_D>=500);
index_0(error_index,:)=[];
place_name(error_index)=[];
%% 数据整理
%计算每个点之间的距离
para_n_dot=size(index_0,1);
para_D=zeros(para_n_dot);
for i=1:para_n_dot
for j=1+i:para_n_dot
para_D(i,j)=Spherediff(index_0(i,4:5),index_0(j,4:5))/1000;
para_D(j,i)=para_D(i,j);
end
end
%% 开始规划
sizepop=100;
%待求解参数的维度
dim=21;
%最大迭代次数
ger=100;
xlimit_max=ones(1,dim)-0.0001;
xlimit_min=zeros(1,dim)+0.0001;
%粒子群点的每次迭代的最大与最小变化速率
vlimit_max=0.2*(xlimit_max-xlimit_min);
vlimit_min =-1*vlimit_max;
%惯性权重，个体与群体权重，属于优化的参数
c_1=0.8;
c_2=0.5;
c_3=0.5;
%初始化种群，得到第一批解
pop_x=zeros(dim,sizepop);
pop_v=zeros(dim,sizepop);
for i=1:dim
for j=1:sizepop
pop_x(i,j)=(xlimit_min(i)+(xlimit_max(i)-xlimit_min(i))*rand);
pop_v(i,j)=(vlimit_min(i)+(vlimit_max(i)-vlimit_min(i))*rand);
end
end
gbest=pop_x;
%对每组解带入目标函数求解不通过限制函数的 设为无穷大
for j=1:sizepop
fitness_gbest(j)=fun_1(pop_x(:,j));
end
%找到最优解
zbest=pop_x(:,1);
fitness_zbest=fitness_gbest(1);
for j=1:sizepop
if fitness_gbest(j)<fitness_zbest
zbest=pop_x(:,j);
fitness_zbest=fitness_gbest(j);
end
end
%迭代次数的计数器
iter=1;
%初始化最优解的储存矩阵
record_2=zeros(ger,1);
%利用初始结果进行迭代
while iter<=ger
%根据最优个体的位置，更新每个个体，逐渐靠拢
for j=1:sizepop
pop_v(:,j)=(c_1*pop_v(:,j)+c_2*rand*(gbest(:,j)-pop_x(:,j))+c_3*rand*(zbest-pop_x(:,j)));
%边界更新，避免变化速度大于限制
for i=1:dim
if pop_v(i,j)>vlimit_max(i)
pop_v(i,j)=vlimit_max(i);
end
if pop_v(i,j)<vlimit_min(i)
pop_v(i,j)=vlimit_min(i);
end
end
pop_x(:,j)=pop_x(:,j)+pop_v(:,j);
%每个个体边界更新，避免个体超出限制
for i=1:dim
if pop_x(i,j)>xlimit_max(i)
pop_x(i,j)=xlimit_max(i);
end
if pop_x(i,j)<xlimit_min(i)
pop_x(i,j)=xlimit_min(i);
end
end
%随机变异，百分40的概率每个点重新生成，避免过早收敛
if rand>0.75
i=ceil(dim*rand);
pop_x(i,j)=(xlimit_min(i)+(xlimit_max(i)-xlimit_min(i))*rand);
end
%对种群进行目标值计算
fitness_pop(j)=fun_1(pop_x(:,j));
%更新最优解
if fitness_pop(j)<fitness_gbest(j)
gbest(:,j)=pop_x(:,j);
fitness_gbest(j)=fitness_pop(j);
end
if fitness_gbest(j)<fitness_zbest
zbest=gbest(:,j);
fitness_zbest=fitness_gbest(j);
end
end
%更新此代的最优值
record_2(iter)=fitness_zbest;
iter=iter+1;
disp(iter)
disp(fitness_zbest)
end
%% 回代
x=zbest;
%聚类数量
X_1=ceil(10*x(1));
%初始质心坐标
X_2=reshape(x(2:1+2*10),[],2);
X_2=X_2(1:X_1,:);
for i=1:2
X_2(:,i)=para_size(i,1)+(para_size(i,2)-para_size(i,1)).*X_2(:,i);
end
%初始聚类
para_n_r=size(index_0,1);
real_select=zeros(para_n_r,1);
for i=1:para_n_r
temp_D=[];
for j=1:X_1
temp_D(j,1)=Spherediff(index_0(i,4:5),X_2(j,:))/1000;
end
[min_D,min_index]=min(temp_D);
if min_D>para_D_max
real_select(i)=0;
else
real_select(i)=min_index;
end
end
%更新质心
for i=1:X_1
temp_index=find(real_select==i);
if length(temp_index)>0
X_2(i,:)=[sum(index_0(temp_index,4).*index_0(temp_index,3)/sum(index_0(temp_index,3))) sum(index_0(temp_index,5).*index_0(temp_index,3)/sum(index_0(temp_index,3)))];
end
end
old_type=real_select;
%重分配
for i=1:para_n_r
temp_D=[];
for j=1:X_1
temp_D(j,1)=Spherediff(index_0(i,4:5),X_2(j,:))/1000;
end
[min_D,min_index]=min(temp_D);
if min_D>para_D_max
real_select(i)=0;
else
real_select(i)=min_index;
end
end
while sum(abs(real_select-old_type))>0
%更新质心
for i=1:X_1
temp_index=find(real_select==i);
if length(temp_index)>0
X_2(i,:)=[sum(index_0(temp_index,4).*index_0(temp_index,3)/sum(index_0(temp_index,3))) sum(index_0(temp_index,5).*index_0(temp_index,3)/sum(index_0(temp_index,3)))];
end
end
old_type=real_select;
%重分配
for i=1:para_n_r
temp_D=[];
for j=1:X_1
temp_D(j,1)=Spherediff(index_0(i,4:5),X_2(j,:))/1000;
end
[min_D,min_index]=min(temp_D);
if min_D>para_D_max
real_select(i)=0;
else
real_select(i)=min_index;
end
end
end
%确定每个点的服务量
y_5=zeros(X_1,1);
y_6=zeros(X_1,1);
for i=1:X_1
temp_index=find(real_select==i);
y_5(i)=sum(index_0(temp_index,3));
y_6(i)=length(temp_index);
end
%% 给予修正
error_index=find(y_6<=10);
X_1=X_1-length(error_index);
X_2(error_index,:)=[];
%重分配
for i=1:para_n_r
temp_D=[];
for j=1:X_1
temp_D(j,1)=Spherediff(index_0(i,4:5),X_2(j,:))/1000;
end
[min_D,min_index]=min(temp_D);
real_select(i)=min_index;
end
for i=1:X_1
temp_index=find(real_select==i);
if length(temp_index)>0
X_2(i,:)=[sum(index_0(temp_index,4).*index_0(temp_index,3)/sum(index_0(temp_index,3))) sum(index_0(temp_index,5).*index_0(temp_index,3)/sum(index_0(temp_index,3)))];
end
end
%确定每个点到目标点的距离
y_1=0;
y_2=0;
y_4=zeros(para_n_r,1);
for i=1:para_n_r
temp_D=[];
for j=1:X_1
temp_D(j,1)=Spherediff(X_2(j,:),index_0(i,4:5))/1000;
end
y_4(i,1)=sum(temp_D<para_D_max);
end
for i=1:X_1
temp_index=find(real_select==i);
temp_D=[];
for j=1:length(temp_index)
temp_D(j,1)=Spherediff(index_0(temp_index(j),4:5),X_2(i,:))/1000;
end
y_1=y_1+sum(temp_D.*index_0(temp_index,3));
if max(temp_D)>y_2
y_2=max(temp_D);
end
end
y=y_1*mean(y_4)*(1+max(y_5)-min(y_5))+y_2-para_D_max+10^9*sum(real_select==0);
figure()
scatter(index_0(:,4),index_0(:,5),50,’blue’,’filled’)
hold on
for i=1:length(place_name)
text(index_0(i,4),index_0(i,5),place_name{i})
end
scatter(X_2(:,1),X_2(:,2),’r’,’filled’)
for i=1:para_n_r
plot([index_0(i,4) X_2(real_select(i),1)],[index_0(i,5) X_2(real_select(i),2)],’LineWidth’,1.5,’Color’,’r’)
end
save data_1 X_2 real_select](https://pic.bokeyuan.net/image/liwu/%E6%A3%92%E6%A3%92%E7%B3%96.jpg)
![[s-8]](https://www.bokeyuan.net/pic/image/emoji/cas/8.png)
