MATLAB数学建模方法与实践（第3版）——读书笔记

1.数学建模的五类问题
2.MATLAB快速入门
- 2.1从问题到发布脚本
- 2.2matlab提高
3.数据建模技术
4.优化技术
5.连续模型的matlab求解方法
6.评价模型的求解
- 6.1线性加权法
- 6.2层次分析法
7.机理建模的matlab实现方法
- 7.1推导机理建模(略)
- 7.2元胞自动机——仿真
8.出租车补贴方案优化问题2015B
9.开放小区对道路通行影响的问题2016

1.数学建模的五类问题

数据型——>拟合、回归、分类、聚类、主成分
离散型——>目标规划、智能算法(神经网络、遗传、模拟退火、蚁群、粒子群)
连续型——>微分、偏微分、差分、函数
评价型——>层次分析法、多因子决策、回归
机理建模——>自主编程

2.MATLAB快速入门

2.1从问题到发布脚本

problem：已知股票的交易数据，包括日期、开盘价、最高价、最低价、收盘价、成交量、换手率，试评价这只股票的价值和风险。

导入数据。
数据可视化(选择比较关心的两个变量：日期、收盘价)：选中工作区的这两个变量：右击，plot。

doc ；help命令可以帮助我们

如何评估股票价格及风险？

3.1 增幅越大价格越好——曲线斜率越大——多项式拟合——得到拟合方程——进而得到斜率。

3.2风险——最大回撤
如何拟合？ployfit。最大回撤的函数？maxdrawdown

P = polyfit(DateNum,Pclose,1)%多项式拟合
value=p(1)%将斜率赋值给value，作为股票价值

MaxDD = maxdrawdown(Pclose);%计算最大回撤
risk=MaxDD%返回给risk，作为股票风险

导入数据的时候选择”生成脚本”，就会自动生成一些代码，在此基础上添加代码，形成脚本。

%链接：https://pan.baidu.com/s/1g0Et-S-RQcNuaDlirSLHkg 
%提取码：y9dy
%% 导入电子表格中的数据
% 用于从以下电子表格导入数据的脚本:
%
%    工作簿: C:\Users\LSdarker\Desktop\建模啊建模\MATLAB数学建模方法与实践（第3版）程序及数据\程序_MATLAB数学建模方法与实践_卓金武等\程序_MATLAB数学建模方法与实践_卓金武等\Cha2\sz000004.xls
%    工作表: Sheet1
%
% 要扩展代码以供其他选定数据或其他电子表格使用，请生成函数来代替脚本。

% 由 MATLAB 自动生成于 2020/09/01 19:15:05

%% 导入数据
[~, ~, raw] = xlsread('C:\Users\LSdarker\Desktop\建模啊建模\MATLAB数学建模方法与实践（第3版）程序及数据\程序_MATLAB数学建模方法与实践_卓金武等\程序_MATLAB数学建模方法与实践_卓金武等\Cha2\sz000004.xls','Sheet1','A2:H99');

%% 创建输出变量
data = reshape([raw{:}],size(raw));

%% 将导入的数组分配给列变量名称
Date1 = data(:,1);
DateNum1 = data(:,2);
Popen1 = data(:,3);
Phigh1 = data(:,4);
Plow1 = data(:,5);
Pclose1 = data(:,6);
Volum1 = data(:,7);
Turn1 = data(:,8);

%% 清除临时变量
clearvars data raw;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上是自动生成的，以下是编写的%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %% 数据探索
figure % 创建一个新的图像窗口
plot(DateNum,Pclose,'k') % 更改图的的颜色的黑色(打印后不失真)
datetick('x','mm');% 更改日期显示类型
xlabel('日期'); % x轴说明
ylabel('收盘价'); % y轴说明
figure
bar(Pclose) % 作为对照图形

%% 股票价值的评估
p = polyfit(DateNum,Pclose,1); % 多项式拟合,
% 分号作用为不在命令窗口显示执行结果
P1 = polyval(p,DateNum); % 得到多项式模型的结果
figure
plot(DateNum,P1,DateNum,Pclose,'*g'); % 模型与原始数据的对照
value = p(1) % 将斜率赋值给value， 作为股票的价值。
 
%% 股票风险的评估
MaxDD = maxdrawdown(Pclose); % 计算最大回撤
risk = MaxDD  % 将最大回撤赋值给risk， 作为股票的风险

发布——编辑发布选项——将HTML改为doc，至此我们便会得到一份详细的运行报告。
改进——这个脚本是求1只股票的风险及价值，如果求1000只，可将该脚本抽象成函数。

2.2matlab提高

Tab：补全命令
:(冒号)：用做多下标援引时，表示对应维度上的全部元素。
{}：元胞数组标记符，也就是广义矩阵，可以存放任何数据类型，每个元素可以有不同的尺寸，每个元素的内容也可以完全不同。这些元素叫做元胞(cell)。
table也是一种数据类型，展示数据和操作数据方面比元胞数组更有优势。

3.数据建模技术

3.1数据的获取、分析、预处理

3.1.1数据读取

%1.从excel读取数据
xlsread('excel所在位置',3,'B1:C5')
%3代表sheet3；B1:C5代表读入数据范围
%写入
xlsread('excel所在位置',a,3,'B1:C5')
%a代表待写入数据；3代表sheet3；B1:C5代表具体位置


%2.从TXT读取数据
load('***.txt')
%保存
save c:\***.txt a -ascii;
%a是变量 以ASCII码存放


%3.如果TXT文件存储了不同类型数据，分类读取就需要textread
[name,type,x,y,answer]=textread('**.txt','format',N,'headerlines',M)
%format 代表格式；N表示读取的次数，每次读一行；headerlines表示从M+1行开始读取。


%4.将信息写入TXT，可以控制写入精度
file=fopen('**.txt','w')
fprintf(file,'%6.2f',3.14)%6.2f表示输出宽度为6，精确到小数点2位


%5.读取图片
A=imread(filename)
%A存放像素矩阵


%6.读入并将图像拼接。此程序需要一定的算法提高拼接正确率
%% 读取图片
%链接：https://pan.baidu.com/s/1Ow6L7a-ArvFllaK4M6j4tw 
%提取码：sqyx
%链接：https://pan.baidu.com/s/1z0NPoyB96CT4FEQUJsZCfg 
%提取码：u3df
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
clc, clear, close all
a1=imread('000.bmp');
[m,n]=size(a1);
%% 批量读取图片
dirname = 'ImageChips';
files = dir(fullfile(dirname, '*.bmp'));
a=zeros(m,n,19);
pic=[];
for ii = 1:length(files)
  filename = fullfile(dirname, files(ii).name);
  a(:,:,ii)=imread(filename);
  pic=[pic,a(:,:,ii)];
end
double(pic);
figure
imshow(pic,[])


%7.读取视频
%% 读取视频数据
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
videoFReader = vision.VideoFileReader('vippedtracking.mp4');%写上自己视频的名字

% 播放视频文件
videoPlayer = vision.VideoPlayer;
while ~isDone(videoFReader)
  videoFrame = step(videoFReader);
  step(videoPlayer, videoFrame);
end
release(videoPlayer);

%% 设置播放方式
% 重置播放器
reset(videoFReader)
% 增加播放器的尺寸
r = groot;
scrPos = r.ScreenSize;
%  Size/position is always a 4-element vector: [x0 y0 dx dy]
dx = scrPos(3); dy = scrPos(4);
videoPlayer = vision.VideoPlayer('Position',[dx/8, dy/8, dx*(3/4), dy*(3/4)]);
while ~isDone(videoFReader)
  videoFrame = step(videoFReader);
  step(videoPlayer, videoFrame);
end
release(videoPlayer);
reset(videoFReader)

%% 获取视频中的图像
videoFrame = step(videoFReader);
n = 0;
while n~=15
  videoFrame = step(videoFReader);
  n = n+1;
end
figure, imshow(videoFrame)
release(videoPlayer);
release(videoFReader)

3.1.2数据预处理

缺失值处理

1.2删除法：删除小部分

1.2插补法：

a. 均值插补：定距用平均值、非定距用众数。

b. 回归插补。

c.极大似然估计(ML)：观测数据的边际分布对未知参数进行极大似然估计。也可以通过期望最大化来参数估计。有效样本的数量足够保证ML估计值是渐进无偏的并服从正态分布。

噪音过滤

noise是数据随机误差，正常的，但会影响变量真值反映。

方法有：

a. 回归法：前提是要有线性趋势。

回归法是用函数拟合数据来光滑数据，线性回归得到两个属性的最佳直线，使得通过一个预测另一个。多元线性回归的属性多于两个。通过回归后的函数值代替原始数据，避免噪声感染。

b.均值平滑法：前提是具有序列特征的变量。

用临近的若干数据均值替换原来数据。

c. 离散点分析：通过聚类的方法检测离散点，删除它。

d.小波过滤法：本质是一个函数逼近问题。

即如何在小波母函数伸缩和平移所展成的函数空间中，根据提出的衡量准则，对原信号的最佳逼近，完成原信号和噪声信号的区分。也就是从实际信号空间到小波函数空间的最佳映射，从而得到原信号的最佳恢复。

数据集成

将分散的数据源中的数据，集成到一个统一的数据集合。(针对数学建模。就是获取所需要的数据)

数据归约

属性选择：删除不相关的属性：通过数据相关性分析、数据统计分析、数据可视化、主成分分析。

样本选择。

数据变换

标准化：转化为无量纲的纯数值，便于比较和加权。

0-1标准化

，也叫离差标准化，对原始数据线性变换，使其落入[0,1]区间

∗

−

x^*=\frac{-x_{min}}{x_{max}-x_{min}}

$x^{*} = \frac{- x _{m i n}}{x _{m a x} - x _{m i n}}$

这种方法的缺陷是有数据新加入时，导致xmax,xmin需要重新定义。

Z标准化

，标准差标准化，处理的数据符合

标准正态分布

∗

−

x^*=\frac{x-μ}{σ}

$x^{*} = \frac{x - μ}{σ}$

μ为所有样本的均值，σ为所有样本标准差。

离散化：消除连续观察点之间的差异

文档：离散化.note

链接：http://note.youdao.com/noteshare?id=a704d2330cc3edf3d39000c6675970bb&sub=A74F96CF9001453EBE68FD27B2E82797

语义转换：将{非常好,好,一般,差,非常差}–>{1,2,3,4,5}

3.1.3数据的统计

了解数据的基本特征。从样本推断总体，推断总体的数据特征。

统计的总体:人们研究对象的全体。

总体中的一个基本单位称为个体，个体的特征用变量x表示。

从总体中随机产生若干个体的集合称为样本，记为{x1,x2,…,xn}，n为样本容量。

基本描述性统计

假设一个容量为n的样本，记为x=(x1,x2…,xn),对它加工后，提出有用的信息。

统计量

是加工出来反映样本数量特征的函数，不包含任何未知量。

常见的统计量

表示位置：算法平均值和中位数

平均值定义：

‾

∑

\overline{x}=\frac{1}{n}\sum_{i=1}^nx_i

$\overline{x} = \frac{1}{n} i = 1 \sum n x_{i}$

中位数是将数据从小到大排序后位于中间位置的数。

mean(x);%x的均值
median(x);%返回中位数

表示数据散度：标准差、方差、极差

标准差定义：

[

−

∑

(

−

‾

)

]

s=[\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\overline{x})^2]^{\frac{1}{2}}

$s = [\frac{1}{n - 1} i = 1 \sum n (x_{i} - \overline{x})^{2}]^{\frac{1}{2}}$

它是各个数据与均值偏离程度的度量，这种偏移可称为

变异

。

被除以n-1是对无偏估计的要求。

若需要被n除，可用matlab中std(x,1)和var(x,1)。

方差：标准差的平方s^2。

极差x=(x1,x2…,xn)的最大值与最小值之差。

std(x)%标准差
var(x)%方差
range(x)%极差

表示分布形态的统计量：偏度、峰度

偏度反映分布的对称性，v1>0称右偏性，位于均值右边的数据比左边的数据多；v1<0相反。v1=0认为分布对称。

峰度是衡量偏离正态分布的尺度之一，正态分布的峰度为3，若v2比3大得多，表示分布有沉重的尾巴，说明样本中有较多远离均值的数据。

skewness(x)%返回x的偏度
kurtosis(x)%峰度

以上matlab计算统计量命令中，若x为矩阵，则作用于x的列返回一个行向量。

分布描述性统计

随机变量的特性完全由它的概率分布函数或概率密度函数来描述。F(x)=P{X≤x},分布函数定义为X≤x的概率。若X是连续型随机变量，则p(x)与F(x)关系为

