+ LFPSO files

Matlab_LFPSO/Calc_refl.m

@@ -0,0 +1,14 @@
+function Chi_Square=Calc_refl(T,Size)
+
+ for i=1:Size
+    layers(:,:)=T(:,:,i);%A three-dimensional array turns into a two-dimensional array
+    Layers=array2table(layers);
+    filename= ['.\Input\input_structure_' num2str(i) '.txt'];
+    writetable(Layers,  filename, 'Delimiter',';','WriteVariableNames',false);
+ end
+system("xrccmd.exe  -f input_structure_#.txt -i exp_data.dat -n "+num2str(Size));
+fileID = fopen("ChiSquare.dat","r");
+formatSpec = '%f';
+A = fscanf(fileID,formatSpec);  
+fclose(fileID);
+Chi_Square = A;

Matlab_LFPSO/Tent.m

@@ -0,0 +1,14 @@
+ %TentÓ³Éä
+ function [x] = Tent(Max_iter,a)
+ x(1)=rand; %³õʼµã
+%  a = 0.5;%²ÎÊýaµÄÖµ
+ for i=1:Max_iter-1
+     if x(i)<a
+         x(i+1)=2*x(i);
+     end
+     if x(i)>=a
+         x(i+1)=2*(1-x(i)) ;
+     end
+ end
+ end
+

Matlab_LFPSO/extract_H.m

@@ -0,0 +1,26 @@
+function [x,y1,y2,y3]=extract_H(~)
+filename = 'current_best_structure.txt';
+fileID = fopen(filename);
+A = textscan(fileID,'%*q %f %*f %*f','Delimiter',';');
+B=cell2mat(A);% whole thickness
+
+%% extract Boron thickness
+B1=B(3:end-2);
+n = length(B1); 
+y1 = B1(1:4:n); 
+x=linspace(1,length(y1),length(y1)); % number of layer
+
+%% extract MoB2 thickness
+B2=B(4:end-2);
+y2= B2(1:4:n); 
+
+%% extract Mo thickness
+B3=B(5:end-2);
+y3= B3(1:4:n); 
+
+%% plot
+% plot(x,y1);hold on
+% plot(x,y2);hold on
+% plot(x,y3);
+% ylim([0 B1(1)+10])
+% xlabel('Layer number'); ylabel('Thickness(A)');

Matlab_LFPSO/initializationNew.m

@@ -0,0 +1,37 @@
+%% »ùÓÚTentÓ³ÉäµÄÖÖȺ³õʼ»¯
+function Positions=initializationNew(SearchAgents_no,dim,ub,lb,a)
+Boundary_no= size(ub,2); % numnber of boundaries
+Positions=zeros(SearchAgents_no,Boundary_no);
+% If the boundaries of all variables are equal and user enter a signle
+% number for both ub and lb
+if Boundary_no==1
+    for i = 1:SearchAgents_no
+        Temp = Tent(dim,a).*(ub-lb)+lb;
+        Temp(Temp>ub) = ub;
+        Temp(Temp<lb) = lb;
+        Positions(i,:)=Temp;
+    end
+end
+
+% If each variable has a different lb and ub
+if Boundary_no>1
+   for j = 1:SearchAgents_no
+       TentValue =  Tent(dim,a);
+        for i=1:dim
+            ub_i=ub(i);
+            lb_i=lb(i);
+            Positions(j,i)=TentValue(i).*(ub_i-lb_i)+lb_i;
+            if(Positions(j,i)>ub_i)
+                Positions(j,i) = ub_i;
+            end
+            if(Positions(j,i)<lb_i)
+                Positions(j,i) = lb_i;
+            end
+        end
+   end
+end
+end
+
+
+
+

Matlab_LFPSO/input_stucture.m

