技术标签： matlab
曲线拟合又叫做曲线，和插值函数不同的是，它只要求拟合的曲线能反映数据的基本趋势，并不要求曲线一定过数据点。
曲线拟合的判别准则有：
数据拟合的方法有：
x=[1 2 3 4 5 6 7 8 9];
y=[9 7 6 3 -1 2 5 7 20];
P=polyfit(x,y,3);%多项式拟合，返回降幂排序的多项式系数，3代表拟合的最高次幂
xi=0:.2:10;
yi=polyval(P,xi);%计算多项式值
plot(xi,yi,x,y,'r*')
syms t
x=[0;0.4;1.2;2;2.8;3.6;4.4;5.2;6;7.2;8;9.2;10.4;11.6;12.4;13.6;14.4;15];
y=[1;0.85;0.29;-0.27;-0.53;-0.4;-0.12;0.17;0.28;0.15;-0.03;-0.15;-0.071;0.059;0.08;0.032;-0.015;-0.02];
f=fittype('a*cos(k*t)*exp(w*t)','independent','t','coefficients',{'a','k','w'});%自定义拟合函数
cfun=fit(x,y,f)%根据f来拟合数据x,y。注意：数据必须为行向量
xi=0:.1:20;
yi=cfun(xi);
plot(x,y,'r*',xi,yi,'b-');
人口预测模型
人口随时间的变化可以用logistic曲线模型拟合
logistic基本形式为
y = 1 a + b e − t y=\frac{1}{a+be^{-t}} y=a+be−t1
令
y ′ = 1 y , x ′ = e − t y'=\frac{1}{y},x'=e^{-t} y′=y1,x′=e−t
转换为直线模型
y ′ = a + b x ′ y'=a+bx' y′=a+bx′
代码如下
%population.m
%% MATLAB Version: 2013a
clc, clear all, close all
% 读入人口数据（1971－2000年）
Y=[33815 33981 34004 34165 34212 34327 34344 34458 34498 34476 34483 34488 34513 34497 34511 34520 34507 34509 34521 34513 34515 34517 34519 34519 34521 34521 34523 34525 34525 34527]
% 读入时间变量数据（t＝年份－1970）
T=[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 28 29 30]
% 线性化处理
for t = 1:30,
x(t)=exp(-t);
y(t)=1/Y(t);
end
% 计算，并输出回归系数B
c=zeros(30,1)+1;
X=[c,x'];
B=inv(X'*X)*X'*y'
for i=1:30,
% 计算回归拟合值
z(i)=B(1,1)+B(2,1)*x(i);
% 计算离差
s(i)=y(i)-sum(y)/30;
% 计算误差
w(i)=z(i)-y(i);
end
% 计算离差平方和S
S=s*s';
% 回归误差平方和Q
Q=w*w';
% 计算回归平方和U
U=S-Q;
% 计算，并输出F检验值
F=28*U/Q
% 计算非线性回归模型的拟合值
for j=1:30,
Y1(j)=1/(B(1,1)+B(2,1)*exp(-j));
end
% 输出非线性回归模型的拟合曲线（Logisic曲线）
plot(T,Y)
%figure
%h1=plot(T,Y,'k*');
%set(gca,'linewidth',2);
%set(h1,'MarkerSize',8);
%xlabel('时间/年', 'fontsize',12);
%ylabel('人口/人', 'fontsize',12);
%figure
%h2=plot( T, Y1,'--k*',...
% 'LineWidth',2)
%set(gca,'linewidth',2);
%set(h2,'MarkerSize',8);
%xlabel('时间/年', 'fontsize',12);
%ylabel('人口/人', 'fontsize',12);
薄膜渗透率测定：http://note.youdao.com/noteshare?id=626fc2468f3f4b14346da42eb70d2370&sub=7250B3DD2E844C06BC3235456C9D8F76
%curvefun.m
function f = curvefun(x,tdata)
f = x(1)+x(2)*exp(-0.02*x(3)*tdata); %其中 x(1) = a;x(2) = b;x(3) = k;
%test1.m
tdata = linspace(100,1000,10);
cdata = 1e-05.*[454 499 535 565 590 610 626 639 650 659];
x0 = [0.2,0.05,0.05];
opts = optimset( 'lsqcurvefit' );
opts = optimset( opts, 'PrecondBandWidth', 0 )
x=lsqcurvefit ('curvefun',x0,tdata,cdata,[],[],opts)
f= curvefun(x,tdata)
plot(tdata,cdata,'o',tdata,f,'r-')
xlabel('时间(秒)')
ylabel('浓度（毫克/立方厘米）')
地形地貌图形绘制：http://note.youdao.com/noteshare?id=5b6cfedcca8915dfbd82c2811eb493f3&sub=DAC9341398874C0E880F69A42A61D024
2002A-车灯光源投影区域的绘制
%Fig1_17.m
p=0.03;x=25.0216;
for y1=-0.002:0.0004:0.002
y0=(-0.036:0.001:0.036)'*ones(1,73);
z0=ones(73,1)*(-0.036:0.001:0.036);
x0=(y0.^2+z0.^2)/(2*p);
xn=(p^3+4*x0*2*p.*x0+p*(-4*y1*y0+3*2*p*x0))./(2*(p^2+2*p*x0));
yn=(2*p*x0.*y0+p^2*(-y1+y0)+y1*(y0.^2-z0.^2))./(p^2+2*p*x0);
zn=(p^2+2*p*x0+2*y1*y0).*z0./(p^2+2*p*x0);
y=y0+(yn-y0).*(x-x0)./(xn-x0);
z=z0+(zn-z0).*(x-x0)./(xn-x0);
plot(y,z,'k.')
hold on
end
xlabel('y/m', 'fontsize',12)
ylabel('z/m', 'fontsize',12)
matlab将评判矩阵转换为因素的权重矩阵
%AHP.m
%% AHP法权重计算MATLAB程序
%% 数据读入
clc
clear all
A=[1 2 6; 1/2 1 4; 1/6 1/4 1];% 评判矩阵
%% 一致性检验和权向量计算
[n,n]=size(A);
[v,d]=eig(A);
r=d(1,1);
CI=(r-n)/(n-1);
RI=[0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.52 1.54 1.56 1.58 1.59];
CR=CI/RI(n);
if CR<0.10
CR_Result='通过';
else
CR_Result='不通过';
end
%% 权向量计算
w=v(:,1)/sum(v(:,1));
w=w';
%% 结果输出
disp('该判断矩阵权向量计算报告：');
disp(['一致性指标:' num2str(CI)]);
disp(['一致性比例:' num2str(CR)]);
disp(['一致性检验结果:' CR_Result]);
disp(['特征值:' num2str(r)]);
disp(['权向量:' num2str(w)]);
%线性规划
[x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
%c行向量--目标函数系数
%A矩阵--不等式约束左端
%b列向量--不等式约束右端
%Aeq矩阵--等式约束左端
%beq列向量--等式约束右端
%LB--变量下界 UB--变量上界
%fval--在向量x处的目标函数值
%x--向量x
%非线性约束
X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)
%FUN--函数f(x)
%x无上界和下界LB=[],UB=[]
%x无下界LB=-inf;x无上界UB=inf
%NONLCON定义非线性约束条件中的函数
%二次规划
[X,FVAL]=quadprog(h,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
%h--目标函数二次系数
%f--目标函数一次系数
%SUMT
[x,y]=fminunc('test',rand(2,1))
%'test'中写外罚函数 -- rand(2,1)中的2是行向量x的维数
文档：非线性规划.note
链接：http://note.youdao.com/noteshare?id=a797a12637b8b216c0ea2e9bc5c12bcf&sub=231E65D41B83495497B46ABDBACD8C17
文档：线性规划.note
链接：http://note.youdao.com/noteshare?id=f9515e3bd2525b927a3df0ebbaf9b97d&sub=05CFB6B4B67F4D29BE53EAF8AD7BCAEC
非线性规划的目标函数为二次函数且约束条件都是线性的，称该规划为二次规划。
文档：二次规划.note
链接：http://note.youdao.com/noteshare?id=67353867ef4fe04be3495a7f435b695a&sub=35AC74753987483DB06F8EE28B3B0EC7
利用罚函数法，可将非线性规划问题转化为求解一系列无约束极值问题，该方法称为序列无约束最小化技术(SUMT)
文档：SUMT.note
链接：http://note.youdao.com/noteshare?id=8efd48376b19e26faa3720d51d204187&sub=DCAA044B23494299B926AC213EDD9320
规划中的变量(部分或者全部)限制为整数时，称为整数规划。
常见的整数规划问题求解算法：
文档：隐枚举法在0-1整数规划中的应用.note
链接：http://note.youdao.com/noteshare?id=294c77eba573d26d0c15c437f43dc2c6&sub=4CD1813C310C43FF995C3F91DFE4CE60
文档：随机取样计算法在非线性规划中的应用.n…
链接：http://note.youdao.com/noteshare?id=e329f09e3052cdb513a10b04555980ae&sub=F7765789BFE94C5984D7E62F7827259E
云模型系我国李德毅院士首创，是一门相对比较新型的学科，属于不确定性人工智能范畴，主要用于定性与定量之间的相互转换。自然界中的不确定性从属性角度来说主要有随机性和模糊性，这跟单色光的“波粒二象性”有点类似。
“云”或者“云滴”(cloud)是云模型的基本单元。所谓“云”是指其在论域U上的一个分布，可以用联合概率的形式(x,μ)类比。
云模型用三个数据来表示其特征：
期望：云滴在论域空间分布的期望，一般用符号Ex表示。
熵：不确定性程度，由离群程度和模糊程度共同决定，一般用符号En表示。
超熵：用来度量熵的不确定性，亦即熵的熵，一般用符号He表示。
云有两种发生器：正向云发生器和逆向云发生器，分别用来生成足够的云滴和计算云滴的云数字特征( Ex，En，He)。
正向云发生器的触发机制：
①生成以En为期望，以He：为方差的正态随机数En‘。
②生成以Ex为期望，以En’’为方差的正态随机数x。
③计算隶属度也就是确定度
μ = e − ( x − E x ) 2 2 E n ′ 2 μ=e^{-\frac{(x-Ex)^2}{2En'^2}} μ=e−2En′2(x−Ex)2
那么(x,μ)便是相对于论域u的一个云滴。这里选择常用的“钟形”函数"为隶属度函数。
μ = e − ( x − a ) 2 2 b 2 μ=e^{-\frac{(x-a)^2}{2b^2}} μ=e−2b2(x−a)2
④重复①②③步骤n至生成足够的云滴。
逆向云发生器触发机制
计算样本均值X，方差S^2，得到
E x = X ; E n = π 2 × 1 n ∑ 1 n ∣ x − E x ∣ ; H e = S 2 − E n 2 Ex=X;En=\sqrt{\frac{π}{2}}×\frac{1}{n}\sum^{n}_{1}|x-Ex|;He=\sqrt{S^2-En^2} Ex=X;En=2π×n11∑n∣x−Ex∣;He=S2−En2
应用：射击比赛
%cloud_main.m
%Input：10次射击ABCD的成绩
%Output：发挥出色的选手是第三位C
%analyzation：从期望、倾向的稳定程度(确定度)、离散程度。
% 以下是主函数cloud_main.m
clc;
clear all;
close all;
% 每幅图生成N个云滴
N = 1500;
% 射击成绩原始数据，这里数据按列存储所以要转置
Y = [9.5 10.3 10.1 8.1
10.3 9.7 10.4 10.1
10.6 8.6 9.2 10.0
10.5 10.4 10.1 10.1
10.9 9.8 10.0 10.1
10.6 9.8 9.7 10.0
10.4 10.5 10.6 10.3
10.1 10.2 10.8 8.4
9.3 10.2 9.6 10.0
10.5 10.0 10.7 9.9]';
for i = 1:size(Y,1)
subplot(size(Y,1)/2,2,i)
% 调用函数
[x,y,Ex,En,He] = cloud_transform(Y(i,:),N);
plot(x,y,'r.');
xlabel('射击成绩分布/环');
ylabel('确定度');
title(strcat('第',num2str(i),'人射击云模型还原图谱'));
% 控制坐标轴的范围
% 统一坐标轴范围才会在云模型形态上具有可比性
axis([8,12,0,1]);
end
%cloud_transform.m
function [x,y,Ex,En,He] = cloud_transform(y_spor,n)
% x 表示云滴
% y 表示隶属度（这里是“钟形”隶属度），意义是度量倾向的稳定程度
% Ex 云模型的数字特征，表示期望
% En 云模型的数字特征，表示熵
% He 云模型的数字特征，表示超熵
% 通过统计数据样本计算云模型的数字特征
Ex = mean(y_spor);
En = mean(abs(y_spor-Ex)).*sqrt(pi./2);
He = sqrt(var(y_spor)-En.^2);
% 通过云模型的数字特征还原更多的“云滴”
for q = 1:n
Enn = randn(1).*He + En ;
x(q) = randn(1).*Enn + Ex ;
y(q) = exp(-(x(q) - Ex).^2./(2 .* Enn.^2));
end
x;
y;
在回归分析中，因变量y可能有两种情形：
①y是一个定量的变量，这时就用通常的regress函数进行回归；
②y是一个定性的变量，比如。y=0或1，这时就不能用通常的
regress函数对y进行回归，而是使用所谓的Logistic回归。
Logistic方法主要应用于研究某些现象发生的概率P，比如股票涨还是跌，公司成功或失败的概率。除此之外，本章还讨论概率P与哪些因素有关。
文档：logistic模型.note
链接：http://note.youdao.com/noteshare?id=3916a6294742bf6df41aa182c25c99d6&sub=970DFA001C9A48B2B764B95049D32803
应用：还款能力评估
文档：还款能力评估题目.note
链接：http://note.youdao.com/noteshare?id=c7f59813f4ef0722b05a3a19494c6a8e&sub=51FB19A546574A4FA42790EF6AE3EA69
%已知20个公司的输入X1 X2 X3 输出Y
%评估5个公司是否能还款（已知他们的X1 X2 X3），根据
%数据链接：https://pan.baidu.com/s/1xQd0JhQu-C1MWy8X5YqEtg
%提取码：p35r
%% 第2章 数学建模方法及Matlab实现
% 程序2-1： logistic方法Matlab实现程序
%--------------------------------------------------------------------------
%% 数据准备
clear all
clc
X0=xlsread('Ch3_logistic_ex1.xlsx', 'A2:C21'); % 回归数据X值
XE=xlsread('Ch3_logistic_ex1.xlsx', 'A2:C26'); % 验证与预测数据
Y0=xlsread('Ch3_logistic_ex1.xlsx', 'D2:D21'); % 回归数据P值
%--------------------------------------------------------------------------
%% 数据转化和参数回归
n=size(Y0,1);
for i=1:n
if Y0(i)==0
Y1(i,1)=0.25;
else
Y1(i,1)=0.75;
end
end
X1=ones(size(X0,1),1); % 构建常数项系数
X=[X1, X0];
Y=log(Y1./(1-Y1));
b=regress(Y,X);
%--------------------------------------------------------------------------
%% 模型验证和应用
for i=1:size(XE,1)
Pai0=exp(b(1)+b(2)*XE(i,1)+b(3)*XE(i,2)+b(4)*XE(i,3))/(1+exp(b(1)+b(2)*XE(i,1)+b(3)*XE(i,2)+b(4)*XE(i,3)));
if Pai0<=0.5
P(i)=0;
else
P(i)=1;
end
end
%--------------------------------------------------------------------------
%% 显示结果
disp(['回归系数:' num2str(b')]);
disp(['评价结果:' num2str(P)]);
背景：当变量个数多且关系复杂时，增加了问题的复杂性。
解决：由于多个变量之间存在相关性，将多个变量综合为少数几个有代表性变量，使他们代表变量大多数信息又互不相关。——PCA——是一种数学降维方法
PCA：
文档：PCA步骤.note
链接：http://note.youdao.com/noteshare?id=fa3ff82e9f1697da892373886e627632&sub=0209224CECC340F6904B8EEFC24EF676
应用：企业综合实力排序
%链接：https://pan.baidu.com/s/198GcoR4dNJUp3vb9RRHHsA
%提取码：2se9
%Input:8项指标
%Output:实力排序9最强，12最弱
%== PCA stepping demonstration program==%
% Read data from a file (e.g. excel) and place it in a matrix.
A=xlsread('Coporation_evaluation.xlsx', 'B2:I16');
% Transfer orginal data to standard data
a=size(A,1); % Get the row number of A
b=size(A,2); % Get the column number of A
for i=1:b
SA(:,i)=(A(:,i)-mean(A(:,i)))/std(A(:,i)); % Matrix normalization
end
% Calculate correlation matrix of A.
CM=corrcoef(SA);
% Calculate eigenvectors and eigenvalues of correlation matrix.
[V, D]=eig(CM);
% Get the eigenvalue sequence according to descending and the corrosponding
% attribution rates and accumulation rates.
for j=1:b
DS(j,1)=D(b+1-j, b+1-j);%特征值降序排序
end
for i=1:b
DS(i,2)=DS(i,1)/sum(DS(:,1));%贡献率(各主成分的权重信息)
DS(i,3)=sum(DS(1:i,1))/sum(DS(:,1));%累计贡献率
end
% Calculate the numvber of principal components.选择主成分及其对应特征向量
T=0.9; % set the threshold value for evaluating information preservation level.主成分信息保留率
for K=1:b
if DS(K,3)>=T
Com_num=K;
break;
end
end
% Get the eigenvectors of the Com_num principal components提取主成分的特征向量(与原始变量的关联关系)
for j=1:Com_num
PV(:,j)=V(:,b+1-j);
end
% Calculate the new socres of the orginal items计算评价对象的主成分得分
new_score=SA*PV;
for i=1:a
total_score(i,2)=sum(new_score(i,:));
total_score(i,1)=i;
end
new_score_s=sortrows(total_score,-2);
% Displays result reports
disp('特征值及贡献率：')
DS
disp('阀值T对应的主成分数与特征向量：')
Com_num
PV
disp('主成分分数：')
new_score
disp('主成分分数排序：')
new_score_s
SVM就是与分类器平行的两个平面，此两个平面能很好的分开两类不同的数据，且穿越两类数据区域集中的点，现在想要寻找最佳超几何分隔平面使之与两个平面间的距离最大，如此便能实现分类总误差最小。
落在边界平面上的(数据集中)的点称为支持向量，这些点至关重用，故称为支撑向量机。
文档：SVM理论基础.note
链接：http://note.youdao.com/noteshare?id=c00525bf6b12f14e0035d6ff847b807b&sub=28D134E5022D4947855BF27B8E9A7BD9
应用：决定是否批准贷款(转化为线性问题来做)
文档：SVM使用案例.note
链接：http://note.youdao.com/noteshare?id=d3bc05599c9733a793b6b247fdb14107&sub=63D8FF72CC0F48E8B62055BD4EE8051A
%链接：https://pan.baidu.com/s/1WDRg_c9QhNsndwTLz8EbAg
%提取码：6xgx
%% 线性规划法求解SVM模型
%% 清空环境变量、读入数据
clc
clear all
X0=xlsread('SVM_ex2.xlsx', 'B2:E19');
for i=1:3
X(:,i)=(X0(:,i)-mean(X0(:,i)))/std(X0(:,i)); % 数据标准化
end
% 数据预处理
[m,n]=size(X);
e=ones(m,1);
D=[X0(:,4)];
B=zeros(m,m);
C=zeros(m,m);
for i=1:m
B(i,i)=1;
C(i,i)=D(i,1);
end
%% 转化成规划模型进行求解
A=[-X(:,1).*D, -X(:,2).*D, -X(:,3).*D, D, -B];
b=-e;
f=[0,0,0,0, ones(1,m)];
lb=[-inf,-inf,-inf,-inf,zeros(1,m)]';
x = linprog(f,A,b,[],[],lb);
%% 模型验证与结果显示
W=[x(1,1), x(2,1), x(3,1)]; % 提取系数
CC=x(4,1); % 提取截距
X1=[X(:,1), X(:,2), X(:,3)];
R1=X1*W'-CC;
R2=sign(R1); %进行分类
disp('程序输出结果：');
disp('超平面方程为：');
disp(['X1:' num2str(x(1,1))]);
disp(['X2:' num2str(x(2,1))]);
disp(['X3:' num2str(x(3,1))]);
disp(['intercept:' num2str(x(4,1))]); % 常数项
disp('超平面分类结果：');
R=[R1, R2]
聚类分析是研究(样品、指标)分类的一种统计分析方法。聚类分析计算方法有
文档：k-means思想.note
链接：http://note.youdao.com/noteshare?id=b00c16993977390c8d5c2c9be27590b9&sub=BFF69556770F4C04B28F5D30126AE7DE
应用：已知20个样本，每个样本有两个特征，要求对这些数据分类
%% 第3章 数学建模方法及MATLAB实现
% K-means 算法MATLAB实现
%--------------------------------------------------------------------------
%% 数据准备和初始化
clc
clear
x=[0 0;1 0; 0 1; 1 1;2 1;1 2; 2 2;3 2; 6 6; 7 6; 8 6; 6 7; 7 7; 8 7; 9 7 ; 7 8; 8 8; 9 8; 8 9 ; 9 9];
z=zeros(2,2);
z1=zeros(2,2);
z=x(1:2, 1:2);
%% 寻找聚类中心
while 1
count=zeros(2,1);
allsum=zeros(2,2);
for i=1:20 % 对每一个样本i，计算到2个聚类中心的距离
temp1=sqrt((z(1,1)-x(i,1)).^2+(z(1,2)-x(i,2)).^2);
temp2=sqrt((z(2,1)-x(i,1)).^2+(z(2,2)-x(i,2)).^2);
if(temp1<temp2)
count(1)=count(1)+1;
allsum(1,1)=allsum(1,1)+x(i,1);
allsum(1,2)=allsum(1,2)+x(i,2);
else
count(2)=count(2)+1;
allsum(2,1)=allsum(2,1)+x(i,1);
allsum(2,2)=allsum(2,2)+x(i,2);
end
end
z1(1,1)=allsum(1,1)/count(1);
z1(1,2)=allsum(1,2)/count(1);
z1(2,1)=allsum(2,1)/count(2);
z1(2,2)=allsum(2,2)/count(2);
if(z==z1)
break;
else
z=z1;
end
end
%% 结果显示
disp(z1);% 输出聚类中心
plot( x(:,1), x(:,2),'b*')
hold on
plot(z1(:,1),z1(:,2),'ro')
title('K-means分类图');
xlabel('特征x1');
ylabel('特征x2');
贝叶斯判别法
文档：朴素贝叶斯使用.note
链接：http://note.youdao.com/noteshare?id=55027131b3743c9b833c0d9f389d6edc&sub=8F5C2D11D14F45EB8A494804CBA95506
%matlab为朴素贝叶斯提供了一个类
ObjBayes=NaiveBayes.fit(training_sample,class);
%training_sample为训练样本变量值，每一行为一个观测值，每一列代表一个变量
%class为分类结果
pred=ObjBayes.predict(test);
%test表示待分类样本
文档：数据建模文献.note
链接：http://note.youdao.com/noteshare?id=4d1c025441c629606db5c1bc8a9055f7&sub=6AACBA5449A7482ABAD642883B7C4A74
Gray Model,GM的形式是在耦合指数e的基础上发展起来的，同样遵循自然界的能量系统。
自然界中的三种不确定性有：
对于不确定的问题，适合用灰色系统的知识来求解。
文档：灰色系统基本理论.note
链接：http://note.youdao.com/noteshare?id=9e4089a881247832e9f6c8790707296b&sub=E68528CE81094FCCAE9B1D08CD0F27D6
灰色关联度矩阵应用：分析投资和收入之间定量关系并给出建议性意见
%5项投资是待比较数列即子因素
%6项收入为参考数列即母因素
%此代码中教会了我如何改在X轴的刻度上写汉字
%结果分析过程：
%文档：结果分析过程.note
%链接：http://note.youdao.com/noteshare?id=da34e76b4e24df497b670168bc5977b7&sub=F40AFFB9935F48568429E956A959B729
clc;
close;
clear all;
% 控制输出结果精度
format short;
% 原始数据
x=[308.58 310 295 346 367
195.4 189.9 187.2 205 222.7
24.6 21 12.2 15.1 14.57
20 25.6 23.3 29.2 30
18.98 19 22.3 23.5 27.655
170 174 197 216.4 235.8
57.55 70.74 76.8 80.7 89.85
88.56 70 85.38 99.83 103.4
11.19 13.28 16.82 18.9 22.8
4.03 4.26 4.34 5.06 5.78
13.7 15.6 13.77 11.98 13.95];
n1=size(x,1);
% 数据标准化处理
for i = 1:n1
x(i,:) = x(i,:)/x(i,1);
end
% 保存中间变量，亦可省略此步，将原始数据赋予变量data
data=x;
% 分离参考数列（母因素）
consult=data(6:n1,:);
m1=size(consult,1);
% 分离比较数列（子因素）
compare=data(1:5,:);
m2=size(compare,1);
for i=1:m1
for j=1:m2
t(j,:)=compare(j,:)-consult(i,:);
end
min_min=min(min(abs(t')));
max_max=max(max(abs(t')));
% 通常分辨率都是取0.5
resolution=0.5;
% 计算关联系数
coefficient=(min_min+resolution*max_max)./(abs(t)+resolution*max_max);
% 计算关联度
corr_degree=sum(coefficient')/size(coefficient,2);
r(i,:)=corr_degree;
end
% 输出关联度值并绘制柱形图
r
bar(r,0.90);
axis tight;
% 从左至右依次是 固定资产投资,工业投资,农业投资,科技投资,交通投资
legend('固定资产投资','工业投资','农业投资','科技投资','交通投资');
% 以下程序是为了给x轴添加汉字标签
% 其基本原理是先去掉x轴上的固有标签，然后用文本标注x轴
% 去掉X轴上默认的标签
set(gca,'XTickLabel','');
% 设定X轴刻度的位置，这里有6个母因素
n=6;
% 这里注意：x_range范围如果是[1 n]会导致部门柱形条不能显示出来
% 所以范围要缩一点
x_value = 1:1:n;
x_range = [0.6 n+.4];
% 获取当前图形的句柄
set(gca,'XTick',x_value,'XLim',x_range);
% 在X轴上标记6个母因素
profits={'国民收入','工业收入','农业收入','商业收入','交通收入','建筑业收入'};
y_range = ylim;
% 用文本标注母因素名称
handle_date = text(x_value,y_range(1)*ones(1,n)+.018,profits(1:1:n));
% 把文本的字体设置合适的格式和大小并旋转一定的角度
set(handle_date,'HorizontalAlignment','right','VerticalAlignment','top', 'Rotation',35);
% y轴标记
ylabel('y');
title('投资对收入的作用');
GM(1,1)应用：客户根据自己资金进行投资，并且有10天犹豫期，期间可以自由退，并且给出了每一天继续持有客户的比例
%output：判断模型精度都是一级的，精度非常高(和小样本有关)
% 在程序的所有节骨眼都进行了详尽的注释
clc;
clear all;
% 以下数据背景：银行各种理财金融衍生品，客户在购买某种理财产品的时候
% 一般会有10天的“犹豫期”，在这10天可以自由退买
% 我们分析客户购买行为数据发现，在当天（姑且称之为第0天）购买之后
% 在第1天没有发生退买（保有量）的比率大概是92.81%，在第2天没有退买的比率大概是97.66%
% 以此类推
% 这里取前8天的数据：
x0 = [92.810 97.660 98.800 99.281 99.537 99.537 99.817 100.000];
n = length(x0);
% 做级比判断，看看是否适合用GM(1，1)建模
lamda = x0(1:n-1)./x0(2:n);
range = minmax(lamda);
% 判定是否适合用一阶灰色模型建模
if range(1,1) < exp(-(2/(n+2))) | range(1,2) > exp(2/(n+2))
error('级比没有落入灰色模型的范围内');
else
% 空行输出
disp(' ');
disp('可用G(1，1)建模');
end
% 做AGO累加处理
x1 = cumsum(x0);
for i = 2:n
% 计算紧邻均值，也就是白化背景值
z(i) = 0.5*(x1(i)+x1(i-1));
end
B = [-z(2:n)',ones(n-1,1)];
Y = x0(2:n)';
% 矩阵做除法，计算发展系数和灰色作用量
% 千万注意：这里是右除，不是左除
u = B\Y;
% 在MATLAB中，用大写做字母D表示导数,dsolve函数用来求解符号常微分方程
x = dsolve('Dx+a*x=b','x(0)=x0');
% subs函数的作用是替换元素，这里是把常量a,b,x0换成具体u(1),u(2),x1(1)数值
x = subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});
forecast1 = subs(x,'t',[0:n-1]);
% digits和vpa函数用来控制计算的有效数字位数
digits(6);
% y值是AGO形式的（还是累加的）
y = vpa(x);
% 把AGO输出值进行累减
% diff用于对符号表达式进行求导
% 但是如果输入的是向量，则表示计算原向量相邻元素的差
forecast11 = double(forecast1);
exchange = diff(forecast11);
% 输出灰色模型预测的值
forecast = [x0(1),exchange]
% 计算残差
epsilon = x0 - forecast;
% 计算相对误差
delta = abs(epsilon./x0);
% 检验模型的误差
% 检验方法一：相对误差Q检验法
Q = mean(delta)
% 检验方法二：方差比C检验法
% 计算标准差函数为std（x，a）
% 如果后面一个参数a取0表示的是除以n－1，如果是1就是最后除以n
C = std(epsilon,1)/std(x0,1)
% 检验方法三：小误差概率P检验法
S1 = std(x0,1);
S1_new = S1*0.6745;
temp_P = find(abs(epsilon-mean(epsilon)) < S1_new);
P = length(temp_P)/n
% 绘制原始数列与灰色模型预测得出的数列差异折线图
plot(1:n,x0,'ro','markersize',11);
hold on
plot(1:n,forecast,'k-','linewidth',2.5);
grid on;
axis tight;
xlabel('x');
ylabel('y');
title('保有量比例与时间序列的关系');
legend('原始数列','模型数列');
x ~ ( 1 ) ( k + 1 ) = ( x ( 0 ) ( 1 ) − b a ) e − a b + b a \tilde{x}^{(1)}(k+1)=(x^{(0)}(1)-\frac{b}{a})e^{-ab}+\frac{b}{a} x~(1)(k+1)=(x(0)(1)−ab)e−ab+ab
% 在程序的所有节骨眼都进行了详尽的注释
clc;
clear all;
% 以下数据背景：银行各种理财金融衍生品，客户在购买某种理财产品的时候
% 一般会有10天的“犹豫期”，在这10天可以自由退买
% 我们分析客户购买行为数据发现，在当天（姑且称之为第0天）购买之后
% 在第1天没有发生退买（保有量）的比率大概是92.81%，在第2天没有退买的比率大概是97.66%
% 以此类推
% 这里取前8天的数据：
x0 = [92.810 97.660 98.800 99.281 99.537 99.537 99.817 100.000];
n = length(x0);
% 做级比判断，看看是否适合用GM(1，1)建模
lamda = x0(1:n-1)./x0(2:n);
range = minmax(lamda);
% 判定是否适合用一阶灰色模型建模
if range(1,1) < exp(-(2/(n+2))) | range(1,2) > exp(2/(n+2))
error('级比没有落入灰色模型的范围内');
else
% 空行输出
disp(' ');
disp('可用G(1，1)建模');
end
%----------以下程序是修改的地方-------%
% 做AGO累加处理
x1 = cumsum(x0);
for i = 2:n
% 计算紧邻均值，也就是白化背景值
z(i) = 0.5*(x1(i)+x1(i-1));
end
B = [-z(2:n)',ones(n-1,1)];
Y = x0(2:n)';
% 用最小二乘法计算发展系数和灰色作用量
u = inv(B'*B)*B'*Y;
% 利用GM(1，1)具体表达式计算原始数列的AGO
forecast1 = (x1(1)-u(2)./u(1)).*exp(-u(1).*([0:n-1]))+u(2)./u(1) ;
% 把AGO输出值进行累减
% diff用于对符号表达式进行求导
% 但是如果输入的是向量，则表示计算原向量相邻元素的差
forecast11 = double(forecast1);
exchange = diff(forecast11);
%----------以上程序是修改的地方-------%
% 输出灰色模型预测的值
forecast = [x0(1),exchange]
% 计算残差
epsilon = x0 - forecast;
% 计算相对误差
delta = abs(epsilon./x0);
% 检验模型的误差
% 检验方法一：相对误差Q检验法
Q = mean(delta)
% 检验方法二：方差比C检验法
% 计算标准差函数为std（x，a）
% 如果后面一个参数a取0表示的是除以n－1，如果是1就是最后除以n
C = std(epsilon,1)/std(x0,1)
% 检验方法三：小误差概率P检验法
S1 = std(x0,1);
S1_new = S1*0.6745;
temp_P = find(abs(epsilon-mean(epsilon)) < S1_new);
P = length(temp_P)/n
% 绘制原始数列与灰色模型预测得出的数列差异折线图
plot(1:n,x0,'ro','markersize',11);
hold on
plot(1:n,forecast,'k-','linewidth',2.5);
grid on;
axis tight;
xlabel('x');
ylabel('y');
title('保有量比例与时间序列的关系');
legend('原始数列','模型数列');
灰色verhulst模型应用：根据光密度值可得知细菌数量，利用该模型预测细菌数量
clc, clear all, close all
% 大肠杆菌测定的光密度值
x1=[0.025 0.023 0.029 0.044 0.084 0.164 0.332 0.521 0.97 1.6 2.45 3.11 3.57 3.76 3.96 4 4.46 4.4 4.49 4.76 5.01];
n=length(x1);
year=0:n-1;
figure(1);
plot(year,x1,'k*');
x0=diff(x1);
x0=[x1(1),x0];
for i=2:n
z1(i)=0.5*(x1(i)+x1(i-1));
end
z1;
B=[-z1(2:end)',z1(2:end)'.^2];
Y=x0(2:end)';
% 当然也可以使用最小二乘法和verhulst灰色模型的具体表达式来求解
% 前面章节已经有很详细的讲解了
abvalue=B\Y;
x=dsolve('Dx+a*x=b*x^2','x(0)=x0');
x=subs(x,{'a','b','x0'},{abvalue(1),abvalue(2),x1(1)});
forecast=subs(x,'t',0:n-1);
digits(6);x=vpa(x);
forecast
% 绘图
hold on;
plot(year,forecast,'k-.','linewidth',4);
xlabel('时间均匀采样/5小时');
ylabel('细菌培养液吸光度/OD600');
legend('实际数量','预测数量');
title('大肠杆菌培养S形增长曲线');
axis tight;
灰色预测步骤：
灰色预测应用：已知公司利润，预测未来几年利润情况
clear
syms a b;
c=[a b]';
A=[89677,99215,109655,120333,135823,159878,182321,209407,246619,300670];
B=cumsum(A); % 原始数据累加
n=length(A);
for i=1:(n-1)
C(i)=(B(i)+B(i+1))/2; % 生成累加矩阵
end
% 计算待定参数的值
D=A;D(1)=[];
D=D';
E=[-C;ones(1,n-1)];
c=inv(E*E')*E*D;
c=c';
a=c(1);b=c(2);
% 预测后续数据
F=[];F(1)=A(1);
for i=2:(n+10)
F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a ;
end
G=[];G(1)=A(1);
for i=2:(n+10)
G(i)=F(i)-F(i-1); %得到预测出来的数据
end
t1=1999:2008;
t2=1999:2018;
G
plot(t1,A,'o',t2,G)%原始数据和预测数据比较
xlabel('年份')%'fontsize',12
ylabel('利润/(元/年)')
%set(gca, 'LineWidth',2);
灰色预测应用：2005A长江水质的预测(数据样本少，污水排放量是不确定的系统)
clear all
syms a b;
c=[a b]';
A=[174 179 183 189 207 234 220.5 256 270 285];
B=cumsum(A); % 原始数据累加
n=length(A);
for i=1:(n-1)
C(i)=(B(i)+B(i+1))/2; % 生成累加矩阵
end
% 计算待定参数的值
D=A;D(1)=[];
D=D';
E=[-C;ones(1,n-1)];
c=inv(E*E')*E*D;
c=c';
a=c(1);b=c(2);
% 预测后续数据
F=[];F(1)=A(1);
for i=2:(n+10)
F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a ;
end
G=[];G(1)=A(1);
for i=2:(n+10)
G(i)=F(i)-F(i-1); %得到预测出来的数据
end
t1=1995:2004;
t2=1995:2014;
G, a, b % 输出预测值，发展系数和灰色作用量
plot(t1,A,'ko')
hold on
plot(t2,G,'k')%, 'LineWidth',2
xlabel('年份')%, 'fontsize',12
ylabel('污水量/亿吨')
%set(gca, 'LineWidth',2);
灰色预测应用：2009D预测与会代表人数
文档：2009D.note
链接：http://note.youdao.com/noteshare?id=89c38a169e26b59fee919655045f432b&sub=BBA2980695224C8ABB9D1FCAE56311E2
小结：学习了灰色关联度矩阵、GM(1,1)、灰色verhulst模型
文档：灰色预测文献.note
链接：http://note.youdao.com/noteshare?id=930fa9e3a9d9608e3cedd3d7c7755f66&sub=5CC8B509727D4C739E17ABCDC4CD6285
数据、信息、知识之间的关系：
智能算法就是通过数据挖掘工具从数据中提取有益的信息和知识。有粒子群算法和蚁群算法、遗传算法和免疫算法、神经网络和模拟退火算法，本质上都是机器学习，其具有明显的特征——机械性。
数学机械化是数学在运算和证明过程中，只要前进一步，都会有确定、必然的下一步，直到到达绝伦，整个过程按照既定的刻板规律进行。
文档：SGA基本步骤.note
链接：http://note.youdao.com/noteshare?id=e30e715b3c32acd9b615cf9ca87227da&sub=7349E321581243CFB8759023596FAB77
%1．遗传算法的程序设计伪代码
BEGIN
t=0； %遗传代数
初始化P(t)； %初始化种群或染色体
计算P(t)的适应值；
while（不满足停止准则）do
begin
t = t+l ：
从P(t -1)中选择P(t)； %选择
重组P(t)； %交叉和变异
计算P(t)的适应值；
end
END
文档：GA工具箱使用.note
链接：http://note.youdao.com/noteshare?id=330963492d9e0d395bc5c93d4c90c3dd&sub=4BA3F4D4D3244591BC8DAF6AAE9E2492
遗传算法应用(源程序，不使用工具箱)：无约束目标函数最大值
m a x f ( x ) = 200 e − 0.05 x s i n x , x ∈ [ − 2 , 2 ] maxf(x)=200e^{-0.05x}sinx,x∈[-2,2] maxf(x)=200e−0.05xsinx,x∈[−2,2]
%结论1：最大适应度值在进化过程中变化幅度不大：因为[-2,2]区间较小，种群为50条染色体较多，让算法最大初始适应度更有机会直接落在理论最优解附近。
%结论2：平均适应度值在进化过程中快速增加：说明遗传算法是进化方向非常明确的机器学习算法。
%如果不用遗传算法，而是将染色体作为蒙特卡洛模拟寻优的随机数产生器，多次产生随机数并保存历史最优值，这种寻优方向是不定向的，与遗传算法有本质区别。
%主程序：用遗传算法求解y=200*exp(-0.05*x).*sin(x)在[-2 2]区间上的最大值
clc;
clear all;
close all;
global BitLength
global boundsbegin
global boundsend
bounds=[-2 2];%一维自变量的取值范围
precision=0.0001; %运算精度
boundsbegin=bounds(:,1);
boundsend=bounds(:,2);
%计算如果满足求解精度至少需要多长的染色体
BitLength=ceil(log2((boundsend-boundsbegin)' ./ precision));
popsize=50; %初始种群大小
Generationnmax=12; %最大代数
pcrossover=0.90; %交配概率
pmutation=0.09; %变异概率
%产生初始种群
population=round(rand(popsize,BitLength));
%计算适应度,返回适应度Fitvalue和累积概率cumsump
[Fitvalue,cumsump]=fitnessfun(population);
Generation=1;
while Generation<Generationnmax+1
for j=1:2:popsize
%选择操作
seln=selection(population,cumsump);
%交叉操作
scro=crossover(population,seln,pcrossover);
scnew(j,:)=scro(1,:);
scnew(j+1,:)=scro(2,:);
%变异操作
smnew(j,:)=mutation(scnew(j,:),pmutation);
smnew(j+1,:)=mutation(scnew(j+1,:),pmutation);
end
population=smnew; %产生了新的种群
%计算新种群的适应度
[Fitvalue,cumsump]=fitnessfun(population);
%记录当前代最好的适应度和平均适应度
[fmax,nmax]=max(Fitvalue);
fmean=mean(Fitvalue);
ymax(Generation)=fmax;
ymean(Generation)=fmean;
%记录当前代的最佳染色体个体
x=transform2to10(population(nmax,:));
%自变量取值范围是[-2 2],需要把经过遗传运算的最佳染色体整合到[-2 2]区间
xx=boundsbegin+x*(boundsend-boundsbegin)/(power((boundsend),BitLength)-1);
xmax(Generation)=xx;
Generation=Generation+1
end
Generation=Generation-1;
Bestpopulation=xx
Besttargetfunvalue=targetfun(xx)
%绘制经过遗传运算后的适应度曲线。一般地，如果进化过程中种群的平均适应度与最大适
%应度在曲线上有相互趋同的形态，表示算法收敛进行得很顺利，没有出现震荡；在这种前
%提下，最大适应度个体连续若干代都没有发生进化表明种群已经成熟。
figure(1);
hand1=plot(1:Generation,ymax);
set(hand1,'linestyle','-','linewidth',1.8,'marker','*','markersize',6)
hold on;
hand2=plot(1:Generation,ymean);
set(hand2,'color','r','linestyle','-','linewidth',1.8,...
'marker','h','markersize',6)
xlabel('进化代数');ylabel('最大/平均适应度');xlim([1 Generationnmax]);
legend('最大适应度','平均适应度');
box off;hold off;
%子程序：新种群交叉操作,函数名称存储为crossover.m
function scro=crossover(population,seln,pc);
BitLength=size(population,2);
pcc=IfCroIfMut(pc); %根据交叉概率决定是否进行交叉操作，1则是，0则否
if pcc==1
chb=round(rand*(BitLength-2))+1; %在[1,BitLength-1]范围内随机产生一个交叉位
scro(1,:)=[population(seln(1),1:chb) population(seln(2),chb+1:BitLength)];
scro(2,:)=[population(seln(2),1:chb) population(seln(1),chb+1:BitLength)];
else
scro(1,:)=population(seln(1),:);
scro(2,:)=population(seln(2),:);
end
%子程序：计算适应度函数, 函数名称存储为fitnessfun
function [Fitvalue,cumsump]=fitnessfun(population);
global BitLength
global boundsbegin
global boundsend
popsize=size(population,1); %有popsize个个体
for i=1:popsize
x=transform2to10(population(i,:)); %将二进制转换为十进制
%转化为[-2,2]区间的实数
xx=boundsbegin+x*(boundsend-boundsbegin)/(power((boundsend),BitLength)-1);
Fitvalue(i)=targetfun(xx); %计算函数值，即适应度
end
%给适应度函数加上一个大小合理的数以便保证种群适应值为正数
Fitvalue=Fitvalue'+230;
%计算选择概率
fsum=sum(Fitvalue);
Pperpopulation=Fitvalue/fsum;
%计算累积概率
cumsump(1)=Pperpopulation(1);
for i=2:popsize
cumsump(i)=cumsump(i-1)+Pperpopulation(i);
end
cumsump=cumsump';
%子程序：新种群变异操作，函数名称存储为mutation.m
function snnew=mutation(snew,pmutation);
BitLength=size(snew,2);
snnew=snew;
pmm=IfCroIfMut(pmutation); %根据变异概率决定是否进行变异操作，1则是，0则否
if pmm==1
chb=round(rand*(BitLength-1))+1; %在[1,BitLength]范围内随机产生一个变异位
snnew(chb)=abs(snew(chb)-1);
end
%子程序：判断遗传运算是否需要进行交叉或变异, 函数名称存储为IfCroIfMut.m
function pcc=IfCroIfMut(mutORcro);
test(1:100)=0;
l=round(100*mutORcro);
test(1:l)=1;
n=round(rand*99)+1;
pcc=test(n);
%子程序：新种群选择操作, 函数名称存储为selection.m
function seln=selection(population,cumsump);
%从种群中选择两个个体
for i=1:2
r=rand; %产生一个随机数
prand=cumsump-r;
j=1;
while prand(j)<0
j=j+1;
end
seln(i)=j; %选中个体的序号
end
%子程序：将2进制数转换为10进制数,函数名称存储为transform2to10.m
function x=transform2to10(Population);
BitLength=size(Population,2);
x=Population(BitLength);
for i=1:BitLength-1
x=x+Population(BitLength-i)*power(2,i);
end
%子程序：对于优化最大值或极大值函数问题，目标函数可以作为适应度函数
%函数名称存储为targetfun.m
function y=targetfun(x); %目标函数
y=200*exp(-0.05*x).*sin(x);
data crusor可以在函数图像上确定具体坐标
遗传算法应用：多约束非线性规划问题的求解
%主程序：本程序采用遗传算法接力进化，
%将上次进化结束后得到的最终种群作为下次输入的初始种群
clc;
close all;
clear all;
%进化的代数
T=100;
optionsOrigin=gaoptimset('Generations',T/2);
[x,fval,reason,output,finnal_pop]=ga(@ch14_2f,2,optionsOrigin);
%进行第二次接力进化
options1=gaoptimset('Generations',T/2,'InitialPopulation',finnal_pop,...
'PlotFcns',@gaplotbestf);
[x,fval,reason,output,finnal_pop]=ga(@ch14_2f,2,options1);
Bestx=x
BestFval=fval
%子函数：适应度函数同时也是目标函数,函数存储名称为ch14_2f.m
function f=ch14_2f(x)
g1=1.5+x(1)*x(2)-x(1)-x(2);
g2=-x(1)*x(2);
if(g1>0|g2>10)
f=100;
else
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
end
文档：电子商务中转化率影响因素研究——GEAT…
链接：http://note.youdao.com/noteshare?id=ec02fdbfa84ab8bb768736831b5597b5&sub=9BA9AEC84ADE496ABBAE10CB7BF94696
文档：遗传算法文献.note
链接：http://note.youdao.com/noteshare?id=1da272bd2d515d57d95e2f6f672a3a8e&sub=854893AC8A1D4514A9766D817BFAADEA
许多优化问题的目标函数是非凸的，存在许多局部最优解。
求解全局优化问题的方法有两类：
而模拟退火算法(SA)是一种通用概率算法，在大的搜寻空间内寻找问题的最优解。能有效解决NP难问题、避免陷入局部最优、对初值没有强依赖关系。
文档：模拟退火算法基本思想.note
链接：http://note.youdao.com/noteshare?id=225ff03a6cb79987cc17f0e74fd2b305&sub=1763A00DCA7A4B54B0B9BFF514EC0C9A
文档：TSP算法设计.note
链接：http://note.youdao.com/noteshare?id=63ee358a3d02b6d50f1d9d904339d288&sub=CE8D0C9FD5B0422EBF7660D8325D5BD6
应用：TSP模拟退火算法
clear
clc
a = 0.99; % 温度衰减函数的参数
t0 = 97; tf = 3; t = t0;
Markov_length = 10000; % Markov链长度
coordinates = [
1 565.0 575.0; 2 25.0 185.0; 3 345.0 750.0;
4 945.0 685.0; 5 845.0 655.0; 6 880.0 660.0;
7 25.0 230.0; 8 525.0 1000.0; 9 580.0 1175.0;
10 650.0 1130.0; 11 1605.0 620.0; 12 1220.0 580.0;
13 1465.0 200.0; 14 1530.0 5.0; 15 845.0 680.0;
16 725.0 370.0; 17 145.0 665.0; 18 415.0 635.0;
19 510.0 875.0; 20 560.0 365.0; 21 300.0 465.0;
22 520.0 585.0; 23 480.0 415.0; 24 835.0 625.0;
25 975.0 580.0; 26 1215.0 245.0; 27 1320.0 315.0;
28 1250.0 400.0; 29 660.0 180.0; 30 410.0 250.0;
31 420.0 555.0; 32 575.0 665.0; 33 1150.0 1160.0;
34 700.0 580.0; 35 685.0 595.0; 36 685.0 610.0;
37 770.0 610.0; 38 795.0 645.0; 39 720.0 635.0;
40 760.0 650.0; 41 475.0 960.0; 42 95.0 260.0;
43 875.0 920.0; 44 700.0 500.0; 45 555.0 815.0;
46 830.0 485.0; 47 1170.0 65.0; 48 830.0 610.0;
49 605.0 625.0; 50 595.0 360.0; 51 1340.0 725.0;
52 1740.0 245.0;
];
coordinates(:,1) = [];
amount = size(coordinates,1); % 城市的数目
% 通过向量化的方法计算距离矩阵
dist_matrix = zeros(amount, amount);
coor_x_tmp1 = coordinates(:,1) * ones(1,amount);
coor_x_tmp2 = coor_x_tmp1';
coor_y_tmp1 = coordinates(:,2) * ones(1,amount);
coor_y_tmp2 = coor_y_tmp1';
dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + ...
(coor_y_tmp1-coor_y_tmp2).^2);
sol_new = 1:amount; % 产生初始解
% sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解；
E_current = inf;E_best = inf; % E_current是当前解对应的回路距离；
% E_new是新解的回路距离；
% E_best是最优解的
sol_current = sol_new; sol_best = sol_new;
p = 1;
while t>=tf
for r=1:Markov_length % Markov链长度
% 产生随机扰动
if (rand < 0.5) % 随机决定是进行两交换还是三交换
% 两交换
ind1 = 0; ind2 = 0;
while (ind1 == ind2)
ind1 = ceil(rand.*amount);
ind2 = ceil(rand.*amount);
end
tmp1 = sol_new(ind1);
sol_new(ind1) = sol_new(ind2);
sol_new(ind2) = tmp1;
else
% 三交换
ind1 = 0; ind2 = 0; ind3 = 0;
while (ind1 == ind2) || (ind1 == ind3) ...
|| (ind2 == ind3) || (abs(ind1-ind2) == 1)
ind1 = ceil(rand.*amount);
ind2 = ceil(rand.*amount);
ind3 = ceil(rand.*amount);
end
tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
% 确保ind1 < ind2 < ind3
if (ind1 < ind2) && (ind2 < ind3)
;
elseif (ind1 < ind3) && (ind3 < ind2)
ind2 = tmp3;ind3 = tmp2;
elseif (ind2 < ind1) && (ind1 < ind3)
ind1 = tmp2;ind2 = tmp1;
elseif (ind2 < ind3) && (ind3 < ind1)
ind1 = tmp2;ind2 = tmp3; ind3 = tmp1;
elseif (ind3 < ind1) && (ind1 < ind2)
ind1 = tmp3;ind2 = tmp1; ind3 = tmp2;
elseif (ind3 < ind2) && (ind2 < ind1)
ind1 = tmp3;ind2 = tmp2; ind3 = tmp1;
end
tmplist1 = sol_new((ind1+1):(ind2-1));
sol_new((ind1+1):(ind1+ind3-ind2+1)) = ...
sol_new((ind2):(ind3));
sol_new((ind1+ind3-ind2+2):ind3) = ...
tmplist1;
end
%检查是否满足约束
% 计算目标函数值（即内能）
E_new = 0;
for i = 1 : (amount-1)
E_new = E_new + ...
dist_matrix(sol_new(i),sol_new(i+1));
end
% 再算上从最后一个城市到第一个城市的距离
E_new = E_new + ...
dist_matrix(sol_new(amount),sol_new(1));
if E_new < E_current
E_current = E_new;
sol_current = sol_new;
if E_new < E_best
% 把冷却过程中最好的解保存下来
E_best = E_new;
sol_best = sol_new;
end
else
% 若新解的目标函数值小于当前解的，
% 则仅以一定概率接受新解
if rand < exp(-(E_new-E_current)./t)
E_current = E_new;
sol_current = sol_new;
else
sol_new = sol_current;
end
end
end
t=t.*a; % 控制参数t（温度）减少为原来的a倍
end
disp('最优解为：')
disp(sol_best)
disp('最短距离：')
disp(E_best)
应用：背包问题模拟退火算法
clear
clc
a = 0.95
k = [5;10;13;4;3;11;13;10;8;16;7;4];%12件物品质量
k = -k; % 模拟退火算法是求解最小值，故取负数
d = [2;5;18;3;2;5;10;4;11;7;14;6];%12件物品对应的价值
restriction = 46;%背包最大载重
num = 12;
sol_new = ones(1,num); % 生成初始解
E_current = inf;E_best = inf;
% E_current是当前解对应的目标函数值（即背包中物品总价值）；
% E_new是新解的目标函数值；
% E_best是最优解的
sol_current = sol_new; sol_best = sol_new;
t0=97; tf=3; t=t0;
p=1;
while t>=tf
for r=1:100
%产生随机扰动
tmp=ceil(rand.*num);
sol_new(1,tmp)=~sol_new(1,tmp);
%检查是否满足约束
while 1
q=(sol_new*d <= restriction);
if ~q
p=~p; %实现交错着逆转头尾的第一个1
tmp=find(sol_new==1);
if p
sol_new(1,tmp)=0;
else
sol_new(1,tmp(end))=0;
end
else
break
end
end
% 计算背包中的物品价值
E_new=sol_new*k;
if E_new<E_current
E_current=E_new;
sol_current=sol_new;
if E_new<E_best
% 把冷却过程中最好的解保存下来
E_best=E_new;
sol_best=sol_new;
end
else
if rand<exp(-(E_new-E_current)./t)
E_current=E_new;
sol_current=sol_new;
else
sol_new=sol_current;
end
end
end
t=t.*a;
end
disp('最优解为：')
sol_best
disp('物品总价值等于：')
val=-E_best;
disp(val)
disp('背包中物品重量是：')
disp(sol_best * d)
文档：ASA模拟退火程序包.note
链接：http://note.youdao.com/noteshare?id=c3c3b39f72cc2ce17466a466958f7cb0&sub=F0766F4B838B4C9CB1F83DCECA89797D
文档：小结.note
链接：http://note.youdao.com/noteshare?id=ed501ea8cd73835bc5aa00f98984dee7&sub=6564713C89664C7B851D548482D4E5F1
文档：模拟退火文献.note
链接：http://note.youdao.com/noteshare?id=bd5f93b397983e4cd299e0bee92fb36b&sub=FFAB2E47FA1E42CB9396E1D098C68CB1
人工神经网络(ANN)由基本元件——神经元相互连接。
前向反馈(BP)网络和径向基(RBF)网络是应用范围广泛的两种网络。
常规BP是梯度学习算法的RBF的特殊形式。
文档：常用激励函数.note
链接：http://note.youdao.com/noteshare?id=f494d2b14013ca9179770fc90a944bb6&sub=1CE8FEDF31CD4769AF01AE58610E3167
文档：BP、RBF数学原理.note
链接：http://note.youdao.com/noteshare?id=d8231fa67ba4cde9898337bfc89459df&sub=4135A4B2C6544C59983DFC962B1A69FB
文档：BP神经网络的结构设计.note
链接：http://note.youdao.com/noteshare?id=9f7f18304ab669aa0b617225bb2c6c1f&sub=C3E6A73395784AD9BE7BA6C722D80BE2
文档：RBF神经网络的结构设计.note
链接：http://note.youdao.com/noteshare?id=27b41254c5e998673173dae7c66f7101&sub=45B899BB6DF74812B4D34DD1C2FBA1C7
应用：基于BP公路运量预测(不依赖于神经网络工具箱，而是源程序)
function main()
clc % 清屏
clear all; %清除内存以便加快运算速度
close all; %关闭当前所有figure图像
SamNum=20; %输入样本数量为20
TestSamNum=20; %测试样本数量也是20
ForcastSamNum=2; %预测样本数量为2
HiddenUnitNum=8; %中间层隐节点数量取8,比工具箱程序多了1个
InDim=3; %网络输入维度为3
OutDim=2; %网络输出维度为2
%一、原始数据的输入
%人数(单位：万人)
sqrs=[20.55 22.44 25.37 27.13 29.45 30.10 30.96 34.06 36.42 38.09 39.13 39.99 ...
41.93 44.59 47.30 52.89 55.73 56.76 59.17 60.63];
%机动车数(单位：万辆)
sqjdcs=[0.6 0.75 0.85 0.9 1.05 1.35 1.45 1.6 1.7 1.85 2.15 2.2 2.25 2.35 2.5 2.6...
2.7 2.85 2.95 3.1];
%公路面积(单位：万平方公里)
sqglmj=[0.09 0.11 0.11 0.14 0.20 0.23 0.23 0.32 0.32 0.34 0.36 0.36 0.38 0.49 ...
0.56 0.59 0.59 0.67 0.69 0.79];
%公路客运量(单位：万人)
glkyl=[5126 6217 7730 9145 10460 11387 12353 15750 18304 19836 21024 19490 20433 ...
22598 25107 33442 36836 40548 42927 43462];
%公路货运量(单位：万吨)
glhyl=[1237 1379 1385 1399 1663 1714 1834 4322 8132 8936 11099 11203 10524 11115 ...
13320 16762 18673 20724 20803 21804];
p=[sqrs;sqjdcs;sqglmj]; %输入数据矩阵
t=[glkyl;glhyl]; %目标数据矩阵
%二、数据归一化
[SamIn,minp,maxp,tn,mint,maxt]=premnmx(p,t); %原始样本对（输入和输出）初始化
rand('state',sum(100*clock)) %依据系统时钟种子产生随机数
NoiseVar=0.01; %噪声强度为0.01（添加噪声的目的是为了防止网络过度拟合）
Noise=NoiseVar*randn(2,SamNum); %生成噪声
SamOut=tn + Noise; %将噪声添加到输出样本上
TestSamIn=SamIn; %这里取输入样本与测试样本相同因为样本容量偏少
TestSamOut=SamOut; %也取输出样本与测试样本相同
MaxEpochs=50000; %最多训练次数为50000
lr=0.035; %学习速率为0.035
E0=0.65*10^(-3); %目标误差为0.65*10^(-3)
W1=0.5*rand(HiddenUnitNum,InDim)-0.1; %初始化输入层与隐含层之间的权值
B1=0.5*rand(HiddenUnitNum,1)-0.1; %初始化输入层与隐含层之间的阈值
W2=0.5*rand(OutDim,HiddenUnitNum)-0.1; %初始化输出层与隐含层之间的权值
B2=0.5*rand(OutDim,1)-0.1; %初始化输出层与隐含层之间的阈值
%三、网络训练
ErrHistory=[]; %给中间变量预先占据内存
for i=1:MaxEpochs
HiddenOut=logsig(W1*SamIn+repmat(B1,1,SamNum)); % 隐含层网络输出
NetworkOut=W2*HiddenOut+repmat(B2,1,SamNum); % 输出层网络输出
Error=SamOut-NetworkOut; % 实际输出与网络输出之差
SSE=sumsqr(Error) %能量函数（误差平方和）
ErrHistory=[ErrHistory SSE];
if SSE<E0,break, end %如果达到误差要求则跳出学习循环
% 以下六行是BP网络最核心的程序
% 他们是权值（阈值）依据能量函数负梯度下降原理所作的每一步动态调整量
Delta2=Error;
Delta1=W2'*Delta2.*HiddenOut.*(1-HiddenOut);
dW2=Delta2*HiddenOut';
dB2=Delta2*ones(SamNum,1);
dW1=Delta1*SamIn';
dB1=Delta1*ones(SamNum,1);
%对输出层与隐含层之间的权值和阈值进行修正
W2=W2+lr*dW2;
B2=B2+lr*dB2;
%对输入层与隐含层之间的权值和阈值进行修正
W1=W1+lr*dW1;
B1=B1+lr*dB1;
end
%四、对原始数据进行仿真
HiddenOut=logsig(W1*SamIn+repmat(B1,1,TestSamNum)); % 隐含层输出最终结果
NetworkOut=W2*HiddenOut+repmat(B2,1,TestSamNum); % 输出层输出最终结果
a=postmnmx(NetworkOut,mint,maxt); % 还原网络输出层的结果
x=1990:2009; % 时间轴刻度
newk=a(1,:); % 网络输出客运量
newh=a(2,:); % 网络输出货运量
%五、将原始数据仿真与已知样本进行对比
figure ;
subplot(2,1,1);plot(x,newk,'r-o',x,glkyl,'b--+') %绘值公路客运量对比图；
legend('网络输出客运量','实际客运量');
xlabel('年份');ylabel('客运量/万人');
subplot(2,1,2);plot(x,newh,'r-o',x,glhyl,'b--+') %绘制公路货运量对比图；
legend('网络输出货运量','实际货运量');
xlabel('年份');ylabel('货运量/万吨');
% 六、利用训练好的网络进行预测
% 当用训练好的网络对新数据pnew进行预测时，也应作相应的处理
pnew=[73.39 75.55
3.9635 4.0975
0.9880 1.0268]; %2010年和2011年的相关数据；
pnewn=tramnmx(pnew,minp,maxp); %利用原始输入数据的归一化参数对新数据进行归一化；
HiddenOut=logsig(W1*pnewn+repmat(B1,1,ForcastSamNum)); % 隐含层输出预测结果
anewn=W2*HiddenOut+repmat(B2,1,ForcastSamNum); % 输出层输出预测结果
%把网络预测得到的数据还原为原始的数量级；
anew=postmnmx(anewn,mint,maxt)
应用：基于BP公路运量预测(神经网络工具箱版本)
clc;
close all;
clear all;
% 原始数据
% 人数(单位：万人)
sqrs=[20.55 22.44 25.37 27.13 29.45 30.10 30.96 34.06 36.42 38.09 39.13 39.99 ...
41.93 44.59 47.30 52.89 55.73 56.76 59.17 60.63];
%机动车数(单位：万辆)
sqjdcs=[0.6 0.75 0.85 0.9 1.05 1.35 1.45 1.6 1.7 1.85 2.15 2.2 2.25 2.35 2.5 2.6...
2.7 2.85 2.95 3.1];
%公路面积(单位：万平方公里)
sqglmj=[0.09 0.11 0.11 0.14 0.20 0.23 0.23 0.32 0.32 0.34 0.36 0.36 0.38 0.49 ...
0.56 0.59 0.59 0.67 0.69 0.79];
%公路客运量(单位：万人)
glkyl=[5126 6217 7730 9145 10460 11387 12353 15750 18304 19836 21024 19490 20433 ...
22598 25107 33442 36836 40548 42927 43462];
%公路货运量(单位：万吨)
glhyl=[1237 1379 1385 1399 1663 1714 1834 4322 8132 8936 11099 11203 10524 11115 ...
13320 16762 18673 20724 20803 21804];
% 输入数据矩阵
p=[sqrs;sqjdcs;sqglmj];
% 目标数据矩阵
t=[glkyl;glhyl];
% 利用mapminmax函数对数据进行归一化
[pn,input_str] = mapminmax(p) ;
[tn,output_str] = mapminmax(t) ;
% 建立BP神经网络，相对旧一点的MATLAB版本，新版本 newff 函数使用更简洁一些
% 但是本质和性能没有区别
net=newff(pn,tn,[3 7 2],{'purelin','logsig','purelin'});
% 10轮回显示一次结果
net.trainParam.show=10;
% 学习速度为0.05
net.trainParam.lr=0.05;
% 最大训练次数为5000次
net.trainParam.epochs=5000;
% 均方误差
net.trainParam.goal=0.65*10^(-3);
% 网络误差如果连续6次迭代都没有变化，训练将会自动终止（系统默认的）
% 为了让程序继续运行，用以下命令取消这条设置
net.divideFcn = '';
% 开始训练，其中pn,tn分别为输入输出样本
net=train(net,pn,tn);
% 利用训练好的网络，基于原始数据对BP网络仿真
an=sim(net,pn);
% 利用函数mapminmax把仿真得到的数据还原为原始的数量级
% 新版本推荐训练样本归一化和反归一化都使用 mapminmax 函数
a = mapminmax('reverse',an,output_str);
% 本例因样本容量有限使用训练数据进行测试，通常必须用新鲜数据进行测试
x=1990:2009;
newk=a(1,:);
newh=a(2,:);
figure (2);
% 绘值公路客运量对比图
subplot(2,1,1);plot(x,newk,'r-o',x,glkyl,'b--+')
legend('网络输出客运量','实际客运量');
xlabel('年份');ylabel('客运量/万人');
title('运用工具箱客运量学习和测试对比图');
% 绘制公路货运量对比图
subplot(2,1,2);plot(x,newh,'r-o',x,glhyl,'b--+')
legend('网络输出货运量','实际货运量');
xlabel('年份');ylabel('货运量/万吨');
title('运用工具箱货运量学习和测试对比图');
% 利用训练好的网络进行预测
% 2010年和2011年的相关数据
% 当用训练好的网络对新数据mapminmax进行预测时，也应作相应的处理：
pnew=[73.39 75.55
3.9635 4.0975
0.9880 1.0268];
% 利用原始输入数据的归一化参数对新数据进行归一化
pnewn = mapminmax('apply',pnew,input_str);
% 利用归一化后的数据进行仿真
anewn=sim(net,pnewn);
% 把仿真得到的数据还原为原始的数量级
anew = mapminmax('reverse',anewn,output_str)
基于BP的应用：2006B艾滋病治疗最佳停药时间的确定
尝试预测继续治疗的效果，确定最佳治疗终止时间(只是样本数据和参数有些变化，其他和上述程序相同)
文档：2006B.note
链接：http://note.youdao.com/noteshare?id=ed4379423eec7c31b6f5c83585a3b32e&sub=CE16EBE337E04CC3BB6AD403E7C57ECB
基于RBF的应用：预测新顾客流失概率
文档：客服流失题目.note
链接：http://note.youdao.com/noteshare?id=1f0949b901a7bf58be7d7781b986333a&sub=DC4AA853D5744D7B89EFA451E8B17506
% 本程序添加了动量因子
% 本程序是基于梯度训练算法的RBF网络
% 所以添加动量因子以便不会轻易陷入局部极小值
function main()
clc;
close all;
clear all;
warning off;
% 初始化参数
% 该网络有三层，输入层和输出层都是线性函数，隐含层由距离函数和激活函数构成
SamNum = 120;
TargetSamNum = 60;
% 样本输入维度
InDim = 1;
% 隐含层神经元数量
UnitNum = 10;
MaxEpoch = 10000;
% 目标误差
E0 = 0.09;
% 输入置于[1,60]区间的随机数
% 理论样本
SamIn = sort(59*rand(1,SamNum)+1);
SamOut = 0.5447*SamIn.^0.1489;
% 实际样本
TargetIn = 1:60;
TargetOut = [0.53173198482933
0.599828865
0.644564773
0.671027441
0.697281167
0.717013297
0.732752613
0.745040151
0.75565936
0.763524144
0.779177473
0.792189854
0.806571209
0.813644571
0.822233807
0.826976013
0.837737352
0.842773177
0.854878049
0.859771055
0.863536819
0.865907219
0.869966906
0.872734818
0.875641915
0.878079332
0.881514601
0.886842845
0.891857506
0.898078292
0.906074968
0.910126947
0.91328894
0.917005814
0.920081668
0.924666569
0.928067079
0.932732111
0.936609264
0.940518784
0.94417839
0.946870779
0.958960328
0.961151737
0.963206107
0.964973998
0.967341306
0.96778647
0.968232044
0.970466082
0.974362934
0.98011496
0.98424337
0.987633062
0.991046183
0.995581505
0.997785861
1
1
1]';
figure;
hold on;
% 添加边框和网格线
box on;
grid on;
plot(SamIn,SamOut,'kO');
plot(TargetIn,TargetOut,'b-');
xlabel('x');
ylabel('y');
title('训练和测试图');
% 原始样本对（输入和输出）初始化
[SamIn,minp,maxp,SamOut,mint,maxt] = premnmx(SamIn,SamOut);
% 利用原始输入数据的归一化参数对新数据进行归一化；
TargetIn = tramnmx(TargetIn,minp,maxp);
TargetOut = tramnmx(TargetOut,mint,maxt);
% 初始化数据中心
Center = 8*rand(InDim,UnitNum)-4;
% 初始化宽度
SP = 0.2*rand(1,UnitNum)+0.1;
% 初始化权值，这里的权值是指隐含层与输出层之间的权值
% 跟BP网络不同的是，这里没有输入层与隐含层之间的权值
% 所以单从网络结构上讲，RBF网其实更简单
W = 0.2*rand(1,UnitNum)-0.1;
% 数据中心学习速度（速率）
lrCent = 0.02;
% 宽度学习速度（速率）
lrSP = 0.001;
% 权值学习速度（速率）
lrW = 0.001;
% 动量因子系数
arf = 0.001;
% 用来存储误差
ErrHistory = [];
for epoch = 1:MaxEpoch
% 计算书中输出样本与数据中心之间的距离（这里是欧式距离）
% 相当于书中的||x-c||表达式，其具体表达式为 sum((x-y).^2).^0.5
AllDist = dist(Center',SamIn);
SPMat = repmat(SP',1,SamNum);
% 以高斯函数作为激活函数，高斯函数表达式为exp(-n^2)
UnitOut = radbas(AllDist./SPMat);
% 隐含层神经元数据经过加权后在输出层输出
% 输出层是线性激活函数所以直接加权输出即可
NetOut = W*UnitOut;
Error = SamOut-NetOut;
SSE = sumsqr(Error)
ErrHistory = [ErrHistory SSE];
if SSE<E0,break,end
% 初始化用来存储前一次调整量的变量，全部用0矩阵填充
CentGrad = zeros(size(Center(:,1)));
SPGrad = zeros(size(SP(1)));
WGrad = zeros(size(W(1)));
for i = 1:UnitNum
% 存储前一次训练的调整量，以下是数据中心前一次调整量赋给变量CenterPre
% 以便后续动量因子会用到
CenterPre = CentGrad ;
% 宽度前一次调整量
SPPre = SPGrad ;
% 隐含层与输出层之间的连接权值矩阵前一次调整量
WPre = WGrad ;
% 数据中心的网络反向调整
% 原理仍然是链式偏导
% 根据求导，数据中心的调整量CentGrad等于 （样本-数据中心）*误差*隐含层网络输出
% *权值/宽度平方
% 以下宽度调整量和权值调整量类似
CentGrad=(SamIn-repmat(Center(:,i),1,SamNum))*(Error.*UnitOut(i,:)*W(i)/(SP(i)^2))';
SPGrad = AllDist(i,:).^2*(Error.*UnitOut(i,:)*W(i)/(SP(i)^3))';
WGrad = Error*UnitOut(i,:)';
% 更新网络的数据中心、宽度以及权值
Center(:,i) = Center(:,i)+lrCent*CentGrad;
SP(i) = SP(i)+lrSP*SPGrad;
W(i) = W(i)+lrW*WGrad;
% 梯度法训练RBF网络同样不能幸免极易陷入局部极小值的缺陷，所以这里添加动量因子
Center(:,i) = Center(:,i)+ arf*CenterPre ;
SP(i) = SP(i)+ arf*SPPre ;
W(i) = W(i) + arf*WPre ;
end
end
% 测试网络的性能如何
TestDistance = dist(Center',TargetIn);
TestSpreadsMat = repmat(SP',1,TargetSamNum);
TestHiddenUnitOut = radbas(TestDistance./TestSpreadsMat);
TestNNOut = W*TestHiddenUnitOut;
% 训练样本进行了归一化，所以需要还原
TestNNOut = postmnmx(TestNNOut,mint,maxt);
TargetIn = postmnmx(TargetIn,minp,maxp);
plot(TargetIn,TestNNOut,'r-');
axis tight;
legend('理论样本','实际样本','测试样本');
figure;
hold on;
grid;
box;
[xx,Num] = size(ErrHistory);
plot(1:Num,ErrHistory,'k-');
xlabel('训练次数');
ylabel('误差大小');
title('训练误差图');
梯度法的RBF神经网络比常规的BP网络训练速度更快。
文档：神经网络文献.note
链接：http://note.youdao.com/noteshare?id=d094c8818bb8c1ba130fb30385f98311&sub=6080D953253E46BEB80BD4A89420DC0C
现代算法分为硬计算和软计算。
硬计算需要建立数学模型，软计算(智能算法)是一种动态自适应求解方式，不需要深入数学模型。
PSO依托群鸟觅食的模型寻求最优解。
文档：PSO算法基本理论.note
链接：http://note.youdao.com/noteshare?id=c47f5e0365e99528c1d1e30ef813cf09&sub=8C5E09CF93004071A78EE8CD804AD2DD文档：PSO算法的约束优化.note
链接：http://note.youdao.com/noteshare?id=0f0ebf18afe5c4d27ac955d3a60bfea5&sub=9C897D7C9E694FBFB09C52993CE1EF0A文档：PSO优缺点.note
链接：http://note.youdao.com/noteshare?id=fae7b4fc2ce0cf2dd2fcd94158b82027&sub=49A2132BB7404B4C88B799747AF6580B文档：PSO算法设计.note
链接：http://note.youdao.com/noteshare?id=faf2db5dace6fa48c028eb59e345f4de&sub=72FFC23F0DA5437A829D083FC52F0ED4
根据PSO思路，计算
m a x f ( x ) = 2.1 ( 1 − x + 2 x 2 ) e x p ( − x 2 2 ) , x ∈ [ − 5 , 5 ] maxf(x)=2.1(1-x+2x^2)exp(-\frac{x^2}{2}),x∈[-5,5] maxf(x)=2.1(1−x+2x2)exp(−2x2),x∈[−5,5]
function main()
clc;clear all;close all;
tic; %程序运行计时
E0=0.001; %允许误差
MaxNum=100; %粒子最大迭代次数
narvs=1; %目标函数的自变量个数
particlesize=30; %粒子群规模
c1=2; %每个粒子的个体学习因子，也称为加速常数
c2=2; %每个粒子的社会学习因子，也称为加速常数
w=0.6; %惯性因子
vmax=0.8; %粒子的最大飞翔速度
x=-5+10*rand(particlesize,narvs); %粒子所在的位置
v=2*rand(particlesize,narvs); %粒子的飞翔速度
%用inline定义适应度函数以便将子函数文件与主程序文件放在一起，
%目标函数是：y=1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2))
%inline命令定义适应度函数如下：
fitness=inline('1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');
%inline定义的适应度函数会使程序运行速度大大降低
for i=1:particlesize
for j=1:narvs
f(i)=fitness(x(i,j));
end
end
personalbest_x=x;
personalbest_faval=f;
[globalbest_faval i]=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
k=1;
while k<=MaxNum
for i=1:particlesize
for j=1:narvs
f(i)=fitness(x(i,j));
end
if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置
personalbest_faval(i)=f(i);
personalbest_x(i,:)=x(i,:);
end
end
[globalbest_faval i]=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
for i=1:particlesize %更新粒子群里每个个体的最新位置
v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...
+c2*rand*(globalbest_x-x(i,:));
for j=1:narvs %判断粒子的飞翔速度是否超过了最大飞翔速度
if v(i,j)>vmax;
v(i,j)=vmax;
elseif v(i,j)<-vmax;
v(i,j)=-vmax;
end
end
x(i,:)=x(i,:)+v(i,:);
end
if abs(globalbest_faval)<E0,break,end
k=k+1;
end
Value1=1/globalbest_faval-1; Value1=num2str(Value1);
% strcat指令可以实现字符的组合输出
disp(strcat('the maximum value','=',Value1));
%输出最大值所在的横坐标位置
Value2=globalbest_x; Value2=num2str(Value2);
disp(strcat('the corresponding coordinate','=',Value2));
x=-5:0.01:5;
y=2.1*(1-x+2*x.^2).*exp(-x.^2/2);
plot(x,y,'m-','linewidth',3);
hold on;
plot(globalbest_x,1/globalbest_faval-1,'kp','linewidth',4);
legend('目标函数','搜索到的最大值');xlabel('x');ylabel('y');grid on;toc;
BP神经网络容易陷入局部极小值，PSO在无约束非线性函数优化方面可以找到全局最优解。因此利用PSO算法和BP算法共同训练神经网络，先将网络进行PSO算法训练(将输出误差函数看做目标函数，进行全局寻找最小值)，然后BP算法接着进行小范围精细搜索。
评价网络的性能
文档：BP能够找到极值的原理.note
链接：http://note.youdao.com/noteshare?id=5542ca38cdcc3cf643706756e6917a61&sub=830245F7A41F49D99A7D3BAB09FA4EFD文档：PSO-BP设计.note
链接：http://note.youdao.com/noteshare?id=91bc791fc7a2705888860417a3419a34&sub=6E638802C7A449538DC6C7DE362F4C39文档：SPOS–PSO优化神经网络结构.note
链接：http://note.youdao.com/noteshare?id=7e5ca3d242b5fdc06f9234f60c04b2f5&sub=7C461C7C6D3F4A7499B419A80658FA2C
PSO优化BP神经网路实现：
文档：PSO、BP、PSO-BP对比.note
链接：http://note.youdao.com/noteshare?id=55e94936768f1949bc1b37594c24605c&sub=223C1E8F80784C01875A95080DBFFBD6
%PSOBP502.m
%耗时10min
function main
clc;clear all;close all;
MaxRunningTime=1; %该参数是为了使网络集成，重复训练MaxRunningTime次
HiddenUnitNum=12;%神经元个数
rand('state',sum(100*clock));
TrainSamIn=-4:0.07:2.5;
TrainSamOut=1.1*(1-TrainSamIn+2*TrainSamIn.^2).*exp(-TrainSamIn.^2/2);
TestSamIn=2:0.04:3;
TestSamOut=1.1*(1-TestSamIn+2*TestSamIn.^2).*exp(-TestSamIn.^2/2);
[xxx,TrainSamNum]=size(TrainSamIn);
[xxx,TestSamNum]=size(TestSamIn);
% for HiddenUnitNum=3:MaxHiddenLayerNode %隐含层神经元的个数可以取逐渐增大的合理整数
fprintf('\n the hidden layer node');HiddenUnitNum
TrainNNOut=[];
TestNNOut=[];
for t=1:MaxRunningTime
fprintf('the current running times is');t
[NewW1,NewB1,NewW2,NewB2]=PSOTrain(TrainSamIn,TrainSamOut,HiddenUnitNum);
disp('PSO算法训练神经网络结束，BP算法接着训练网络……');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%BP算法参数初始化，注意与上面PSO参数一致
SamInNum=length(TrainSamIn);
TestSamNum=length(TestSamIn);
InDim=1;
OutDim=1;
%学习样本添加噪声
rand('state',sum(100*clock))
NoiseVar=0.01;
Noise=NoiseVar*randn(1,SamInNum);
SamIn=TrainSamIn;
SamOutNoNoise=TrainSamOut;
SamOut=SamOutNoNoise + Noise;
%BP算法参数
MaxEpochs=300;%训练次数
lr=0.003;%学习速率
E0=0.0001;
W1=NewW1;
B1=NewB1;
W2=NewW2';
B2=NewB2;
W1Ex=[W1 B1];
W2Ex=[W2 B2];
SamInEx=[SamIn' ones(SamInNum,1)]';
ErrHistory=[];
%网络参数初始化完毕
%给动画初始画图和构建动画框架和背景
HiddenOut=logsig(W1Ex*SamInEx);
HiddenOutEx=[HiddenOut' ones(SamInNum,1)]';
NetworkOut=W2Ex*HiddenOutEx;
Error=SamOut-NetworkOut;
%给误差的动画显示提供空矩阵和其维数
SSEINIT=zeros(1,MaxEpochs);
%这仅限于输出是一维的情况
SSE=sumsqr(Error);
%绘画动画显示的图形框架
figure(1);
rangecolour=linspace(0,1,MaxEpochs);
%采用分区画图，把两幅动画在一个figure中显示
%先画第一幅图形
subplot(2,1,1);
hold on
axis([1 SamInNum min(SamOut) max(SamOut)]);
hflash1=line(1:SamInNum,SamOut,'color',[rangecolour(1) 0 1-rangecolour(1)],...
'linestyle','-','linewidth',2,'marker','d',...
'markersize',2,'erasemode','none');
hflash2=line(1:SamInNum,NetworkOut,'color',[rangecolour(1) 0 1-rangecolour(1)],...
'linestyle','-','linewidth',2.5,'marker','h',...
'markersize',2.3,'erasemode','xor');
xlabel('训练样本的数目');
ylabel('样本的输出值或网络的输出值');
title('样本的输出值与网络的输出值动画显示','fontsize',13);
legend('样本的输出值','网络的输出值');
hold off
%再画第二幅图形
subplot(2,1,2);
hold on
axis([1 MaxEpochs -0.2*SSE SSE]);
hflash3=line(1:MaxEpochs,E0*ones(1,MaxEpochs),'color',...
[rangecolour(1) 0 1-rangecolour(1)],...
'linestyle','--','linewidth',2,'marker','h',...
'markersize',2,'erasemode','none');
hflash4=line(1,SSE,'color',...
[rangecolour(1) 0 1-rangecolour(1)],...
'linestyle','-','linewidth',2,'marker','*',...
'markersize',2,'erasemode','xor');
xlabel('网络训练次数');
ylabel('目标误差或网络输出误差');
title('目标误差与网络输出误差动画显示','fontsize',13);
legend('目标误差','网络输出误差');
hold off
for i=2:MaxEpochs
HiddenOut=logsig(W1Ex*SamInEx);
HiddenOutEx=[HiddenOut' ones(SamInNum,1)]';
NetworkOut=W2Ex*HiddenOutEx;
Error=SamOut-NetworkOut;
SSE=sumsqr(Error)
%让第二幅动画逐点显示
SSEINIT(:,i)=SSE;
%对于神经网络训练过程中发生震荡的瞬间图像将其显示出来
ErrHistory=[ErrHistory SSE];
SSEINIT(:,1);
SSEINIT(:,2);
SSEINIT2=SSEINIT(:,i);
SSEINIT1=SSEINIT(:,i-1);
if SSE<E0,break, end
Delta2=Error;
Delta1=W2'*Delta2.*HiddenOut.*(1-HiddenOut);
dW2Ex=Delta2*HiddenOutEx';
dW1Ex=Delta1*SamInEx';
W1Ex=W1Ex+lr*dW1Ex;
W2Ex=W2Ex+lr*dW2Ex;
W2=W2Ex(:,1:HiddenUnitNum);
if SSEINIT2>SSEINIT1
%如果网络学习时发生震荡，每10步显示一次
if mod(i,10)==0
Counter(i)=SSEINIT(:,i);
Len=size(Counter);
figure(Len(1,2));
subplot(2,1,1);
hold on
axis([1 SamInNum min(SamOut) max(SamOut)]);
hflash5=line(1:SamInNum,SamOut,'color',[rangecolour(1) 0 1-rangecolour(1)],...
'linestyle','-','linewidth',2,'marker','d',...
'markersize',2,'erasemode','none');
hflash6=line(1:SamInNum,NetworkOut,'color',[rangecolour(MaxEpochs) 0 1-rangecolour(MaxEpochs)],...
'linestyle','-','linewidth',2.5,'marker','h',...
'markersize',2.3,'erasemode','xor');
xlabel('训练样本的数目');
ylabel('样本的输出值或网络的输出值');
title('神经网络学习震荡时拟合曲线','fontsize',13);
legend('样本的输出值','网络的输出值');
hold off
%再画第二幅图形
subplot(2,1,2);
hold on
axis([1 MaxEpochs -2*SSEINIT(:,2) 2*SSEINIT(:,2)]);
hflash7=line(1:MaxEpochs,E0*ones(1,MaxEpochs),'color',...
[rangecolour(1) 0 1-rangecolour(1)],...
'linestyle','--','linewidth',2,'marker','h',...
'markersize',2,'erasemode','none');
hflash8=line(1:i,SSEINIT(:,1:i),'color',...
[rangecolour(1) 0 1-rangecolour(1)],...
'linestyle','-','linewidth',2,'marker','*',...
'markersize',2,'erasemode','xor');
xlabel('网络训练次数');
ylabel('目标误差或网络输出误差');
title('神经网络学习震荡时误差','fontsize',13);
legend('目标误差','网络输出误差');
hold off
end
end
%动画开始放映
set(hflash2,'XData',1:SamInNum,'YData',NetworkOut,'color',...
[rangecolour(MaxEpochs) 0 1-rangecolour(MaxEpochs)]);
set(hflash4,'XData',1:i,'YData',SSEINIT(:,1:i),...
'color',[rangecolour(MaxEpochs) 0 1-rangecolour(MaxEpochs)]);
drawnow;
end
W2=W2Ex(:,1:HiddenUnitNum);
W1=W1Ex(:,1:InDim);
B1=W1Ex(:,InDim+1);
B2=W2Ex(:,1+HiddenUnitNum);
TrainHiddenOut=logsig(W1*SamIn+repmat(B1,1,SamInNum));
TrainNNOut=W2*TrainHiddenOut+repmat(B2,1,SamInNum);
TestHiddenOut=logsig(W1*TestSamIn+repmat(B1,1,TestSamNum));
TestNNOut=W2*TestHiddenOut+repmat(B2,1,TestSamNum);
figure(MaxEpochs+1);
hold on;
grid;
h1=plot(SamIn,SamOut);
set(h1,'color','r','linestyle','-',...
'linewidth',2.5,'marker','p','markersize',5);
hold on
h2=plot(TestSamIn,TestSamOut);
set(h2,'color','g','linestyle','--',...
'linewidth',2.5,'marker','^','markersize',7);
h3=plot(SamIn,TrainNNOut);
set(h3,'color','c','linestyle','-.',...
'linewidth',2.5,'marker','o','markersize',5);
h4=plot(TestSamIn,TestNNOut);
set(h4,'color','m','linestyle',':',...
'linewidth',2.5,'marker','s','markersize',5);
xlabel('Input x','fontsize',13);ylabel('Output y','fontsize',13);
box on;axis tight;
%title('PSO-BP神经网络误差测试图');
legend('网络学习实际样本值','网络测试实际样本值',...
'网络学习网络输出值','网络测试网络输出值');
hold off;
end
% end
fidW1=fopen('dW1.txt','a+');fidB1=fopen('B1.txt','a+');
fidW2=fopen('W2.txt','a+');fidB2=fopen('B2.txt','a+');
for i=1:length(W1)
fprintf(fidW1,'\n %6.5f',W1(i));
end
for i=1:length(B1)
fprintf(fidB1,'\n %6.5f',B1(i));
end
for i=1:length(W2)
fprintf(fidW2,'\n %6.5f',W2(i));
end
for i=1:length(B2)
fprintf(fidB2,'\n %6.5f',B2(i));
end
fclose(fidW1);fclose(fidB1);fclose(fidW2);fclose(fidB2);
%子函数PSOTrain.m
function [NewW1,NewB1,NewW2,NewB2]=PSOTrain(SamIn,SamOut,HiddenUnitNum);
%PSO算法参数取值
Maxgeneration=700;%最大迭代次数
E0=0.0001;%最小允许误差
%粒子飞翔幅度
Xmin=-10;
Xmax=10;
%粒子速度范围
Vmin=-5;
Vmax=5;
M=100;%粒子规模
c1=2.7;%加速常数
c2=1.3;%加速常数
w=0.9;%惯性因子
[R,SamNum]=size(SamIn);
[S2,SamNum]=size(SamOut);
generation=1;
Done=0;
Pb1=zeros(HiddenUnitNum,R+S2+1,M);
Pb2=zeros(S2,M);
Pg1=zeros(HiddenUnitNum,R+S2+1);
Pg2=zeros(S2,1);
E=zeros(size(SamOut));
rand('state',sum(100*clock));
startP1=rand(HiddenUnitNum,R+S2+1,M)-0.5;
startP2=rand(S2,M)-0.5;
startV1=rand(HiddenUnitNum,R+S2+1,M)-0.5;
startV2=rand(S2,M)-0.5;
endP1=zeros(HiddenUnitNum,R+S2+1,M);
endP2=zeros(S2,M);
endV1=zeros(HiddenUnitNum,R+S2+1,M);
endV2=zeros(S2,M);
startE=zeros(1,M);
endE=zeros(1,M);
for i=1:M
W1=startP1(1:HiddenUnitNum,1:R,i);
W2=startP1(1:HiddenUnitNum,R+1:R+S2,i);
B1=startP1(1:HiddenUnitNum,R+S2+1,i);
B2=startP2(1:S2,i);
for q=1:SamNum
TempOut=logsig(W1*SamIn(:,q)+B1);
NetworkOut(1,q)=W2'*TempOut+B2;
end
E=NetworkOut-SamOut;
startE(1,i)=sumsqr(E)/(SamNum*S2);
Pb1(:,:,i)=startP1(:,:,i);
Pb2(:,i)=startP2(:,i);
end
[val,position]=min(startE(1,:));
Pg1=startP1(:,:,position);
Pg2=startP2(:,position);
Pgvalue=val;
Pgvalue_last=Pgvalue;
while(~Done)
for num=1:M
endV1(:,:,num)=w*startV1(:,:,num)+c1*rand*(Pb1(:,:,num)-startP1(:,:,num))+c2*rand*(Pg1-startP1(:,:,num));
endV2(:,num)=w*startV2(:,num)+c1*rand*(Pb2(:,num)-startP2(:,num))+c2*rand*(Pg2-startP2(:,num));
for i=1:HiddenUnitNum
for j=1:(R+S2+1)
endV1(i,j,num)=endV1(i,j,num);
if endV1(i,j,num)>Vmax
endV1(i,j,num)=Vmax;
elseif endV1(i,j,num)<Vmin
endV1(i,j,num)=Vmin;
end
end
end
for s2=1:S2
endV2(s2,num)=endV2(s2,num);
if endV2(s2,num)>Vmax
endV2(s2,num)=Vmax;
elseif endV2(s2,num)<Vmin
endV2(s2,num)=Vmin;
end
end
endP1(:,:,num)=startP1(:,:,num)+endV1(:,:,num);
endP2(:,num)=startP2(:,num)+endV2(:,num);
for i=1:HiddenUnitNum
for j=1:(R+S2+1)
if endP1(i,j,num)>Xmax
endP1(i,j,num)=Xmax;
elseif endP1(i,j,num)<Xmin
endP1(i,j,num)=Xmin;
end
end
end
for s2=1:S2
if endP2(s2,num)>Xmax
endP2(s2,num)=Xmax;
elseif endP2(s2,num)<Xmin
endP2(s2,num)=Xmin;
end
end
W1=endP1(1:HiddenUnitNum,1:R,num);
W2=endP1(1:HiddenUnitNum,R+1:R+S2,num);
B1=endP1(1:HiddenUnitNum,R+S2+1,num);
B2=endP2(1:S2,num);
for q=1:SamNum
TempOut=logsig(W1*SamIn(:,q)+B1);
NetworkOut(1,q)=W2'*TempOut+B2;
end
E=NetworkOut-SamOut;
SSE=sumsqr(E) %便于在命令窗口观察网络误差的变化情况
endE(1,num)=sumsqr(E)/(SamNum*S2);
if endE(1,num)<startE(1,num)
Pb1(:,:,num)=endP1(:,:,num);
Pb2(:,num)=endP2(:,num);
startE(1,num)=endE(1,num);
end
end
w=0.9-(0.5/Maxgeneration)*generation;
[value,position]=min(startE(1,:));
if value<Pgvalue
Pg1=Pb1(:,:,position);
Pg2=Pb2(:,position);
Pgvalue=value;
end
if (generation>=Maxgeneration)
Done=1;
end
if Pgvalue<E0
Done=1;
end
startP1=endP1;
startP2=endP2;
startV1=endV1;
startV2=endV2;
startE=endE;
generation=generation+1;
end
W1=Pg1(1:HiddenUnitNum,1:R);
W2=Pg1(1:HiddenUnitNum,R+1:R+S2);
B1=Pg1(1:HiddenUnitNum,R+S2+1);
B2=Pg2(:,1);
NewW1=W1;
NewW2=W2;
NewB1=B1;
NewB2=B2;
%B1.txt
%链接：https://pan.baidu.com/s/1s3Umy5uw-zuQ6YP2NBwyZQ
%提取码：656i
%B2.txt
%链接：https://pan.baidu.com/s/1jmGVgPcW3_tEza1YcUyuGQ
%提取码：7p4p
%dW1.txt
%链接：https://pan.baidu.com/s/1FYtTSA4uPEudRco3v-DD7g
%提取码：4sef
%W1.txt
%链接：https://pan.baidu.com/s/1Mq0vr1RhFyZp3LGiF7jOng
%提取码：x7fs
%W2.txt
%链接：https://pan.baidu.com/s/1Lb2faQHuw2uY4r5L95m9uA
%提取码：3f00
文档：PSO参考文献.note
链接：http://note.youdao.com/noteshare?id=88aad5e95830b97a6c9763c33f149b96&sub=91744A67372A4B8C8C168B4E1E594709
通用型的随机优化方法。适合解决组合优化问题。在寻找路径的过程中表现出一种正反馈的过程。
文档：TSP蚁群算法数学模型.note
链接：http://note.youdao.com/noteshare?id=2237bb48bdba110b8439698b2243d0e9&sub=095674755E934C519852D9EB66355D22
TSP蚁群算法实现：
%链接：https://pan.baidu.com/s/16f0v5CCvzz5WVkdUQ3vFbA
%提取码：flcs
%% 第8章 蚁群算法及Matlab实现——TSP问题
% 程序8-1
%--------------------------------------------------------------------------
%% 一、数据准备
% 清空环境变量
clear all
clc
% 程序运行计时开始
t0 = clock;
%导入数据
citys=xlsread('Chap9_citys_data.xlsx', 'B2:C53');
%--------------------------------------------------------------------------
%% 二、计算城市间相互距离
n = size(citys,1);
D = zeros(n,n);
for i = 1:n
for j = 1:n
if i ~= j
D(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2));
else
D(i,j) = 1e-4; %设定的对角矩阵修正值,因为如果矩阵对角线上是0，会导致启发函数分母为0，因此设置一个足够小的数
end
end
end
%--------------------------------------------------------------------------
%% 三、初始化参数
m = 75; % 蚂蚁数量
alpha = 1; % 信息素重要程度因子
beta = 5; % 启发函数重要程度因子
vol = 0.2; % 信息素挥发(volatilization)因子
Q = 10; % 常系数
Heu_F = 1./D; % 启发函数(heuristic function)
Tau = ones(n,n); % 信息素矩阵
Table = zeros(m,n); % 路径记录表
iter = 1; % 迭代次数初值
iter_max = 100; % 最大迭代次数
Route_best = zeros(iter_max,n); % 各代最佳路径
Length_best = zeros(iter_max,1); % 各代最佳路径的长度
Length_ave = zeros(iter_max,1); % 各代路径的平均长度
Limit_iter = 0; % 程序收敛时迭代次数
%-------------------------------------------------------------------------
%% 四、迭代寻找最佳路径
while iter <= iter_max
% 随机产生各个蚂蚁的起点城市
start = zeros(m,1);
for i = 1:m
temp = randperm(n);
start(i) = temp(1);
end
Table(:,1) = start;
% 4.1构建解空间
citys_index = 1:n;
% 4.2逐个蚂蚁路径选择
for i = 1:m
% 4.2逐个城市路径选择
for j = 2:n
tabu = Table(i,1:(j - 1)); % 已访问的城市集合(禁忌表)
allow_index = ~ismember(citys_index,tabu); % 参加说明1（程序底部）
allow = citys_index(allow_index); % 待访问的城市集合
P = allow;
% 计算城市间转移概率
for k = 1:length(allow)
P(k) = Tau(tabu(end),allow(k))^alpha * Heu_F(tabu(end),allow(k))^beta;
end
P = P/sum(P);
% 轮盘赌法选择下一个访问城市
Pc = cumsum(P); %参加说明2(程序底部)
target_index = find(Pc >= rand);
target = allow(target_index(1));
Table(i,j) = target;
end
end
% 4.3计算各个蚂蚁的路径距离
Length = zeros(m,1);
for i = 1:m
Route = Table(i,:);
for j = 1:(n - 1)
Length(i) = Length(i) + D(Route(j),Route(j + 1));
end
Length(i) = Length(i) + D(Route(n),Route(1));
end
% 4.4计算最短路径距离及平均距离
if iter == 1
[min_Length,min_index] = min(Length);
Length_best(iter) = min_Length;
Length_ave(iter) = mean(Length);
Route_best(iter,:) = Table(min_index,:);
Limit_iter = 1;
else
[min_Length,min_index] = min(Length);
Length_best(iter) = min(Length_best(iter - 1),min_Length);
Length_ave(iter) = mean(Length);
if Length_best(iter) == min_Length
Route_best(iter,:) = Table(min_index,:);
Limit_iter = iter;
else
Route_best(iter,:) = Route_best((iter-1),:);
end
end
% 4.5更新信息素
Delta_Tau = zeros(n,n);
% 逐个蚂蚁计算
for i = 1:m
% 逐个城市计算
for j = 1:(n - 1)
Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i);
end
Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i);
end
Tau = (1-vol) * Tau + Delta_Tau;
% 迭代次数加1，清空路径记录表
iter = iter + 1;
Table = zeros(m,n);
end
%--------------------------------------------------------------------------
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost=etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost) '秒']);
%--------------------------------------------------------------------------
%% 绘图
figure(1)
plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],... %三点省略符为Matlab续行符
[citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-');
grid on
for i = 1:size(citys,1)
text(citys(i,1),citys(i,2),[' ' num2str(i)]);
end
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),' 起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 终点');
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')'])
figure(2)
plot(1:iter_max,Length_best,'b')
legend('最短距离')
xlabel('迭代次数')
ylabel('距离')
title('算法收敛轨迹')
%--------------------------------------------------------------------------
%% 程序解释或说明
% 1. ismember函数判断一个变量中的元素是否在另一个变量中出现，返回0-1矩阵；
% 2. cumsum函数用于求变量中累加元素的和，如A=[1, 2, 3, 4, 5], 那么cumsum(A)=[1, 3, 6, 10, 15]。
文档：蚁群算法参数设定.note
链接：http://note.youdao.com/noteshare?id=bea29b34364be899f9c02c89bed5c209&sub=D63C08501CA444EDBE369E20DD4F6437
%% 绘制蚁群算法中蚂蚁数目与最短路程、收敛迭代次数、程序执行时间关系图
%% 准备环境与数据
clc
clear all
% 输入数据
R=[0.6 0.8 1 1.2 1.4 1.6 1.8 2];
Y1=[7652 7718 7710 7760 7674 7696 7675 7681];
Y2=[64 59 67 60 48 66 55 81];
Y3=[21 27 33 39 46 52 58 66];
%% 绘图
set(gca,'linewidth',2)
% 与最短路程关系图
subplot(3,1,1);
plot(R,Y1, '-r*', 'LineWidth', 2);
set(gca,'linewidth',2); % 设置坐标轴线宽
xlabel('蚂蚁数与城市数之比');
ylabel('最短路程');
title('蚂蚁数与最短路程关系图','fontsize',12);
% legend('最短路程');
% 与收敛迭代次数关系
subplot(3,1,2);
plot(R,Y2, '--bs', 'LineWidth', 2);
set(gca,'linewidth',2) ;
xlabel('蚂蚁数与城市数之比','LineWidth', 2);
ylabel('收敛迭代次数', 'LineWidth', 2);
title('蚂蚁数与收敛迭代次数关系图','fontsize',12);
% legend('收敛迭代次数');
% 与程序执行时间关系
subplot(3,1,3);
plot(R,Y3, '-ko', 'LineWidth', 2);
set(gca,'linewidth',2) ;
xlabel('蚂蚁数与城市数之比');
ylabel('执行时间');
title('蚂蚁数与执行时间关系图','fontsize',12);
% legend('执行时间');
应用：最佳旅游方案(苏北赛2011B)
文档：苏北赛.note
链接：http://note.youdao.com/noteshare?id=868967eb3825da026c4dfcacbe77246a&sub=88E34AED702841C497376414DC549627
%链接：https://pan.baidu.com/s/1WVz_QXZkDyk5WZGOmO7Q-Q
%提取码：ix6z
%% 第8章 蚁群算法及Matlab实现——TSP问题
% 程序8-1
%--------------------------------------------------------------------------
%% 数据准备
% 清空环境变量
clear all
clc
% 程序运行计时开始
t0 = clock;
%导入数据
citys=xlsread('Ch9_spots_data.xlsx', 'B2:L12');
%--------------------------------------------------------------------------
%% 计算城市间相互距离
n = size(citys,1);
D = zeros(n,n);
for i = 1:n
for j = 1:n
if i ~= j
D(i,j) = citys(i,j);
else
D(i,j) = 1e-4; %设定的对角矩阵修正值
end
end
end
%--------------------------------------------------------------------------
%% 初始化参数
m = 15; % 蚂蚁数量
alpha = 1; % 信息素重要程度因子
beta = 5; % 启发函数重要程度因子
vol = 0.2; % 信息素挥发(volatilization)因子
Q = 10; % 常系数
Heu_F = 1./D; % 启发函数(heuristic function)
Tau = ones(n,n); % 信息素矩阵
Table = zeros(m,n); % 路径记录表
iter = 1; % 迭代次数初值
iter_max = 100; % 最大迭代次数
Route_best = zeros(iter_max,n); % 各代最佳路径
Length_best = zeros(iter_max,1); % 各代最佳路径的长度
Length_ave = zeros(iter_max,1); % 各代路径的平均长度
Limit_iter = 0; % 程序收敛时迭代次数
%-------------------------------------------------------------------------
%% 迭代寻找最佳路径
while iter <= iter_max
% 随机产生各个蚂蚁的起点城市
start = zeros(m,1);
for i = 1:m
temp = randperm(n);
start(i) = temp(1);
end
Table(:,1) = start;
% 构建解空间
citys_index = 1:n;
% 逐个蚂蚁路径选择
for i = 1:m
% 逐个城市路径选择
for j = 2:n
tabu = Table(i,1:(j - 1)); % 已访问的城市集合(禁忌表)
allow_index = ~ismember(citys_index,tabu); % 参加说明1（程序底部）
allow = citys_index(allow_index); % 待访问的城市集合
P = allow;
% 计算城市间转移概率
for k = 1:length(allow)
P(k) = Tau(tabu(end),allow(k))^alpha * Heu_F(tabu(end),allow(k))^beta;
end
P = P/sum(P);
% 轮盘赌法选择下一个访问城市
Pc = cumsum(P); %参加说明2(程序底部)
target_index = find(Pc >= rand);
target = allow(target_index(1));
Table(i,j) = target;
end
end
% 计算各个蚂蚁的路径距离
Length = zeros(m,1);
for i = 1:m
Route = Table(i,:);
for j = 1:(n - 1)
Length(i) = Length(i) + D(Route(j),Route(j + 1));
end
Length(i) = Length(i) + D(Route(n),Route(1));
end
% 计算最短路径距离及平均距离
if iter == 1
[min_Length,min_index] = min(Length);
Length_best(iter) = min_Length;
Length_ave(iter) = mean(Length);
Route_best(iter,:) = Table(min_index,:);
Limit_iter = 1;
else
[min_Length,min_index] = min(Length);
Length_best(iter) = min(Length_best(iter - 1),min_Length);
Length_ave(iter) = mean(Length);
if Length_best(iter) == min_Length
Route_best(iter,:) = Table(min_index,:);
Limit_iter = iter;
else
Route_best(iter,:) = Route_best((iter-1),:);
end
end
% 更新信息素
Delta_Tau = zeros(n,n);
% 逐个蚂蚁计算
for i = 1:m
% 逐个城市计算
for j = 1:(n - 1)
Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i);
end
Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i);
end
Tau = (1-vol) * Tau + Delta_Tau;
% 迭代次数加1，清空路径记录表
iter = iter + 1;
Table = zeros(m,n);
end
%--------------------------------------------------------------------------
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost=etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost) '秒']);
%--------------------------------------------------------------------------
%% 绘图
figure(1)
plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],... %三点省略符为Matlab续行符
[citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-');
grid on
for i = 1:size(citys,1)
text(citys(i,1),citys(i,2),[' ' num2str(i)]);
end
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),' 起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 终点');
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')'])
figure(2)
plot(1:iter_max,Length_best,'b')
legend('最少费用')
xlabel('迭代次数')
ylabel('距离')
title('算法收敛轨迹')
%--------------------------------------------------------------------------
%% 程序解释或说明
% 1. ismember函数判断一个变量中的元素是否在另一个变量中出现，返回0-1矩阵；
% 2. cumsum函数用于求变量中累加元素的和，如A=[1, 2, 3, 4, 5], 那么cumsum(A)=[1, 3, 6, 10, 15]。
小结：
文档：蚁群算法参考文献.note
链接：http://note.youdao.com/noteshare?id=a16825f77f6b0aa669fc0e113390a300&sub=CFC318CCAFDB484AACF9837D205EF786
文档：傅里叶变换局限性.note
链接：http://note.youdao.com/noteshare?id=51bdcdfa6af2192636f8dac46945d5ea&sub=B39D3850493C460D8546BDF5130CB8D0
小波特点：
文档：伸缩平移和小波变换.note
链接：http://note.youdao.com/noteshare?id=27c81f238e525edc4a68a696103bad56&sub=DBAE52886A7540E4A82D7D62C7B9FD9E文档：多尺度分析.note
链接：http://note.youdao.com/noteshare?id=b2772ca864c3a462ae7b0f67e219ce98&sub=48938523FFF349FA9E84B108EA9E008D文档：自适应分析.note
链接：http://note.youdao.com/noteshare?id=9dcbb2bf9ad19c72a04d1114552f3dcb&sub=61714C203ED049ADA3BD3698B62CDBEC文档：小波分析工具箱.note
链接：http://note.youdao.com/noteshare?id=a239e489be3feb1ef528b02f4abc1c76&sub=3EF121397C054B0192D527C91CC88248
文档：小波分析案例：对RGB图像进行多尺度分?..
链接：http://note.youdao.com/noteshare?id=352279e5fb3e5de4caaa92dbe1fd6ffb&sub=5F9AC4D6F616452FBA1CB753D1ACAA11
应用1：
文档：融合拓扑结构的小波神经网络.note
链接：http://note.youdao.com/noteshare?id=f0dce239bc7d47ed90c32914e37c75ad&sub=8CE34F4BBD704BCC97BB0CAC9461F8D1
function main()
clc;clear all;close all;
%用Mexihat函数作为样本输入和输出
x=0:0.03:3; %样本输入值
c=2/(sqrt(3).*pi.^(1/4));
d=1/sqrt(2);
u=x/2-1;
targ=d.*c.*exp(-u.^2/2).*(1-u.^2); % 目标函数的样本输出值
eta=0.02;aerfa=0.735; %赋予网络学习速率和动量因子初始值
%初始化输出层和隐层的连接权wjh和隐层和输出层的连接权.
%假设小波函数节点数为H,样本数为P,输出节点数为J,输入节点数为I.
H=15;P=2;I=length(x);J=length(targ);
b=rand(H,1);a=rand(H,1); %初始化小波参数
whi=rand(I,H);wjh=rand(H,J); %初始化权系数；
b1=rand(H,1);b2=rand(J,1);%阈值初始化;
p=0;
Err_NetOut=[];%保存的误差；
flag=1;count=0;
while flag>0
flag=0;
count=count+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xhp1=0;
for h=1:H
for i=1:I
xhp1=xhp1+whi(i,h)*x(i);
end
ixhp(h)=xhp1+b1(h);
xhp1=0;
end
for h=1:H
oxhp(h)=fai((ixhp(h)-b(h))/a(h));
end
ixjp1=0;
for j=1:J
for h=1:H
ixjp1=ixjp1+wjh(h,j)*oxhp(h);
end
ixjp(j)=ixjp1+b2(j);
ixjp1=0;
end
for i=1:J
oxjp(i)=fnn(ixjp(i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wuchayy=1/2*sumsqr(oxjp-targ);
Err_NetOut=[Err_NetOut wuchayy];%保存每次的误差;
%求解小波网络运用BP算法，各参数每次学习的调整量
for j=1:J
detaj(j)=-(oxjp(j)-targ(j))*oxjp(j)*(1-oxjp(j));
end
for j=1:J
for h=1:H
detawjh(h,j)=eta*detaj(j)*oxhp(h);
end
end
detab2=eta*detaj;
sum=0;
for h=1:H
for j=1:J
sum=detaj(j)*wjh(h,j)*diffai((ixhp(h)-b(h))/a(h))/a(h)+sum;
end
detah(h)=sum;
sum=0;
end
for h=1:H
for i=1:I
detawhi(i,h)=eta*detah(h)*x(i);
end
end
detab1=eta*detah;
detab=-eta*detah;
for h=1:H
detaa(h)=-eta*detah(h)*((ixhp(h)-b(h))/a(h));
end
%引入动量因子aerfa，加快收敛速度和阻碍陷入局部极小值.
wjh=wjh+(1+aerfa)*detawjh;
whi=whi+(1+aerfa)*detawhi;
a=a+(1+aerfa)*detaa';
b=b+(1+aerfa)*detab';
b1=b1+(1+aerfa)*detab1';
b2=b2+(1+aerfa)*detab2';
%本算法采用的是样本逐个处理而不是数据批处理
p=p+1;
if p~=P
flag=flag+1;
else
if Err_NetOut(end)>0.008
flag=flag+1;
else
figure;
plot(Err_NetOut);
xlabel('网络学习的次数');ylabel('网络输出的误差');
title('网络学习误差曲线','fontsize',20,'color',[0 1 1],'fontname','隶书');
end
end
if count>6000
figure(1);
subplot(1,2,1)
plot(Err_NetOut,'color','b','linestyle','-','linewidth',2.2,...
'marker','^','markersize',3.5);
xlabel('网络学习的次数');ylabel('网络输出的误差');
title('误差曲线','fontsize',20,'color',[1 1 1],'fontname','隶书');
subplot(1,2,2)
handle1=plot(x,targ,'color','r','linestyle','--','linewidth',2.2,...
'marker','p','markersize',3.5);
hold on
handle1=plot(x,oxjp,'color','g','linestyle','-.','linewidth',2.2,...
'marker','d','markersize',3.5);
xlabel('样本输入值');ylabel('样本目标值与网络输出值');
title('目标值与网络输出值比较','fontsize',20,'color',[1 1 1],'fontname','隶书');
legend('样本目标值','网络仿真值');
break;
end
end
function y3=diffai(x) %子程序
y3=-1.75*sin(1.75*x).*exp(-x.^2/2)-cos(1.75*x).*exp(-x.^2/2).*x;
function yl=fai(x) %子程序
yl=cos(1.75.*x).*exp(-x.^2/2);
function y2=fnn(x) %子程序
y2=1/(1+exp(-x));
结果分析：
应用2：血管重建引出的图像数字水印
文档：血管重建引出的图像数字水印.note
链接：http://note.youdao.com/noteshare?id=1d9733bcfe436a9749966d59a01a9011&sub=E48C9FDAFB9641AB9B9ED6C02C331310
%链接：https://pan.baidu.com/s/1KK_qggS84sxy7KgVpOZ30A
%提取码：xoat
%链接：https://pan.baidu.com/s/1Q_sQHwRO17wrU-FNLcXi6Q
%提取码：gsln
%此程序实现数字图像的水印嵌入
function dwtgl
clc;clear all;close all;
disp('计算机正在准备输出数字水印图像，请耐心等待……');
input=imread('image.jpg');%从D盘里读出原始图像
input_add=input;
water=imread('water1.jpg');%读出水印图像
water_add=water;
%将原始图像和水印的RGB三色分离
input=double(input);water=double(water);inputr=input(:,:,1);waterr=water(:,:,1);
inputg=input(:,:,2);waterg=water(:,:,2);inputb=input(:,:,3);waterb=water(:,:,3);
%系数r大，增加鲁棒性，不易随几何变形和物理变形失去水印,但是水印的透明性较差
r=0.06;
[Cwr,Swr]=wavedec2(waterr,1,'haar');%水印R分解
[Cr,Sr]=wavedec2(inputr,2,'haar');%图像R分解
%水印的嵌入
Cr(1:size(Cwr,2)/16)=Cr(1:size(Cwr,2)/16)+r*Cwr(1:size(Cwr,2)/16);
k=0;
while k<=size(Cr,2)/size(Cwr,2)-1
Cr(1+size(Cr,2)/4+k*size(Cwr,2)/4:size(Cr,2)/4+...
(k+1)*size(Cwr,2)/4)=Cr(1+size(Cr,2)/4+...
k*size(Cwr,2)/4:size(Cr,2)/4+(k+1)*size(Cwr,2)/4)+...
r*Cwr(1+size(Cwr,2)/4:size(Cwr,2)/2);
Cr(1+size(Cr,2)/2+k*size(Cwr,2)/4:size(Cr,2)/2+...
(k+1)*size(Cwr,2)/4)=Cr(1+size(Cr,2)/2+...
k*size(Cwr,2)/4:size(Cr,2)/2+(k+1)*size(Cwr,2)/4)+...
r*Cwr(1+size(Cwr,2)/2:3*size(Cwr,2)/4);
Cr(1+3*size(Cr,2)/4+k*size(Cwr,2)/4:3*size(Cr,2)/4+...
(k+1)*size(Cwr,2)/4)=Cr(1+3*size(Cr,2)/4+...
k*size(Cwr,2)/4:3*size(Cr,2)/4+(k+1)*size(Cwr,2)/4)+...
r*Cwr(1+3*size(Cwr,2)/4:size(Cwr,2));
k=k+1;
end;
Cr(1:size(Cwr,2)/4)=Cr(1:size(Cwr,2)/4)+r*Cwr(1:size(Cwr,2)/4);
g=0.03;
[Cwg,Swg]=wavedec2(waterg,1,'haar');%水印G分解
[Cg,Sg]=wavedec2(inputg,2,'haar');%图像G分解
%水印的嵌入
Cg(1:size(Cwg,2)/16)= Cg(1:size(Cwg,2)/16)+r*Cwg(1:size(Cwg,2)/16);
k=0;
while k<=size(Cg,2)/size(Cwg,2)-1
Cg(1+size(Cg,2)/4+k*size(Cwg,2)/4:size(Cg,2)/4+...
(k+1)*size(Cwg,2)/4)=Cg(1+size(Cg,2)/4+...
k*size(Cwg,2)/4:size(Cg,2)/4+(k+1)*size(Cwg,2)/4)+...
r*Cwg(1+size(Cwg,2)/4:size(Cwg,2)/2);
Cg(1+size(Cg,2)/2+k*size(Cwg,2)/4:size(Cg,2)/2+...
(k+1)*size(Cwg,2)/4)=Cg(1+size(Cg,2)/2+...
k*size(Cwg,2)/4:size(Cg,2)/2+(k+1)*size(Cwg,2)/4)+...
r*Cwg(1+size(Cwg,2)/2:3*size(Cwg,2)/4);
Cg(1+3*size(Cg,2)/4+k*size(Cwg,2)/4:3*size(Cg,2)/4+...
(k+1)*size(Cwg,2)/4)=Cg(1+3*size(Cg,2)/4+...
k*size(Cwg,2)/4:3*size(Cg,2)/4+(k+1)*size(Cwg,2)/4)+...
r*Cwg(1+3*size(Cwg,2)/4:size(Cwg,2));
k=k+1;
end;
Cg(1:size(Cwg,2)/4)=Cg(1:size(Cwg,2)/4)+g*Cwg(1:size(Cwg,2)/4);
b=0.12;
[Cwb,Swb]=wavedec2(waterb,1,'haar');%水印B分解
[Cb,Sb]=wavedec2(inputb,2,'haar');%图像B分解
%水印的嵌入
Cb(1:size(Cwb,2)/16)=Cb(1:size(Cwb,2)/16)+r*Cwb(1:size(Cwb,2)/16);
k=0;
while k<=size(Cb,2)/size(Cwb,2)-1
Cb(1+size(Cb,2)/4+k*size(Cwb,2)/4:size(Cb,2)/4+...
(k+1)*size(Cwb,2)/4)=Cb(1+size(Cb,2)/4+...
k*size(Cwb,2)/4:size(Cb,2)/4+(k+1)*size(Cwb,2)/4)+...
r*Cwb(1+size(Cwb,2)/4:size(Cwb,2)/2);
Cb(1+size(Cb,2)/2+k*size(Cwb,2)/4:size(Cb,2)/2+...
(k+1)*size(Cwb,2)/4)=Cb(1+size(Cb,2)/2+...
k*size(Cwb,2)/4:size(Cb,2)/2+(k+1)*size(Cwb,2)/4)+...
r*Cwb(1+size(Cwb,2)/2:3*size(Cwb,2)/4);
Cb(1+3*size(Cb,2)/4+k*size(Cwb,2)/4:3*size(Cb,2)/4+...
(k+1)*size(Cwb,2)/4)=Cb(1+3*size(Cb,2)/4+...
k*size(Cwb,2)/4:3*size(Cb,2)/4+(k+1)*size(Cwb,2)/4)+...
r*Cwb(1+3*size(Cwb,2)/4:size(Cwb,2));
k=k+1;
end;
Cb(1:size(Cwb,2)/4)=Cb(1:size(Cwb,2)/4)+b*Cwb(1:size(Cwb,2)/4);
%用哈尔小波重构图像，亦可用mexh等小波重构
inputr=waverec2(Cr,Sr,'haar');inputg=waverec2(Cg,Sg,'haar');inputb=waverec2(Cb,Sb,'haar');
%三色的叠加
temp=size(inputr);pic=zeros(temp(1),temp(2),3);
for i=1:temp(1);
for j=1:temp(2);
pic(i,j,1)=inputr(i,j);
pic(i,j,2)=inputg(i,j);
pic(i,j,3)=inputb(i,j);
end
end
output=uint8(round(pic));
watermarked_image_uint8=uint8(output);%转化为uint8
imwrite(watermarked_image_uint8,'watermarked.jpg','jpg');
%画出水印图像
figure(2);subplot(3,1,1);image(input_add);
t1=title('原始图像');
set(t1,'fontsize',15,'color',[0 0 0],'fontname','隶书');
subplot(3,1,2);image(water_add);t2=title('水印图像');
set(t2,'fontsize',15,'color',[0 0 0],'fontname','隶书');
subplot(3,1,3);image(watermarked_image_uint8);t3=title('嵌入水印后的图像');
set(t3,'fontsize',15,'color',[0 0 0],'fontname','隶书');
文档：小波分析参考文献.note
链接：http://note.youdao.com/noteshare?id=7a8e0479b63702a82558d27dbe62a337&sub=93566E7B513D4E4593B5E422282F690E
数值模拟：一般的递推方程组和微分方程组(1996B洗衣机节水问题、2003A SARS传播问题)，用于难以求出解析解、理论确定但是初始状态随机的模型。
文档：1996B.note
链接：http://note.youdao.com/noteshare?id=3343198f6fa0f00ab88acdb2bf521f72&sub=D2344B960C7A44EE8043FB2D8A83357A
1.1 微分方程组模拟
文档：微分方程组模拟题目.note
链接：http://note.youdao.com/noteshare?id=045acec2e4cd11b952c61ae0070fe5be&sub=E62D408DA04C463ABC33AF25522AF01E
%微分方程组求解主程序
clc;clear all;clf;close all;
%Windows时钟自动计时
T1=clock;%Clock函数返回的值是clock = [year month day hour minute seconds]
disp('计算机正在准备输出湖泊有机物新陈代谢结果，请耐心等待……');
[tt,y]=ode45('lbwfun',[0:10:2020],[95.9,0.83,0.003,0.0001,0.0,0.0]);
t=tt(191:end,:)
ys=y(191:end,1)
yp=y(191:end,2)
yh=y(191:end,3)
yr=y(191:end,4)
yo=y(191:end,5)
ye=y(191:end,6)
T2=clock;
API_elapsed_time=T2-T1;
if API_elapsed_time(6)<0
API_elapsed_time(6)=API_elapsed_time(6)+60;
API_elapsed_time(5)=API_elapsed_time(5)-1;
end
if API_elapsed_time(5)<0
API_elapsed_time(5)=API_elapsed_time(5)+60;
API_elapsed_time(4)=API_elapsed_time(4)-1;
end
if API_elapsed_time(4)<0
API_elapsed_time(4)=API_elapsed_time(4)+60;
API_elapsed_time(3)=API_elapsed_time(4)-1;
end
str=sprintf('湖泊新陈代谢模拟程序共运行 %d 小时 %d 分钟 %.4f 秒',...
API_elapsed_time(4),API_elapsed_time(5),API_elapsed_time(6));
disp(str);
%子函数：微分方程组odefile文件并命名为lbwfun.m
function ydot=lbwfun(t,y);
ydot=[121.793*pi*cos(2*pi*t);
y(1)-4.03*y(2);
0.48*y(2)-17.87*y(3);
4.85*y(3)-4.65*y(4);
2.55*y(2)+6.12*y(3)+1.95*y(4);
1.10*y(2)+6.90*y(3)+2.70*y(4)];
1.2 服从概率分布的随机模拟
文档：服从概率分布的随机模拟.note
链接：http://note.youdao.com/noteshare?id=ea3f5d09b0938bd7a51ad0db01f901e1&sub=BD617DEC79D74DAA99957B6C5FE5B595
clc;close all;clear all;
disp('计算机大概需要两分钟的运行时间，请等待……');
p=zeros(10,100);avert=zeros(10,100);
%分别在顾客人数为10、100、500等情况时，模拟系统工作强度和顾客平均逗留时间
nn=[10 100 500 1000 5000 10000 20000 50000 100000 500000];
for d=1:10 %length(nn)=10
for s=1:100 %每种情况重复模拟100次以便消除随机因素
n=nn(d);%模拟顾客数目
dt=exprnd(10,1,n);%到达时间间隔
st=normrnd(6.5,1.2,1,n);%服务台服务时间
a=zeros(1,n);%每个顾客到达时刻
b=zeros(1,n);%每个顾客开始接受服务时刻
c=zeros(1,n);%每个顾客离开时刻
a(1)=0;
for i=2:n
a(i)=a(i-1)+dt(i-1);%第i顾客到达时刻
end
b(1)=0;%第一个顾客开始接受服务的时刻是其到达的时刻
c(1)=b(1)+st(1);%第一个顾客的离开时刻为其接受的服务时间加上开始接受服务的时刻
for i=2:n
%如果第i个顾客到达时间比前一个顾客离开时间早，则接受服务时间为前一人离开时间
if(a(i)<=c(i-1))b(i)=c(i-1);
%如果第i个顾客到达时间比前一个顾客离开时间晚，则接受服务时间为其到达时间
else b(i)=a(i);
end
%第i个顾客离开时间为其开始接受服务的时刻+接受服务的时间长度
c(i)=b(i)+st(i);
end
cost=zeros(1,n);%记录每个顾客在系统逗留时间
for i=1:n
cost(i)=c(i)-a(i);%第i个顾客在系统逗留时间
end
T=c(n);%总时间
p(d,s)=sum(st)/T;
avert(d,s)=sum(cost)/n;
end
end
pc=sum(p')/100;avertc=sum(avert')/100;
figure(1);subplot(2,1,1);%分区画图
plot(pc,'color','g','linestyle','-','linewidth',2.5,'marker','*','markersize',5);
text(1,pc(:,1),texlabel('10人'),'fontsize',11);
text(2,pc(:,2),texlabel('10^2 人'),'fontsize',11);
text(3-0.15,pc(:,3)-0.004,texlabel('5x10^2 人'),'fontsize',11);
text(4-0.15,pc(:,4)-0.004,texlabel('10^3 人'),'fontsize',11);
text(5-0.15,pc(:,5)-0.004,texlabel('5x10^3 人'),'fontsize',11);
text(6-0.15,pc(:,6)-0.004,texlabel('10^4 人'),'fontsize',11);
text(7-0.15,pc(:,7)-0.004,texlabel('2x10^4 人'),'fontsize',11);
text(8-0.15,pc(:,8)-0.004,texlabel('5x10^4 人'),'fontsize',11);
text(9-0.15,pc(:,9)-0.004,texlabel('10^5 人'),'fontsize',11);
text(10-0.15,pc(:,10)-0.004,texlabel('5x10^5 人'),'fontsize',11);
xlim([1 11]);xlabel('顾客数量/个','fontsize',11);ylabel('系统工作强度','fontsize',11);
subplot(2,1,2);
plot(avertc,'color','r','linestyle','-','linewidth',2.5,'marker','s','markersize',5);
text(1,avertc(:,1)+0.4,texlabel('10人'),'fontsize',11);
text(2-0.15,avertc(:,2)-0.4,texlabel('10^2 人'),'fontsize',11);
text(3-0.15,avertc(:,3)-0.4,texlabel('5x10^2 人'),'fontsize',11);
text(4-0.15,avertc(:,4)-0.4,texlabel('10^3 人'),'fontsize',11);
text(5-0.15,avertc(:,5)-0.4,texlabel('5x10^3 人'),'fontsize',11);
text(6-0.15,avertc(:,6)-0.4,texlabel('10^4 人'),'fontsize',11);
text(7-0.15,avertc(:,7)-0.4,texlabel('2x10^4 人'),'fontsize',11);
text(8-0.15,avertc(:,8)-0.4,texlabel('5x10^4 人'),'fontsize',11);
text(9-0.15,avertc(:,9)-0.4,texlabel('10^5 人'),'fontsize',11);
text(10-0.15,avertc(:,10)-0.4,texlabel('5x10^5 人'),'fontsize',11);
xlim([1 11]);xlabel('顾客数量/个','fontsize',11);ylabel('顾客逗留时间/秒','fontsize',11);
1.3 蒙特卡洛模拟
文档：蒙特卡洛.note
链接：http://note.youdao.com/noteshare?id=d215bc6819827ea36811c271d072bf91&sub=C72991CCBA0645959020C7567E7CF662
%P11-3.求π
%本程序画出在正方体内嵌套1/4圆的组合图形
%程序虽小，但是蕴含着用结构体设置图形属性和图形对象及坐标轴属性底层设置
clc;
x=0:0.01:1;
y=sqrt(1-x.^2);
figure;
h=plot(x,y);
active.linestyle='-';
active.linewidth=3;
active.color='k';
set(h,active);
set(gcf,'color','y');
set(gca,'color','b');
xlabel('x');
ylabel('y');
grid on;
axis square;
%P11-4.
%本例没有使用循环语句而是将相关量矢量化
%程序运行前在“file”菜单下“preferences”将输出格式调成“long”形式
clc;
clear all;
n=1000000;
a=rand(n,1);
%必须用rand指令而不能用randn指令因为产生的随机数必须是均匀的
b=rand(n,1);
c=find(a.^2+b.^2<=1);
d=length(c);
pi_value=d/n*4
蒙特卡洛求解多约束非线性规划
m a x f ( x ) = x 1 x 2 x 3 maxf(x)=x_1x_2x_3 maxf(x)=x1x2x3
s . t . { − x 1 + 2 x 2 + 2 x 3 ≥ 0 x 1 + 2 x 2 + 2 x 3 ≤ 72 10 ≤ x 2 ≤ 20 x 1 − x 2 = 10 s.t.\begin{cases} -x_1+2x_2+2x_3≥0\\ x_1+2x_2+2x_3≤72\\ 10≤x_2≤20\\ x_1-x_2=10 \end{cases} s.t.⎩⎪⎪⎪⎨⎪⎪⎪⎧−x1+2x2+2x3≥0x1+2x2+2x3≤7210≤x2≤20x1−x2=10
可知20≤x1≤30,10≤x2≤20,-10≤x3≤16
文档：非线性规划出现inf情况的改法.note
链接：http://note.youdao.com/noteshare?id=b6d2584938622653c5f20c519a985a83&sub=F7607168B3CE4B71B8A745B8E18FF384
%wrong:
%程序运行前在“file”菜单下“preferences”将输出格式调成“long”形式
clc;clear all;
T1=clock;
N=1000;
x10=[];x20=[];x30=[];vmax=-inf;
x1=unifrnd(20,30,N,1);
x2=unifrnd(10,20,N,1);
x3=unifrnd(-10,16,N,1);
for i=1:N
for j=1:N
for k=1:N
if -x1(i)+2*x2(j)+2*x3(k)>=0&...
x1(i)+2*x2(j)+2*x3(k)<=72&...
x1(i)-x2(j)<=0.1,%x1(i)-x2(j)==10,
v=x1(i)*x2(j)*x3(k);
if v>vmax,
vmax=v;x10=x1(i);x20=x2(j);x30=x3(k);
end
end
end
end
end
x=[x10,x20,x30],vmax
T2=clock;
API_elapsed_time=T2-T1;
if API_elapsed_time(6)<0
API_elapsed_time(6)=API_elapsed_time(6)+60;
API_elapsed_time(5)=API_elapsed_time(5)-1;
end
if API_elapsed_time(5)<0
API_elapsed_time(5)=API_elapsed_time(5)+60;
API_elapsed_time(4)=API_elapsed_time(4)-1;
end
if API_elapsed_time(4)<0
API_elapsed_time(4)=API_elapsed_time(4)+60;
API_elapsed_time(3)=API_elapsed_time(4)-1;
end
str=sprintf('MC程序共运行 %d 小时 %d 分钟 %.4f 秒',...
API_elapsed_time(4),API_elapsed_time(5),API_elapsed_time(6));
disp(str);
动画仿真：预测的功能。
图像静态模拟：将生成的伪随机数作为初始值进行模拟，以figure形式输出。(2002A wind and waterspray 2002A 车灯光源的优化设计问题)
文档：matlab音频处理.note
链接：http://note.youdao.com/noteshare?id=4b2f348ba25a7d1bd0c3a4217383a939&sub=38AD5DFB88D34608B352552789BF8263文档：常规动画实现.note
链接：http://note.youdao.com/noteshare?id=8fecfc960fc246ebe1ec104a75ccba3d&sub=1FB56FE9E48C42558C85489A93B4CA41
《功夫熊猫》动画：
%链接：https://pan.baidu.com/s/1axGg5KG5MCslqAgrUdtuBw
%提取码：dad3
clc;close all;clear all;figure(1);
%制作电影动画
m=moviein(28); %分配存储帧的内存
a0=int2str(0);
for i=1:9 % 依次读入1-9张数字图像
a=int2str(i);
a1=strcat('KonghuPanda\KonghuPanda',a0,a0,a,'.jpg');
X1=imread(a1);imshow(X1);m(i)=getframe;
end
for i=10:28 % 依次读入10-28张数字图像
a2=int2str(i);
a3=strcat('KonghuPanda\KonghuPanda',a0,a2,'.jpg');
X2=imread(a3);imshow(X2);m(i)=getframe;
end
movie(m,8); %播放动画
%输出各个电影帧
figure(2);
for i=1:9
subplot(7,4,i);
a=int2str(i);
a1=strcat('KonghuPanda\KonghuPanda',a0,a0,a,'.jpg');
X1=imread(a1);
image(X1);axis off;box off;
a4=int2str(i);
a5=strcat('第',a4,'帧');
title(a5,'fontsize',9,'fontname','隶书');
end
for i=10:28
subplot(7,4,i);
a2=int2str(i);
a3=strcat('KonghuPanda\KonghuPanda',a0,a2,'.jpg');
X2=imread(a3);
image(X2);axis off;box off;
a6=int2str(i);
a7=strcat('第',a6,'帧');
title(a7,'fontsize',9,'fontname','隶书');
end
文档：实时动画.note
链接：http://note.youdao.com/noteshare?id=3fcd0c92e90693ea96e087a279408c62&sub=9FA832034B4845FE80612EF703C3B6A7文档：avi视频、照相机函数、C语言动画.note
链接：http://note.youdao.com/noteshare?id=974f99409238ea81a1ba17cf5ef41e01&sub=85BA89EBEAFB4CD58311FE48608D2D58
应用：四维水质模型(2005A 长江水质的评价和预测——一维水质模型)
文档：2005A问题分析.note
链接：http://note.youdao.com/noteshare?id=6c030504733fe2197a73a72f91493174&sub=A9BAB7CD0300482ABBAD00E70737D90C
四维水质模型广泛应用于各类复杂的水污染和大气污染。
文档：四维水质模型.note
链接：http://note.youdao.com/noteshare?id=b20d822b10aba8581ff64f377ca9e880&sub=57872391D3724339A1B2CA94A20342A4
%问题一：需要用到污染点源四维扩散膜模型
%这是主程序,命名为Main1206.m
clc;clear all;close all;
%本程序研究的对象是瞬时点源污染扩散问题
subplot(1,2,1);
daspect([1 1 1]);
xmin=-500;dx=10;xmax=500;ymin=-100;dy=5;ymax=100;zmin=-10;dz=1;zmax=10;
[cx,cy,ct]=meshgrid(xmin:dx:xmax,ymin:dy:ymax,10:10:50);
[grad_cx,grad_cy,grad_cz,grad_ct]=grad(cx,cy,0,ct);
hcone=coneplot(-500:10:500,-100:5:100,10:10:50,grad_cx,grad_cy,grad_ct,cx,cy,ct,0.9);
xlim([-60 70]);ylim([-50 25]);zlim([0 150]);
grid on;camlight right;lighting phong;camproj perspective;
xlabel('横向/(mg/L*s)','fontsize',18);ylabel('纵向/(mg/L*s)','fontsize',18);
zlabel('时间/s','fontsize',18);view(3);
title('10秒至50秒XOY切面扩散梯度流锥图','fontsize',18);
subplot(1,2,2);
[x,y,z]=meshgrid(xmin:dx:xmax,ymin:dy:ymax,zmin:dz:zmax);
[grad_x,grad_y,grad_z,grad_t]=grad(x,y,z,20);
quiver3(x,y,z,grad_x,grad_y,grad_z,1.1,'r-');
xlabel('横向/(mg/L*s)','fontsize',18);ylabel('纵向/(mg/L*s)','fontsize',18);
zlabel('河流竖直方向(mg/L*s)','fontsize',18);
axis tight;
title('时间T=20秒时3D扩散梯度向量图','fontsize',18);
% 这是子程序，命名为grad.m
function [grad_x,grad_y,grad_z,grad_t]=grad(x,y,z,t)
%以下各梯度grad_x,grad_y,grad_z,grad_t的表达式可以由此命令用计算机求出：
%syms t M Dx Dy Dz x ux y uy z uz K
%grad_maple=maple('grad(M/(8*(pi*t)^(3*2)*sqrt(Dx*Dy*Dz))*exp(-(x-ux*t)^2/(
%4*Dx*t)-(y-uy*t)^2/(4*Dy*t)-(z-uz*t)^2/(4*Dz*t))*exp(-K*t),vector([x,y,z,t]));')
%先用simplify和collect命令对各变量梯度表达式进行化简和合并
M=200400;K=4.2/(24*60*6);
Dx=50;Dy=5;Dz=2;ux=1.5;uy=0.2;uz=0.1;
grad_x=-1./16.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*(x-ux.*t)./Dx.*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t);
grad_y=-1./16.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*(y-uy.*t)./Dy.*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t);
grad_z=-1./16.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*(z-uz.*t)./Dz.*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t);
grad_t=-3./4.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t)+1./8.*M./pi.^6./t.^6./(Dx.*Dy.*Dz).^(1./2).*(1./2.*(x-ux.*t)./Dx./t.*ux+1./4.*(x-ux.*t).^2./Dx./t.^2+1./2.*(y-uy.*t)./Dy./t.*uy+1./4.*(y-uy.*t).^2./Dy./t.^2+1./2.*(z-uz.*t)./Dz./t.*uz+1./4.*(z-uz.*t).^2./Dz./t.^2).*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t)-1./8.*M./pi.^6./t.^6./(Dx.*Dy.*Dz).^(1./2).*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*K.*exp(-K.*t);
%问题二：
%这是主程序,命名为lianxuwuran.m
%分别输入20 180
function Value=lianxuwuran(t)
%现在开始编写连续污染源模型
global xx yy zz;
xmin=-500;dx=10;xmax=500;ymin=-100;dy=5;ymax=100;zmin=-10;dz=1;zmax=10;
Cxyz_t((xmax-xmin)/dx+1,(ymax-ymin)/dy+1,(zmax-zmin)/dz+1)=0;
ii=0;jj=0;kk=0;
for zz=zmin:dz:zmax;kk=kk+1;
for yy=ymin:dy:ymax;jj=jj+1;
for xx=xmin:dx:xmax;ii=ii+1;
Cxyz_t(ii,jj,kk)=quadl(@fun3D,1,t);
end
ii=0;
end
jj=0;
end
Cxyz_tt2=Cxyz_t;
%现在开始编写瞬时点源污染模型
M=200400;K=4.2/(24*60*6);
Dx=50;Dy=5;Dz=2;ux=1.5;uy=0.2;uz=0.1;
[x,y,z]=meshgrid(xmin:dx:xmax,ymin:dy:ymax,zmin:dz:zmax);
C_Point=M./(8.*(pi.*t).^(3./2).*sqrt(Dx.*Dy.*Dz)).*exp(-(x-ux.*t).^2./(4.*Dx.*...
t)-(y-uy.*t).^2./(4.*Dy.*t)-(z-uz.*t).^2./(4.*Dz.*t)).*exp(-K.*t);
if (t<=30&t>0)
figure(1);
subplot(4,1,1);[c,h]=contour(Cxyz_tt2(:,:,11),200);grid on;
axis tight;title('连续污染源中心扩散示意图(a)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
subplot(4,1,2);[c,h]=contour(C_Point(:,:,11),200);grid on;
axis tight;title('瞬时点源中心扩散示意图(b)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
subplot(4,1,3);[c,h]=contour(Cxyz_tt2(:,:,1),200);grid on;
axis tight;title('连续污染源水下10米扩散示意图(c)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
subplot(4,1,4);[c,h]=contour(C_Point(:,:,1),200);grid on;
axis tight;title('瞬时点源水下10米扩散示意图(d)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
else t>30
figure(1);
subplot(4,1,1);[c,h]=contour(Cxyz_tt2(:,:,11),8);grid on;
clabel(c,h);axis tight;title('连续污染源中心扩散示意图(a)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
subplot(4,1,2);[c,h]=contour(C_Point(:,:,11),8);grid on;
clabel(c,h);axis tight;title('瞬时点源中心扩散示意图(b)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
subplot(4,1,3);[c,h]=contour(Cxyz_tt2(:,:,1),8);grid on;
clabel(c,h);axis tight;title('连续污染源水下10米扩散示意图(c)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
subplot(4,1,4);[c,h]=contour(C_Point(:,:,1),8);grid on;
clabel(c,h);axis tight;title('瞬时点源水下10米扩散示意图(d)');
xlabel('x方向污染物浓度/(mg/L)'); ylabel('y方向污染物浓度/(mg/L)');
end
%这是子程序，命名为fun3D.m
function CPoint=fun3D(t)
global xx yy zz;
Q=30000;Cq=6.68;K=4.2/(24*60*6);
Dx=50;Dy=5;Dz=2;ux=1.5;uy=0.2;uz=0.1;
CPoint=Cq*Q./(8.*(pi.*t).^(3./2).*sqrt(Dx.*Dy.*Dz)).*exp(-(xx-ux.*t).^2./(4.*Dx.*...
t)-(yy-uy.*t).^2./(4.*Dy.*t)-(zz-uz.*t).^2./(4.*Dz.*t)).*exp(-K.*t);
文档：四维水质模型.note
链接：http://note.youdao.com/noteshare?id=b20d822b10aba8581ff64f377ca9e880&sub=57872391D3724339A1B2CA94A20342A4文档：计算机虚拟参考文献.note
链接：http://note.youdao.com/noteshare?id=4cd28db8e362cc1239456a50762e301c&sub=14E1A9C519934ED2BD9554C3D563B9D5
文章浏览阅读1w次，点赞2次，收藏27次。来源：机器人小妹 很多时候企业拥有重复，乏味且困难的工作流程，这些流程往往会减慢生产速度并增加运营成本。为了降低生产成本，企业别无选择，只能自动化某些功能以降低生产成本。 通过数字化..._人工智能平台
文章浏览阅读2.2k次。热加载能够在每次保存修改的代码后自动刷新 electron 应用界面，而不必每次去手动操作重新运行，这极大的提升了开发效率。安装 electron 热加载插件热加载虽然很方便，但是不是每个 electron 项目必须的，所以想要舒服的开发 electron 就只能给 electron 项目单独的安装热加载插件[electron-reloader]:// 在项目的根目录下安装 electron-reloader，国内建议使用 cnpm 代替 npmnpm install electron-relo._electron-reloader
文章浏览阅读942次。在11.0 进行定制化开发，会根据需要去掉recovery模式的一些选项 就是在device.cpp去掉一些选项就可以了。_android recovery 删除 部分菜单
文章浏览阅读3.7k次。https://www.yuque.com/mnn/cn/cvrt_linux_mac基础依赖这些依赖是无关编译选项的基础编译依赖• cmake（3.10 以上）• protobuf (3.0 以上)• 指protobuf库以及protobuf编译器。版本号使用 protoc --version 打印出来。• 在某些Linux发行版上这两个包是分开发布的，需要手动安装• Ubuntu需要分别安装 libprotobuf-dev 以及 protobuf-compiler 两个包•..._mnn 编译linux
文章浏览阅读1.8k次。CSS3新增动画属性“@-webkit-keyframes”，从字面就可以看出其含义——关键帧，这与Flash中的含义一致。利用CSS3制作动画效果其原理与Flash一样，我们需要定义关键帧处的状态效果，由CSS3来驱动产生动画效果。下面讲解一下如何利用CSS3制作淡入淡出的动画效果。具体实例可参考刚进入本站时的淡入效果。1. 定义动画，名称为fadeIn@-webkit-keyf_css3入场效果淡入淡出
文章浏览阅读2.8k次。计算机系统应包括硬件和软件两个子系统，硬件和软件又必须依次分别包括中央处理器和系统软件。按人的要求接收和存储信息，自动进行数据处理和计算，并输出结果信息的机器系统。计算机是脑力的延伸和扩充，是近代科学的重大成就之一。计算机系统由硬件(子)系统和软件(子)系统组成。前者是借助电、磁、光、机械等原理构成的各种物理部件的有机组合，是系统赖以工作的实体。后者是各种程序和文件，用于指挥全系统按指定的要求进行..._计算机系统包括硬件系统和软件系统 软件又必须包括
文章浏览阅读7.9k次，点赞3次，收藏22次。一 定义这是最早出现的置换算法。该算法总是淘汰最先进入内存的页面，即选择在内存中驻留时间最久的页面予以淘汰。该算法实现简单，只需把一个进程已调入内存的页面，按先后次序链接成一个队列，并设置一个指针，称为替换指针，使它总是指向最老的页面。但该算法与进程实际运行的规律不相适应，因为在进程中，有些页面经常被访问，比如，含有全局变量、常用函数、例程等的页面，FIFO 算法并不能保证这些页面不被淘汰。这里，我_进程调度fifo算法代码
文章浏览阅读133次。rownum是oracle才有的写法，rownum在oracle中可以用于取第一条数据，或者批量写数据时限定批量写的数量等mysql取第一条数据写法SELECT * FROM t order by id LIMIT 1;oracle取第一条数据写法SELECT * FROM t where rownum =1 order by id;ok，上面是mysql和oracle取第一条数据的写法对比，不过..._mysql 替换@rownum的写法
文章浏览阅读790次，点赞3次，收藏4次。官网下载下载链接：http://www.eclipse.org/downloads/点击Download下载完成后双击运行我选择第2个，看自己需要（我选择企业级应用，如果只是单纯学习java选第一个就行）进入下一步后选择jre和安装路径修改jvm/jre的时候也可以选择本地的（点后面的文件夹进去），但是我们没有11版本的，所以还是用他的吧选择接受安装中安装过程中如果有其他界面弹出就点accept就行..._ecjelm
文章浏览阅读245次。原文链接：https://linux.cn/article-7801-1.htmlifconfigping <IP地址>:发送ICMP echo消息到某个主机traceroute <IP地址>：用于跟踪IP包的路由路由：netstat -r: 打印路由表route add ：添加静态路由路径routed：控制动态路由的BSD守护程序。运行RIP路由协议gat..._ifconfig 删除vlan
文章浏览阅读224次。reduxredux里要求把数据都放在公共的存储区域叫store里面，组件中尽量少放数据，假如绿色的组件要给很多灰色的组件传值，绿色的组件只需要改变store里面对应的数据就行了，接着灰色的组件会自动感知到store里的数据发生了改变，store只要有变化，灰色的组件就会自动从store里重新取数据，这样绿色组件的数据就很方便的传到其它灰色组件里了。redux就是把公用的数据放在公共的区域去存..._redux redis
文章浏览阅读2.2k次，点赞3次，收藏6次。unzip版本不支持4G以上的压缩包所以要使用p7zip：Linux一个高压缩率软件wget http://sourceforge.net/projects/p7zip/files/p7zip/9.20.1/p7zip_9.20.1_src_all.tar.bz2tar jxvf p7zip_9.20.1_src_all.tar.bz2cd p7zip_9.20.1make && make install 如果安装失败，看一下报错是不是因为没有下载gcc 和 gcc ++(p7_linux 7za解压中文乱码