(

)

∫

−

∞

(

)

F(x)=\int_{-\infty}^{x}{p(x)dx}

$F (x) = \int_{- \infty}^{x} p (x) d x$

分位数

定义为：对于0<α<1，使某分布函数F(x)=α的x称为这个分布的α分位数，记做x_α。

直方图是频数分布图。频数/样本容量n=频率。当n->∞，频率是概率的近似，因此直方图可看做密度函数图形的离散化近似，

3.1.4数据可视化

数据分布形态、中心分布、关联情况。

%链接：https://pan.baidu.com/s/1mvusLgomW1i6bUoS7z3DKA 
%提取码：pdws
% 数据可视化——基本绘图
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
% 读取数据
clc, clear al, close all
X=xlsread('dataTableA2.xlsx');
% 绘制变量dv1的基本分布
N=size(X,1);
id=1:N;
figure
plot( id', X(:,2),'LineWidth',1)
set(gca,'linewidth',2);
xlabel('编号','fontsize',12);
ylabel('dv1', 'fontsize',12);
title('变量dv1分布图','fontsize',12);
% 同时绘制变量dv1-dv4的柱状分布图
figure
subplot(2,2,1);
hist(X(:,2));
title('dv1柱状分布图','fontsize',12)
subplot(2,2,2);
hist(X(:,3));
title('dv2柱状分布图','fontsize',12)
subplot(2,2,3);
hist(X(:,4));
title('dv3柱状分布图','fontsize',12)
subplot(2,2,4);
hist(X(:,5));
title('dv4柱状分布图','fontsize',12)

%统计量绘图
% 数据可视化——数据分布形状图
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
% 读取数据
clc, clear al, close all
X=xlsread('dataTableA2.xlsx');
dv1=X(:,2);

% 绘制变量dv1的柱状分布图
h = -5:0.5:5;
n = hist(dv1,h);
figure
bar(h, n)

% 计算常用的形状度量指标
mn = mean(dv1); % 均值
sdev = std(dv1); % 标准差
mdsprd = iqr(dv1); % 四分位数
mnad = mad(dv1); % 中位数
rng = range(dv1); % 极差

% 标识度量数值
x = round(quantile(dv1,[0.25,0.5,0.75]));
y = (n(h==x(1)) + n(h==x(3)))/2;
line(x,[y,y,y],'marker','x','color','r')

x = round(mn + sdev*[-1,0,1]);
y = (n(h==x(1)) + n(h==x(3)))/2;
line(x,[y,y,y],'marker','o','color',[0 0.5 0])

x = round(mn + mnad*[-1,0,1]);
y = (n(h==x(1)) + n(h==x(3)))/2;
line(x,[y,y,y],'marker','*','color',[0.75 0 0.75])

x = round([min(dv1),max(dv1)]);
line(x,[1,1],'marker','.','color',[0 0.75 0.75])

legend('Data','Midspread','Std Dev','Mean Abs Dev','Range')

%数据关联可视化
% 数据可视化——变量想相关性
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
% 读取数据
clc, clear al, close all
X=xlsread('dataTableA2.xlsx');
Vars = X(:,7:12);
%  绘制变量间相关性关联图
figure
plotmatrix(Vars)
%  绘制变量间相关性强度图
covmat = corrcoef(Vars);
figure
imagesc(covmat);
grid;
colorbar;
%figure1绘制的变量间相关性关联图，可以看出任意两个变量的数据关联趋向
%figure2是变量间相关性强度图，用于筛选变量

%数据分组可视化
%箱体图按照不同的分位数将数据进行分组
% 数据可视化——数据分组
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
% 读取数据
clc, clear al, close all
X=xlsread('dataTableA2.xlsx');
dv1=X(:,2);
eva=X(:,12);
% Boxplot
figure
boxplot(X(:,2:12))
figure
boxplot(dv1, eva)
figure
boxplot(X(:,5))

文档：箱体图含义.note

链接：http://note.youdao.com/noteshare?id=0728298a7a1f716377feb14c32e13c2c&sub=96220FAE8BBE47E1A111648CC26CAD0F

3.1.5数据降维

PCA：使得众多相关的变成较少的无关的变量。

使用PCA分析案例中的影响因素

文档：PCA.note

链接：http://note.youdao.com/noteshare?id=c29d8e16615ccdd07ad20f865207a6f1&sub=B5D76F7ACA004ACDA5A9BAA3585776BD

PCA应用：企业综合实力排序

%链接：https://pan.baidu.com/s/198GcoR4dNJUp3vb9RRHHsA 
%提取码：2se9 
%% PCA数据降维实例
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 读取数据
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);

%% 显示结果
disp('特征值及贡献率：')
DS
disp('阀值T对应的主成分数与特征向量：')
Com_num
PV
disp('主成分分数：')
new_score
disp('主成分分数排序：')
new_score_s

%分析：第9综合实力最强 12综合实力最弱。
%各成分的权重信息：贡献率
%原始变量的关联关系：特征向量

3.2常用的数据建模方法

3.2.1回归

使用回归时，首先判断自变量个数，超过两个用多元回归。判断线性还是非线性，对于多元回归，保持其它变量不变，将多元转化为一元去判断。

根据回归方法因变量的个数和回归函数的类型(线性回归、非线性回归)，可将回归方法分为以下。

一元线性回归

%% 一元线性回归实例
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 输入数据
clc, clear all, close all
x=[23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40];
y=[41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0];	
%% 采用最小二乘回归
% 绘制散点图，判断是否具有线性关系
figure
plot(x,y,'r*')                         %作散点图
xlabel('x（职工工资总额）','fontsize', 12)           %横坐标名
ylabel('y（商品零售总额）', 'fontsize',12)           %纵坐标名
set(gca,'linewidth',2);
% 采用最小二乘拟合
Lxx=sum((x-mean(x)).^2);
Lxy=sum((x-mean(x)).*(y-mean(y)));
b1=Lxy/Lxx;
b0=mean(y)-b1*mean(x);
y1=b1*x+b0;
hold on
plot(x, y1,'linewidth',2);

%% 采用LinearModel.fit函数进行回归
m2 = LinearModel.fit(x,y)

