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=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 ccpcpoint=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)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 =T_min iter_num=1; s_num=1; plot(T,totaldis1,r.) hold on while iter_numiter_max&s_nums_max; order2=exhgpath(order1); %随机交换两个城市位置 totaldis2=distance(address,order2); R=rand; DeltaDis=totaldis2-totaldis1; %新的距离-原来的距离 if DeltaDisR)%本算法最核心的思想：以一定概率接受坏的结果，防止局部最优 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; % 由蚂蚁产生的路径数量 M_vehicle=500; t=0; %最佳路径数组的元素数为0 while sumload=QV; for k=1:length(rn); if sumload+g(rn(k)maxp maxp=p(j); index_max=j; end; end; old_pos=cur_pos(i); if rand(1)q0 cur_pos(i)=A(index_max); else krnd=randperm(na); cur_pos(i)=A(krnd(1); bbb=old_pos cur_pos(i); ccc=1 1; if bbb=ccc; cur_pos(i)=A(krnd(2); end; end; tao(old_pos,cur_pos(i)=taolocalupdate(tao(old_pos,cur_pos(i),rou,tao0);%对所经弧进行局部更新 sumload=sumload+g(cur_pos(i); nn=nn+1; part_sol(nn)=cur_pos(i); temp_load=sumload; if cur_pos(i)=1; rn=setdiff(rn,cur_pos(i); n=1; A=; end; if cur_pos(i)=1; % 如果当前点为车场,将当前路径中的已访问用户去掉后,开始产生新路径 if setdiff(part_sol,1)=; n_sol=n_sol+1; % 表示产生的路径数,n_sol=1,2,3,.5,6.,超过5条对其费用加上车辆的派遣费用 fid=fopen(out_solution.txt,a+); fprintf(fid,%s%i%s,NO.,n_sol,条路径是:); fprintf(fid,%i ,part_sol); fprintf(fid,n); fprintf(fid,%s,当前的用户需求量是:); fprintf(fid,%in,temp_load); fprintf(fid,-n); fclose(fid); % 对所得路径进行路径内3-opt优化 final_sol=exchange(part_sol); for nt=1:length(final_sol); % 将所有产生的路径传给一个数组 temp(t+nt)=final_sol(nt); end; t=t+length(final_sol)-1; sumload=0; final