MATLAB数学建模14个范例

1.整数规划的蒙特卡洛解法2015-06-10 (2)
2. 罚函数法 2015-06-11 (3)
3. 层次分析 2015-06-12 (4)
4. 粒子群优化算法的寻优算法--非线性函数极值寻优 2015-06-13 (5)
5有约束函数极值APSO寻优 2015-06-14 (12)
6.模拟退火算法 TSP问题2015-06-15 (17)
7. 右端步连续微分方程求解2015-06-16 (19)
8. 多元方差分析 2015-06-17 (22)
9. 基于MIV的神经网络变量筛选 2015-06-18 (25)
10. RBF网络的回归--非线性函数回归的实现 2015-06-19 (29)
11. 极限学习机在回归拟合中的应用 2015-06-20 (32)
12. 极限学习机在分类中的应用 2015-06-21 (34)
13. 基于PSO改进策略 2015-06-22 (37)
14. 神经网络遗传算法函数极值寻优 2015-06-23 (46)
1.
1.整数规划的蒙特卡洛解法2015-06-10 已知非线性整数规划为：
⎪⎪⎪⎩⎪
⎪⎪
⎨⎧≤++≤++≤++++≤++++=≤≤-----++++=200
5200
62800622400)5,....,1(9902328243max 54233
21
64321543215
43212
524232221x x x x x x x x x x x x x x x x i x x x x x x x x x x x z i
如果用显枚举试探，共需要计算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 。

罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。

主要有两种形式，一种叫外函数法，另一种叫内函数法，下面介绍外函数法。

考虑问题：
⎪
⎩⎪
⎨⎧===≥=≤t m x k s j x h r
i x g t s x f m j i .....
1,0)(,....1,0)(......1,0)(..)
(min 取一个充分大的数M>0,构造函数： ∑∑∑===+-+=t
i i
s
i i
r
i i
x k M x h M x g M
x f M x P 1
1
1
|)(|)0),((min )0),(max ()(),(
增广目标函数P(x,M)为目标函数的无约束极值问题
minP(x,M)
的最优解x 也是原问题的最优解。

Eg:求解下列非线性规划问题
⎪⎩⎪⎨⎧≥=+--≥-++=0,020..8
)(min 2
12212212
22
1x x x x x x t s x x x f
解 编写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最合适。

准则B1 B2 B3 B4 B5 B6 总权值
准则层权值0.1600 0.1600 0.1600 0.0400 0.1600 0.3200 工作1 0.1429 0.1000 0.2308 0.2381 0.4667 0.7975 0.4152 工作2 0.5714 0.4000 0.0769 0.7143 0.4667 0.1139 0.3074 工作3 0.2857 0.5000 0.6923 0.0476 0.0667 0.0886 0.2774
权重的计算
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);
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 位置。

假设在一个D 维的搜索空间中有n 个粒子组成的种群）（n 1,.......,X X X =，其中第i 个粒子
表示为一个D 维的向量T iD i x x X ].......[1i
=,代表第i 个粒子在D 为搜索空间中的位置，
亦代表问题的一个潜在解。

根据目标函数即可以计算出每个粒子位置T iD i x x X ].......[1i =对应的适应度
值。

第i 个粒子的速度T iD i V V V ].......[1i
=，其中个体极值为T iD i P P P ].......[1i =,种群的全局
极值为T iD i P P P ] (1)
=。

在每一次迭代过程中，粒子通过个体极值和全局极值更新自身的速度和位置，更新公式如下： 速度：)()(22111
k gd k gd k id k id k id k id
X P r c X P r c V V -+-+=+ω 位置：1
id
id 1
id
+++=k k k V X X 式中，ω为惯性权重；21c c 和为非负的常数，称为加速因子；21r r 和为分布于[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,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:);
pop(j,find(pop(j,:)>popmax))=popmax;
pop(j,find(pop(j,:)<popmin))=popmin;
%适应度值
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,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:);
pop(j,find(pop(j,:)>popmax))=popmax;
pop(j,find(pop(j,:)<popmin))=popmin;
if rand>0.98
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,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:);
pop(j,find(pop(j,:)>popmax))=popmax;
pop(j,find(pop(j,:)<popmin))=popmin;
if rand>0.98
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
目标函数：
约束条件：
Bestsolution =
0.2124 3.3917 8.8909 0.2125
fmin =
1.7501 ()()
2112234214f x a x x a x x x =++3430.25 2.1952
x x
≥()36.14231010.02826000x ⨯-≥()()
2112342114 5.0a x x -+≤140.125x x ≤≤20.04811a =1 1.10471a =10.1 2.0x ≤≤20.110.0x ≤≤30.110.0x ≤≤40.1 2.0x ≤≤
使用者只需要改前三部分
% 改进的快速粒子群优化算法 (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<Lb;
ns_tmp(I)=Lb(I);
% 上边界约束
J=ns_tmp>Ub;
ns_tmp(J)=Ub(J);
%更新粒子
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;
6.模拟退火算法 TSP 问题2015-06-15
已知敌方100个目标的经度、纬度如下表所示：
我方有一个基地，经度和纬度为（70，40）。

假设我方飞机的速度为1000公里/小时。

我方派出一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。

在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间。

这是一个旅行商问题。

我们依次给出基地编号为1，敌方目标依次编号为2，3，...,101,最后我方基地再重复编号102.距离矩阵102102ij )d (⨯=D ，其中ij d 表示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 两点的实际距离为
)|
|||arccos(d →→→
→
•=OB OA OB OA R 化简得：
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(s j(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<long
path=path0; long=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<e
break; 退火温度
end
end
path, long % 输出巡航路径及路径长度
xx=sj(path,1);yy=sj(path,2);
plot(xx,yy,'-*') %画出巡航路径
7.右端步连续微分方程求解2015-06-16
考虑下面的具体例子：
，⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧-=-=-++-=-+--=)
649.0)((0.8)4008.0)((0.8)()25.005.0(')()1.0100011('2'21'128.018.0y y g E E x x h E E y y g E x y x y y x x h E x y y x x x ⎩⎨⎧≥<=⎩⎨⎧≥<=70
,170,0)(492
,1492,0)(y y y g x x x h
程序如下：
dvp2.m
function dy=dvp2(t,y)
dy=zeros(4,1);
p1=.8;
c1=400;
p2=.9;
c2=64;
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水平的正交实验，通过试验找出饲料
根据这些数据分析因素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.8889] 'A=2,B=2,C=2,D=3,E=1' [634.3333] 'A=1,B=2,C=3,D=2,E=1' [634.8889] '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' [640.8889] 'A=3,B=2,C=2,D=2,E=3' [643.1111] 'A=3,B=1,C=2,D=3,E=1' [643.6667] 'A=3,B=1,C=2,D=2,E=1' [645.8889] 'A=3,B=2,C=2,D=3,E=1' [649.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 本例产生网络训练数据的方法如下：
)2cos(10)2cos(102022
2121x x x x F ππ-+-+=
随机产生的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 网络进行回归分析，并评价拟合效果。

)2cos(10)2cos(1020y 22
2121x x x x ππ-+-+=①
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; ③计算输出层权值'T H +=ββ：
值得一提的是，相关研究结果表明，在ELM 中不仅许多非线性激活函数都可以使用，还可以使用步可微函数，甚至可以使用不连续的函数作为激活函数。

程序：
ELM 训练函数----elmtrain
e lmtrain()函数作为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。