%% 采用regress函数进行回归
Y=y';
X=[ones(size(x,2),1),x'];
[b, bint, r, rint, s] = regress(Y, X)
%注意Y为列向量，因变量   X为自变量组成的矩阵
%regress后缺省alpha(显著性水平)时默认是0.05
%b为β1，β2
%bint为b的置信区间
%r为残差
%rint为残差的置信区间


%s包含四个统计量：决定系数R^2(R为相关系数)、F值
%F(1,n-2)分布大于F值的概率p、剩余方差s^2
%s^2也可以由sum(r.^2)/(n-2)计算

%作用：1.R^2越接近1，变量线性相关越强，模型越有效
%2.(在1-α水平下)F(1,n-2)<F，认为y与x显著地有线性关系。
%它的值可直接通过finv(1-α,1,n-2)计算得到
%3.p<α表示线性模型可用
%4.s^2比较模型是否有所改进，值越小说明模型精度越高

一元非线性回归

%% 一元非线性回归实例
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 输入数据
clc, clear all, close all
x=[1.5, 4.5, 7.5,10.5,13.5,16.5,19.5,22.5,25.5];
y=[7.0,4.8,3.6,3.1,2.7,2.5,2.4,2.3,2.2];  
plot(x,y,'*','linewidth',2);
set(gca,'linewidth',2);
xlabel('销售额x/万元','fontsize', 12)           
ylabel('流通费率y/%', 'fontsize',12)           

%% 对数形式
m1 = @(b,x) b(1) + b(2)*log(x);
nonlinfit1 = fitnlm(x,y,m1,[0.01;0.01])
b=nonlinfit1.Coefficients.Estimate;
Y1=b(1,1)+b(2,1)*log(x);
hold on
plot(x,Y1,'--k','linewidth',2)

%% 指数形式拟合
m2 = 'y ~ b1*x^b2';
nonlinfit2 = fitnlm(x,y,m2,[1;1])
b1=nonlinfit2.Coefficients.Estimate(1,1);
b2=nonlinfit2.Coefficients.Estimate(2,1);
Y2=b1*x.^b2;
hold on
plot(x,Y2,'r','linewidth',2)
legend('原始数据','a+b*lnx','a*x^b')
%分析结果，对数的决定系数为0.973
%指数的决定系数为0.993
%选择指数形式的函数

多元回归

包括多元线性回归、多元多项式回归。

x1=[3.5 5.3 5.1 5.8 4.2 6.0 6.8 5.5 3.1 7.2 4.5 4.9 8.0 6.5 6.6 3.7 6.2 7.0 4.0 4.5 5.9 5.6 4.8 3.9];
x2=[9 20 18 33 31 13 25 30 5 47 25 11 23 35 39 21 7 40 35 23 33 27 34 15];
x3=[6.1 6.4 7.4 6.7 7.5 5.9 6.0 4.0 5.8 8.3 5.0 6.4 7.6 7.0 5.0 4.4 5.5 7.0 6.0 3.5 4.9 4.3 8.0 5.8];
Y=[33.2 40.3 38.7 46.8 41.4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.6 34.2 48.0 38.0 35.9 40.4 36.8 45.2 35.1];
%% 做出Y与各自变量的散点图
subplot(1,3,1),plot(x1,Y,'g*'),
subplot(1,3,2),plot(x2,Y,'k+'),
subplot(1,3,3),plot(x3,Y,'ro'),
%这些点大致分布在一条直线旁边，有比较好的线性相关，做线性回归
%% 多元线性回归
n=24;m=3;
X=[ones(n,1),x1',x2',x3'];
[b,bint,r,rint,s]=regress(Y',X,0.05);
b%回归系数
bint%置信区间
s%R^2、F、p、s^2

%%模型检验
%1.相关系数R^2接近1，相关性较强
%2.F检验：当F>F(m,n-m-1)时，认为y为x1,...,xm之间有显著的线性相关。(在1-α的水平下)
finv(0.95,3,n-4)%F>3.098
%3.p检验:p<α(预订的显著水平),那么y与x1,...,xm之间有显著的线性相关
%s^2越小越好，在改进模型时做参考

还有两种特殊的回归方式，一种是在回归过程中可以调整变量数量的回归方法，称为逐步回归。另一种是以指数结构函数作为回归模型的回归方法，称为logistic回归。

逐步回归

是多元线性回归、多元多项式回归的结合，使得方程因子最合理

X=[7,26,6,60
    1,29,15,52
    11,56,8,20
    11,31,8,47
    7,52,6,33
    11 55 9 22
    3 71 17 6
    1 31 22 44
    2 54 18 22
    21 47 4 26
    1 40 23 34
    11 66 9 12];%自变量  一行是一个样本，每个样本四个因素
Y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3];
stepwise(X,Y,[1,2,3,4],0.05,0.10)%[1,2,3,4]表示4个自变量均保存在模型中
%一直单击next step直至变灰，便为最终结果

logistic回归(预测)

%链接：https://pan.baidu.com/s/1ljksQlTi53wCML_SfHm4CA 
%提取码：wny1 
% 程序P7-1： logistic方法Matlab实现程序
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 数据准备
clear all
clc
X0=xlsread('logistic_ex1.xlsx', 'A2:C21'); % 回归数据X值
X1=xlsread('logistic_ex1.xlsx', 'A2:C26'); % 验证与预测数据
Y0=xlsread('logistic_ex1.xlsx', 'D2:D21'); % 回归数据P值
%--------------------------------------------------------------------------
%% logistic函数
GM=fitglm(X0,Y0,'Distribution','binomial');
Y1=predict(GM,X1);

%% 模型评估
N0=1:size(Y0,1);N1=1:size(Y1,1);
plot(N0',Y0,'-kd');
hold on;scatter(N1',Y1,'b')
xlabel('数据点编号');ylabel('输出值');

3.2.2机器学习

利用历史数据”训练”,”学习”到某种规律、模式,建立预测未来的模型。

有两种学习方法：

有监督学习：用于决策支持。包括分类、回归。
无监督学习：用于知识发现。包括聚类(隐藏的模式)。

分类方法：

K-NN(K-近邻)
贝叶斯
神经网络
logistic
判别分析
SVM
决策树

聚类方法：

K-mens
层次聚类
神经网络
高斯混合
模糊C均值

分类

文档：K-NN.note

链接：http://note.youdao.com/noteshare?id=fedb9e76f023e0db060fdc01093b1ab8&sub=74C32A6EBAF44C9FA67CEDE19EA94FB8

文档：贝叶斯分类.note

链接：http://note.youdao.com/noteshare?id=eb1fe28035ea3efd72ce6b4699236614&sub=76A05961496045E58B698991F4E18AB7

文档：支持向量机分类.note

链接：http://note.youdao.com/noteshare?id=39f394075d8d4460c57397a2eae7b956&sub=EF8A15FDAE46411EB49AB2CD4B1C03D3

K-NN实现分类器、朴素贝叶斯算法实现分类器、神经网络、Logistic、判别分析、SVM、决策树的matlab代码实现

一家银行通过电话调查客户是否愿意买理财产品，并记录调查结果为y。另外银行有这些客户的资料X，包括十六个属性(age、job…)。希望建立一个分类器，预测一个新客服是够愿意购买。

%链接：https://pan.baidu.com/s/1RQdaIx6_EYTgKuf4DJtP2w 
%提取码：ue1o 

%% 分类方法示例程序
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 1.
clc, clear all, close all
%% 2.导入数据及数据预处理 
load bank.mat
% 将分类变量转换成分类数组
names = bank.Properties.VariableNames;
category = varfun(@iscellstr, bank, 'Output', 'uniform');
for i = find(category)
    bank.(names{i}) = categorical(bank.(names{i}));
end
% 跟踪分类变量
catPred = category(1:end-1);
% 设置默认随机数生成方式确保该脚本中的结果是可以重现的
rng('default');
% 数据探索----数据可视化
figure(1)
gscatter(bank.balance,bank.duration,bank.y,'kk','xo')
xlabel('年平均余额/万元', 'fontsize',12)
ylabel('上次接触时间/秒', 'fontsize',12)
title('数据可视化效果图', 'fontsize',12)
set(gca,'linewidth',2);

% 设置响应变量和预测变量
X = table2array(varfun(@double, bank(:,1:end-1)));  % 预测变量
Y = bank.y;   % 响应变量
disp('数据中Yes & No的统计结果：')
tabulate(Y)
%将分类数组进一步转换成二进制数组以便于某些算法对分类变量的处理
XNum = [X(:,~catPred) dummyvar(X(:,catPred))];
YNum = double(Y)-1;

%% 3.设置交叉验证方式
% 随机选择40%的样本作为测试样本
cv = cvpartition(height(bank),'holdout',0.40);
% 训练集
Xtrain = X(training(cv),:);
Ytrain = Y(training(cv),:);
XtrainNum = XNum(training(cv),:);
YtrainNum = YNum(training(cv),:);
% 测试集
Xtest = X(test(cv),:);
Ytest = Y(test(cv),:);
XtestNum = XNum(test(cv),:);
YtestNum = YNum(test(cv),:);
disp('训练集：')
tabulate(Ytrain)
disp('测试集：')
tabulate(Ytest)

%% 最近邻
% 4.训练分类器
knn = ClassificationKNN.fit(Xtrain,Ytrain,'Distance','seuclidean',...
                            'NumNeighbors',5);
% 进行预测
[Y_knn, Yscore_knn] = knn.predict(Xtest);
Yscore_knn = Yscore_knn(:,2);
% 计算混淆矩阵
disp('最近邻方法分类结果：')
C_knn = confusionmat(Ytest,Y_knn)






%% 贝叶斯
% 设置分布类型
dist = repmat({'normal'},1,width(bank)-1);
dist(catPred) = {'mvmn'};
% 训练分类器
Nb = NaiveBayes.fit(Xtrain,Ytrain,'Distribution',dist);
% 进行预测
Y_Nb = Nb.predict(Xtest);
Yscore_Nb = Nb.posterior(Xtest);
Yscore_Nb = Yscore_Nb(:,2);
% 计算混淆矩阵
disp('贝叶斯方法分类结果：')
C_nb = confusionmat(Ytest,Y_Nb)




%% 神经网络
% 设置神经网络模式及参数
hiddenLayerSize = 5;
net = patternnet(hiddenLayerSize);
% 设置训练集、验证机和测试集
net.divideParam.trainRatio = 70/100;
net.divideParam.valRatio = 15/100;
net.divideParam.testRatio = 15/100;
% 训练网络
net.trainParam.showWindow = false;
inputs = XtrainNum';
targets = YtrainNum';
[net,~] = train(net,inputs,targets);
% 用测试集数据进行预测
Yscore_nn = net(XtestNum')';
Y_nn = round(Yscore_nn);
% 计算混淆矩阵
disp('神经网络方法分类结果：')
C_nn = confusionmat(YtestNum,Y_nn)

%% Logistic
% 训练分类器
glm = fitglm(Xtrain,YtrainNum,'linear', 'Distribution','binomial',...
    'link','logit','CategoricalVars',catPred, 'VarNames', names);
% 用测试集数据进行预测
Yscore_glm = glm.predict(Xtest);
Y_glm = round(Yscore_glm);
% 计算混淆矩阵
disp('Logistic方法分类结果：')
C_glm = confusionmat(YtestNum,Y_glm)

%% 判别分析
% 训练分类器
da = ClassificationDiscriminant.fit(XtrainNum,Ytrain);
% 进行预测
[Y_da, Yscore_da] = da.predict(XtestNum); 
Yscore_da = Yscore_da(:,2);
% 计算混淆矩阵
disp('判别方法分类结果：')
C_da = confusionmat(Ytest,Y_da)



%% 支持向量机(SVM)
% 设置最大迭代次数
opts = statset('MaxIter',45000);
% 训练分类器
svmStruct = svmtrain(Xtrain,Ytrain,'kernel_function','linear','kktviolationlevel',0.2,'options',opts);
% 进行预测
Y_svm = svmclassify(svmStruct,Xtest);
Yscore_svm = svmscore(svmStruct, Xtest);
Yscore_svm = (Yscore_svm - min(Yscore_svm))/range(Yscore_svm);
% 计算混淆矩阵
disp('SVM方法分类结果：')
C_svm = confusionmat(Ytest,Y_svm)



%% 决策树
% 训练分类器
t = ClassificationTree.fit(Xtrain,Ytrain,'CategoricalPredictors',catPred);
% 进行预测
Y_t = t.predict(Xtest);
% 计算混淆矩阵
disp('决策树方法分类结果：')
C_t = confusionmat(Ytest,Y_t)

%% 通过ROC曲线来比较方法
methods = {'KNN','NBayes','NNet', 'GLM',  'LDA', 'SVM'};
scores = [Yscore_knn, Yscore_Nb, Yscore_nn, Yscore_glm, Yscore_da,  Yscore_svm];
%绘制ROC曲线
figure
auc= zeros(6); hCurve = zeros(1,6);
for ii=1:6;
  [rocx, rocy, ~, auc(ii)] = perfcurve(Ytest, scores(:,ii), 'yes');
  hCurve(ii,:) = plot(rocx, rocy, 'k','LineWidth',2); hold on;
end
legend(hCurve(:,1), methods)
set(gca,'linewidth',2);
grid on;
title('各方法ROC曲线', 'fontsize',12); 
xlabel('假阳率 [ = FP/(TN+FP)]', 'fontsize',12); 
ylabel('真阳率 [ = TP/(TP+FN)]', 'fontsize',12);
% 绘制各方法分类正确率
figure;
bar(auc); set(gca,'YGrid', 'on','XTickLabel',methods); 
xlabel('方法简称', 'fontsize',12); 
ylabel('分类正确率', 'fontsize',12);
title('各方法分类正确率','fontsize',12);
set(gca,'linewidth',2);



%svmscore.m
function f = svmscore(svm,X)
% Compute SVM score for classes -1 and +1

% Copyright 2014 The MathWorks, Inc.

shift = svm.ScaleData.shift;
scale = svm.ScaleData.scaleFactor;

X = bsxfun(@plus,X,shift);
X = bsxfun(@times,X,scale);

sv = svm.SupportVectors;
alphaHat = svm.Alpha;
bias = svm.Bias;
kfun = svm.KernelFunction;
kfunargs = svm.KernelFunctionArgs;

f = kfun(sv,X,kfunargs{:})'*alphaHat(:) + bias;
f = -f; % flip the sign to get the score for the +1 class

end

聚类

文档：k-means聚类.note

链接：http://note.youdao.com/noteshare?id=34c9772b906cca6bda8c1c37dc3b1f7d&sub=F37E921C3837492F8689DDC3FCAA1993

已知20个样本，每个样本两种特征，对这些数据分类

%% K-means方法的MATLAB实现
% 《MATLAB数学建模方法与实践》(《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),'k*',...
    'LineWidth',2,...
    'MarkerSize',10,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.5,0.5,0.5])
hold on
plot(z1(:,1),z1(:,2),'ko',...
    'LineWidth',2,...
    'MarkerSize',10,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor',[0.5,0.5,0.5])
set(gca,'linewidth',2) ;
xlabel('特征x1','fontsize',12);
ylabel('特征x2', 'fontsize',12);
title('K-means分类图','fontsize',12);

一家银行希望对债券分类，但不知道分成几类，已知债券一些属性(Type、name)，请确定分成几类。

k-means聚类(使用kmeans函数)、层次聚类、模糊C-均值聚类、高斯混合聚类 (GMM)、神经网络matlab实现分类。

文档：层次聚类.note

链接：http://note.youdao.com/noteshare?id=04f1ec00ee41e0ebc4b178d4a509b4bc&sub=37626CCFD47D4170BD4FBEAE3B4C43DD

文档：模糊C-均值聚类(FCM).note

链接：http://note.youdao.com/noteshare?id=e8e959eebc264887ee101d0ec3ed1c3d&sub=3512CD8F115F4551B199F2C6687950CB

%链接：https://pan.baidu.com/s/14D7BeBF2PXjyWoO5MF1BrA 
%提取码：ob7o
%% 聚类方法
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 导入数据和预处理数据
clc, clear all, close all
load BondData
settle = floor(date);
%数据预处理
bondData.MaturityN = datenum(bondData.Maturity, 'dd-mmm-yyyy');
bondData.SettleN = settle * ones(height(bondData),1);
% 筛选数据
corp = bondData(bondData.MaturityN > settle & ...
             bondData.Type == 'Corp' & ...
             bondData.Rating >= 'CC' & ...
             bondData.YTM < 30 & ...
             bondData.YTM >= 0, :);
% 设置随机数生成方式保证结果可重现
rng('default');

%% 探索数据
% 数据可视化
figure
gscatter(corp.Coupon,corp.YTM,corp.Rating)
set(gca,'linewidth',2);
xlabel('票面利率')
ylabel('到期收益率')

% 选择聚类变量
corp.RatingNum = double(corp.Rating);
bonds = corp{:,{'Coupon','YTM','CurrentYield','RatingNum'}};

% 设置类别数量
numClust = 3;

% 设置用于可视化聚类效果的变量
VX=[corp.Coupon, double(corp.Rating), corp.YTM];

%% K-Means 聚类
dist_k = 'cosine';
kidx = kmeans(bonds, numClust, 'distance', dist_k);

%绘制聚类效果图
figure
F1 = plot3(VX(kidx==1,1), VX(kidx==1,2),VX(kidx==1,3),'r*', ...
           VX(kidx==2,1), VX(kidx==2,2),VX(kidx==2,3), 'bo', ...
           VX(kidx==3,1), VX(kidx==3,2),VX(kidx==3,3), 'kd');
set(gca,'linewidth',2);
grid on;
set(F1,'linewidth',2, 'MarkerSize',8);
xlabel('票面利率','fontsize',12);
ylabel('评级得分','fontsize',12);
ylabel('到期收益率','fontsize',12);
title('Kmeans方法聚类结果')

% 评估各类别的相关程度
dist_metric_k = pdist(bonds,dist_k);
dd_k = squareform(dist_metric_k);
[~,idx] = sort(kidx);
dd_k = dd_k(idx,idx);
figure
imagesc(dd_k)
set(gca,'linewidth',2);
xlabel('数据点','fontsize',12)
ylabel('数据点', 'fontsize',12)
title('k-Means聚类结果相关程度图', 'fontsize',12)
ylabel(colorbar,['距离矩阵:', dist_k])
axis square

%% 层次聚类
dist_h = 'spearman';
link = 'weighted';
hidx = clusterdata(bonds, 'maxclust', numClust, 'distance' , dist_h, 'linkage', link);

%绘制聚类效果图
figure
F2 = plot3(VX(hidx==1,1), VX(hidx==1,2),VX(hidx==1,3),'r*', ...
           VX(hidx==2,1), VX(hidx==2,2),VX(hidx==2,3), 'bo', ...
           VX(hidx==3,1), VX(hidx==3,2),VX(hidx==3,3), 'kd');
set(gca,'linewidth',2);
grid on
set(F2,'linewidth',2, 'MarkerSize',8);
set(gca,'linewidth',2);
xlabel('票面利率','fontsize',12);
ylabel('评级得分','fontsize',12);
ylabel('到期收益率','fontsize',12);
title('层次聚类方法聚类结果')

% 评估各类别的相关程度
dist_metric_h = pdist(bonds,dist_h);
dd_h = squareform(dist_metric_h);
[~,idx] = sort(hidx);
dd_h = dd_h(idx,idx);
figure
imagesc(dd_h)
set(gca,'linewidth',2);
xlabel('数据点', 'fontsize',12)
ylabel('数据点', 'fontsize',12)
title('层次聚类结果相关程度图')
ylabel(colorbar,['距离矩阵:', dist_h])
axis square

% 计算同型相关系数
Z = linkage(dist_metric_h,link);
cpcc = cophenet(Z,dist_metric_h);
disp('同型相关系数: ')
disp(cpcc)

% 层次结构图
set(0,'RecursionLimit',5000)
figure
dendrogram(Z)
set(gca,'linewidth',2);
set(0,'RecursionLimit',500)
xlabel('数据点', 'fontsize',12)
ylabel ('标准距离', 'fontsize',12)
 title('层次聚类法层次结构图')

%% 神经网络
%设置网络
dimension1 = 3;
dimension2 = 1;
net = selforgmap([dimension1 dimension2]);
net.trainParam.showWindow = 0;

%训练网络
[net,tr] = train(net,bonds');
nidx = net(bonds');
nidx = vec2ind(nidx)';

%绘制聚类效果图
figure
F3 = plot3(VX(nidx==1,1), VX(nidx==1,2),VX(nidx==1,3),'r*', ...
           VX(nidx==2,1), VX(nidx==2,2),VX(nidx==2,3), 'bo', ...
           VX(nidx==3,1), VX(nidx==3,2),VX(nidx==3,3), 'kd');
set(gca,'linewidth',2);
grid on
set(F3,'linewidth',2, 'MarkerSize',8);
xlabel('票面利率','fontsize',12);
ylabel('评级得分','fontsize',12);
ylabel('到期收益率','fontsize',12);
title('神经网络方法聚类结果')

%% 模糊C-Means聚类
options = nan(4,1);
options(4) = 0;
[centres,U] = fcm(bonds,numClust, options);
[~, fidx] = max(U);
fidx = fidx';

%绘制聚类效果图
figure
F4 = plot3(VX(fidx==1,1), VX(fidx==1,2),VX(fidx==1,3),'r*', ...
           VX(fidx==2,1), VX(fidx==2,2),VX(fidx==2,3), 'bo', ...
           VX(fidx==3,1), VX(fidx==3,2),VX(fidx==3,3), 'kd');
set(gca,'linewidth',2);
grid on
set(F4,'linewidth',2, 'MarkerSize',8);
xlabel('票面利率','fontsize',12);
ylabel('评级得分','fontsize',12);
ylabel('到期收益率','fontsize',12);
title('模糊C-Means方法聚类结果')

%% 高斯混合聚类 (GMM)
gmobj = gmdistribution.fit(bonds,numClust);
gidx = cluster(gmobj,bonds);

%绘制聚类效果图
figure
F5 = plot3(VX(fidx==1,1), VX(fidx==1,2),VX(fidx==1,3),'r*', ...
           VX(fidx==2,1), VX(fidx==2,2),VX(fidx==2,3), 'bo', ...
           VX(fidx==3,1), VX(fidx==3,2),VX(fidx==3,3), 'kd');
set(gca,'linewidth',2);
grid on
set(F5,'linewidth',2, 'MarkerSize',8);
xlabel('票面利率','fontsize',12);
ylabel('评级得分','fontsize',12);
ylabel('到期收益率','fontsize',12);
title('高斯混合方法聚类结果')

%% k-Means方法确定最佳的聚类类别数
% 绘制几个典型类别数情况下的平均轮廓值图
figure
for i=2:4
    kidx = kmeans(bonds,i,'distance',dist_k);
    subplot(3,1,i-1)
    [~,F6] = silhouette(bonds,kidx,dist_k);
    xlabel('轮廓值','fontsize',12);
    ylabel('类别数','fontsize',12);
    set(gca,'linewidth',2);
    title([num2str(i) '类对应的轮廓值图 ' ])
    snapnow
end


% 计算平均轮廓值
numC = 15;
silh_m = nan(numC,1);
for i=1:numC
    kidx = kmeans(bonds,i,'distance',dist_k,'MaxIter',500);
    silh = silhouette(bonds,kidx,dist_k);
    silh_m(i) = mean(silh);
end

%绘制各类别数对应的平均轮廓值图
figure
F7 = plot(1:numC,silh_m,'o-');
set(gca,'linewidth',2);
set(F7, 'linewidth',2, 'MarkerSize',8);
xlabel('类别数', 'fontsize',12)
ylabel('平均轮廓值','fontsize',12)
title('平均轮廓值vs.类别数')
%% 聚类方案案例

小结

机器学习主要还是研究分类和聚类。分类为主，聚类为辅，往往先通过聚类确定分类的最佳类别。

在需要选择特征的一类建模问题中，深度学习就有用了。

3.2.3深度学习

文档：深度学习.note

链接：http://note.youdao.com/noteshare?id=a7eb6815a95db0160302154dc38d822c&sub=F301FCC29CBB46479F32A2ED52C86747

研究的问题：人类活动的分类，人类活动传感器数据来自于人们不同活动时手机传感器测量的观测值。目的是建立一个分类器，自动识别给定传感器测量的活动类型。

数据来自https://archive.ics.uci.edu

%注意matlab版本要在2018，否则运行出错
%% 深度学习案例：人类行为的分类
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 

%链接：https://pan.baidu.com/s/1MY6e0EPRvEKjTiUhS9tpmQ 
%提取码：p1fm

%Human_Activity_Learning_DL.m
%% 下载数据
if false %~exist('UCI HAR Dataset','file')
    downloadSensorData;
end
if ~exist('rawSensorData_train.mat','file') && ~exist('rawSensorData_test.mat','file')
    LoadSensorData;
end
load rawSensorData_train

%% 定义深度学习结构
allRawDataDL = cat(3, body_gyro_x_train, body_gyro_y_train, body_gyro_z_train, total_acc_x_train, total_acc_y_train, total_acc_z_train);
C = num2cell(allRawDataDL, [2 3]);
C = cellfun(@squeeze, C, 'UniformOutput', false);
trainingData = table(C);
trainingData.activity = categorical(trainActivity);
% class(trainingData{:,1}) % should be cell

layers = [imageInputLayer([128 6])
            convolution2dLayer(3, 2)
            reluLayer
            maxPooling2dLayer([12 2], 'Stride', 1)
            fullyConnectedLayer(5)
            softmaxLayer
            classificationLayer()];
        
options = trainingOptions('sgdm','MaxEpochs',15, ...
	'InitialLearnRate',0.005);

convnet = trainNetwork(trainingData, layers, options);

%% 训练深度网络
load rawSensorData_test
%
allRawDataTestDL = cat(3, body_gyro_x_test, body_gyro_y_test, body_gyro_z_test, total_acc_x_test, total_acc_y_test, total_acc_z_test);
Ctest = num2cell(allRawDataTestDL, [2 3]);
Ctest = cellfun(@squeeze, Ctest, 'UniformOutput', false);
testData = table(Ctest);
testData.activity = categorical(testActivity);
Y_test = classify(convnet, testData(:,1));
accuracy_test = sum(Y_test == testActivity)/numel(testActivity) %#ok<*NOPTS>
cm = confusionmat(testActivity, Y_test);

% Display in a table
test_results = array2table(cm, ...
    'RowNames', {'Walking', 'ClibmingStairs', 'Sitting', 'Standing', 'Laying'}, ...  
    'VariableNames', {'Walking', 'ClibmingStairs', 'Sitting', 'Standing', 'Laying'})


%downloadSensorData.m
function downloadSensorData
%% Download and extract data if data folder does not exist
if ~(exist('UCI HAR Dataset','file') == 7)
    downloadlink = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00240/UCI%20HAR%20Dataset.zip';
    toc
    disp('UCI HAR Dataset does not exist, downloading dataset...')
    fname = websave('UCI_HAR_Dataset.zip',downloadlink);
    toc
    disp('Done downloading, file downloaded to:')
    disp(fname)
    disp(' ')
    yn = input('It may be faster to extract the file manually, do you want MATLAB to extract the file (y/n)? ','s');
    if strcmp(yn,'y')
        tic
        disp('Extracting file, this may take a while...')
        foldername = unzip(fname);
        disp('Done extracting, extracted to:')
        disp(foldername)
        toc
    else 
        disp(' ')
        disp('OK, You must manually extract the file contents to the current folder before proceeding. Don''t change the default folder names.')
    end
end


%LoadSensorData.m
function LoadSensorData

if exist('rawSensorData_train.mat','file') && exist('rawSensorData_test.mat','file')
    fprintf(1,'rawSensorData_train.mat and rawSensorData_test.mat already exists at location:\n');
    disp(['* ', which('rawSensorData_train.mat')]);
    disp(['* ', which('rawSensorData_test.mat')]);
    disp(' ')
else
    %% Load training data from files
    activity_labels = {'Walking','WalkingUpstairs','WalkingDownstairs','Sitting','Standing','Laying'};
    trainActivity = categorical(importdata('UCI HAR Dataset\train\y_train.txt'),1:6,activity_labels);
    trainActivity = mergecats(trainActivity,{'WalkingUpstairs','WalkingDownstairs'},'ClimbingStairs');
    % Choose this if you want to only load the total acca and gyro data
    filestoload = strcat('UCI HAR Dataset\train\Inertial Signals\',{'*gyro*','total*'});
	% Choose this if you want to load all files
    % filestoload = strcat('UCI HAR Dataset\train\Inertial Signals\'); 

    disp('Loading training data from files:')
    try
        dstrain = datastore(filestoload,'TextscanFormats',repmat({'%f'},1,128), 'ReadVariableNames',false);
    catch err
        if strcmp(err.identifier,'MATLAB:datastoreio:pathlookup:fileNotFound')
            error('File not found. Please make sure that you download and extract the data first using ''downloadSensorData'' function')
        end
    end
    dstrain.ReadSize = 'file';
    [~,fnames] = cellfun(@fileparts,dstrain.Files,'UniformOutput',false);
    iter = 1;
%     rawSensorDataTrain = table;
    while hasdata(dstrain)
        fprintf('Importing: %16s ...',fnames{iter})
        M = table2array(read(dstrain));
        rawSensorDataTrain.(fnames{iter}) = M;
        iter = iter + 1;
        fprintf('Done\n')
    end
    rawSensorDataTrain.trainActivity = trainActivity;
    disp(' ')
    %% Load test data from files
    testActivity = categorical(importdata('UCI HAR Dataset\test\y_test.txt'),1:6,activity_labels);
    testActivity = mergecats(testActivity,{'WalkingUpstairs','WalkingDownstairs'},'ClimbingStairs');
	% Choose this if you want to only load the total acca and gyro data
    filestoload = strcat('UCI HAR Dataset\test\Inertial Signals\',{'*gyro*','total*'});
	% Choose this if you want to load all files
%     filestoload = strcat('UCI HAR Dataset\test\Inertial Signals\'); 
    disp('Loading test data from files:')
    dstest = datastore(filestoload,'TextscanFormats',repmat({'%f'},1,128),'DatastoreType','tabulartext',...
        'ReadVariableNames',false);
    [~,fnames] = cellfun(@fileparts,dstest.Files,'UniformOutput',false);
    dstest.ReadSize = 'file';
    iter = 1;
%     rawSensorDataTest = table;
    while hasdata(dstest)
        fprintf('Importing: %16s ...',fnames{iter})
        M = table2array(read(dstest));
        rawSensorDataTest.(fnames{iter}) = M;
        iter = iter + 1;
        fprintf('Done\n')
    end
	rawSensorDataTest.testActivity = testActivity;
disp(' ')
    %% Saving MAT file with raw data
    fprintf('Saving MAT files: rawSensorData_train.mat ...')
    save rawSensorData_train.mat -struct rawSensorDataTrain
    disp('Done')
    fprintf('Saving MAT files: rawSensorData_test.mat ...')
    save rawSensorData_test.mat -struct rawSensorDataTest
    disp('Done')
end

3.2.4灰色预测

文档：灰色预测.note

链接：http://note.youdao.com/noteshare?id=ff4969e2b9bb9ef22c53dde351c2191c&sub=ED51AD0AB4DB4B399698543351A2A545

% 灰色预测
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%根据已知的利润，预测未来几年的利润
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,'ko', 'LineWidth',2)
hold on
plot(t2,G,'k', 'LineWidth',2)
xlabel('年份', 'fontsize',12)
ylabel('利润/(元/年)','fontsize',12)
set(gca,  'LineWidth',2);

长江水质预测2005A

文档：2009D预测与会代表人数.note

链接：http://note.youdao.com/noteshare?id=bb58edae81c8c534415bccb79bda35e8&sub=F8AC603142904C118AA7BEEA79444729

% 灰色预测
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
clear
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', 'LineWidth',2)
hold on
plot(t2,G,'k', 'LineWidth',2)
xlabel('年份', 'fontsize',12)
ylabel('污水量/亿吨','fontsize',12)
set(gca,  'LineWidth',2);

3.2.5神经网络

文档：神经网络.note

链接：http://note.youdao.com/noteshare?id=402c9bce01ad986fefb6bc8d76fde8c1&sub=D1ABB335BC144500BFC83761EDFDCBC9

银行市场调查的分类器：一家银行通过电话调查客户是否愿意买理财产品，并记录调查结果为y。另外银行有这些客户的资料X，包括十六个属性(age、job…)。希望建立一个分类器，预测一个新客服是够愿意购买。

代码在机器学习分类的地方。

3.2.6小波分析

文档：小波分析.note

链接：http://note.youdao.com/noteshare?id=3b06893f971f1540a53243dbc4e290ca&sub=689B92DD8CBF4B8E9FEF0788DAB88DE5

小波图像处理：给定一个图像信号，由于图像经过二维小波分解后，图像的轮廓主要体现在低频部分，细节部分在高频部分，因此可以通过对低频分解系数进行增强处理，对高频分解系数进行衰弱处理，从而到达图像增强的效果。

%wmandril.mat
%链接：https://pan.baidu.com/s/1wDV7p4a7Q9dNBfeWmduOOQ 
%提取码：rcci

load wmandril
%图像增强处理
%用小波函数sym4对X进行2层小波分解
[c,s]=wavedec2(X,2,'sym4');
sizec=size(c);c1=c;
%对分解系数进行处理突出轮廓，弱化细节
for i=1:sizec(2)
	if(c(i)>350)
        c1(i)=2*c(i);
    else
        c1(i)=0.5*c(i);
	end
end
%对处理后的系数重构
xx=waverec2(c1,s,'sym4');
%画图
colormap(map);
subplot(121);image(X);title('原始图像');axis square
subplot(122);image(xx);title('增强图像');axis square

小波数据去噪：


% 小波去噪
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
clc, clear all, close all,
load nelec.mat;
sig = nelec;
denPAR = {[1 94 5.9 ; 94 1110 19.5 ; 1110 2000 4.5]};
wname = 'sym4';
level = 5;
sorh  = 's'; % type of thresholding
thr = 4.5;
[sigden_1,~,~,perf0,perfl2] = wdencmp('gbl',sig,wname,level,thr,sorh,1);
res = sig-sigden_1;
subplot(3,1,1);plot(sig,'r');       axis tight
title('Original Signal')
subplot(3,1,2);plot(sigden_1,'b');  axis tight
title('Denoised Signal');
subplot(3,1,3);plot(res,'k');       axis tight
title('Residual');
% perf0,perfl2

小结

灰色和神经网络一般用于预测，灰色系统适合小样本数据，神经网络更适合大样本数据。神经网络也适合多输入多输出的复杂预测问题。

小波方法主要用于数据的预处理，比如去噪、提取数据特征、图像的增强。

3.2.7统计、主成分分析、时间序列分析略

4.优化技术

离散系统的优化问题一般都可以通过规划模型求解。

4.1标准规划问题(基础)

4.1.1线性规划

文档：线性规划.note

链接：http://note.youdao.com/noteshare?id=d72a55966042c60799ff3a4a928986af&sub=06A486777B484A38B1698CE6414CDF61

4.1.2非线性规划

罚函数法将求解非线性规划转化为

无约束极值问题

。

文档：非线性规划、二次规划.note

链接：http://note.youdao.com/noteshare?id=1b312c603bd99e9d48b638c48fd5012a&sub=B1035D918BEA4B888B9D109EE4602A89

4.1.3整数规划

文档：整数规划.note

链接：http://note.youdao.com/noteshare?id=8efd37cc9c498f2a07a290473077fe74&sub=9E9F24B605CA4E61BAB5FB904A6CE6ED

% 0-1整数规划
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
clc, clear, close all
f=[-3; 2; -5];
intcon=3;
A=[1 2 -1; 1 4 1; 1 1 0; 0 4 1];
b=[2; 4; 3; 6];
lb=[0, 0 , 0];
ub=[1,1,1];
Aeq=[0,0,0];
beq=0;
x = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)

4.2全局优化算法(组合优化)

包括灾情巡视、公交车调度、彩票、露天卡车调度、交巡警服务平台、太阳影子定位。

用于求解组合类，不适合大规模遍历。遗传算法和模拟退火算法、蚁群算法最常用。

4.2.1全局优化工具箱

global optimization toolbox(2018a及以后才有)

全局搜索：globalsearch 寻找全局最小值
多起点搜索：multistart 寻找多个局部最小值
模式搜索：patternsearch 用模式搜索方法寻找函数最小值
遗传：Ga 用遗传算法寻找函数最小值
粒子群：particlewarm 用粒子群算法寻找函数最小值
模拟退火：simulannealbnd 用模拟退火算法寻找函数最小值

https://cn.mathworks.com/help/gads/index.html

中文：

https://cn.mathworks.com/products/global-optimization.html

4.2.2遗传算法

文档：遗传算法.note

链接：http://note.youdao.com/noteshare?id=6bb173ca2bd09021b96d6bc61f5b4577&sub=B6AD34532DB8461DBB4DED5EA4A68244

遗传算法求解TSP

%traveling_salesman_main.m
%% GA算法求解旅行商问题
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
 %% 加载问题的数据
load('usborder.mat','x','y','xx','yy');

cities = 40;
locations = zeros(cities,2);
n = 1;
while (n <= cities)
    xp = rand*1.5;
    yp = rand;
    if inpolygon(xp,yp,xx,yy)
        locations(n,1) = xp;
        locations(n,2) = yp;
        n = n+1;
    end
end
plot(locations(:,1),locations(:,2),'bo');
xlabel('城市的横坐标x'); ylabel('城市的纵坐标y'); 
grid on

%% 计算城市距离

distances = zeros(cities);
for count1=1:cities
    for count2=1:count1
        x1 = locations(count1,1);
        y1 = locations(count1,2);
        x2 = locations(count2,1);
        y2 = locations(count2,2);
        distances(count1,count2)=sqrt((x1-x2)^2+(y1-y2)^2);
        distances(count2,count1)=distances(count1,count2);
    end
end


%% 定义目标函数
FitnessFcn = @(x) traveling_salesman_fitness(x,distances);

my_plot = @(options,state,flag) traveling_salesman_plot(options, ...
    state,flag,locations);

%% 设置优化属性并执行遗传算法求解

options = optimoptions(@ga, 'PopulationType', 'custom','InitialPopulationRange', ...
                            [1;cities]);

options = optimoptions(options,'CreationFcn',@create_permutations, ...
                        'CrossoverFcn',@crossover_permutation, ...
                        'MutationFcn',@mutate_permutation, ...
                        'PlotFcn', my_plot, ...
                        'MaxGenerations',500,'PopulationSize',60, ...
                        'MaxStallGenerations',200,'UseVectorized',true);

numberOfVariables = cities;
[x,fval,reason,output] = ...
    ga(FitnessFcn,numberOfVariables,[],[],[],[],[],[],[],options);




%traveling_salesman_fitness.m
function scores = traveling_salesman_fitness(x,distances)
%TRAVELING_SALESMAN_FITNESS  Custom fitness function for TSP. 
%   SCORES = TRAVELING_SALESMAN_FITNESS(X,DISTANCES) Calculate the fitness 
%   of an individual. The fitness is the total distance traveled for an
%   ordered set of cities in X. DISTANCE(A,B) is the distance from the city
%   A to the city B.

%   Copyright 2004-2007 The MathWorks, Inc.

scores = zeros(size(x,1),1);
for j = 1:size(x,1)
    % here is where the special knowledge that the population is a cell
    % array is used. Normally, this would be pop(j,:);
    p = x{j}; 
    f = distances(p(end),p(1));
    for i = 2:length(p)
        f = f + distances(p(i-1),p(i));
    end
    scores(j) = f;
end




%traveling_salesman_plot.m
function state = traveling_salesman_plot(options,state,flag,locations)
%   TRAVELING_SALESMAN_PLOT Custom plot function for traveling salesman.
%   STATE = TRAVELING_SALESMAN_PLOT(OPTIONS,STATE,FLAG,LOCATIONS) Plot city
%   LOCATIONS and connecting route between them. This function is specific
%   to the traveling salesman problem.

%   Copyright 2004-2006 The MathWorks, Inc.
persistent x y xx yy
if strcmpi(flag,'init')
  load('usborder.mat','x','y','xx','yy');
end
% plot(x,y,'Color','red');
% hold on;
[unused,i] = min(state.Score);
genotype = state.Population{i};

plot(locations(:,1),locations(:,2),'bo');
axis([-0.1 1.5 -0.2 1.2]);
hold on
plot(locations(genotype,1),locations(genotype,2));
xlabel('城市的横坐标x'); ylabel('城市的纵坐标y'); 
grid on
hold off

4.2.3模拟退火算法

文档：模拟退火算法.note

链接：http://note.youdao.com/noteshare?id=07e5c2a4685beae2a4d89c2219e28708&sub=58F3CAAB22964113AE23EA0E7F819CDA

TSP模拟退火原理代码：

% 模拟退火算法
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
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)

TSP模拟退火全局优化工具箱工具箱中SA函数simulannealbnd代码：

使用2016a测试的代码，在2018a之后应该不需要那两个子文件

%SA1_PeaksExample.m
%% 模拟退火算法求解器
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% Objective Function
clear, close all, clc
peaks

%% Nonlinear Constraint Function
type circularConstraint

%% Define Optimization Problem
problem = createOptimProblem('fmincon',...
                             'objective',@(x) peaks(x(1),x(2)), ...
                             'nonlcon',@circularConstraint,...
                             'x0',[-1 -1],...
                             'lb',[-3 -3],...
                             'ub',[3 3],...
                             'options',optimset('OutputFcn',...
                                                @peaksPlotIterates))
                             
%% Run the solver |fmincon| from the inital point
[x,f] = fmincon(problem)                                                    

%% Use Simmulated Annealing to Find the Global Minimum

problem.solver  = 'simulannealbnd';
problem.objective = @(x) peaks(x(1),x(2)) + (x(1)^2 + x(2)^2 - 9);
problem.options = saoptimset('OutputFcn',@peaksPlotIterates,...
                             'Display','iter',...
                             'InitialTemperature',10,...
                             'MaxIter',300)

[x,f] = simulannealbnd(problem) 




%circularConstraint.m
function [c,ceq] = circularConstraint(x)
% Nonlinear constraint definition

%  Copyright (c) 2010, The MathWorks, Inc.
%  All rights reserved.

% Define nonlinear equality constraint (none)
ceq = [];

% Define nonlinear inequality constraint
% circular region with radius 3: x1^2 + x^2 -3^2 <= 0 
c = x(:,1).^2 + x(:,2).^2 - 9;





%peaksPlotIterates.m
function varargout = peaksPlotIterates(varargin)
% Output function that plots the iterates of the optimization algorithm.

%  Copyright (c) 2010, The MathWorks, Inc.
%  All rights reserved.

% Check if caller is from global or optimization toolbox
optimValues = varargin{2};
state = varargin{3};
if nargout > 1
    if isfield(optimValues,'x') % simulated annealing options,optimvalues,flag
        x = optimValues.x;
        varargout{1} = false;
        varargout{2} = x; % options field
        varargout{3} = false;
    else % pattern search optimvalues,options,flag,interval
        optimValues = varargin{1};
        state = varargin{3};
        if isfield(optimValues,'x')
            x = optimValues.x;
            varargout{1} = false;
            varargout{2} = x; % options field
            varargout{3} = false;
        else % gentic algorithm options,state,flag,interval
            x = varargin{2}.Population;
            optimValues.iteration = -1;
            varargout{1} = varargin{2};
            varargout{2} = varargin{1};
            varargout{3} = false;
        end
    end
else
    x = varargin{1};
    varargout{1} = false;
end

% Check for state
switch state
    case 'init'
        % Plot objective function surface
        PlotSurface(x,peaks(x(:,1),x(:,2)));
    case 'iter'
        if ~(optimValues.iteration == 0)
            % Update surface plot to show current solution
            PlotUpdate(x,peaks(x(:,1),x(:,2)));
        end
    case 'done'
        if ~(optimValues.iteration == 0)
            % After optimization, display solution in plot title
            DisplayTitle(x,peaks(x(:,1),x(:,2)))
        end
end

% -------------------------------------------------------------------------
% helper function PlotSurface
% -------------------------------------------------------------------------
function PlotSurface(x,z,varargin)

% Check to see if figure exists, if not create it
h = findobj('Tag','PlotIterates');
if isempty(h)
    h = figure('Tag','PlotIterates','Name','Plot of Iterates', ...
        'NumberTitle','off');
    
    % Plot the objective function
    [X,Y,Z] = peaks(100);
    zlower = -15;
    axis([-3 3 -3 3 zlower 10]);
    hold on
    surfc(gca,X,Y,Z,'EdgeColor','None','FaceColor','interp')
    xlabel('X'), ylabel('Y'), zlabel('Z')
    view([-45 30])
    shading interp
    lightangle(-45,30)
    set(findobj(gca,'type','surface'),...
        'FaceLighting','phong',...
        'AmbientStrength',.3,'DiffuseStrength',.8,...
        'SpecularStrength',.9,'SpecularExponent',25,...
        'BackFaceLighting','unlit');
    
    % Plot constraint on lower contour plot
    hc=0; k=0; r=3; N=256; % circle parameters
    t=(0:N)*2*pi/N;
    xc = r*cos(t)+hc;
    yc = r*sin(t)+k;
    
    % bounds
    ax = axis;%.*[1.1 1.1 1.1 1.1 1 1]; 
    xbound = ( ax(1):(ax(2)-ax(1))/N*4:ax(2) )';
    ybound = ( ax(3):(ax(4)-ax(3))/N*4:ax(4) )';
    len = length(xbound);
    xbox = [xbound;  xbound(end)*ones(len-1,1);
            xbound(end-1:-1:1); xbound(1)*ones(len-2,1)]; 
    ybox = [ybound(1)*ones(len,1); ybound(2:end);
            ybound(end)*ones(len-1,1); ybound(end-1:-1:2)];
    
    boxCon = [(1:length(xbox)-1)' (2:length(ybox))'; length(xbox) 1];
    cirCon = [(1:length(xc)-1)' (2:length(yc))'; length(x) 1] + length(xbox);
    
    warning off
    DT = DelaunayTri([[xbox(:); xc(:)] [ybox(:); yc(:)]], [boxCon; cirCon]);
    warning on
    inside = inOutStatus(DT);
    cx = caxis;
    trisurf(DT(inside,:),DT.X(:,1),DT.X(:,2),...
           zlower*ones(size(DT.X(:,1))),'EdgeColor','none',...
           'FaceColor',[0.9 0.9 0.9]);
    caxis(cx)
    hold off
    
    % colors to use for multiple staring points
    ms.index = 1;
    ms.Colors = ['rgbcmyk'];
    set(h,'UserData',ms);
end

PlotUpdate(x,z)
if nargin > 2
    DisplayTitle(x,z,varargin{1})
else
    DisplayTitle(x,z,'Initial')
end

% -------------------------------------------------------------------------
% helper function PlotUpdate
% -------------------------------------------------------------------------
function PlotUpdate(x,z)

% Check to see if figure exists, if not, create
h = findobj('Tag','PlotIterates');
if isempty(h)
    PlotSurface(x,z,'Current')
    h = gcf;
end
% Update Plot with New Point
figure(h)
ms = get(h,'UserData');
hold on
spts = findobj('Tag','SurfacePoints');
if isempty(spts)
    spts = plot3(x(:,1),x(:,2),z*1.02,'MarkerFaceColor',ms.Colors(ms.index),...
        'MarkerSize',10,...
        'Marker','diamond',...
        'LineStyle','none',...
        'Color',ms.Colors(ms.index));
    set(spts,'Tag','SurfacePoints');
else
    set(spts, 'XData', [get(spts,'XData'),x(:,1)']);
    set(spts, 'YData', [get(spts,'YData'),x(:,2)']);
    set(spts, 'ZData', [get(spts,'ZData'),z']);
end

ax1 = findobj('Tag','LowerContour');
if isempty(ax1)
    if isvector(x)
        mk = '.-';
    else
        mk = '.';
    end
    ax1 = plot3(x(:,1),x(:,2),min(get(gca,'ZLim'))*ones(size(x(:,1))),...
        [ms.Colors(ms.index),mk],'MarkerSize',16);
    set(ax1,'Tag','LowerContour');
else
    set(ax1, 'XData', [get(ax1,'XData'),x(:,1)']);
    set(ax1, 'YData', [get(ax1,'YData'),x(:,2)']);
    set(ax1, 'ZData', [get(ax1,'ZData'),...
                       min(get(gca,'ZLim'))*ones(size(x(:,1)'))]);
end

DisplayTitle(x,z,'Current')

% -------------------------------------------------------------------------
% helper function DisplayTitle
% -------------------------------------------------------------------------
function DisplayTitle(x,z,varargin)
% colors to use for iterates

% Check to see if figure exists, if not, create
h = findobj('Tag','PlotIterates');
if isempty(h)
    PlotUpdate(x,z)
    h = gcf;
end

% Update Plot Title
if nargin < 3
    varargin{1} = 'Final';
end

ms = get(h,'UserData');
[mz,indx] = min(z);
switch lower(varargin{1})
    case 'current'
        str = 'Current';
    case 'initial'
        str = 'Initial';
        text(x(indx,1),x(indx,2),z(indx)*1.1,'\bf Start','Color',...
            ms.Colors(ms.index))
    case 'final'
        str = 'Final';
        text(x(indx,1),x(indx,2),z(indx)*1.1,'\bf End','Color',...
            ms.Colors(ms.index))
        ms.index = ms.index+1;
        set(h,'UserData',ms);
        ax1 = findobj('Tag','LowerContour');
        set(ax1,'Tag',['LowerContour',num2str(ms.index)]);
        spts = findobj('Tag','SurfacePoints');
        set(spts,'Tag',['SurfacePoints',num2str(ms.index)]);
end

str = sprintf('%s x = [%6.4f %6.4f]',str, x(indx,1),x(indx,2));
figure(h), title(str)
drawnow;

4.2.4蚁群算法(局部最优)

文档：蚁群算法.note

链接：http://note.youdao.com/noteshare?id=debe2b99081cdec49d9c694b2d3679a1&sub=3E729AEDEBCD4147AF366015B26641B1

所需数据：

链接：https://pan.baidu.com/s/10PprCHy3_nlhzdZtNzquKg

提取码：go9r

蚁群算法求解TSP

%数据chapter9——Chap9_citys_data.xlsx

%% 蚁群算法及Matlab实现——TSP问题
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 数据准备
% 清空环境变量
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;      %设定的对角矩阵修正值
        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; 
      % 构建解空间
      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]。

直观分析蚂蚁数量对蚁群算法的影响

我们希望三个指标越小越好，判断出m/M约为1.5，即蚂蚁数量是城市的1.5倍

%% 绘制蚁群算法中蚂蚁数目与最短路程、收敛迭代次数、程序执行时间关系图
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 准备环境与数据
clc
clear 

% 输入数据
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('执行时间');

最佳旅游方案

文档：最佳旅游方案题目.note

链接：http://note.youdao.com/noteshare?id=d55a86c28fc00b7856ffc298a53e05c0&sub=B55C4814162D452594CE51DFFC8EE68D

%Ch9_spots_data.xlsx
%% 第8章 蚁群算法及Matlab实现——TSP问题
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 数据准备
% 清空环境变量
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]。

5.连续模型的matlab求解方法

指模型是连续函数的模型，建模方法是微分方程建模。

确定自变量、未知函数、参数。并确定坐标系
找出这些量满足的基本规律(物理、几何)
运用规律列出方程和定解条件。

matlab用于求解析解，以及数值模拟(给出比变量之间的函数形式)，常见的方法有：

dsolve求解常见的微分方程解析解
ODE家族求解器求数值解
专用求解器

5.1常规微分方程求解

文档：常规.note

链接：http://note.youdao.com/noteshare?id=b1c10869ac1a5333eb3ddaa38c8a79f7&sub=85B656E31C3544F5AA15A19EA97B4CC3

5.2ODE家族求解器

文档：ode家族求解器.note

链接：http://note.youdao.com/noteshare?id=dd7e9613d552a2d9156e0e12f115fc9d&sub=B474242961C441F08478E4943E9B479E

5.3专用求解器

文档：专用求解器.note

链接：http://note.youdao.com/noteshare?id=5cbbc575aa94966d13e6352123b5b9e2&sub=4E96856A74574FE0A2FB43CD66E9D113

%版本问题，未能复现

%% Wave Equation on a Square Domain
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% Problem Definition设置参数
c = 1;
a = 0;
f = 0;
m = 1;

%% Geometry定义波的空间位置
numberOfPDE = 1;
model = createpde(numberOfPDE);
geometryFromEdges(model,@squareg);
pdegplot(model,'EdgeLabels','on'); 
ylim([-1.1 1.1]);
axis equal
title 'Geometry With Edge Labels Displayed';
xlabel x
ylabel y

%% Specify PDE Coefficients定义微分方程的系数和边界条件
specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);
applyBoundaryCondition(model,'dirichlet','Edge',[2,4],'u',0);
applyBoundaryCondition(model,'neumann','Edge',([1 3]),'g',0);

%% Generate Mesh%定义有限元网格
generateMesh(model);
figure
pdemesh(model);
ylim([-1.1 1.1]);
axis equal
xlabel x
ylabel y

%% Create Initial Conditions定义初始条件
u0 = @(location) atan(cos(pi/2*location.x));
ut0 = @(location) 3*sin(pi*location.x).*exp(sin(pi/2*location.y));
setInitialConditions(model,u0,ut0);

%% Define Solution Times方程求解
n = 31;%求解次数
tlist = linspace(0,5,n);
model.SolverOptions.ReportStatistics ='on';
result = solvepde(model,tlist);
u = result.NodalSolution;

%% Animate the Solution模型的数值仿真
figure
umax = max(max(u));
umin = min(min(u));
for i = 1:n
    pdeplot(model,'XYData',u(:,i),'ZData',u(:,i),'ZStyle','continuous',...
                  'Mesh','off','XYGrid','on','ColorBar','off');
    axis([-1 1 -1 1 umin umax]); 
    caxis([umin umax]);
    xlabel x
    ylabel y
    zlabel u
    M(i) = getframe;
end

6.评价模型的求解

6.1线性加权法

适用于指标相互独立。

最重要的是确定权重。

文档：matlab评价模型.note

链接：http://note.youdao.com/noteshare?id=80dc635082abefd5396d3559a20ee5a7&sub=A1B359AA15E446EC9D10A883F8C1AB05

数据：

链接：https://pan.baidu.com/s/1v-aeHPp2ZcRDAnqK9eHIQQ

提取码：ym92

评价的对象是股票，已知股票的指标及其表现。1代表上涨、0代表一般、-1代表下跌。利用模型评价新的股票

%% 多因子选股模型
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%% 导入数据
clc, clear all, close all
s = dataset('xlsfile', 'SampleA1.xlsx');

%% 多元线性回归(对所有变量回归)
myFit = LinearModel.fit(s);
disp(myFit) 
sx=s(:,1:10);
sy=s(:,11);
n=1:size(s,1);
sy1= predict(myFit,sx);
figure
plot(n,sy, 'ob', n, sy1,'*r')
xlabel('样本编号', 'fontsize',12)
ylabel('综合得分', 'fontsize',12)
title('多元线性回归模型', 'fontsize',12)
set(gca, 'linewidth',2)

%% 逐步回归(因子筛选，得到优选因子)
myFit2 = LinearModel.stepwise(s);
disp(myFit2)
sy2= predict(myFit2,sx);
figure
plot(n,sy, 'ob', n, sy2,'*r')
xlabel('样本编号', 'fontsize',12)
ylabel('综合得分', 'fontsize',12)
title('逐步回归模型', 'fontsize',12)
set(gca, 'linewidth',2)
%少了5个单一因子，多了5个组合因子

6.2层次分析法

%% AHP法权重计算MATLAB程序
% 《MATLAB数学建模方法与实践》(《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)]);

7.机理建模的matlab实现方法

文档：matlab机理建模.note

链接：http://note.youdao.com/noteshare?id=840ba253406837dab35c37d19ce2c9a5&sub=773CAAAC7AA34C5181A3F680B89032AE

7.1推导机理建模(略)

7.2元胞自动机——仿真

CA,也称为细胞自动机。所研究的问题是一个系统问题，系统由若干一个或几个不同类的对象组成。如滴滴打车2015、开发小区2016。

实现步骤：

定义元胞的初始状态
定义系统内元胞的变化规则
设置仿真时间，输出仿真结果

比如2015打车问题中，

元胞

就是打车人和出租车；

更新规则

是当打车人发出打车信号，周边出租车响应规则；

系统输出

则是评价指标。

以下给出一个元胞自动机的框架，需要建模时灵活定义元胞，更新规则，系统输出。

%% 元胞自动机（CA）MATLAB实现程序
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
clc, clf,clear 

%% 界面设计（环境的定义）
plotbutton=uicontrol('style','pushbutton',...
   'string','Run', ...
   'fontsize',12, ...
   'position',[100,400,50,20], ...
   'callback', 'run=1;');

% 定义 stop button
erasebutton=uicontrol('style','pushbutton',...
   'string','Stop', ...
   'fontsize',12, ...
   'position',[200,400,50,20], ...
   'callback','freeze=1;');

% 定义 Quit button
quitbutton=uicontrol('style','pushbutton',...
   'string','Quit', ...
   'fontsize',12, ...
   'position',[300,400,50,20], ...
   'callback','stop=1;close;');

number = uicontrol('style','text', ...
    'string','1', ...
   'fontsize',12, ...
   'position',[20,400,50,20]);

%%  元胞自动机的设置
n=128;
z = zeros(n,n);
cells = z;
sum = z;
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;
cells = (rand(n,n))<.5 ;
imh = image(cat(3,cells,z,z));
axis equal
axis tight

% 元胞索引更新的定义
x = 2:n-1;
y = 2:n-1;

% 元胞更新的规则定义
stop= 0; % wait for a quit button push
run = 0; % wait for a draw 
freeze = 0; % wait for a freeze
while (stop==0) 
    if (run==1)
        %nearest neighbor sum
        sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
            cells(x-1, y) + cells(x+1,y) + ...
            cells(x-1,y-1) + cells(x-1,y+1) + ...
            cells(3:n,y-1) + cells(x+1,y+1);
        % The CA rule
        cells = (sum==3) | (sum==2 & cells);       
        % draw the new image
        set(imh, 'cdata', cat(3,cells,z,z) )
        % update the step number diaplay
        stepnumber = 1 + str2num(get(number,'string'));
        set(number,'string',num2str(stepnumber))
    end
    if (freeze==1)
        run = 0;
        freeze = 0;
    end
    drawnow  %need this in the loop for controls to work
end

8.出租车补贴方案优化问题2015B

融合了数据分析、优化、机理建模、评价。

文档：2015B出租车补贴方案优化.note

链接：http://note.youdao.com/noteshare?id=53dddb4848f03e26ee0d9264e621e700&sub=49AF6D20D97443EC8975D659486DA032

链接：https://pan.baidu.com/s/1e2Zt6Ji2Do3-dSOIIk5v9w

提取码：s2es

%运行有问题line89
%P1_taxi.m
%% 出租车补贴方案仿真程序
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
clc, clear, close all
%% 数据结构设计
% passengers:
% [出发点横坐标，出发点纵坐标，目的地横坐标，目的地纵坐标，出行里程]
% 即[xs,ys,xd,yd,l]
% taxis
% [出租车位置横坐标，出租车位置纵坐标，出租车被占用里程]
% 即[x_taxi,y_taxi,lo]
r_valid = 2/10;%出租车有效覆盖半径
xmax = 111*cos(pi*34/180)*1.4;
ymax = 0.7*111;
xmax = xmax/10;
ymax = ymax/10;
psnger_total = 80;
taxi_total = 152;
%先生成5000个出发点
for i = 1:psnger_total
    passengers(i,:) = gen_passenger();
end
for i = 1:taxi_total
    taxis(i,:) = gen_taxi();
end
figure
scatter(taxis(:,1)*10,taxis(:,2)*10)
xlabel('x(km)')
ylabel('y(km)')
all_B = [];
all_K = [];
for i = 1:200
    %% 首先更新出租车状态
    lc = taxis(:,3) - 0.01;%出租车被占用里程
    lc(lc < 0) = 0;
    taxis(:,3) = lc;
    %空车随机一个方向前进0.01
    valid_lines = find( lc == 0 );
    all_K = [all_K,1-length(valid_lines)/taxi_total]; 
    for m = 1:length(valid_lines)
        k = valid_lines(m);
        while(1)
            degree = 2*pi*rand();%出行方向
            new_x = taxis(k,1) + 0.01.*cos(degree);
            new_y = taxis(k,2) + 0.01.*sin(degree);
            if(new_x>=0 && new_x<=xmax && new_y>=0 && new_y<=ymax)
                taxis(k,1:2) = [new_x,new_y];
                break
            end
        end
    end
    
    %% 乘客加入系统
    add_passengers_total = 4;%round(normrnd(10,3));
    add_passengers = zeros(add_passengers_total,5);
    for n = 1:add_passengers_total
        add_passengers(n,:) = gen_passenger();
    end
    passengers = [passengers;add_passengers];    
    %% 计算各乘客视野内出租车数目  
    for j = 1:length(passengers)  
        p = passengers(j,:);
        if isnan(p(1))
            continue
        end
        temp_taxis = taxis; 
        %被占用的出租车不参与打车
        invalid_lines = find(temp_taxis(:,3)>0);
        temp_taxis(invalid_lines,:) = nan;
       %% 然后是乘客乘车
        r = sqrt((temp_taxis(:,1)-p(1)).^2+(temp_taxis(:,2)-p(2)).^2);
        taxi_num = find(r<r_valid);%视野范围内的车辆 
        if isempty(taxi_num)%视野内没有车，下一位乘客
            continue;
        else
            %随机选一辆乘坐
            index = round(rand()*(length(taxi_num)-1))+1;
            taxi_num = taxi_num(index);
            taxis(taxi_num,3) = p(5) + sqrt((taxis(taxi_num,1)-p(1))^2+(taxis(taxi_num,2)-p(2))^2);%此乘客p乘坐的出租车被占用
            taxis(taxi_num,1) = p(3);taxis(taxi_num,2) = p(4);%将其更新到目的地
            passengers(j,:) = nan;%更新乘客状态，上车的乘客变为nan，移出系统
        end
    end
        all_B = [all_B,calcu_b(passengers,taxis)];
end
%% 结果可视化
figure
hold on
scatter(taxis(:,1)*10,taxis(:,2)*10,'g*')
legend('初始位置','一段时间后位置')

%20次演化后才得到平时状态，故只保留20次之后的数据
pos = (20:length(all_B));
figure
plot((pos-20)/4.5,all_B(pos));
xlabel('时间(分钟)')
ylabel('数目不平衡度')
figure
plot((pos-20)/4.5,all_K(pos));
xlabel('时间(分钟)')
ylabel('里程利用率')

figure
res = sqrt( ((all_K-0.66)./0.66).^2 + (all_B-1).^2 );
plot((pos-20)/4.5,res(pos));
xlabel('时间(分钟)')
ylabel('供需不平衡度')



%calcu_b.m计算供求比率
function [ B ] = calcu_b( passengers,taxis )
invalid_lines = find(isnan(passengers(:,1)));
passengers(invalid_lines,:) = [];
%划分网格
area_d = zeros(9,14);%需
[h,w] = size(passengers);
for i = 1:h
    y = floor(passengers(i,1))+1;
    x = floor(passengers(i,2))+1;
    area_d(x,y) = area_d(x,y) + 1;
end
area_s = zeros(9,14);%供
[h,w] = size(taxis);
for i = 1:h
    y = floor(taxis(i,1))+1;
    x = floor(taxis(i,2))+1;
    area_s(x,y) = area_s(x,y) + 1;
end
%供需相等
s1 = sum(area_s(area_s == area_d));
d1 = sum(area_d(area_s == area_d));
%供大于求
s2 = sum(area_s(area_s > area_d));
d2 = sum(area_d(area_s > area_d));
%供小于求
s3 = sum(area_s(area_s < area_d));
d3 = sum(area_d(area_s < area_d));
d = d1+d2+d3;
if d3 == s3
    B = s1/d + s2/d;
    return
end
B = s1/d + s2/d + d3/d*d3/s3;
end


%gen_passenger.m计算乘客位置参数
function [ ret ] = gen_passenger()
% 计算乘客位置参数
% xmax = 111*cos(pi*34/180)*1.4 = 129
% ymax = 0.7*111 = 78

xmax = 111*cos(pi*34/180)*1.4/10;
ymax = 0.7*111/10;
ux = xmax/2;uy = ymax/2;
sigmax = xmax/2/3;sigmay = ymax/2/3;
    th = 6.5/((pi/2)^0.5);
    while(1)
        xs = normrnd(ux,sigmax);%出发点横坐标
        ys = normrnd(uy,sigmay);%出发点纵坐标
        if xs<0 || xs>xmax || ys<0 || ys > ymax
            continue
        end
        d_go = sqrt(-2*th^2*log(1-rand()))/10;%出行距离
        degree = 2*pi*rand();%出行角度
        xd = xs + d_go.*cos(degree);
        yd = ys + d_go.*sin(degree);
        if(xd>=0 && xd<=xmax && yd>=0 && yd<=ymax)
            ret = [xs,ys,xd,yd,d_go];
            break
        end

    end
end


%gen_taxi.m计算出租车位置参数
function [ ret ] = gen_taxi( input_args )
% 计算出租车位置参数
%   Detailed explanation goes here
    xmax = 111*cos(pi*34/180)*1.4/10;
    ymax = 0.7*111/10;
    ux = xmax/2;uy = ymax/2;
    sigmax = xmax/2/3;sigmay = ymax/2/3;
    while(1)
        x = normrnd(ux,sigmax);%出发点横坐标
        y = normrnd(uy,sigmay);%出发点纵坐标
        if x>=0 && x<=xmax && y>=0 && y <= ymax
            break
        end    
    end
    ret = [x,y,0];
end

9.开放小区对道路通行影响的问题2016

涉及数据分析、评价、机理分析、仿真。

文档：开放小区对道路通行影响的问题2016.not…

链接：http://note.youdao.com/noteshare?id=fdf0d3da3677db039b806cd3bea6d09f&sub=551C58934D194753B33340DE9E404703

链接：https://pan.baidu.com/s/1VcphutEb7BdacVKLw5JQWw

提取码：kjwv

matlab实现车辆通行模型，仿真如下

%% 2016.09  全国大学生数学建模竞赛 B 题 小区开放对道路通行的影响
% 《MATLAB数学建模方法与实践》(《MATLAB在数学建模中的应用》升级版)，北航出版社，卓金武、王鸿钧编著. 
%P19_1.m
%% 参数设置
clc;
clear;
Tmax = 200;        %考虑的时间的上限
carnum = 0;         %carnum  是当前通过的车的总数目
fieldCarNum = 0;    %fieldCarNum 进入小区的车辆数
beta = 5;           %1-alpha / 10,为出现车的概率，体现车流量 
fieldCapure = 800;  %小区内道路的承载量
fieldDistance = 80; %小区内道路的逻辑长度
fRpb = 0.8;         %行人自行车修正系数
R = 0.5;            %主路进入小区的比例
L1 = 0.005;         %车身长度,km
h = 0.0025;         %标准饱和车头时距
alpha = 2*pi / 3;   %转弯角度
mu = 0.18;          %横向力系数
v = 50;             %车辆速度
t_avg = L1 / v ;    %畅通情况下车辆直行通过计算截面的平均耗时
ER = (L1 + alpha * R)/(h * sqrt(127 * R * mu));
                    %右转车转换系数
fR = 100/(100 + R * (ER - 1));
                    %右转修正系数
theta = 0;          %车辆入口延时的影响因子
delay = zeros(Tmax, 1); %记录入口延时
T = 20; %路灯周期
Tg = 5; %绿灯时间
Car = cell(carnum,1);  %Car 是所有的车的集合
Car0 = cell(carnum,1); %Car 是所有的车的集合
car = struct('road',0,'distance',0,'state',0); %road 是当前车所在的道路，distance表示在这条路上的位置,state表示是否在区域里，1在里面，0在外面
t = 0; %当前时间
Dt = cell(t,1); %每个时刻的平均延误
Dt0 = cell(t,1); %每个时刻的平均延误
Light = 0; % Light 表示当前灯是红灯0，还是绿灯1

roadnum = 4;  %roadnum 是所有的道路数目
Outroad = []; % 与出口相连的路
Outroad(1) = 4;
Inroad = [2,3,]; %与入口相连的路
FieldRoad = []; %小区内的路
FieldRoad(1) = 3;

Troad = zeros(roadnum,1);%每条路上的路灯周期
Tgroad = zeros(roadnum,1);%每条路上绿灯时间
RoadMap = zeros(roadnum,roadnum); %路的可达性矩阵
RoadMap0 = zeros(roadnum,roadnum); %路的可达性矩阵
Roadcapture = zeros(roadnum,1); %Roadcapture 路的设计承载量
Roadcarnum = zeros(roadnum,1); %Roadcarnum 路当前有的车数量
Roadcarnum0 = zeros(roadnum,1); %Roadcarnum 路当前有的车数量
Roaddistance = zeros(roadnum,1);%Roaddistance 路的距离
Roaddt = zeros(roadnum,1); %Roaddt 每条路上的平均延迟时间


Roadcapture = [1000,1000,1000,1];
Roadcapture(FieldRoad) = fieldCapure;
Roadcapture = Roadcapture';
Roaddistance = [100,100,100,1];
Roaddistance(FieldRoad) = fieldDistance;

Troad(:) = 20;
Tgroad(:) = 5;

RoadMap = [0,0,0,0;
           1,0,0,0;
           1,0,0,0;
           0,1,1,0]; %%小区内部交通可达性矩阵
RoadMap0 = [0,0,0,0;
            1,0,0,0;
            0,0,0,0;
            0,1,0,0];

myres = zeros(500,4);
Flag = [0,1,0,0]'; 
mainflow = 3600 * (1 - beta / 10);         %主路的车流量
%% 仿真过程
for t = 1:Tmax
    %%%%%%判断当前是否有车到来0没有，1有
    Dt{t} = 0;
    Dt0{t} = 0;
    ra = rand();
    if ra >= beta / 10
        ra = 1;
    else
        ra = 0;
    end
    %%%%%%增加车的数量,更新车的情况
    if ra == 1
        carnum = carnum+1;
        car.road = 1;
        car.distance = Roaddistance(1);
        car.state = 1;
        Car{carnum} = car;
        Car0{carnum} = car;
        Roadcarnum(1) = Roadcarnum(1)+1;
        Roadcarnum0(1) = Roadcarnum0(1)+1;
    end
   
    %%%%%%判断当前红绿灯情况
    Light = mod(t,T);
    if Light <= Tg
        Light = 1;
    else
        Light = 0;
    end
    
    for cari = 1:carnum
        if(Car0{cari}.state == 1)
        %%不考虑小区开放时  
        if Light == 1  %%%绿灯时
            [nextroad,nextdistance,nextstate] = nextdir(cari,RoadMap0,Car0,Roaddistance,Outroad,Roadcarnum0); %%%函数nextdir返回cari这辆车下一次所在的路的在路上的距离
        else %%%%红灯时
            [nextroad,nextdistance,nextstate] = nextdir_red(cari,RoadMap0,Car0,Roaddistance,Outroad,Roadcarnum0);
        end
        Roadcarnum0(Car0{cari}.road) =  Roadcarnum0(Car0{cari}.road)-1;
        Car0{cari}.road = nextroad;
        Car0{cari}.distance = nextdistance;
        Car0{cari}.state = nextstate;
        if nextstate == 1
            Roadcarnum0(nextroad) =  Roadcarnum0(nextroad)+1;  %%%%%%当前这条路上的车的数量
        end
        end
        %%考虑小区开放时
    %   disp(strcat('car', num2str(cari)));
        if(Car{cari}.state == 1)
        if Light == 1  %%%绿灯时
            [nextroad,nextdistance,nextstate] = nextdir(cari,RoadMap,Car,Roaddistance,Outroad,Roadcarnum); %%%函数nextdir返回cari这辆车下一次所在的路的在路上的距离
        else %%%%红灯时
            [nextroad,nextdistance,nextstate] = nextdir_red(cari,RoadMap,Car,Roaddistance,Outroad,Roadcarnum);
        end
        myres(t, 1) = nextroad;
        myres(t, 2) = nextdistance;
        myres(t, 3) = nextstate;
        myres(t, 4) = cari;
        Roadcarnum(Car{cari}.road) =  Roadcarnum(Car{cari}.road)-1;
        Car{cari}.road = nextroad;
        Car{cari}.distance = nextdistance;
        Car{cari}.state = nextstate;
        if nextstate == 1
            Roadcarnum(nextroad) =  Roadcarnum(nextroad)+1;  %%%%%%当前这条路上的车的数量
        end
        end
    end
    if( ra == 1 && ismember(nextroad, FieldRoad) )
            fieldCarNum = fieldCarNum + 1;
            theta = theta + 1;
    else
        if(theta >= 0.3)
            theta = theta - 0.3;
        else
            theta = 0;
        end
    end
    if(carnum > 0)
        R = R * 0.5 + 0.5 * fieldCarNum / carnum;
        ER = (L1 + alpha * R)/(h * sqrt(127 * R * mu));%右转车转换系数
        fR = 100/(100 + R * (ER - 1));%右转修正系数
        x = ( mainflow) / Roadcapture(1); %道路饱和度
        C1 = x * (fRpb + fR);               %主流向通行能力
        NDT = (1/C1 - t_avg) * (1/t_avg + 1) * t_avg / 2; %信号交叉路口平均延误
        Dt{t} = theta * NDT;
        delay(t) = Dt{t};
    end
    %%%%%%计算每条路的平均延误时间
    dt = ((0.5*Troad).*(1-Tgroad./Troad))./(1-min(1,Roadcarnum.*(720./Roaddistance')./Roadcapture).*(Tgroad./Troad)).*Flag;
    Dt{t} = Dt{t} + sum(dt);
        
    dt0 = ((0.5*Troad).*(1-Tgroad./Troad))./(1-min(1,Roadcarnum0.*(720./Roaddistance')./Roadcapture).*(Tgroad./Troad)).*Flag;
    Dt0{t} = Dt0{t} + sum(dt0);
end
%% 画图
matDt = cell2mat(Dt);
matDt0 = cell2mat(Dt0);
plot([1:Tmax],matDt,'-','markersize',5);
hold on;
plot([1:Tmax],matDt0,'-','markersize',5);
legend('小区开放','小区未开放')
grid on;
dtDelta = cell2mat(Dt0) - cell2mat(Dt);
save(strcat(num2str(alpha), 'dtDelta_7.mat'), 'dtDelta');
save(strcat(num2str(alpha), 'Dt_7.mat'), 'Dt');
save(strcat(num2str(alpha), 'Dt0_7.mat'), 'Dt0');


%chooseRoad.m
function [nextroad] = chooseRoad(cari,currentRoad, RoadMap,Car,Roaddistance, RoadCarNum)
%% function [nextroad] = chooseRoad(cari,currentRoad, RoadMap,Car,Roaddistance, RoadCarNum)
%Input arguments：
%               cari：1 x 1 int  车的序号
%               currentRoad:1 x 1 int 当前所在的路的标号
%               RoadMap: n x n array 路的可达性矩阵
%               Car：n x 1 cell 当前各辆车的状态
%               Roaddistance：n x 1 int 每条路的长度
%               RoadCarNum: n x 1 int 当前每条路上有多少车
%Output arguments：
%               nextroad: 1 x 1 int 选择的下一条路的标号
%%
k1 = 1;%两个参数
k2 = 1;
nextroad = 0;
min = intmax();
roadNum = length(RoadMap);
% disp(strcat('currentroadNum', num2str(currentRoad)));
myNeighbor = zeros(roadNum);
myNeighborNum = 0;
for i = 1 : roadNum
    if(RoadMap(i, currentRoad) == 0) %道路不通的情况
        continue;
    end
    myNeighborNum = myNeighborNum + 1;
    myNeighbor(myNeighborNum) = i;
    % disp(strcat('roadNum', num2str(i)));
    
    % 根据最优选择
    myvalue = 0;
    for j = 1 : (cari - 1)
        if(Car{j}.state == 1 && Car{j}.road == i && Car{j}.distance == 1)
            %判断走这条路是否需要等
            myvalue = intmax() - 100000;
            break;
        end
    end
    %计算每条路的指标
    myvalue = myvalue + k1 * Roaddistance(i) + k2 * RoadCarNum(i);
    
    if(myvalue < min)
        min = myvalue;
        nextroad = i;
    end
    
    %     disp(i);
    %    disp(myvalue);
    
end
% 随即选择
% nextroad = myNeighbor(ceil(rand() * (myNeighborNum - 1)) + 1);


%nextdir.m
function [nextroad,nextdistance,nextstate]=nextdir(cari,RoadMap,Car,Roaddistance,Outroad, RoadCarNum)
%% function [nextroad,nextdistance,nextstate]=nextdir(cari,RoadMap,Car,Roaddistance,Outroad, RoadCarNum)
% 函数nextdir返回cari这辆车这一秒之后的状态
%Input arguments:
%              cari：1 x 1 int 车的序号
%              RoadMap: n x n array 路的可达性矩阵
%              Car：n x 1 cell 当前各辆车的状态
%              Roaddistance：n x 1 int 每条路的长度
%              Outroad：k x 1 array 所有可能出去的路 k 是出口数目
%              RoadCarNum: n x 1 int 当前每条路上有多少车
%Output arguments:
%              nextroad: 1 x 1 int 1s后所在的路
%              nextstate: 1 x 1 int dual value(0,1) 1s后，该车是否还在区域内
%              nextdistance：1 x 1 int 1s后在路上距路入口的逻辑距离
%% 
    car = Car{cari}; % 当前正在考虑的车
    nextroad = car.road;
    nextdistance = car.distance;
     if car.distance >= Roaddistance(Car{cari}.road) && ismember(car.road,Outroad);
        nextstate = 0;
        return;
    else
        nextstate = 1;
    end
    nextdistance = nextdistance + 1;
    if(Roaddistance(nextroad) < nextdistance)
        nextroad = chooseRoad(cari,nextroad,RoadMap,Car,Roaddistance, RoadCarNum); %%%交叉路口选择路的方向
        nextdistance = 1;
    end
    for i = 1 : (cari - 1)
           if(Car{i}.state == 1 && Car{i}.road == nextroad && Car{i}.distance == nextdistance)
               % 车辆不能走，
               nextroad = car.road;
               nextdistance = car.distance;
           end
    end
    

%nextdir_red.m
function [nextroad,nextdistance,nextstate]=nextdir_red(cari,RoadMap,Car,Roaddistance,Outroad, RoadCarNum)
%% 函数nextdir返回cari这辆车这一秒之后的状态
%Input arguments：
%                cari：1 x 1 int 车的序号
%                RoadMap: n x n array 路的可达性矩阵
%                Car：n x 1 cell 当前各辆车的状态
%                Roaddistance：n x 1 int 每条路的长度
%                Outroad：k x 1 int 所有可能出去的路 k是出口数
%                RoadCarNum: n x 1 int 当前每条路上有多少车
%Output arguments:
%                nextroad: 1 x 1 int 1s后所在的路
%                nextdistance：1 x 1 int 1s后在路上距路入口的逻辑距离
%                nextstate: 1 x 1 int dual value(0,1)1s后，该车是否还在区域内
%% 
    car = Car{cari}; % 当前正在考虑的车
    nextstate = 1;
    nextroad = car.road;
    nextdistance = car.distance;
    if car.distance >= Roaddistance(Car{cari}.road) && ismember(car.road,Outroad);
        return;
    end
    nextdistance = nextdistance + 1;
    if(Roaddistance(nextroad) < nextdistance)
        nextroad = chooseRoad(cari,nextroad,RoadMap,Car,Roaddistance, RoadCarNum); %%%交叉路口选择路的方向
        nextdistance = 1;
    end
    for i = 1 : (cari - 1)
           if(Car{i}.state == 1 && Car{i}.road == nextroad && Car{i}.distance == nextdistance)
               % 车辆不能走，
               nextroad = car.road;
               nextdistance = car.distance;
           end
    end

原文链接：https://blog.csdn.net/qq_44809326/article/details/108380242