数学建模必备程序doc.docx
《数学建模必备程序doc.docx》由会员分享，可在线阅读，更多相关《数学建模必备程序doc.docx（22页珍藏版）》请在冰豆网上搜索。
数学建模必备程序doc
一基于均值生成函数时间序列预测算法程序
1.predict」un.m为主程序；
2.timeseries.m和serie**pan.m为调川的子程序
functionima_pre=predict_fun（b,step）
%mainprograminvokestimeseries.mandserie**pan.m%inputparameters:
%bthetrainingdata（vector）;
%step-mberofpredictiondata;
%outputparameters:
%ima_pre--thepredictiondata（vector）;old_b=b;
mean_b=sum（old_b）/length（old_b）;std_b=std（old_b）;
old_b=（oId_b-mean_b）/std_b;
[f,x]=timeseries（old_b）;
old_f2=serie**pan（old_b,step）;
%f（f<0.0001&f>-0.0001）=f（f<0.0001&f>-0.0001）+eps;R=corrcoef（f）;
[eigvectoreigroot]=eig（R）;
eigroot=diag（eigroot）;
a=eigroot（end:
-1:
1）;
vector=eigvector（:
end:
-1:
1）;
Devote=a./sum（a）;
Devotem=cumsum（Devote）;m=find（Devotem>=0.995）;
m=m
（1）;
V1=f*eigvector';
V=V1（:
1:
m）;
%old_b=old_b;
old」ai=inv（V,*V）*V,*old_b;
eigvector=eigvector（1:
m,1:
m）;
fai=eigvector*old_fai;
f2=oldI_f2（:
1:
m）;
predictvalue=f2*fai;
ima_pre=std_b*predictvalue+mean_b;
1•子函数：
timeseries.m
%timeseriesprogram%
%thisprogramisusedtogeneratemeanvaluematrixf;
function[f,x]=timeseries（data）
%datatheinputsequenee（vector）;
%f-meanvaluematrixf;
n=length（data）;
forL=1:
n/2
nL=floor（n/L）;
fori=1:
L
sum=0;
forj=1:
nLsum=sum+data（i+（j-1）*L）;
end
x{L,i}=sum/nL;
end
end
L=n/2;
f=zeros（n3L）;
fori=1:
L
rep=floor（n/i）;
res=mod（n,i）;
b=[x{i,1:
i}];b=b';
f（1:
rep*i,i）=repmat（b,rep,1）;
ifres~=0
c=rep*i+1:
n;
f（rep*i+1:
end,i）=b（1:
length（c））;
end
end
%serie**pan.m
%theprogramisusedtogeneratethepredictionmatrixf;functionf=serieipan（data,step）;
%data-•-theinputsequenee（vector）
%setp-—thepredictionnumber;
n=length（data）;
forL=1:
n/2
nL=floor（n/L）;
fori=1:
L
sum=0;
forj=1:
nLsum=sum+data（i+（j-1）*L）;
end
x{L,i}=sum/nL;
end
end
L=n/2;
f=zeros（n+step丄）；
fori=1:
L
rep=floor（（n+step）/i）;
res=mod（n+step3i）;
b=[x{i,1:
i}];b=b,;
f（1:
rep*ij）=repmat（b,rep,1）;
ifres~=0
c=rep*i+1:
n+step;
f（rep*i+1:
end,i）=b（1:
length（c））;
end
end
二最短路Dijkstra算法
%dijkstraalgorithmcodeprogram%
%theshortestpathlengthalgorithm
function[path,short_distanee]=ShortF^th_Dijkstra（Input_weight,start,endpoint）
%Inputparameters:
%Input_weighttheinputnodeweight!
%startthestartnodenumber;
%endpointtheendnodenumber;
%Outputparameters:
%paththeshortestlenghtpathfromthestartnodetoendnode;
%short_distance・……thedistaneeoftheshortestlenghtpathfromthe
%startnodetoendnode.
[row,col]=size（lnput_weight）;
%inputdetection
ifrow~=col
error（1nputmatrixisnotasquarematrix,inputerror1）;
end
ifendpoint>row
error（*inputparameterendpointexceedthemaximalpointnumberf）;
end
%initialization
s_path=[start];
distance=inf*ones（1,row）;distance（start）=O;
flag（start）=start;temp=start;
whilelength（s_path）vrow
pos=find（lnput_weight（temp,:
）~=inf）;
fori=1:
length（pos）
if（length（find（s_path==pos（i）））==0）&
（distance（pos（i））>（distance（temp）+lnput_weight（temp,pos（i））））
distance（pos（i））=distance（temp）+lnput_weight（temp,pos（i））;
什ag（pos（i））=temp;
end
end
k=inf;
fori=1:
row
if（length（find（s_path==i））==O）&（k>distance（i））k=distance（i）;
temp_2=i;
end
end
s_path=[s_path,temp_2];
temp=temp_2;
end
%outputtheresult
path
（1）=endpoint;
i=1;
whilepath（i）^=start
path（i+1）=flag（path（i））;
i=i+1;
end
path（i）=start;
path=path（end>1:
1）;
shortdistance=distance（endpoint）;
三绘制差分方程的映射分叉图
functionfork1（a）;
%绘制x_（n+1）=1-a*xA2_n映射的分叉图
%Example:
%fork1（[0,2]）;
N=300;%取样点数
A=linspace（a
（1）,a
（2）,N）;
starx=0.9;
Z=[];
h=waitbar（O,'pleasewaitf）;m=1;
forap=A;
x=starx;
fork=1:
50;
x=1-ap*xA2;
end
fork=1:
201;
x=1-ap*xA2;
Z=[Z,ap-x*i];
end
waitbar（m/N,h,['completed',num2str（round（100*m/N））,*%']^）;
m=m+1;
end
delete（h）;
plot（Z；.\,markersize\2）
xlim（a）;
四最短路算法——floyd算法
functionShortF^th_floyd（w,start,terminal）
%w-—adjoinmatrix,w=[050infinfinf;inf0infinf80;
%inf30020inf;infinfinf070;65inf100inf0];
%startthestartnode;
%terminal-theendnode;
n=size（w,1）;
[D,path]=floyd1（w）;%调用floyd算法程序
%找出任意两点之间的最短路径，并输出
fori=1:
n
forj=1:
n
MinjDath（iJ）.distance=D（i,j）;
%将匚到j的最短路程赋值Min_path（i,j）.distance
%将i至I」j所经路径赋给Min_path（i,j）.path
MinjDath（i,j）.path
（1）=i;
k=1;
whileMin_path（i,j）.path（k）~=j
k=k+1;
Min_path（i,j）.path（k）=path（Min_path（i,j）.path（k-1）,j）;
end
end
end
s=sprintf（'任意两点之间的最短路径如下：
’）；
disp（s）;
fori=1:
n
forj=1:
n
s=sprintf（'从％d到％（1的最短路径长度为：
%d\n所经路径为:
J,j5Min_path（iJ）.distanee）;
disp（s）;
disp（Min_path（i,j）.path）;
end
end
%找出在指定从start点到terminal点的最短路径，并输出
str1=sprintf（%%d到％（1的最短路径长度为:
%d\n所经路径为:
'，…start,terminal,Min_path（start,terminal）.distanee）;
disp（strl）;
disp（MinjDath（start,terminal）.path）;
%Foldy'sAlgorithm算法程序function[D,path]=floyd1（a）n=size（a,1）;
D=a;path=zeros（n,n）;%设置D和path的初值fori=1:
n
forj=1:
n
ifD（i,j）~=inf
path（i,j）=j;%j是i的后点
end
end
end
%做n次迭代，每次迭代都更新D（i,j）和path（i,j）fork=1:
n
fori=1:
n
forj=1:
n
ifD（i,k）+D（k,j）vD（i,j）D（i,j）=D（i,k）+D（k,j）;%修改长度path（i,j）=path（i,k）;%修改路径end
end
end
end
五模拟退火算法源程序
function[MinD,BestF^th]=MainAneal（QtyRosition3pn）
function[MinD,BestF^th]=MainAneal2（QtyFbsition,pn）
%此题以屮国31省会城市的最短旅行路径为例，给出TSP问题的模拟退火程序
%QtyFbsition_31=[13042312;36391315;41772244;37121399;3488
1556;..・
%
1756；...
%
2370；...
%
2376;...
%
2826；23702975]；
1535;3326
3238
2788
3780
3394
1229;4196
1044;4312
790;4386
570;3007
1970;2562
1491;2381
1676；1332
695；3715
1678;3918
2179;4061
2212;3676
2578;4029
2838;4263
2931;3429
1908;3507
2643;3439
3201;2935
3240;3140
3550;2545
2357;2778
%T0=clock
globalpathp2D;
[m,n]=size（QtyFbsition）;
%牛成初始解空间，这样可以比逐步分配空间运行快i些
TraceF^th=zeros（1e3,m）;
Distance=inf*zeros（1,1e3）;
D=sqrt（（QtyFbsition（:
ones（1,m））-OtyFbsition（:
ones（1,m））'）.A2+...（QtyFbsition（:
2*ones（1,m））・QtyFbsition（:
2*ones（1,m））*）.A2）;
%将城市的坐标矩阵转换为邻接矩阵（城市间距离矩阵）
fori=1:
pn
path（i,:
）=randperm（m）;%构造一个初始可行解
end
t=zeros（1,pn）;
p2=zeros（1,m）;
iter_max=100;%input（'ih输入固定温度下最人迭代次数iter_max=*）;m_max=5;%inputC请输入El定温度下日标函数值允许的最大连续未改进次数m_nax=*）;
%如果考虑到降温初期新解被吸收概率较人，容易陷入局部最优
%而随着降温的进行新解被吸收的概率逐渐减少，乂难以跳出局限
%人为的使初期iter_max,m_max较小，然后使Z随温度降低而逐步增大，可能
%会收到到比较好的效果
T=1e5;
N=1;
tau=1e-5;%input（*iR输入最低温度tau='）;
%nn=ceil（log10（tau/T）/log10（0.9））;
whileT>=tau%&m_numvm_max
iter_num=1;%某固定温度下迭代计数器m_num=1;%某固定温度下目标函数值连续未改进次数计算器%iter_max=100;
%m_max=10;%ceil（10+0.5*nn-0.3*N）;
whilem_num%MRF0T（Metropolis,Rosenbluth,Rosenbluth,Teller,Teller）过程：
%川任意启发式算法在path的领域叫path）中找出新的更优解
fori=1:
pn
Len1（i）=sum（[D（path（i,1:
m-1）+m*（path（i,2:
m）-1））D（path（i,m）+m*（path（i,1）-1））]）;
%计算一次行遍所有城市的总路程
[path2（i,:
）]=ChangePath2（path（i,:
）,m）;%更新路线
Len2（i）=sum（[D（path2（i,1:
m-1）+m*（path2（i,2:
m）-1））D（path2（i,m）+m*（path2（i,1）-1））]）;
end
%Len1
%Len2
%ifLenyLenlAo-exp（（LenrLen2）、（T））Vrand
R=rand?
pn）_
%Len2,Len1A二exp（（LenrLen2）二T））VRiffind（（Len2,Lerdci-exp（（LenrLen2）、（T））vRTHO）path（find（（Len2-LerdAi-exp（（LenrLen2）、（T））VR）〜H0L二Hpaih2（find（（Len2,Lenlci_ex
p（（Len「Len2）、（T}vR）；0r）_
Lenl（find（（Len2,LerHct_exp（（LenrLen2LenrLen2）、（T））VR）〜H0））八
（TempMinpRmp-ndex」nmin（Lerd）-
%TempMinD
TraceF^ihz..Tpaih（Temp-ndex..「
Disiance（N..TTempMinp
NHN++
%TUT69
m—numHp
e-se
m—numHm—num+Jt
end
iiernumHiiernum+t
end
THT69%m—nummer—num>end
【MinDJndex」nmin（Distance）八BesiF^hHTracepaihundex..）八disp（MinD）
%Tfc-ock
%M連察function（p2Tchangepa52（plQiyNum）
g-oba-p2;
whi-e（l）
RHunidmd（CKyNum=2）八
ifabs（R（vR
（2））vl
break」
end
R==unidmd（cityNum・12）-
I=R
（1）;J=R
（2）;
%len1=D（p（l）,p（J））+D（p（l+1）,p（J+1））;
%len2=D（p（l）,p（l+1））+D（p（J）,p（J+1））;
iflp2（1:
l）=p1（1:
l）;
p2（l+1:
J）=p1（J:
-1:
l+1）;
p2（J+1:
QtyNum）=p1（J+1:
QtyNum）;
else
p2（1：
J）=p1（1:
J）；
p2（J+1：
l）=p1（l：
-1：
J+1）；
p2（l+1:
QtyNum）=p1（l+1:
QtyNum）;
end
六遗传算
法程序：
说明：
为遗传算法的主程序；采用二进制&ay编码，采用基于轮盘赌法的非线性排名
选择，均匀交义，变异操作，而且还引入了倒位操作！
function
[BestFbp?
Trace]=fga（FUN,LB,UB,eranum,popsize,pCross,pMutation,plnversion,options）
%[BestFbp,Trace]=fmaxga（FUN,LB,UB,eranum5popsize,pcross,pmutation）
%Findsamaximumofafunctionofseveralvariables.
%fmaxgasolvesproblemsoftheform:
%
maxF（X）
subjectto:
l_B<=X<=UB
%
BestFtp
・最优的群体即为授优的染色体群
%
Trace
・最佳染色体所对应的口标惭数值
%
FUN
・目标函数
%
LB
・自变量下限
%
UB
・自变量上限
%
eranum
-种群的代数,JR100-1000（默认200）
%
popsize
・每一•代种群的规模；此可取50-200（默认100）
%
pcross
・交义概率，一般取0.5-0.85之间较好（默认0.8）
%
pmutation
-初始变异概率，一般取0.05-0.2Z间较好（默认0.1）
%
pinversion
•倒位概率,一般取0.05—0.3Z间较好（默认0.2）
%
options
-T2矩阵‘options⑴=0二进制编码（默认0）,option⑴~=0十进制编
%码,option
（2）设定求解精度（默认1e-4）
%
%
T1=clock;
ifnargin<3,error（fFMAXGArequiresatleastthreeinputarguments'）;end
ifnargin==3,
eranum=200;popsize=100;pCross=0.8;pMutation=0.1;plnversion=0.15;options=[01e-4];end
ifnargin==4,popsize=100;pCross=0.8;pMutation=0.1;plnversion=0.15;options=[0
1e-4];end
ifnargin==5,pCross=0.8;pMutation=0.1;plnversion=0.15;options=[01e-4];end
ifnargin==6,pMutation=0.1;plnversion=0.15;options=[01e-4];end
ifnargin==7,pinversion=0.15;options=[01e-4];end
iffind（（l_B-UB）>0）
errorC数据输入错误，请重新输入（l_B*）;
end
s=sprin廿（'程序运行需要约％.4f秒钟时间，请稍等……',（eranum*popsize/1000））;
disp（s）;
globalmnNewFbpchildren1children2VarNum
bounds=[LB;UB]f;bits=[];VarNum=size（bounds,1）;precision=options
（2）;%dL]求解精度确定二进制编码长度bits=ceil（log2（（bounds（:
2）-bounds（:
1））*./precision））;%由设定精度划分区间[Fbp]=InitFbpGray（popsize,bits）;%初始化种群
[m,n]=size（Fbp）;
NewFtop=zeros（m,n）;
children1=zeros（15n）;
children2=zeros（1,n）;
pmO=pMutation;
BestFbp=zeros（eranum,n）;%分配初始解空间BestFbp,Trace
Trace=zeros（eranumjength（bits）+1）;
i=1；
whilei<=eranum
升级会员 