@@ -0,0 +1,27 @@
+function  T=input_stucture(xi)
+%Step 1: exp. data
+ExpFileName = "exp_data.dat";
+opts = detectImportOptions(ExpFileName);
+opts.PreserveVariableNames = true;
+% ExpData = readtable(ExpFileName,opts);
+% Step 2: structure
+sz = [1204 4];
+varTypes = ["string","single","single","single"];
+varNames = ["Name","H","Sigma","rho"];
+Layers = table('Size',sz,'VariableTypes',varTypes,'VariableNames',varNames);
+Layers(1,:) = {"C", xi(1), xi(2), xi(3)};
+Layers(2,:) = {"B", xi(4), xi(5), xi(6)};
+for i=3:4:1199
+Layers(i,:) = {"B",  xi(7+3*fix(i/4)), xi(7+3*fix(i/4)+1), xi(7+3*fix(i/4)+2)};
+Layers(i+1,:) = {"MoB2", xi(3*300+7+3*fix(i/4)), xi(3*300+7+3*fix(i/4)+1), xi(3*300+7+3*fix(i/4)+2)};
+% Layers(i+1+1,:) = {"Mo", (fix(i/4)*0.37/300+3.7), xi(3*300+3*300+7+3*fix(i/4)+1), xi(3*300+3*300+7+3*fix(i/4)+2)};
+Layers(i+1+1,:) = {"Mo", xi(3*300+3*300+7+3*fix(i/4)), xi(3*300+3*300+7+3*fix(i/4)+1), xi(3*300+3*300+7+3*fix(i/4)+2)};
+Layers(i+1+1+1,:) = {"MoB2", xi(3*300+3*300+3*300+7+3*fix(i/4)), xi(3*300+3*300+3*300+7+3*fix(i/4)+1), xi(3*300+3*300+3*300+7+3*fix(i/4)+2)};
+end
+
+Layers(end-1,:) = {"Mo", xi(end-2), xi(end-1), xi(end)};
+Layers(end,:) = {"Si", 100000, 5, 2.33};
+
+%save Size=100 stuctures.txt
+T=table2array(Layers);
+

Matlab_LFPSO/levy.m

@@ -0,0 +1,35 @@
+function [z] = levy(n,m,beta)
+% This function implements Levy's flight. 
+
+% For more information see 
+%'Multiobjective cuckoo search for design optimization Xin-She Yang, Suash Deb'. 
+
+
+% Input parameters
+% n     -> Number of steps 
+% m     -> Number of Dimensions 
+% beta  -> Power law index  % Note: 1 < beta < 2
+
+% Output 
+% z     -> 'n' levy steps in 'm' dimension
+
+    num = gamma(1+beta)*sin(pi*beta/2); % used for Numerator 
+
+    den = gamma((1+beta)/2)*beta*2^((beta-1)/2); % used for Denominator
+
+    sigma_u = (num/den)^(1/beta);% Standard deviation
+% 
+%     u = random('Normal',0,sigma_u^2,n,m); 
+%     
+%     v = random('Normal',0,1,n,m);
+    u = normrnd(0,sigma_u,n,m);
+
+    v = normrnd(0,1,n,m);
+
+    z = u./(abs(v).^(1/beta));
+
+
+end
+
+
+

Matlab_LFPSO/main_CPSO_v2_vp1.m

