matlab 常用算法大全.doc

Matlab 高级算法程序代码汇总一、灰色预测模型matlab程序% renkou1=renkou(:,1);%年末常住人口数 % renkou2=renkou(:,2);%户籍人口% renkou3=renkou(:,3);%非户籍人口% shjian=1979:2010; %以上数据自己给x0=renkou2'n=length(x0);lamda=x0(1:n-1)./x0(2:n)range=minmax(lamda)x1=cumsum(x0)for i=2:nz(i)=0.5*(x1(i)+x1(i-1);endB=-z(2:n)',ones(n-1,1);Y=x0(2:n)'u=BYx=dsolve('Dx+a*x=b','x(0)=x0');x=subs(x,'a','b','x0',u(1),u(2),x1(1);yuce1=subs(x,'t',0:n-1);digits(6),y=vpa(x) %为提高预测精度，先计算预测值，再显示微分方程的解yuce=x0(1),diff(yuce1)epsilon=x0-yuce %计算残差delta=abs(epsilon./x0) %计算相对误差rho=1-(1-0.5*u(1)/(1+0.5*u(1)*lamda %计算级比偏差值%以深圳人口数据得到预测模型及预测误差相关数据 lamda = Columns 1 through 8 0.9741 0.9611 0.9419 0.8749 0.9311 0.9093 0.9302 0.9254 Columns 9 through 16 0.9245 0.9278 0.9442 0.9376 0.9127 0.9148 0.9332 0.9477 Columns 17 through 24 0.9592 0.9445 0.9551 0.9562 0.9594 0.9461 0.9469 0.9239 Columns 25 through 31 0.9140 0.9077 0.9243 0.9268 0.9312 0.9446 0.9618 range = 0.8749 0.9741 x1 = 1.0e+003 * Columns 1 through 8 0.0313 0.0634 0.0967 0.1322 0.1727 0.2162 0.2641 0.3155 Columns 9 through 16 0.3711 0.4313 0.4961 0.5647 0.6380 0.7182 0.8059 0.8999 Columns 17 through 24 0.9990 1.1024 1.2119 1.3265 1.4463 1.5712 1.7033 1.8427 Columns 25 through 32 1.9936 2.1588 2.3407 2.5375 2.7499 2.9780 3.2194 3.4705 u = -0.0665 31.3737 y = -472.117+503.377*exp(.664533e-1*t)yuce = Columns 1 through 8 31.2600 34.5876 36.9641 39.5040 42.2183 45.1192 48.2194 51.5326 Columns 9 through 16 55.0734 58.8576 62.9017 67.2238 71.8428 76.7792 82.0548 87.6928 Columns 17 through 24 93.7183 100.1578 107.0397 114.3945 122.2547 130.6550 139.6324 149.2267 Columns 25 through 32 159.4802 170.4382 182.1492 194.6649 208.0405 222.3352 237.6121 253.9386 epsilon = Columns 1 through 8 0 -2.4976 -3.5741 -4.0540 -1.6983 -1.5992 -0.3594 -0.0826 Columns 9 through 16 0.5266 1.2824 1.9183 1.4262 1.3772 3.4408 5.6352 6.2772 Columns 17 through 24 5.4417 3.2222 2.4203 0.2055 -2.4047 -5.7350 -7.5924 -9.7767 Columns 25 through 32 -8.5502 -5.3082 -0.2192 2.1651 4.3395 5.7348 3.8379 -2.9086 delta = Columns 1 through 8 0 0.0778 0.1070 0.1144 0.0419 0.0367 0.0075 0.0016 Columns 9 through 16 0.0095 0.0213 0.0296 0.0208 0.0188 0.0429 0.0643 0.0668 Columns 17 through 24 0.0549 0.0312 0.0221 0.0018 0.0201 0.0459 0.0575 0.0701 Columns 25 through 32 0.0567 0.0321 0.0012 0.0110 0.0204 0.0251 0.0159 0.0116 rho = Columns 1 through 8 -0.0411 -0.0271 -0.0066 0.0650 0.0049 0.0282 0.0058 0.0110 Columns 9 through 16 0.0119 0.0084 -0.0091 -0.0020 0.0245 0.0223 0.0027 -0.0128 Columns 17 through 24 -0.0251 -0.0094 -0.0208 -0.0219 -0.0254 -0.0111 -0.0119 0.0126 Columns 25 through 31 0.0232 0.0300 0.0122 0.0095 0.0048 -0.0095 -0.0280二、遗传算法程序代码% Optimizing a function  using Simple Genetic Algorithm with elitist preserved%Max f(x1,x2)=100*(x1*x1-x2).2+(1-x1).2; -2.0480<=x1,x2<=2.0480% Author: Wang Yonglin (wylin77)clc;clear all;format long;%设定数据显示格式%初始化参数T=100;%仿真代数N=80;% 群体规模pm=0.05;pc=0.8;%交叉变异概率umax=2.048;umin=-2.048;%参数取值范围L=10;%单个参数字串长度，总编码长度2Lbval=round(rand(N,2*L);%初始种群bestv=-inf;%最优适应度初值%迭代开始for ii=1:T%解码，计算适应度for i=1:Ny1=0;y2=0;for j=1:1:Ly1=y1+bval(i,L-j+1)*2(j-1);endx1=(umax-umin)*y1/(2L-1)+umin;for j=1:1:Ly2=y2+bval(i,2*L-j+1)*2(j-1);endx2=(umax-umin)*y2/(2L-1)+umin;obj(i)=100*(x1*x1-x2).2+(1-x1).2; %目标函数 xx(i,:)=x1,x2;endfunc=obj;%目标函数转换为适应度函数p=func./sum(func);q=cumsum(p);%累加fmax,indmax=max(func);%求当代最佳个体if fmax>=bestvbestv=fmax;%到目前为止最优适应度值bvalxx=bval(indmax,:);%到目前为止最佳位串optxx=xx(indmax,:);%到目前为止最优参数end   Bfit1(ii)=bestv; % 存储每代的最优适应度%遗传操作开始%轮盘赌选择for i=1:(N-1)r=rand;tmp=find(r<=q);newbval(i,:)=bval(tmp(1),:);end newbval(N,:)=bvalxx;%最优保留bval=newbval;%单点交叉for i=1:2:(N-1)cc=rand;if cc<pcpoint=ceil(rand*(2*L-1);%取得一个1到2L-1的整数ch=bval(i,:);bval(i,point+1:2*L)=bval(i+1,point+1:2*L);bval(i+1,point+1:2*L)=ch(1,point+1:2*L);endend   bval(N,:)=bvalxx;%最优保留%位点变异mm=rand(N,2*L)<pm;%N行mm(N,:)=zeros(1,2*L);%最后一行不变异，强制赋0bval(mm)=1-bval(mm); end%输出plot(Bfit1);% 绘制最优适应度进化曲线bestv   %输出最优适应度值optxx    %输出最优参数三、种子群算法程序代码%declare the parameters of the optimizationmax_iterations = 1000;no_of_particles = 50;dimensions = 1;delta_min = -0.003;delta_max = 0.003;c1 = 1.3;c2 = 1.3;%initialise the particles and teir velocity componentsfor count_x = 1:no_of_particlesfor count_y = 1:dimensionsparticle_position(count_x,count_y) = rand*10;particle_velocity(count_x,count_y) = rand;p_best(count_x,count_y) = particle_position(count_x,count_y);endend%initialize the p_best_fitness arrayfor count = 1:no_of_particlesp_best_fitness(count) = -1000;end%particle_position%particle_velocity%main particle swrm routinefor count = 1:max_iterations%find the fitness of each particle%change fitness function as per equation requiresd and dimensionsfor count_x = 1:no_of_particles%x = particle_position(count_x,1);%y = particle_position(count_x,2);%z = particle_position(count_x,3);%soln = x2 - 3*y*x + z;%x = particle_position(count_x);%soln = x2-2*x+1;x = particle_position(count_x);soln = x-7;if soln=0current_fitness(count_x) = 1/abs(soln);elsecurrent_fitness =1000;endend%decide on p_best etc for each particlefor count_x = 1:no_of_particlesif current_fitness(count_x) >p_best_fitness(count_x)p_best_fitness(count_x) = current_fitness(count_x);for count_y = 1:dimensionsp_best(count_x,count_y) = particle_position(count_x,count_y);endendend%decide on the global best among all the particlesg_best_val,g_best_index = max(current_fitness);%g_best contains the position of teh global bestfor count_y = 1:dimensionsg_best(count_y) = particle_position(g_best_index,count_y);end%update the position and velocity compponentsfor count_x = 1:no_of_particlesfor count_y = 1:dimensionsp_current(count_y) = particle_position(count_x,count_y);endfor count_y = 1:dimensionsparticle_velocity(count_y) = particle_velocity(count_y) + c1*rand*(p_best(count_y)-p_current(count_y) + c2*rand*(g_best(count_y)-p_current(count_y);particle_positon(count_x,count_y) = p_current(count_y) +particle_velocity(count_y);endendendg_bestcurrent_fitness(g_best_index)clear all, clc % pso exampleiter = 1000; % number of algorithm iterationsnp = 2; % number of model parametersns = 10; % number of sets of model parametersWmax = 0.9; % maximum inertial weightWmin = 0.4; % minimum inertial weightc1 = 2.0; % parameter in PSO methodologyc2 = 2.0; % parameter in PSO methodologyPmax = 10 10; % maximum model parameter valuePmin = -10 -10; % minimum model parameter valueVmax = 1 1; % maximum change in model parameterVmin = -1 -1; % minimum change in model parametermodelparameters(1:np,1:ns) = 0; % set all model parameter estimates for all model parameter sets to zeromodelparameterchanges(1:np,1:ns) = 0; % set all change in model parameter estimates for all model parameter sets to zerobestmodelparameters(1:np,1:ns) = 0; % set best model parameter estimates for all model parameter sets to zerosetbestcostfunction(1:ns) = 1e6; % set best cost function of each model parameter set to a large numberglobalbestparameters(1:np) = 0; % set best model parameter values for all model parameter sets to zerobestparameters = globalbestparameters' % best model parameter values for all model parameter sets (to plot)globalbestcostfunction = 1e6; % set best cost function for all model parameter sets to a large numberi = 0; % indicates ith algorithm iterationj = 0; % indicates jth set of model parametersk = 0; % indicates kth model parameterfor k = 1:np % initializationfor j = 1:nsmodelparameters(k,j) = (Pmax(k)-Pmin(k)*rand(1) + Pmin(k); % randomly distribute model parametersmodelparameterchanges(k,j) = (Vmax(k)-Vmin(k)*rand(1) + Vmin(k); % randomly distribute change in model parametersendendfor i = 2:iterfor j = 1:nsx = modelparameters(:,j);% calculate cost functioncostfunction = 105*(x(2)-x(1)2)2 + (1-x(1)2;if costfunction <setbestcostfunction(j) % best cost function for jth set of model parametersbestmodelparameters(:,j) = modelparameters(:,j);setbestcostfunction(j) = costfunction;end四、模拟退火算法% for d=1:50 %循环10次发现最小路径为4.115，循环50次有3次出现4.115 T_max=80; %input('please input the start temprature'); T_min=0.001; %input('please input the end temprature'); iter_max=100;%input('please input the most interp steps on the fit temp'); s_max=100; %input('please input the most steady steps ont the fit temp'); T=T_max; load .address.txt; order1=randperm(size(address,1)'%生成初始解。 figure(1); plot(address(order1,1),address(order1,2),'*r-') title('随机产生的路径'); totaldis1=distance(address,order1); for n=1:size(address,1) text(address(n,1)+0.01,address(n,2),num2str(n)%标号 end text(0.9,0.9,num2str(totaldis1) figure(2); while T>=T_min iter_num=1; s_num=1; plot(T,totaldis1,'r.') hold on while iter_num<iter_max&s_num<s_max; order2=exhgpath(order1); %随机交换两个城市位置 totaldis2=distance(address,order2); R=rand; DeltaDis=totaldis2-totaldis1; %新的距离-原来的距离 if DeltaDis<0; order1=order2; totaldis1=totaldis2; %新距离小，无条件接受 elseif (exp(totaldis1-totaldis2)/(T)>R)%本算法最核心的思想：以一定概率接受坏的结果，防止局部最优 order1=order2; totaldis1=totaldis2; else s_num=s_num+1; end iter_num=iter_num+1; end T=T*0.99; end set(gca,'xscale','log');%或者使用semilogx，有相同效果 xlabel('退火温度');ylabel('总距离'); order1 totaldis1 figure(3) plot(address(order1,1),address(order1,2),'*b-') title('最终路径'); for n=1:size(address,1) text(address(n,1)+0.01,address(n,2),num2str(n)%标号 end text(0.9,0.9,num2str(totaldis1) dstc(d)=totaldis1;% endfunction y=exhgpath(order) while 1 b=size(order,1); r=unidrnd(b,1,2); if r(1)-r(2)=0 break end end b=order(r(2); order(r(2)=order(r(1); order(r(1)=b; y=order;-function y=distance(address,order)nmb=size(address,1);y=0;for i=1:nmb-1 y=y+sqrt(address(order(i+1),1)-address(order(i),1)2+(address(order(i+1),2)-address(order(i),2)2);endy=y+sqrt(address(order(i+1),1)-address(order(1),1)2+(address(order(i+1),2)-address(order(1),2)2);人工神经网络程序代码%产生指定类别的样本点，并在图中绘出 X = 0 1; 0 1;    % 限制类中心的范围 clusters = 5;      % 指定类别数目 points = 10;        % 指定每一类的点的数目 std_dev = 0.05;    % 每一类的标准差 P = nngenc(X,clusters,points,std_dev); plot(P(1,:),P(2,:),'+r'); title('输入样本向量'); xlabel('p(1)'); ylabel('p(2)'); %建立网络 net=newc(0 1;0 1,5,0.1); %设置神经元数目为5 %得到网络权值，并在图上绘出 figure; plot(P(1,:),P(2,:),'+r'); w=net.iw1 hold on; plot(w(:,1),w(:,2),'ob'); hold off; title('输入样本向量及初始权值'); xlabel('p(1)'); ylabel('p(2)'); figure; plot(P(1,:),P(2,:),'+r'); hold on; %训练网络 net.trainParam.epochs=7; net=init(net); net=train(net,P); %得到训练后的网络权值，并在图上绘出 w=net.iw1 plot(w(:,1),w(:,2),'ob'); hold off; title('输入样本向量及更新后的权值'); xlabel('p(1)'); ylabel('p(2)'); a=0; p = 0.6 ;0.8; a=sim(net,p) example8_2 %随机生成1000个二维向量，作为样本，并绘出其分布 P = rands(2,1000); plot(P(1,:),P(2,:),'+r') title('初始随机样本点分布'); xlabel('P(1)'); ylabel('P(2)'); %建立网络，得到初始权值 net=newsom(0 1; 0 1,5 6); w1_init=net.iw1,1 %绘出初始权值分布图 figure; plotsom(w1_init,net.layers1.distances) %分别对不同的步长，训练网络，绘出相应的权值分布图 for i=10:30:100     net.trainParam.epochs=i;     net=train(net,P);     figure;     plotsom(net.iw1,1,net.layers1.distances) end %对于训练好的网络，选择特定的输入向量，得到网络的输出结果 p=0.5;0.3; a=0; a = sim(net,p) example8_3 %指定输入二维向量及其类别 P = -3 -2 -2   0   0   0   0 +2 +2 +3;         0 +1 -1 +2 +1 -1 -2 +1 -1   0; C = 1 1 1 2 2 2 2 1 1 1; %将这些类别转换成学习向量量化网络使用的目标向量 T = ind2vec(C) %用不同的颜色，绘出这些输入向量 plotvec(P,C), title('输入二维向量'); xlabel('P(1)'); ylabel('P(2)'); %建立网络 net = newlvq(minmax(P),4,.6 .4,0.1); %在同一幅图上绘出输入向量及初始权重向量 figure; plotvec(P,C) hold on W1=net.iw1; plot(W1(1,1),W1(1,2),'ow') title('输入以及权重向量'); xlabel('P(1), W(1)'); ylabel('P(2), W(2)'); hold off; %训练网络，并再次绘出权重向量 figure; plotvec(P,C); hold on; net.trainParam.epochs=150; net.trainParam.show=Inf; net=train(net,P,T); plotvec(net.iw1',vec2ind(net.lw2),'o'); %对于一个特定的点，得到网络的输出 p = 0.8; 0.3; a = vec2ind(sim(net,p)五、蚁群算法程序代码% the procedure of ant colony algorithm for VRP% % % % % % % % % % %initialize the parameters of ant colony algorithmsload data.txt;d=data(:,2:3);g=data(:,4);m=31; % 蚂蚁数alpha=1;belta=4;% 决定tao和miu重要性的参数lmda=0;rou=0.9;%衰减系数q0=0.95;% 概率tao0=1/(31*841.04);%初始信息素Q=1;%蚂蚁循环一周所释放的信息素defined_phrm=15.0; % initial pheromone level value QV=100; % 车辆容量vehicle_best=round(sum(g)/QV)+1;%所完成任务所需的最少车数V=40;% 计算两点的距离 for i=1:32; for j=1:32; dist(i,j)=sqrt(d(i,1)-d(j,1)2+(d(i,2)-d(j,2)2); end;end;%给tao miu赋初值for i=1:32; for j=1:32; if i=j; %s(i,j)=dist(i,1)+dist(1,j)-dist(i,j); tao(i,j)=defined_phrm; miu(i,j)=1/dist(i,j); end; end;end;for k=1:32; for k=1:32; deltao(i,j)=0; end;end; best_cost=10000; for n_gen=1:50; print_head(n_gen); for i=1:m; %best_solution=; print_head2(i); sumload=0; cur_pos(i)=1; rn=randperm(32); n=1; nn=1; part_sol(nn)=1; %cost(n_gen,i)=0.0; n_sol=0; % 由