1.整数规划的蒙特卡洛解法2015-06-10 已知非线性整数规划为：
2222maxzx12x23x34x42x58x12x23x3x42x50xi99(i1,....,5)x1x2x3x4x5400x12x22x3x46x68002xx6x2002312x3x45x5200
如果用显枚举试探，共需要计算100^5=10^10个点，其计算量非常大。然而应用蒙特卡洛
去随机模拟计算10^6个点，便可以找到满意解，那么这种方法的可信度究竟怎么样呢？ 下面就分析随机采样10^6个点计算时，应用概率理论估计下可信度。 不是一般性，假设一个整数规划的最优点不是孤立的奇点。
假设目标函数落在高值区的概率分别为0.01，0.00001，则当计算10^6个点后，有任一个点落在高值区的概率分别为：
1-0.99^1000000=0.99...99(100多位) 1-0.99999^1000000=0.999954602
解 （1）首先编写M文件 mengte.m定义目标函数f和约束向量g,程序如下：
function [f,g]=mengte(x);
f=x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)-8*x(1)-2*x(2)-3*x(3)-... x(4)-2*x(5); g=[sum(x)-400
x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800 2*x(1)+x(2)+6*x(3)-200 x(3)*x(3)+x(4)+5*x(5)-200];
（2）编写M文件mainint.m如下求问题的解： rand('state',sum(clock)); p0=0; tic
for i=1:10^5
x=99*rand(5,1);
x1=floor(x);%向下取整 x2=ceil(x);%向上取整 [f,g]=mengte(x1); if sum(g<=0)==4 if p0<=f x0=x1; p0=f; end end
[f,g]=mengte(x2); if sum(g<=0)==4 if p0<=f
x0=x2; p0=f; end end end x0,p0
Matlab求解整数规划祥见第二章（优秀教材）
2.罚函数法 2015-06-11
利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题，因而也称这种方法为系列无约束最小化技术，简记为SUMT。
罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外函数法，另一种叫内函数法，下面介绍外函数法。 考虑问题：
minf(x)gi(x)0,i1......r s.t.hj(x)0,j1,....stkm(x)0,m1.....取一个充分大的数M>0,构造函数： P(x,M)f(x)Mrstmax(g(x),0)Mmin(h(x),0)M|k(x)|
iiii1i1i1增广目标函数P(x,M)为目标函数的无约束极值问题
minP(x,M)
的最优解x也是原问题的最优解。 Eg:求解下列非线性规划问题
2minf(x)x1x282x12x202s.t.x1x220x,x012解 编写M文件test1.m function g=test1(x); M=50000;
f=x(1)^2+x(2)^2+8;
g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)+M*abs(-x(1)-x(2)^2+2);
然后命令窗口输入
[x,y]=fminunc('test1',rand(2,1))
扩展： 将有约束的极值问题转化为无约束的极值问题，这样便可以很好的扩展matlab神经网络43个案例分析里面求无约束的函数极值问题。
Matlab求解非线性规划祥见第三章（优秀教材）
3.层次分析 2015-06-12 层次分析祥见第八章
准则层b=[1 1 1 4 1 1/2 ];
方案层cb1=[1 1/4 1/2]; cb2=[1 1/4 1/5];cb3=[1 3 1/3]; cb4=[1 1/3 5];cb5=[1 1 7]; cb6=[1 7 9];
由总权值可以知道工作1最合适。
准则 准则层权值 工作1 工作2 工作3 权重的计算
weight1.m 传入各因素重要性 例如b=[1 3 9]就可以计算相应的权重了weigth1(b)。 Weight.m是近似算法和weight1.m 用法相同。 主要还是用 weight1.m weight1.m
function [w,cr]=weight1(b) %b=[1 1/3 1/9];%模糊比较 n1=length(b);
ri=[0 0 0.58 0.90 1.12 1.24 1.32 1.41 .145];%一致性指标 a=zeros(n1,n1);
B1 0.1600 0.1429 0.5714 0.2857 B2 0.1600 0.1000 0.4000 0.5000 B3 0.1600 0.2308 0.0769 0.6923 B4 0.0400 0.2381 0.7143 0.0476 B5 0.1600 0.4667 0.4667 0.0667 B6 0.3200 0.7975 0.1139 0.0886 总权值 0.4152 0.3074 0.2774 for i=1:n1 for j=1:n1 a(i,j)=b(j)/b(i); end end
[x,y]=eig(a); lamda=max(diag(y)); num=find(diag(y)==lamda); w=x(:,num)/sum(x(:,num));
cr=(lamda-n1)/(n1-1)/ri(n1);%一致性检验 若cr<0.1,认为矩阵的一致性可以接受
主程序 score3.m
b=[1 1 1 4 1 1/2 ]; cb1=[1 1/4 1/2]; cb2=[1 1/4 1/5]; cb3=[1 3 1/3]; cb4=[1 1/3 5]; cb5=[1 1 7]; cb6=[1 7 9]; w0=weight1(b); w1=weight1(cb1); w2=weight1(cb2); w3=weight1(cb3); w4=weight1(cb4); w5=weight1(cb5); w6=weight1(cb6); ww=[w1 w2 w3 w4 w5 w6]; score3=ww*w0 %计算总的权值
4.粒子群优化算法的寻优算法--非线性函数极值寻优 2015-06-13
PSO算法首先在可解空间中初始化一群粒子，每个粒子都代表极值最优问题的一个潜在最优解。用位置、速度和适应度三项指标表示该粒子的特征，适应度值由适应度函数计算得到，其值的好坏表示例子的优劣。粒子在解空间中运动，通过跟踪个体极值Pbest和群体极值Gbest跟新个体位置；个体极值Pbest是指个体所经历位置中计算得到的适应度最优解的位置，群体的极值Gbest是指种群中的所有粒子搜索到的适应度最优位置。粒子每更新一次位置，就计算一次适应度值，并通过比较粒子的适应度值和个体极值、群体极值的适应度值来跟新个体极值Pbest和群体极值Gbest位置。
（X1,.......,Xn） 假设在一个D维的搜索空间中有n个粒子组成的种群X，其中第i个粒子
表示为一个D维的向量Xi亦代表[xi1.......xiD]T,代表第i个粒子在D为搜索空间中的位置，
问题的一个潜在解。根据目标函数即可以计算出每个粒子位置Xi值。第i个粒子的速度Vi[xi1.......xiD]T对应的适应度
[Vi1.......ViD]T，其中个体极值为Pi[Pi1.......PiD]T,种群的全局
极值为Pg[Pi1.......PiD]T。
在每一次迭代过程中，粒子通过个体极值和全局极值更新自身的速度和位置，更新公式如下： 速度：Vidk1kkkkkVidc1r1(PidXid)c2r2(PgdXgd) kk1 XidVid位置：Xidk1式中，为惯性权重；c1和c2为非负的常数，称为加速因子；r1和r2为分布于[0,1]之间的随机数。为了防止粒子的盲目搜索，一般建议将奇位置和速度在一定的区间。
例子：求非线性函数的最小值 Min f(x)
如果是带约束的可以通过罚函数法转化为无约束的。 程序如下
所求函数的最小值，函数如下： function y = fun(x) %函数用于计算粒子适应度值
%x input 输入粒子 %y output 粒子适应度值
y=-20*exp(-0.2*sqrt((x(1)^2+x(2)^2)/2))-exp((cos(2*pi*x(1))+cos(2*pi*x(2)))/2)+20+exp(1)*(x(1)<10)*(x(2)<10)+x(3); %y=x(1)^2+x(2)^2+7+sin(x(1)+x(2));
%y=x(1)^2+x(2)^2+7+sin(x(1)+x(2))*((x(1)+x(2))>10);
%y=x(1)^2-10*cos(2*pi*x(1))+10+x(2)^2-10*cos(2*pi*x(2))+10; %y=-x(1)^2-x(2)+x(3)*x(2)*x(1);
主程序： PSO.m clc clear
%% 参数初始化
%粒子群算法中的两个参数 c1 = 1.49445; c2 = 1.49445; %c1=20; %c2=20;
maxgen=500; % 进化次数 sizepop=50; %种群规模
Vmax=1;%速度的范围 Vmin=-1;
popmax=5;%变量的取值范围
popmin=-5;
%% 产生初始粒子和速度 for i=1:sizepop %随机产生一个种群
pop(i,:)=5*rands(1,3); %初始种群 变量的维数3 V(i,:)=rands(1,3); %初始化速度 速度的维数和变量的维数相同 %计算适应度
fitness(i)=fun(pop(i,:)); %染色体的适应度 end
%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness); zbest=pop(bestindex,:); %全局最佳 gbest=pop; %个体最佳
fitnessgbest=fitness; %个体最佳适应度值 fitnesszbest=bestfitness; %全局最佳适应度值
%% 迭代寻优
for i=1:maxgen
for j=1:sizepop
%速度更新 V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
V(j,find(V(j,:)>Vmax))=Vmax; V(j,find(V(j,:) pop(j,:)=pop(j,:)+0.5*V(j,:); pop(j,find(pop(j,:)>popmax))=popmax; pop(j,find(pop(j,:) fitness(j)=fun(pop(j,:)); end for j=1:sizepop %个体最优更新 if fitness(j) < fitnessgbest(j) gbest(j,:) = pop(j,:); fitnessgbest(j) = fitness(j); end %群体最优更新 if fitness(j) < fitnesszbest zbest = pop(j,:);%最优值的位置 fitnesszbest = fitness(j);%最小的函数值 end end yy(i)=fitnesszbest; end %% 结果分析 plot(yy) title('最优个体适应度','fontsize',12); xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12); axis([0 500 0 0.4]) 增加自适应变异 变异操作扩展了在迭代中不断缩小的种群空间，使粒子能够跳出先前搜索到的最优值位置，在更大的空间中开展搜索，同时保持了种群多样性，提高了算法寻找到更优值的可能性。 代码： PSOMutation.m %% 该代码为基于变异粒子群算法的函数极值寻优算法 %% 清空环境 clc clear %% 参数初始化 %粒子群算法中的两个参数 %c1 = 1.49445; %c2 = 1.49445; c1=200; c2=100; maxgen=500; % 进化次数 sizepop=50; %种群规模 Vmax=1; Vmin=-1; popmax=5; popmin=-5; %% 产生初始粒子和速度 for i=1:sizepop %随机产生一个种群 pop(i,:)=5*rands(1,3); %初始种群 V(i,:)=rands(1,3); %初始化速度 %计算适应度 fitness(i)=fun(pop(i,:)); %染色体的适应度 %目标函数 end %% 个体极值和群体极值 [bestfitness bestindex]=min(fitness); zbest=pop(bestindex,:); %全局最佳 gbest=pop; %个体最佳 fitnessgbest=fitness; %个体最佳适应度值 fitnesszbest=bestfitness; %全局最佳适应度值 %% 迭代寻优 for i=1:maxgen %maxgen最大的迭代次数 for j=1:sizepop omg=1;%惯性权重 %速度更新 V(j,:) =omg* V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:)); V(j,find(V(j,:)>Vmax))=Vmax; V(j,find(V(j,:) pop(j,:)=pop(j,:)+0.5*V(j,:); pop(j,find(pop(j,:)>popmax))=popmax; pop(j,find(pop(j,:) pop(j,:)=rands(1,2); %变异操作 end %适应度值 fitness(j)=fun(pop(j,:)); end for j=1:sizepop %个体最优更新 if fitness(j) < fitnessgbest(j) gbest(j,:) = pop(j,:); fitnessgbest(j) = fitness(j); end %群体最优更新 if fitness(j) < fitnesszbest zbest = pop(j,:); fitnesszbest = fitness(j); end end yy(i)=fitnesszbest; end %% 结果分析 plot(yy) title('最优个体适应度','fontsize',12); xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12); axis([0 500 0 0.4]) 惯性权重的选择： PSOMutationOmg.m %% 该代码为基于变异粒子群算法的函数极值寻优算法 %% 清空环境 clc clear %% 参数初始化 %粒子群算法中的两个参数 %c1 = 1.49445; %c2 = 1.49445; c1=200; c2=100; maxgen=500; % 进化次数 sizepop=50; %种群规模 omgstar=0.9;%初始惯性权重 omgend=0.4;%迭代最大次数时的惯性权重 Vmax=1; Vmin=-1; popmax=5; popmin=-5; %% 产生初始粒子和速度 for i=1:sizepop %随机产生一个种群 pop(i,:)=5*rands(1,3); %初始种群 V(i,:)=rands(1,3); %初始化速度 %计算适应度 fitness(i)=fun(pop(i,:)); %染色体的适应度 %目标函数 end %% 个体极值和群体极值 [bestfitness bestindex]=min(fitness); zbest=pop(bestindex,:); %全局最佳 gbest=pop; %个体最佳 fitnessgbest=fitness; %个体最佳适应度值 fitnesszbest=bestfitness; %全局最佳适应度值 %% 迭代寻优 for i=1:maxgen %maxgen最大的迭代次数 for j=1:sizepop omg=omgstar-(omgstar-omgend)*i/maxgen;%惯性权重 %速度更新 V(j,:) =omg* V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:)); V(j,find(V(j,:)>Vmax))=Vmax; V(j,find(V(j,:) pop(j,:)=pop(j,:)+0.5*V(j,:); pop(j,find(pop(j,:)>popmax))=popmax; pop(j,find(pop(j,:) pop(j,:)=rands(1,2); %变异操作 end %适应度值 fitness(j)=fun(pop(j,:)); end for j=1:sizepop %个体最优更新 if fitness(j) < fitnessgbest(j) gbest(j,:) = pop(j,:); fitnessgbest(j) = fitness(j); end %群体最优更新 if fitness(j) < fitnesszbest zbest = pop(j,:); fitnesszbest = fitness(j); end end yy(i)=fitnesszbest; end %% 结果分析 plot(yy) title('最优个体适应度','fontsize',12); xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12); axis([0 500 0 0.4]) 主要的文字说明见MATLAB神经网络43个案例分析35章 5有约束函数极值APSO寻优 2015-06-14 目标函数： fxa1x12x2a2x3x4x214约束条件： 3 0.25x4x32.1952 26x3x46 6.14231010.0282x3600036 2a11x1a2x3x4x2145.0 0.125x1x4 a11.10471a20.04811 0.1x12.0 0.1x210.0 0.1x310.0Bestsolution = 0.1x42.0 0.2124 3.3917 8.09 0.2125 fmin = 1.7501 使用者只需要改前三部分 % 改进的快速粒子群优化算法 (APSO): function apso Lb=[0.1 0.1 0.1 0.1]; %下边界 Ub=[2.0 10.0 10.0 2.0]; %上边界 % 默认参数 para=[25 150 0.95]; %[粒子数，迭代次数，gama参数] % APSO 优化求解函数 [gbest,fmin]=pso_mincon(@cost,@constraint,Lb,Ub,para); % 输出结果 Bestsolution=gbest % 全局最优个体 fmin %% 目标函数 function f=cost(x) f=1.10471*x(1)^2*x(2)+0.04811*x(3)*x(4)*(14.0+x(2)); % 非线性约束 function [g,geq]=constraint(x) % 不等式条件 Q=6000*(14+x(2)/2); D=sqrt(x(2)^2/4+(x(1)+x(3))^2/4); J=2*(x(1)*x(2)*sqrt(2)*(x(2)^2/12+(x(1)+x(3))^2/4)); alpha=6000/(sqrt(2)*x(1)*x(2)); beta=Q*D/J; tau=sqrt(alpha^2+2*alpha*beta*x(2)/(2*D)+beta^2); sigma=504000/(x(4)*x(3)^2); delta=65856000/(30*10^6*x(4)*x(3)^3); F=4.013*(30*10^6)/196*sqrt(x(3)^2*x(4)^6/36)*(1-x(3)*sqrt(30/48)/28); g(1)=tau-13600; g(2)=sigma-30000; g(3)=x(1)-x(4); g(4)=0.10471*x(1)^2+0.04811*x(3)*x(4)*(14+x(2))-5.0; g(5)=0.125-x(1); g(6)=delta-0.25; g(7)=6000-F; % 如果没有等式约束，则置geq=[]； geq=[]; %% APSO Solver function [gbest,fbest]=pso_mincon(fhandle,fnonlin,Lb,Ub,para) if nargin<=4, para=[20 150 0.95]; end n=para(1);% 粒子种群大小 time=para(2); %时间步长，迭代次数 gamma=para(3); %gama参数 scale=abs(Ub-Lb); %取值区间 % 验证约束条件是否合乎条件 if abs(length(Lb)-length(Ub))>0, disp('Constraints must have equal size'); return end alpha=0.2; % alpha=[0,1]粒子随机衰减因子 beta=0.5; % 收敛速度(0->1)=(slow->fast); % 初始化粒子群 best=init_pso(n,Lb,Ub); fbest=1.0e+100; % 迭代开始 for t=1:time, %寻找全局最优个体 for i=1:n, fval=Fun(fhandle,fnonlin,best(i,:)); % 更新最有个体 if fval<=fbest, gbest=best(i,:); fbest=fval; end end % 随机性衰减因子 alpha=newPara(alpha,gamma); % 更新粒子位置 best=pso_move(best,gbest,alpha,beta,Lb,Ub); % 结果显示 str=strcat('Best estimates: gbest=',num2str(gbest)); str=strcat(str,' iteration='); str=strcat(str,num2str(t)); disp(str); fitness1(t)=fbest; plot(fitness1,'r','Linewidth',2) grid on hold on title('适应度') end % 初始化粒子函数 function [guess]=init_pso(n,Lb,Ub) ndim=length(Lb); for i=1:n, guess(i,1:ndim)=Lb+rand(1,ndim).*(Ub-Lb); end %更新所有的粒子 toward (xo,yo) function ns=pso_move(best,gbest,alpha,beta,Lb,Ub) % 增加粒子在上下边界区间内的随机性 n=size(best,1); ndim=size(best,2); scale=(Ub-Lb); for i=1:n, ns(i,:)=best(i,:)+beta*(gbest-best(i,:))+alpha.*randn(1,ndim).*scale; end ns=findrange(ns,Lb,Ub); % 边界函数 function ns=findrange(ns,Lb,Ub) n=length(ns); for i=1:n, % 下边界约束 ns_tmp=ns(i,:); I=ns_tmp %更新粒子 ns(i,:)=ns_tmp; end % 随机性衰减因子 function alpha=newPara(alpha,gamma); alpha=alpha*gamma; % 带约束的d维目标函数的求解 function z=Fun(fhandle,fnonlin,u) % 目标 z=fhandle(u); z=z+getconstraints(fnonlin,u); % 非线性约束 function Z=getconstraints(fnonlin,u) % 罚常数 >> 1 PEN=10^15; lam=PEN; lameq=PEN; Z=0; % 非线性约束 [g,geq]=fnonlin(u); %通过不等式约束建立罚函数 for k=1:length(g), Z=Z+ lam*g(k)^2*getH(g(k)); end % 等式条件约束 for k=1:length(geq), Z=Z+lameq*geq(k)^2*geteqH(geq(k)); end % Test if inequalities function H=getH(g) if g<=0, H=0; else H=1; end % Test if equalities hold function H=geteqH(g) if g==0, H=0; else H=1; end 6.模拟退火算法 TSP问题2015-06-15 已知敌方100个目标的经度、纬度如下表所示： 经度和纬度数据表 53.7121 15.3046 51.1758 0.0322 46.3253 28.2753 30.3313 6.9348 56.5432 21.4188 10.8198 16.2529 22.71 23.1045 10.1584 12.4819 20.1050 15.4562 1.9451 0.2057 26.4951 22.1221 31.4847 8.90 26.2418 18.1760 44.0356 13.5401 28.9836 25.9879 38.4722 20.1731 28.2694 29.0011 32.1910 5.8699 36.4863 29.7284 0.9718 28.1477 8.9586 24.6635 16.5618 23.6143 10.5597 15.1178 50.2111 10.2944 8.1519 9.5325 22.1075 18.5569 0.1215 18.8726 48.2077 16.88 31.9499 17.6309 0.7732 0.4656 47.4134 23.7783 41.8671 3.5667 43.5474 3.9061 53.3524 26.7256 30.8165 13.4595 27.7133 5.0706 23.9222 7.6306 51.9612 22.8511 12.7938 15.7307 4.9568 8.3669 21.5051 24.0909 15.2548 27.2111 6.2070 5.1442 49.2430 16.7044 17.1168 20.0354 34.1688 22.7571 9.4402 3.9200 11.5812 14.5677 52.1181 0.4088 9.5559 11.4219 24.4509 6.5634 26.7213 28.5667 37.5848 16.8474 35.6619 9.9333 24.4654 3.14 0.7775 6.9576 14.4703 13.6368 19.8660 15.1224 3.1616 4.2428 18.5245 14.3598 58.6849 27.1485 39.5168 16.9371 56.50 13.7090 52.5211 15.7957 38.4300 8.48 51.8181 23.0159 8.9983 23.40 50.1156 23.7816 13.7909 1.9510 34.0574 23.3960 23.0624 8.4319 19.9857 5.7902 40.8801 14.2978 58.82 14.5229 18.6635 6.7436 52.8423 27.2880 39.9494 29.5114 47.5099 24.06 10.1121 27.2662 28.7812 27.6659 8.0831 27.6705 9.1556 14.1304 53.79 0.2199 33.90 0.3980 1.3496 16.8359 49.9816 6.0828 19.3635 17.6622 36.9545 23.0265 15.7320 19.5697 11.5118 17.3884 44.0398 16.2635 39.7139 28.4203 6.9909 23.1804 38.3392 19.9950 24.6543 19.6057 36.9980 24.3992 4.1591 3.1853 40.1400 20.3030 23.9876 9.4030 41.1084 27.7149 我方有一个基地，经度和纬度为（70，40）。假设我方飞机的速度为1000公里/小时。 我方派出一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间。 这是一个旅行商问题。我们依次给出基地编号为1，敌方目标依次编号为2，3，...,101,最后我方基地再重复编号102.距离矩阵D(dij)102102，其中dij表示i,j两点的距离，i,j=1,2...,102，这里D为实对称矩阵。则问题是求从一个点1出发，走遍所有中间点，达到点102的一个最短路径。 上面问题中给出的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设A,B两点的地理坐标分别为（x1,y1）,(x2,y2)，过A,B两点的大圆的劣弧即为两点的实际距离。以地心为坐标原点O,以赤道平面为XOY平面，以0度经线圈所在的平面为XOZ平面建立三维直角坐标系。则A,B两点的直角坐标分别为： A（Rcosx1cosy1,Rsinx1siny1） B（Rcosx2cosy2,Rsinx2siny2） 其中R=6370为地球半径。 A,B两点的实际距离为 OA•OBdRarccos() |OA||OB|化简得： d=Rarccos[cos(x1-x2)cosy1cosy2+siny1siny2] 求解模拟退火的步骤见第二十三章 现代优化算法 求解程序 monituihuo.m clc, clear sj0=load('sj.txt'); %加载100个目标的数据，数据按照表格中的位置保存在纯文本文件sj.txt中 x=sj0(:,[1:2:8]);x=x(:); y=sj0(:,[2:2:8]);y=y(:); sj=[x y]; d1=[70,40]; sj=[d1;sj;d1]; sj=sj*pi/180; %角度化成弧度 d=zeros(102); %距离矩阵d初始化 %计算距离矩阵 for i=1:101 for j=i+1:102 d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2))); end end d=d+d'; path=[];long=inf; %巡航路径及长度初始化 rand('state',sum(clock)); %初始化随机数发生器 %求较好的初始解 for j=1:1000 path0=[1 1+randperm(100),102]; %随机产生一个路径 %起始点出发走一圈回到起始点 故路径为102个点 temp=0; for i=1:101 temp=temp+d(path0(i),path0(i+1)); end if temp end end e=0.1^30;%终止温度 L=20000; at=0.999;%降温系数 T=1;%初始温度 for k=1:L %退火过程 c=2+floor(100*rand(1,2)); %产生新解 c=sort(c); c1=c(1);c2=c(2); %计算代价函数值的增量 df=d(path(c1-1),path(c2))+d(path(c1),path(c2+1))-d(path(c1-1),path(c1))-d(path(c2),path(c2+1)); if df<0 %接受准则 path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df; elseif exp(-df/T)>rand %以概率接受新解 path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df; end T=T*at; if T path, long % 输出巡航路径及路径长度 xx=sj(path,1);yy=sj(path,2); plot(xx,yy,'-*') %画出巡航路径 0.70.60.50.40.30.20.1000.20.40.60.811.21.4 7.右端步连续微分方程求解 2015-06-16 考虑下面的具体例子： 10.1yx'x(1x)E1h(x)x1000y0.8xy'y(0.050.25x)Eg(y)y，2 y0.8xE1'0.8E1h(x)(0.8x400)'E20.8E2g(y)(0.9y)0,x492h(x)1,x492 0,y70g(y)1,y70程序如下： dvp2.m function dy=dvp2(t,y) dy=zeros(4,1); p1=.8; c1=400; p2=.9; c2=; r1=1; r2=0.05; b=1/1000; a1=0.1; a2=0.25; m=1; gam=0.8; dy(1)=y(1)*(r1-b*y(1)-a1*y(2)/(m*y(2)^gam+y(1)))-y(3)*y(1)*(y(1)>492); dy(2)=y(2)*(-r2+a2*y(1)/(m*y(2)^gam+y(1)))-y(4)*y(2)*(y(2)>70); dy(3)=0.8*y(3)*(p1*y(1)-c1)*(y(1)>492); dy(4)=0.8*y(4)*(p2*y(2)-c2)*(y(2)>70); 主函数Solve2.m [t,y]=ode45('dvp2',[0,20],[100 50 0.1 0.1]); ca=input('请输入画:x y E1 E2用1 2 3 4 表示:'); if ca==1 plot(t,y(:,1),'r') xlabel('t'); ylabel('x(t)'); title('种群x随时间的变化图') elseif ca==2 plot(t,y(:,2),'r') xlabel('t'); ylabel('y(t)'); title('种群y随时间的变化图') elseif ca==3 plot(t,y(:,3)-0.1,'r') xlabel('t'); ylabel('E1(t)'); title('E1随时间的变化图') else plot(t,y(:,4),'r') xlabel('t'); ylabel('E2(t)'); title('E2随时间的变化图') End 8.多元方差分析 2015-06-17 某养鸡场的鸡蛋育成期的配合饲料主要由5种成分（玉米、麦皮、豆饼、鱼粉和食盐）组成，分别记为A,B,C,D,E，为了研究饲料配方对鸡产蛋效果的影响，对哥成分均选取3个水平，进行了5因素3水平的正交实验，通过试验找出饲料的最佳配方，试验要求考虑交互作用AB、AC和AE。5个因素的水平表如下： 因素 1 2 A 61.5 66 B 6.5 8 C 6 9 D 3 5 E 0 0.1 3 70.6 14 15 9 0.25 对这样的5个因素3水平试验的正交实验表如下： 试验号 A B C E D 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 2 3 1 3 1 2 1 2 3 2 3 1 3 1 2 1 2 3 2 3 1 3 1 2 1 2 3 3 1 2 2 3 1 1 2 3 3 1 2 2 3 1 1 2 3 3 1 2 2 3 1 产蛋量y 569 554 637 566 565 8 581 568 535 593 615 620 586 597 617 599 613 580 569 615 591 586 616 630 566 638 573 根据这些数据分析因素A,B,C,D,E以及交互作用AB,AC和AE对产量Y是否有显著影响，并找出饲料的最佳配方。显著水平为0.05。 程序 siliao.m ydata = xlsread('siliao.xls'); y = ydata(:,7); % 提取ydata的第7列数据，即产蛋量y A = ydata(:,2); % 提取ydata的第2列数据，即因素A的水平列表 B = ydata(:,3); % 提取ydata的第3列数据，即因素B的水平列表 C = ydata(:,4); % 提取ydata的第4列数据，即因素C的水平列表 D = ydata(:,6); % 提取ydata的第6列数据，即因素D的水平列表 E = ydata(:,5); % 提取ydata的第5列数据，即因素E的水平列表 varnames = {'A','B','C','D','E'}; % 定义因素名称 % 定义模型的效应项矩阵，考虑主效应：A,B,C,D,E，交互效应：AB,AC,AE model = [eye(5);1 1 0 0 0;1 0 1 0 0;1 0 0 0 1] % 调用anovan函数作多因素一元方差分析 [p,table] = anovan(y,{A,B,C,D,E},'model',model,'varnames',varnames) p = 0.0237 0.0546 0.0183 0.0124 0.0148 0.1873 0.0229 0.1094 从上面返回的p值向量p和方差分析table可以看出，在5个主效应和3个交互效应中，主效应B和交互效应AB、AE的p值均大于0.05，即在显著水平0.05下，B、AB和AE3个效应是不显著的。下面将两个最不显著的交互效应项AB和AE从模型中去除，重新调用anovan函数作方差分析。 %********************************重新作方差分析***************************** % 定义模型的效应项矩阵，考虑主效应：A,B,C,D,E，交互效应：AC model = [eye(5);1 0 1 0 0] % 调用anovan函数作多因素一元方差分析 [p,table,stats] = anovan(y,{A,B,C,D,E},'model',model,'varnames',varnames); p % 查看p的值 table % 查看table的值 p = 0.0366 0.1128 0.0246 0.0129 0.0173 0.0263 重新进行的检验中只有B的p值大于0.05，即主效应B是步显著的。 下面对5个因素的各水平进行多重比较。 ans = '处理' '均值' 'A=2,B=1,C=3,D=2,E=1' [628.5556] 'A=2,B=1,C=2,D=3,E=1' [ 629] 'A=1,B=1,C=3,D=2,E=1' [629.5556] 'A=2,B=1,C=2,D=2,E=1' [631.2222] 'A=3,B=3,C=2,D=3,E=1' [631.4444] 'A=2,B=2,C=3,D=3,E=1' [631.6667] 'A=1,B=2,C=3,D=3,E=1' [632.6667] 'A=3,B=3,C=2,D=2,E=1' [633.6667] 'A=2,B=2,C=3,D=2,E=1' [633.88] 'A=2,B=2,C=2,D=3,E=1' [634.3333] 'A=1,B=2,C=3,D=2,E=1' [634.88] 'A=3,B=1,C=2,D=3,E=3' [635.5556] 'A=2,B=2,C=2,D=2,E=1' [636.5556] 'A=3,B=1,C=2,D=2,E=3' [637.7778] 'A=3,B=2,C=2,D=3,E=3' [0.88] 'A=3,B=2,C=2,D=2,E=3' [3.1111] 'A=3,B=1,C=2,D=3,E=1' [3.6667] 'A=3,B=1,C=2,D=2,E=1' [5.88] 'A=3,B=2,C=2,D=3,E=1' [9.0000] 'A=3,B=2,C=2,D=2,E=1' [651.2222] 上面调用multcompare函数对5个因素A,B,C,D,E的全部水平3^5=243进行了多重比较，由于返回的结果都是比较长的，无法做到全部显示，上面通过对各处理的均值从小到大进行排序只显示了产蛋量均值最大的后20个处理的名称及相应的均值。这20个处理之间没有显著差异，可以结合本文从中选择饲料的最佳配方。 9.基于MIV的神经网络变量筛选 2015-06-18 本例产生网络训练数据的方法如下： 2F20x1210cos(2x1)x210cos(2x2) 随机产生的x1和x2和由它们决定的F作为BP神经网络的训练样本，同时加入 x3,x4噪声，通过MIV方法，筛选对网络结果主要的影响变量。 MIV方法具体的计算过程：在网络训练终止后，将训练样本P中每一个变量特征在其原值的基础上分别加和减10%构成新的训练样本P1和P2，将P1,P2作为仿真的样本利用已经构建的网络进行仿真，得到仿真结果A1和A2，求出A1和A2的差值，即为变动该自变量后对输出产生影响变化值（impact value），最后讲IV按观测例数平均得出该自变量对于应变量--网络输出的MIV。按照上面步骤依次算出各个自变量的MIV值，最后根据MIV绝对值的大小为自变量排序， 得到各自变量对网络输出影响相对重要性的位次表，从而判断输出特征对于网络结果的影响程度，即实现了变量筛选。 MATLAB实现：bianliangselect.m %% Matlab神经网络43个案例分析 % 神经网络变量筛选—基于MIV的神经网络变量筛选 %% 清空环境变量 clc clear %% 产生输入 输出数据 % 设置步长 interval=0.01; % 产生x1 x2 x1=-1.5:interval:1.5; x2=-1.5:interval:1.5; % 产生x3 x4（噪声） x=rand(1,301); x3=(x-0.5)*1.5*2; x4=(x-0.5)*1.5*2; % 按照函数先求得相应的函数值，作为网络的输出。 F =20+x1.^2-10*cos(2*pi*x1)+x2.^2-10*cos(2*pi*x2); %设置网络输入输出值 p=[x1;x2;x3;x4]; t=F; %% 变量筛选 MIV算法的初步实现（增加或者减少自变量） p=p'; [m,n]=size(p); yy_temp=p; % p_increase为增加10%的矩阵 p_decrease为减少10%的矩阵 for i=1:n p=yy_temp; pX=p(:,i); pa=pX*1.1; p(:,i)=pa; aa=['p_increase' int2str(i) '=p;'];%‘p_increasei=p’;字符表达式 eval(aa);%p_increasei=p数值表达式 即将p矩阵的数值赋值给p_increasei（p_increase1 p_increase2 p_increase3 p_increase4） end for i=1:n p=yy_temp; pX=p(:,i); pa=pX*0.9; p(:,i)=pa; aa=['p_decrease' int2str(i) '=p;'];%‘p_decreasei=p’； eval(aa);%eval将字符串转化为数值。 end % aa=['p_decrease' int2str(i) '=p;']; eval(aa); %命名有规则的变量，并进行赋值。 %% 利用原始数据训练一个正确的神经网络 nntwarn off; p=yy_temp; p=p'; % bp网络建立 net=newff(minmax(p),[8,1],{'tansig','purelin'},'traingdm'); % 初始化bp网络 net=init(net); % 网络训练参数设置 net.trainParam.show=50; net.trainParam.lr=0.05; net.trainParam.mc=0.9; net.trainParam.epochs=2000;%迭代次数 % bp网络训练 net=train(net,p,t); %% 变量筛选 MIV算法的后续实现（差值计算） % 转置后sim %转置操作 for i=1:n eval(['p_increase',num2str(i),'=transpose(p_increase',num2str(i),');']) end for i=1:n eval(['p_decrease',num2str(i),'=transpose(p_decrease',num2str(i),');']) end %网络输出 % result_in为增加10%后的输出 result_de为减少10%后的输出 for i=1:n eval(['result_in',num2str(i),'=sim(net,','p_increase',num2str(i),');' ]) end for i=1:n eval(['result_de',num2str(i),'=sim(net,','p_decrease',num2str(i),');']) end %输出结果转置 for i=1:n eval(['result_in',num2str(i),'=transpose(result_in',num2str(i),');']) end for i=1:n eval(['result_de',num2str(i),'=transpose(result_de',num2str(i),');']) end %% MIV的值为各个项网络输出的MIV值 MIV被认为是在神经网络中评价变量相关的最好指标之一，其符号代表相关的方向，绝对值大小代表影响的相对重要性。 for i=1:n IV= ['result_in',num2str(i), '-result_de',num2str(i)]; eval(['MIV_',num2str(i) ,'=mean(',IV,')']) end 结果： MIV_1 = 1.2030 MIV_2 = 1.0120 MIV_3 = -0.0376 MIV_4 = 0.0773 MIV_n的值为各项网络输出的MIV值，MIV被认为是在神经网络应用中评价变量对结果影响大小的最好指标之一，其符号代表相关的方向，绝对值大小代表影响的相对重要性。 由此可见：第一、二个变量得出的MIV值较大；因为F值是依靠想，x1,x2计算出来的，与x3,x4无关，所以MIV筛选出的对结果有重要影响的自变量和真实情况一致。 10.RBF网络的回归--非线性函数回归的实现 2015-06-19 RBF网络的基本思想：用RBF作为隐单元的“基”构成隐藏层空间，隐含层对输入矢量进行变换，将低维的模式输入数据变换到高维的空间内，使得在低维空间内的线性不可分的问题在高维空间内线性可分。 RBF神经网络结构模型，祥见MATLAB神经网络43个案例分析 RBF神经网络的学习算法，祥见MATLAB神经网络43个案例分析 模型的建立 本例子用RBF网络拟合未知函数，预先设定一个非线性函数，假定函数解析式不清楚的情况下，随机产生x1,x2和由这两个变量按式①得出的y。将x1,x2作为RBF网络的输入，将y作为RBF网络的输出数据，分别建立近似和精确RBF网络进行回归分析，并评价拟合效果。 2y20x1210cos(2x1)x210cos(2x2)① RBF网络的相关函数： （1）newrb() 该函数可以用来设计一个近似径向网络。其调用格式为 [net,tr]=newrb(P,T,GOAL,SPREAD,MN,DF) 其中，P为输入矩阵；T为目标矩阵；GOAL为均方误差目标，默认0.0；SPREAD为径向基函数的扩展速度，默认为1；MN为神经元的最大数目；DF为两次显示之间所添加的神经元数目，默认为25；net为返回值，一个RBF网络；tr为返回值，训练记录。 （2）newrbe() 该函数用于设计一个精确径向基网络。其调用格式为 Net=newrbe(P,T,SPREAD) 例子的MATLAB实现 使用newrbe实现非线性的函数回归 精确回归 %% Matlab神经网络43个案例分析 % RBF网络的回归--非线性函数回归的实现 %% 清空环境变量 clc clear %% 产生输入 输出数据 % 设置步长 interval=0.01; % 产生x1 x2 x1=-1.5:interval:1.5; x2=-1.5:interval:1.5; % 按照函数先求得相应的函数值，作为网络的输出。 F =20+x1.^2-10*cos(2*pi*x1)+x2.^2-10*cos(2*pi*x2); %% 网络建立和训练 % 网络建立 输入为[x1;x2],输出为F。Spread使用默认。 net=newrbe([x1;x2],F) %% 网络的效果验证 % 我们将原数据回带，测试网络效果： ty=sim(net,[x1;x2]); % 我们使用图像来看网络对非线性函数的拟合效果 figure plot3(x1,x2,F,'rd'); hold on; plot3(x1,x2,ty,'b-.'); view(113,36) title('可视化的方法观察准确RBF神经网络的拟合效果') xlabel('x1') ylabel('x2') zlabel('F') grid on 使用newrb实现非线性的函数回归 近似回归 %% Matlab神经网络43个案例分析 % RBF网络的回归--非线性函数拟合的实现 精确拟合 %% 清空环境变量 clc clear %% 产生训练样本（训练输入，训练输出） % ld为样本例数 ld=400; % 产生2*ld的矩阵 x=rand(2,ld); % 将x转换到[-1.5 1.5]之间 x=(x-0.5)*1.5*2; % x的第一行为x1，第二行为x2. x1=x(1,:); x2=x(2,:); % 计算网络输出F值 F=20+x1.^2-10*cos(2*pi*x1)+x2.^2-10*cos(2*pi*x2); %% 建立RBF神经网络 % 采用approximate RBF神经网络。spread为默认值 net=newrb(x,F); %% 建立测试样本 % generate the testing data interval=0.1; [i, j]=meshgrid(-1.5:interval:1.5); row=size(i); tx1=i(:); tx1=tx1'; tx2=j(:); tx2=tx2'; tx=[tx1;tx2]; %% 使用建立的RBF网络进行模拟，得出网络输出 ty=sim(net,tx); %% 使用图像，画出3维图 % 真正的函数图像 interval=0.1; [x1, x2]=meshgrid(-1.5:interval:1.5); F = 20+x1.^2-10*cos(2*pi*x1)+x2.^2-10*cos(2*pi*x2); subplot(1,3,1) mesh(x1,x2,F); zlim([0,60]) title('真正的函数图像') % 网络得出的函数图像 v=reshape(ty,row);%正方形 subplot(1,3,2) mesh(i,j,v); zlim([0,60]) title('RBF神经网络结果') % 误差图像 subplot(1,3,3) mesh(x1,x2,F-v); zlim([0,60]) title('误差图像') set(gcf,'position',[300 ,250,900,400]) 11.极限学习机在回归拟合中的应用 2015-06-20 详见matlab神经网络43个案例分析 训练的函数可以作为神经网络，粒子群算法的适应度函数，进而求最优。这个算法也可以与MIV算法结合，分析变量的影响大小。 ELM的学习算法 ELM在训练之前就随机产生w和b，只需确定隐含层神经元个数及隐含层神经元的激活函数（无限可微），即可以计算出β。具体地，ELM的学习算法主要有一下几个步骤： ①确定隐含层神经元个数，随机设定输入层和隐含层间的链接权值w和隐含层神经元的偏置b； ②选择一个无限可微的函数作为隐含层神经元激活函数，进而计算隐含层输出矩阵H; ③计算输出层权值：HT' 值得一提的是，相关研究结果表明，在ELM中不仅许多非线性激活函数都可以使用，还可以使用步可微函数，甚至可以使用不连续的函数作为激活函数。 程序： ELM训练函数----elmtrain elmtrain()函数作为ELM的创建、训练函数，调用格式如下： [IW,B,LW,TF,TYPE]=elmtrain(P,T,N,TF,TYPE) 其中,P为训练集的输入矩阵；T为训练集的输出矩阵；N为隐含层神经元的个数（默认为训练集的样本数）；TF为隐含层神经元的激活函数；TYPE的应用类型，其取值可以为0（默认，表示回归，拟合）和1（表示分类）；IW为输入层与隐含层间的连接权值；B为隐含层神经元的阀值；LW为隐含层与输出层的连接权值。 ELM预测函数----elmpredict() elmpredict()函数作为ELM的预测函数，其调用格式为： y=elmpredict(P,IW,B,LW,TF,TYPE) 其中,P为测试集的输入矩阵；IW为elmtrain（）函数返回的输入层与隐含层间的连接权值；B为elmtrain（）函数返回的隐含层神经元的阀值；LW为elmtrain（）函数返回的隐含层与输出层的连接权值；TF为与elmtrain()中一致的激活函数；TYPE为与elmtrain()函数中一致的ELM应用类型；y为测试集对应的输出预测值矩阵。 主程序：elmmain.m %% 极限学习机在回归拟合问题中的应用研究 %% 清空环境变量 clear all clc %% 导入数据 %load data input=rand(2000,2);%列数代表变量的个数 output1=input(:,1).^2+input(:,2).^2+sin(input(:,1)); %拟合函数 output=output1'; % 随机生成训练集、测试集 k = randperm(size(input,1)); % 训练集——1900个样本 P_train=input(k(1:1900),:)'; T_train=output(k(1:1900)); % 测试集——100个样本 P_test=input(k(1901:2000),:)'; T_test=output(k(1901:2000)); %% 归一化 % 训练集 [Pn_train,inputps] = mapminmax(P_train,-1,1); Pn_test = mapminmax('apply',P_test,inputps); % 测试集 [Tn_train,outputps] = mapminmax(T_train,-1,1); Tn_test = mapminmax('apply',T_test,outputps); tic %% ELM创建/训练 [IW,B,LW,TF,TYPE] = elmtrain(Pn_train,Tn_train,20,'sig',0); %% ELM仿真测试 Tn_sim = elmpredict(Pn_test,IW,B,LW,TF,TYPE); % 反归一化 T_sim = mapminmax('reverse',Tn_sim,outputps); toc %% 结果对比 result = [T_test' T_sim']; % 均方误差 E = mse(T_sim - T_test) % 决定系数 N = length(T_test); R2 = (N*sum(T_sim.*T_test)-sum(T_sim)*sum(T_test))^2/((N*sum((T_sim).^2)-(sum(T_sim))^2)*(N*sum((T_test).^2)-(sum(T_test))^2)) %% 绘图 figure plot(1:length(T_test),T_test,'r*') hold on plot(1:length(T_sim),T_sim,'b:o') xlabel('测试集样本编号') ylabel('测试集输出') title('ELM测试集输出') legend('期望输出','预测输出') figure plot(1:length(T_test),T_test-T_sim,'r-*') xlabel('测试集样本编号') ylabel('绝对误差') title('ELM测试集预测误差') 12.极限学习机在分类中的应用 2015-06-21 程序 ruxianaimain.m %% 极限学习机在分类问题中的应用研究 %% 清空环境变量 clear all clc warning off %% 导入数据 load ruxianaidata.mat % 随机产生训练集/测试集 a = randperm(569); Train = data(a(1:500),:); Test = data(a(501:end),:); % 训练数据 P_train = Train(:,3:end)'; T_train = Train(:,2)'; % 测试数据 P_test = Test(:,3:end)'; T_test = Test(:,2)'; tic %% ELM创建/训练 [IW,B,LW,TF,TYPE] = elmtrain(P_train,T_train,300,'sig',1); %% ELM仿真测试 T_sim_1 = elmpredict(P_train,IW,B,LW,TF,TYPE); T_sim_2 = elmpredict(P_test,IW,B,LW,TF,TYPE); toc %% 结果对比 result_1 = [T_train' T_sim_1']; result_2 = [T_test' T_sim_2']; % 训练集正确率 k1 = length(find(T_train == T_sim_1)); n1 = length(T_train); Accuracy_1 = k1 / n1 * 100; disp(['训练集正确率Accuracy = ' num2str(Accuracy_1) '%(' num2str(k1) '/' num2str(n1) ')']) % 测试集正确率 k2 = length(find(T_test == T_sim_2)); n2 = length(T_test); Accuracy_2 = k2 / n2 * 100; disp(['测试集正确率Accuracy = ' num2str(Accuracy_2) '%(' num2str(k2) '/' num2str(n2) ')']) %% 显示 count_B = length(find(T_train == 1)); count_M = length(find(T_train == 2)); rate_B = count_B / 500; rate_M = count_M / 500; total_B = length(find(data(:,2) == 1)); total_M = length(find(data(:,2) == 2)); number_B = length(find(T_test == 1)); number_M = length(find(T_test == 2)); number_B_sim = length(find(T_sim_2 == 1 & T_test == 1)); number_M_sim = length(find(T_sim_2 == 2 & T_test == 2)); disp(['病例总数：' num2str(569)... ' 良性：' num2str(total_B)... ' 恶性：' num2str(total_M)]); disp(['训练集病例总数：' num2str(500)... ' 良性：' num2str(count_B)... ' 恶性：' num2str(count_M)]); disp(['测试集病例总数：' num2str(69)... ' 良性：' num2str(number_B)... ' 恶性：' num2str(number_M)]); disp(['良性乳腺肿瘤确诊：' num2str(number_B_sim)... ' 误诊：' num2str(number_B - number_B_sim)... ' 确诊率p1=' num2str(number_B_sim/number_B*100) '%']); disp(['恶性乳腺肿瘤确诊：' num2str(number_M_sim)... ' 误诊：' num2str(number_M - number_M_sim)... ' 确诊率p2=' num2str(number_M_sim/number_M*100) '%']); R = []; %寻找合适的隐含层节点数目 for i = 50:50:800 %% ELM创建/训练 [IW,B,LW,TF,TYPE] = elmtrain(P_train,T_train,i,'sig',1); %% ELM仿真测试 T_sim_1 = elmpredict(P_train,IW,B,LW,TF,TYPE); T_sim_2 = elmpredict(P_test,IW,B,LW,TF,TYPE); %% 结果对比 result_1 = [T_train' T_sim_1']; result_2 = [T_test' T_sim_2']; % 训练集正确率 k1 = length(find(T_train == T_sim_1)); n1 = length(T_train); Accuracy_1 = k1 / n1 * 100; % disp(['训练集正确率Accuracy = ' num2str(Accuracy_1) '%(' num2str(k1) '/' num2str(n1) ')']) % 测试集正确率 k2 = length(find(T_test == T_sim_2)); n2 = length(T_test); Accuracy_2 = k2 / n2 * 100; % disp(['测试集正确率Accuracy = ' num2str(Accuracy_2) '%(' num2str(k2) '/' num2str(n2) ')']) R = [R;Accuracy_1 Accuracy_2]; end figure plot(50:50:800,R(:,2),'b:o') xlabel('隐含层神经元个数') ylabel('测试集预测正确率（%）') title('隐含层神经元个数对ELM性能的影响') 13. 基于PSO改进策略 2015-06-22 基本PSO算法 PSO算法中的微粒的飞行的行为规则类似于鸟类运动，从而使整个微粒群的运动表现出与鸟类觅食类似的特性，进而用于求解复杂的优化问题。该算法的基本思想是通过群体中个体之间的协作和信息共享来寻找最优解。 相比于进化算法，PSO算法保留了基于种群的全局搜索策略，采用简单的速度—位移模型，避免了复杂的遗传操作，同时特有的记忆使算法可动态跟踪当前的搜索情况来调整其搜索策略。 微粒群算法中，每个优化问题的解都是搜索空间中的一只鸟，被抽象为没有质量和体积的微粒，并将其延伸到 N 维空间。微粒 i 在 N 维空间里的位置和速度都表示为一个矢量。PSO算法首先在可行解空间和速度空间随机初始化微粒群，即确定微粒的初始位置和速度。 基本的PSO流程 粒子群算法改进 PSO算法的改进研究可以归纳为两方面：一方面的研究是将各种先进理论引入到PSO算法，研究各种改进和PSO算法；另一方面是将PSO算法和其它智能优化算法相结合，研究各种混合优化算法，达到取长补短、改善算法某方面性能的效果。近时期粒子群改进策略主要体现以下几个方面。 （1）PSO算法的惯性权重模型，通过惯性权重的引入，提高了算法的全局搜索能力； （2）带邻域操作的PSO模型，克服了PSO模型在优化搜索后期随迭代次数增加搜索结果无明显改进的缺点； （3）将拉伸技术用于PSO最小化问题的求解，避免了PSO算法易陷于局部最小值的缺点； （4）用适应度定标的方法对PSO算法进行改进，在算法收敛的前提下能够提高粒子间适应度的差异；在每次迭代中，依据杂交概率选取指定数量的粒子放入一个池中的粒子随机地两两杂交，产生同样数目的孩子粒子，并用孩子粒子代替父母粒子，以保持种群的粒子数目不变； （5）协同PSO算法，其基本思想是用K个相互的粒子群分别在D维的目标搜索空问中的不同维度方向上进行搜索。 粒子群算法主要保证PSO算法的收敛性，采用收敛因子能够确保算法的收敛，PSO收敛因子模型如下： vikvic1rand()(Pixi)c2rand()(gxi)k2224c1c2,4,xixivi, 粒子群算法对于求解极值最优问题应用较为广泛，粒子群算法根据约束范围内的所有可能的粒子，计算适应度值，通过不断的迭代，最终得到相应的极值最优解，粒子群算法能够求解很多复杂的工程问题，避免了问题无可行解的情况。总的来说，粒子群算法算法简单，需要调节的参数不多，许多情况可以按经验值设置参数就能获得较好的收敛结果。其次，算法采用实数编码，可直接取目标函数本身作为适应度函数，根据目标函数值进行迭代搜索。另外，算法的各微粒具有记忆性，算法搜索速度快，大多数情况下，微粒能收敛于最优解。 带惯性权重的PSO算法 vi.j(t1)vi.j(t)c1r1pi.jxi.j(t)c2r2pg.jxi.j(t)惯性权重 w 可以影响微粒的局部最优能力和全局最优能力。 权重线性递减的PSO算法 tmaxminmaxtmaxmax0.9min0.4自适应权重的PSO算法 maxminffmin,ff minavg favgfmin ,ffavgmaxf表示微粒当前的目标函数值。 当各微粒的目标值趋于一致或趋于局部最优时，将使惯性权重增大，而各微粒的目标值比较分散时，使惯性权重减小，同时对于目标函数值优于平均目标值的微粒，其对应的惯性权重因子较小，从而保留了该微粒，反之对于目标函数值差于平均目标值的微粒，其对应的惯性权重因子较大，使得该微粒向较好的搜索区域靠拢。 随机权重策略的PSO算法 N0,1minmaxminrand0,1 其中，N(0，l)表示标准正态分布的随机数，rand(0，l)表示0到1之间的随机数。经研究表明，随机权重策略的PSO算法对于多峰函数，能在一定程度上避免陷入局部最优。该方法多用于动态系统中。 随机地选取 值，使得微粒历史速度对当前速度的影响是随机的。为服从某种随机分布的随机数，从一定程度上可从两方面克服 w的线性递减所带来的不足。 带压缩权重PSO function F=fitness(x) F=0; for i=1:2 F=F+x(i)^2-cos(18*x(i)); end function YSPSO1 N=20;%种群规模 M=100;%最大迭代次数 D=2;%自变量维数 c1=2;%学习因子 c2=3; phi = c1 + c2; if phi <= 4 disp('c1 与 c2 的 和 必 须 大 于 4 ！'); xm = NaN; fv = NaN; return; end format long; %------初始化种群的个体------------ for i=1:N for j=1:D x(i,j)=randn; %随机初始化位置 v(i,j)=randn; %随机初始化速度 end end %------先计算各个粒子的适应度，并初始化Pi和Pg---------------------- for i=1:N p(i)=fitness(x(i,:)); y(i,:)=x(i,:); end pg = x(N,:); %Pg为全局最优 for i=1:(N-1) if fitness(x(i,:)) for t=1:M for i=1:N ksi = 2 / abs(2 - phi - sqrt(phi^2 - 4*phi)); v(i,:) = v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:)); v(i,:) = ksi*v(i,:); x(i,:)=x(i,:)+v(i,:); if fitness(x(i,:)) r=[1:1:100]; plot(r,Pbest,'r--','linewidth',2); xlabel('迭代次数') ylabel('适应度值') title('改进PSO算法收敛曲线') legend('带压缩因子的PSO算法') grid on hold on xm = pg'; fv = fitness(pg); 惯性权重线性递减PSO function F = fitness(x) F=100*(x(1)^2-x(2))^2+(1-x(1))^2; function LinWPSO1 N=20; M=100; D=2; c1=1.5; c2=1.5; wmax=0.9; wmin=0.4; %------初始化种群的个体------------ %N种群规模 %M最大迭代次数 %D自变量的维数 %------初始化种群的个体------------ for i=1:N for j=1:D x(i,j)=randn; %随机初始化位置 v(i,j)=randn; %随机初始化速度 end end %------先计算各个粒子的适应度，初始化Pi和Pg---------------------- for i=1:N p(i)=fitness(x(i,:)); y(i,:)=x(i,:); end pg=x(N,:); %Pg为全局最优 for i=1:(N-1) if fitness(x(i,:)) for t=1:M for i=1:N w = wmax - (t-1)*(wmax-wmin)/(M-1); v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:)); x(i,:)=x(i,:)+v(i,:); if fitness(x(i,:)) if p(i) r=[1:1:100]; plot(r,Pbest,'r--','linewidth',2); xlabel('迭代次数') ylabel('适应度值') title('改进PSO算法收敛曲线') legend('惯性权重线性递减PSO算法') grid on hold on xm = pg'; fv = fitness(pg); 程序\\权重自适应PSO算法 function F = fitness(x) F=100*(x(1)^2-x(2))^2+(1-x(1))^2; function SAPSO1 N=20; M=100; D=2; c1=1.5; c2=1.5; wmax=0.9; wmin=0.9; %------初始化种群的个体------------ %N种群规模 %M最大迭代次数 %D自变量的维数 for i=1:N for j=1:D x(i,j)=randn; %随机初始化位置 v(i,j)=randn; %随机初始化速度 end end %------先计算各个粒子的适应度---------------------- for i=1:N p(i)=fitness(x(i,:)); y(i,:)=x(i,:); end pg=x(N,:); %Pg为全局最优 for i=1:(N-1) if fitness(x(i,:)) for t=1:M for j=1:N fv(j) = fitness(x(j,:)); end fvag = sum(fv)/N; fmin = min(fv); for i=1:N if fv(i) <= fvag w = wmin + (fv(i)-fmin)*(wmax-wmin)/(fvag-fmin); else w = wmax; end v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:)); x(i,:)=x(i,:)+v(i,:); if fitness(x(i,:)) r=[1:1:100]; plot(r,Pbest,'r--','linewidth',2); xlabel('迭代次数') ylabel('适应度值') title('改进PSO算法收敛曲线') legend('权重自适应PSO算法') grid on hold on xm = pg'; fv = fitness(pg); 随机权重PSO算法 function F =fitness(x) F=4*x(1)^2-2.1*x(1)^4+x(1)^6/3+x(1)*x(2)-4*x(2)^2+4*x(2)^4; function RandWPSO1 format long; N=20; M=100; D=2; c1=1.5; c2=1.5; mean_max=0.9; mean_min=0.4; sigma=1; %------初始化种群的个体------------ %N种群规模 %M最大迭代次数 %D自变量的维数 for i=1:N for j=1:D x(i,j)=randn; %随机初始化位置 v(i,j)=randn; %随机初始化速度 end end %------先计算各个粒子的适应度，并初始化Pi和Pg---------------------- for i=1:N p(i)=fitness(x(i,:)); y(i,:)=x(i,:); end pg = x(N,:); %Pg为全局最优 for i=1:(N-1) if fitness(x(i,:)) end end %------进入主要循环，按照公式依次迭代------------ for t=1:M for i=1:N miu = mean_min + (mean_max - mean_min)*rand(); w = miu + sigma*randn(); v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:)); x(i,:)=x(i,:)+v(i,:); if fitness(x(i,:)) y(i,:)=x(i,:); end if p(i) end end Pbest(t)=fitness(pg); end r=[1:1:100]; plot(r,Pbest,'r--','linewidth',2); xlabel('迭代次数') ylabel('适应度值') title('改进PSO算法收敛曲线') legend('随机权重PSO算法') grid on hold on xm = pg'; fv = fitness(pg); 14.神经网络遗传算法函数极值寻优 2015-06-23 步骤： 1.神经网络训练 2.训练好的神经网络作为适应度函数。（如果函数已知第一步骤可以省略） 3.遗传算法寻最优值 BP神经网络训练 BP.m %% 清空环境变量 clc clear tic %% 训练数据预测数据提取及归一化 %下载输入输出数据 load data1 input output %从1到2000间随机排序 k=rand(1,4000); [m,n]=sort(k); %找出训练数据和预测数据 input_train=input(n(1:3900),:)'; output_train=output(n(1:3900),:)'; input_test=input(n(3901:4000),:)'; output_test=output(n(3901:4000),:)'; %选连样本输入输出数据归一化 [inputn,inputps]=mapminmax(input_train); [outputn,outputps]=mapminmax(output_train); %% BP网络训练 % %初始化网络结构 net=newff(inputn,outputn,5); net.trainParam.epochs=100; net.trainParam.lr=0.1; net.trainParam.goal=0.0000004; %网络训练 net=train(net,inputn,outputn); %% BP网络预测 %预测数据归一化 inputn_test=mapminmax('apply',input_test,inputps); %网络预测输出 an=sim(net,inputn_test); %网络输出反归一化 BPoutput=mapminmax('reverse',an,outputps); %% 结果分析 figure(1) plot(BPoutput,':og') hold on plot(output_test,'-*'); legend('预测输出','期望输出','fontsize',12) title('BP网络预测输出','fontsize',12) xlabel('样本','fontsize',12) ylabel('输出','fontsize',12) %预测误差 error=BPoutput-output_test; figure(2) plot(error,'-*') title('神经网络预测误差') figure(3) plot((output_test-BPoutput)./BPoutput,'-*'); title('神经网络预测误差百分比') errorsum=sum(abs(error)) toc save data net inputps outputps 适应度函数 fun.m function fitness = fun(x) % 函数功能：计算该个体对应适应度值 % x input 个体 % fitness output 个体适应度值 % load data net inputps outputps %数据归一化 x=x'; inputn_test=mapminmax('apply',x,inputps); %网络预测输出 an=sim(net,inputn_test); %网络输出反归一化 fitness=mapminmax('reverse',an,outputps); 遗传算法寻优 Genetic.m %% 清空环境变量 clc clear %% 初始化遗传算法参数 %初始化参数 maxgen=100; %进化代数，即迭代次数 sizepop=20; %种群规模 pcross=[0.4]; %交叉概率选择，0和1之间 pmutation=[0.2]; %变异概率选择，0和1之间 lenchrom=[1 1]; %每个变量的字串长度，如果是浮点变量，则长度都为1 bound=[-5 5;-5 5]; %数据范围 individuals=struct('fitness',zeros(1,sizepop), 'chrom',[]); %将种群信息定义为一个结构体 avgfitness=[]; %每一代种群的平均适应度 bestfitness=[]; %每一代种群的最佳适应度 bestchrom=[]; %适应度最好的染色体 %% 初始化种群计算适应度值 % 初始化种群 for i=1:sizepop %随机产生一个种群 individuals.chrom(i,:)=Code(lenchrom,bound); x=individuals.chrom(i,:); %计算适应度 individuals.fitness(i)=fun(x); %染色体的适应度 end %找最好的染色体 [bestfitness bestindex]=min(individuals.fitness); bestchrom=individuals.chrom(bestindex,:); %最好的染色体 avgfitness=sum(individuals.fitness)/sizepop; %染色体的平均适应度 % 记录每一代进化中最好的适应度和平均适应度 trace=[avgfitness bestfitness]; %% 迭代寻优 % 进化开始 for i=1:maxgen % 选择 individuals=Select(individuals,sizepop); avgfitness=sum(individuals.fitness)/sizepop; %交叉 individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop,bound); % 变异 individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,[i maxgen],bound); % 计算适应度 for j=1:sizepop x=individuals.chrom(j,:); %解码 individuals.fitness(j)=fun(x); end %找到最小和最大适应度的染色体及它们在种群中的位置 [newbestfitness,newbestindex]=min(individuals.fitness); [worestfitness,worestindex]=max(individuals.fitness); % 代替上一次进化中最好的染色体 if bestfitness>newbestfitness bestfitness=newbestfitness; bestchrom=individuals.chrom(newbestindex,:); end individuals.chrom(worestindex,:)=bestchrom; individuals.fitness(worestindex)=bestfitness; avgfitness=sum(individuals.fitness)/sizepop; trace=[trace;avgfitness bestfitness]; %记录每一代进化中最好的适应度和平均适应度 end %进化结束 %% 结果分析 [r c]=size(trace); plot([1:r]',trace(:,2),'r-'); title('适应度曲线','fontsize',12); xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12); axis([0,100,0,1]) disp(' 适应度 变量'); x=bestchrom; % 窗口显示 disp([bestfitness x]);
因篇幅问题不能全部显示，请点此查看更多更全内容
Copyright © 2019- baoaiwan.cn 版权所有 赣ICP备2024042794号-3
违法及侵权请联系：TEL:199 18 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务