@@ -0,0 +1,157 @@
+function [best_x , f_gg]=main_CPSO_v2_vp1(Stucture,Xmax,Xmin,Size,Tmax)
+%% PSO parameters
+dimension=length(Xmax);%dimension
+% Size=4;% Size
+% Tmax=3;%%iteration  Tmax
+c1=1.49445;
+c2=1.49445;%study factor
+ w_max = 0.8;                        
+ w_min = 0.1;
+k=0.2;%Proportional relationship between position and velocity
+Vmax= Xmax*k;
+Vmin=-Vmax;   
+u=3.999;
+MaxC=10;
+a=0.5;
+%  X=zeros(Size,dimension);
+%  V=zeros(Size,dimension);
+% f_x=zeros(1,Size);
+% f_xbest=zeros(1,Size);
+%% initialization population
+X = rand(Size,dimension).*(Xmax-Xmin)+Xmin;
+V = rand(Size,dimension).*(Vmax-Vmin)+Vmin;
+Pi=X;  
+for j=1:Size
+   structure_3D(:,:,j)=Stucture(X(j,:));% used for generating Size structures      
+end
+f_x=Calc_refl(structure_3D,Size);% used for generating Chi-square  
+f_xbest=f_x;  
+[f_g,I]=min(f_xbest);     
+Pg=X(I,:); 
+time=zeros(1,Tmax);
+f_gg=zeros(1,Tmax);
+%% iteration
+for t=1:Tmax
+    time(t)=t;
+    fprintf('Current evolution number：%d\n',t)
+     w = w_max-(w_max-w_min)*t/Tmax; %Adaptive inertia weights (linearly decreasing)
+    r1=rand(1);r2=rand(1);
+    for j=1:Size
+        V(j,:) = w*V(j,:) + c1*r1*(Pi(j,:)-X(j,:)) + c2*r2*(Pg-X(j,:));  %Speed Update
+        index1 = find(V(j,:)>Vmax);
+        V(j,index1) =  rand*(Vmax(index1)-Vmin(index1))+Vmin(index1);
+          index1=  find(V(j,:)<Vmin);
+          V(j,index1) = rand*(Vmax(index1)-Vmin(index1))+Vmin(index1);
+
+        X(j,:)= X(j,:)+V(j,:);  %position update
+
+      index1 = find(X(j,:)>Xmax);
+      X(j,index1) = rand*(Xmax(index1)-Xmin(index1))+Xmin(index1);
+      index1=  find(X(j,:)<Xmin);
+      X(j,index1) = rand*(Xmax(index1)-Xmin(index1))+Xmin(index1);
+    end
+
+    for j=1:Size
+    structure_3D(:,:,j)=Stucture(X(j,:));%  used for generating  updated Size 3 dimension structures   structure_3D 
+    end
+    f_x=Calc_refl(structure_3D,Size); % used for generating  updated Chi-square f_x
+    for j=1:Size
+        if f_x(j)<f_xbest(j)
+            f_xbest(j)=f_x(j);
+            Pi(j,:)=X(j,:);
+        end
+        if f_xbest(j)<f_g
+            Pg= Pi(j,:);
+            f_g=f_xbest(j);
+        end
+    end
+
+       %chaotic mapping
+%         y=zeros(Size,dimension);
+%         y(1,:)=(Pg-Xmin)/(Xmax-Xmin);
+%         structure_3D(:,:,1)=Stucture(y(1,:));
+%         f_x(1)=Calc_refl(structure_3D(:,:,1),1);
+%         for tt=1:Size-1 
+%             for e=1:dimension
+%                 
+%                  y(tt+1,e)=u*y(tt,e).*(1-y(tt,e)); 
+%                  
+%             end
+%             y(tt+1,:)=Xmin+(Xmax-Xmin).*y(tt+1,:);   
+%         structure_3D(:,:,tt+1)=Stucture(y(tt+1,:));
+%        f_x(tt+1)=Calc_refl(structure_3D(:,:,tt+1),1); 
+%         end
+%          
+%            
+%         [~,ybestindex]=min(f_x);%寻找最优混沌可行解矢量
+%       ybest=y(ybestindex,:);
+%      ran=1+fix(rand()*Size);%产生一随机数1~sizepop之间
+%      X(ran,:)=ybest;
+
+%     f_gg(t)=f_g;
+% %     PGG{t}=Pg;
+%       [sort_f_x,index]=sort(f_xbest);
+%       Nbest=round(Size*0.2);
+% %       if f_g<sort_f_x(1)
+% %           X(index(1),:)=Pg;
+% %           sort_f_x(1) =f_g;
+% %       end
+%       cx=zeros(1,dimension);
+%      for n=1: Nbest
+%          tmpx=X(index(n),:);
+%          for k=1:MaxC 
+%                   
+%              for dim=1:dimension
+%                    cx(dim)=(tmpx(1,dim)-Xmin(dim))/(Xmax(dim)-Xmin(dim));
+%                    cx(dim)=4*cx(dim)*(1-cx(dim));
+%              
+%                 tmpx(1,dim)=Xmin(dim)+cx(dim)*(Xmax(dim)-Xmin(dim));
+%              end
+%               structure_3D_1=Stucture(tmpx);
+%                fcs=Calc_refl(structure_3D_1,1);
+% %               fcs=Fitness_Function(tmpx,deg00,rdata,refl_inv);
+%               if fcs<sort_f_x(n)
+%                   X(index(n),:)=tmpx;
+%                   break;
+%               end
+%          end
+% %          X(index(n),:)=tmpx;
+%      end
+%       X(1:Nbest,:)=X(index(1:Nbest),:);
+%       Pg=X(index(1),:);
+%         for s=1:dimension
+%         XXmin(s)=max(Xmin(s),Pg(s)-rand*(Xmax(s)-Xmin(s)));
+%         XXmax(s)=min(Xmax(s),Pg(s)+rand*(Xmax(s)-Xmin(s)));
+%         end
+%           VVVmin=k* XXmin;
+%           VVVmax=k* XXmax;
+%          for i=Nbest+1:Size
+%             V(i,:) = VVVmin+(VVVmax-VVVmin).*rand(1,dimension);   
+%             X(i,:)=XXmin+(XXmax-XXmin).*rand(1,dimension); 
+%          end
+%          for j=1:Size
+%     structure_3D(:,:,j)=Stucture(X(j,:)); 
+%          end
+%     f_x=Calc_refl(structure_3D,Size);
+%       for j=1:Size
+%         if f_x(j)<f_xbest(j)
+%             f_xbest(j)=f_x(j);
+%             Pi(j,:)=X(j,:);
+%         end
+%         if f_xbest(j)<f_g
+%             Pg= Pi(j,:);
+%             f_g=f_xbest(j);
+%         end
+%      end
+
+    f_gg(t)=f_g;
+    PGG{t}=Pg;
+end
+%% Obtain the optimal individual
+[~,index]=min(f_gg);
+best_x=PGG{index(1)};
+%% Plot adaptation curves
+figure
+plot(time,f_gg);
+xlabel('Number of iterations');ylabel('Adaption value f_g'); 
+end