数学建模

作者: 恰似一碗咸鱼粥 | 来源:发表于2019-07-27 20:58 被阅读0次

1.启发式算法

它主要包括禁忌搜索，模拟退火，遗传算法，神经网络，蚁群算法

模拟退火算法

流程图

Metropolis准则：对于目标函数E，若 $E_{i+1}<E_i$ ，则接收E_{i+1}，否则依 $e^{-(E_{i+1}-E_i)/KT}$ 概率接受。
模拟退火与物理退火的对应关系：

解--粒子状态
目标函数--能量
最优解--能量最低态
设定初温--加热过程
扰动--热涨落
Metropolis采样--热平衡
控制参数的下降--冷却

关于初始解：初始解不宜太好，否则很难从邻域跳出，并且邻域要尽可能小。

降温方式：
1.经典模拟退火： $T(t)=\frac{T_0}{log(1+t)}$
2.快速模拟退火： $T(t)=\frac{T_0}{1+t}$
3.其他： $T(t+\Delta t)=\alpha T(t)$

花费函数COST：
1.一般直接用目标函数构造，或目标函数的倒数/相反数。
2.终止条件：理论上温度要降为0时才终止退火算法，但实际温度较低时接受的概率就接近0了。

旅行商问题(TSP)

有一个商人要拜访n个城市， $d_{ij}$ 表示两城市间的距离， $x_{ij}$ 为0，1变量，表示拜访路径是否包含 $d_{ij}$
限制：每个城市必须且只能拜访一次，且最后要回到原来的城市，即出度与入度都为1。
目标：路径之和最小，即：
$min\sum_{i!=j}d_{ij}x_{ij}$

解：
1.设置初始解S0
[1,...,i,..,j...,n,1]
2.产生邻解S1
[1,...,j...,i,...,n,1]//将i与j位置互换
3.定义COST函数
$\sum_{i=1}^n dist(S_{(i)},S_{(i+1)})$
4.设置温度，降温
$T_0=1000,T_{t+dt}=\alpha T_t$

计算各个城市两两之间距离的函数：

function dis = distancematrix(city)
% DISTANCEMATRIX
% dis = DISTANCEMATRIX(city) return the distance matrix, dis(i,j) is the 
% distance between city_i and city_j

numberofcities = length(city);
R = 6378.137; % The radius of the Earth
for i = 1:numberofcities
    for j = i+1:numberofcities
        dis(i,j) = distance(city(i).lat, city(i).long, ...
                            city(j).lat, city(j).long, R);
        dis(j,i) = dis(i,j);
    end
end

制造扰动的函数：

function route = perturb(route_old, method)
% PERTURB
% route = PERTURB(route_old, method) generate randomly a neighbouring route by
% perturb old route. perturb methods:
%                        ___________            ___________         
%     1. reverse:   [1 2 3 4 5 6 7 8 9] -> [1 2 8 7 6 5 4 3 9]
%                        _         _            _         _
%     2. swap:      [1 2 3 4 5 6 7 8 9] -> [1 2 8 4 5 6 7 3 9]

route = route_old;
numbercities = length(route);
city1 = ceil(numbercities*rand);
city2 = ceil(numbercities*rand);
switch method
    case 'reverse'
        citymin = min(city1,city2);
        citymax = max(city1,city2);
        route(citymin:citymax) = route(citymax:-1:citymin);
    case 'swap'
        route([city1, city2]) = route([city2, city1]);
end

计算总距离的函数：

function d = totaldistance(route, dis)
% TOTALDISTANCE
% d = TOTALDISTANCE(route, dis) calculate total distance of a route with
% the distance matrix dis.

d = dis(route(end),route(1)); % closed path
for k = 1:length(route)-1
    i = route(k);
    j = route(k+1);
    d = d + dis(i,j);
end

main主函数

clear;clc;

load china; % geographic information
plotcities(province, border, city); % draw the map of China

numberofcities = length(city);      % number of cities
% distance matrix: dis(i,j) is the distance between city i and j.
dis = distancematrix(city);   


temperature = 1000;                 % Initialize the temperature.
cooling_rate = 0.94;                % cooling rate
iterations = 1;                     % Initialize the iteration number.

% Initialize random number generator with "seed". 
rand('seed',0);                    

% Initialize the route by generate a sequence of random
route = randperm(numberofcities);
% This is objective function, the total distance for the routes.
previous_distance = totaldistance(route,dis);

% This is a flag used to cool the current temperature after 100 iterations
temperature_iterations = 1;
% This is a flag used to plot the current route after 200 iterations
plot_iterations = 1;

% plot the current route
plotroute(city, route, previous_distance, temperature);

while 1.0 < temperature
    % generate randomly a neighbouring solution
    temp_route = perturb(route,'reverse');
    % compute total distance of the temp_route
    current_distance = totaldistance(temp_route, dis);
    % compute change of distance
    diff = current_distance - previous_distance;
    
    % Metropolis Algorithm
    if (diff < 0) || (rand < exp(-diff/(temperature)))
        route = temp_route;         %accept new route
        previous_distance = current_distance;
        
        % update iterations
        temperature_iterations = temperature_iterations + 1;
        plot_iterations = plot_iterations + 1;
        iterations = iterations + 1;
    end
    
    % reduce the temperature every 100 iterations
    if temperature_iterations >= 100
       temperature = cooling_rate*temperature;
       temperature_iterations = 0;
    end
    
    %  plot the current route every 200 iterations
    if plot_iterations >= 200
       plotroute(city, route, previous_distance,temperature);
       plot_iterations = 0;
    end
end

% plot the final solution
plotroute(city, route, previous_distance,temperature);

遗传算法

遗传算法的流程：
初始阶段：
随机生成一组可行解，也就是第一代染色体。
采用适应度函数计算每一条染色体的适应度，并根据适应度计算每一条染色体在下一代被选中的概率。
进化过程：
通过交叉生成N-M条染色体
对交叉后的染色体进行变异
最后选出这M条染色体，完成这次进化并进行下一次进化。

基因：染色体上的一个单元，解中的一个参数
染色体：由一组基因构成，问题可能的一个解
种群：由一系列染色体组成的一个集合
伪代码：

set initial generation k=0//初代为0
probability of mutation = alpha//变异概率为 $\alpha$
probability of performing crossover = beta//杂交概率
construct a population of n individuals Pk//设置一个含有n个个体的种群
while not termination do
  evalute:compute fitness(i) for each individuals in Pk//评价函数，得到每一个个体的适应度
  select:select m of Pk into Pk+1//选择函数，挑出m个到下一代当中  
  crossover:add alpha*m into Pk+1//杂交
  mutate:add beta*m into Pk+1//变异
  k=k+1
end while
return the fittest individual from P_last//返回一个最优秀的个体

如何对解进行编码

编码方法有二进制编码，格雷编码，实数编码，符号编码
编码的方法取决于变异杂交等等。

适应度函数

一般由目标函数直接或间接改造得到，且一般是非负的

选择

轮盘赌(比例选择)： $p_i=f_i/\sum_{j=1}^Nf_j$
两两竞争:
从父代中随机选择两个个体，比较适应值，保留优秀个体，淘汰较差个体。
排序选择：
根据个体适应度大小进行排序，然后根据序号进行选择。

杂交

有单点交叉和两点交叉，根据交叉点互换部分基因。

变异

分为单点变异和换位变异

TSP问题的遗传算法求解

对问题进行编码：
[1,...,i...j...,n,1]
构造适应度函数：
由于距离越短越好，所以构造适应度函数为目标函数的倒数。
$1/\sum_{i=1}^n dist(S_{(i)},S_{(i+1)})$
选择运算：
两两竞争、轮盘赌

2.元胞自动机

元胞分布于一维线性网格上，且元胞仅具有占和空两种状态，一个元胞的状态由两个邻居决定(某种规则，总共8种状态)。
对于二维的元胞，则其相当于生命游戏，元胞状态由周围八个邻居决定。
将元胞自动机抽象为数学：
A=(L,d,S,N,f)

L：元胞网格空间
d：元胞空间的维数
S：有限离散的状态集合
N：某邻域内所有元胞的集合
f：局部映射或局部规划

常用的二维网格：正方形网格，六边形网格，三角形网格。
对于边界条件：邻居周围是边界，一般有定值型、周期型、反射型、吸收型。
规则：实际上是一种状态转移函数，类型为总和型与合法型。

交通问题

总路程：L
车距：d
车长：l
车的数目：N
满足 $d=\frac{L-Nl}{N}=\frac{1}{\rho}-l$
其中 $\rho=\frac{N}{L}$
流量方程：单位时间内通过某路段的车辆数 $J=\rho v$
对于一段长x的路段，通过求导，有如下方程：
$J_x+\rho_t=0$ ，即流量对x的导数加上密度对t的导数和等于0，
又 $J=\rho v$ ，所以方程变为： $(\rho v)_x+\rho_t=0$
速度函数：
$v(\rho )=v_{max}$

3.预测与评价

评价问题

（1）加权平均

$p=\sum_{i=1}^n w_i p_i$

（2）层次分析

一般来讲分为三个层：目标层、准则层、备选层。
准则可以分为：C1颜值、C2身材....，然后将各个准则两两比较， $C_i$ 相对于 $C_j$ 的重要程度， $a_{ij}=\frac{C_i}{C_j}$ 。
将所有的因素两两比较，可以得到一个判断矩阵A，对角线上的元素都为1。若比较结果前后完全一致 $a_{ij}a_{jk}=a_{ik}$ 。
一致性检验： $CI=\frac{\lambda -n}{n-1}$ ， $CR=\frac{CI}{RI(n)}<0.1$ ，n为判断矩阵的维数， $\lambda$ 为最大特征值，RI(n)的取值可以查表获得。
一致性检验程序：

>> A=[1/1 2/1 5/1 3/1
1/2 1/1 3/1 1/2
1/5 1/3 1/1 1/4
1/3 2/1 4/1 1/1];
>> [V,D]=eig(A);%计算特征向量V与特征值D
>> [lamda,i]=max(diag(D))
>> CI=(lamda-4)/(4-1);
>> CR=CI/0.9

层次单排序，可以求出权重w：

>> W=V(:,i);
>> w=W/sum(W)

w =

    0.4816
    0.1864
    0.0713
    0.2608

对于层次的总排序：
得到的P即为最终的评分

function ahpactor

A = [1/1  2/1  5/1  3/1 
     1/2  1/1  3/1  1/2 
     1/5  1/3  1/1  1/4
     1/3  2/1  4/1  1/1];
[w, CR] = AHP(A);

% face
A1 = [1/1  1/2  3/1
      2/1  1/1  5/1
      1/3  1/5  1/1];
[w1, CR1] = AHP(A1);

% body
A2 = [1/1  1/3  2/1
      3/1  1/1  5/1
      1/2  1/5  1/1];
[w2, CR2] = AHP(A2);

% voice
A3 = [1/1  2/1  1/5
      1/2  1/1  1/7
      5/1  7/1  1/1];
[w3, CR3] = AHP(A3);

% acting
A4 = [1/1  2/1  1/3
      1/2  1/1  1/5
      3/1  5/1  1/1];
[w4, CR4] = AHP(A4);


CRs = [CR1 CR2 CR3 CR4]
P = [w1 w2 w3 w4] * w

 % ------------------------------------------------------------------------
 
function [w, CR] = AHP(A)
% n= [ 1    2    3    4    5    6    7    8    9
RI = [ 0.00 0.00 0.58 0.90 1.12 1.24 1.32 1.41 1.45];

n = size(A,1);
[V, D] = eig(A);

[lamda, i] = max(diag(D));
CI=(lamda-n)/(n-1);
CR = CI/RI(n);

W = V(:,i);
w = W/sum(W);

总结：
通过目标层与准则层两两比较，得到准则层的权重，再通过准则层与备选层两两比较，得到各项的加权权重。同样，如果有多层，也是这样层层分析。

模糊总和评价

模糊总和评价要素：因素集，评语集。
主要由权重以及投票(%)R来表示。投票是对每一个因素进行评语。
假设W是权重，R是投票情况，则模糊合成： $B=max(R.*W')$ ，算出来n个值，n也为评语数目，第几个最大，就对应第几个评语。

预测问题

（1）拟合

多次线性拟合：polyfit/fit用法
具体内容见机器学习

（2）时间序列

时间序列：预测对象按照时间顺序排列而成的序列。
时间预测：根据时序过去的变化规律，推测今后趋势。
时间序列变化形式：

长期变动趋势 $T_t$
季节变动 $S_t$
循环变动 $C_t$
不规则变动 $R_t$

模型：

加法模型
乘法模型
混合模型

移动平均法：
即t+1时刻的值为前N个时刻的值的算术平均
类似还有一次指数平滑法、一次差分指数平滑法等等。

（2）灰色预测

不使用原始数据，使用生成数据，且不需要很多数据，只适用于中短期的预测，只适合指数增长的预测。
最常用GM(1,1)预测模型，表示一阶微分方程，且只含一个变量。
对于原始序列 $X^{(0)}=x^{(0)}(1),...,x^{(0)}(n)$
可行性检验条件 $\lambda (k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)}$ 在 $(e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}})$ 范围内。如果不满足条件，可在原始序列上加上常数c组成新的原始序列。
一次累加生成序列，对于序列 $X^{(1)}$ ， $x^{(1)}(k)=\sum_{i=1}^kx^{(0)}(i)$
均值生成序列，在一次累加序列的基础上 $z^{(1)}(k)=x^{(1)}(k)/2+x^{(1)}(k-1)/2$ 得到。
灰微分方程， $x^{(0)}(k)+az^{(1)}(k)=b,k=2,3...n$ ，a与b是要求的两个参数。
白化微分方程: $\frac{dx^{(1)}}{dt}+ax^{(1)}(t)=b$ ，若求出a与b，则可以得到序列与时间的关系式。

a与b的估计值最小二乘估计自行上网查找。然后白化微分方程求解，然后将模型还原，还原为1阶。
例子：

t0 = [1999:2003]';
X0 = [89, 99, 109, 120, 135]';
n = length(X0);
lambda = X0(1:n-1)./X0(2:n);
%range = minmax(lambda')
%exp([-2/(n+1), 2/(n+2)])
X1 = cumsum(X0);
Z1 = (X1(1:n-1)+X1(2:n))/2
B = [-Z1, ones(n-1,1)];
Y = X0(2:n);
u = B\Y; a = u(1); b = u(2);
k = 0:n+4;
xhat1 = (X0(1) - b/a).*exp(-a*k) + b/a;%代入公式
xhat0 = [X0(1) diff(xhat1)]%还原
plot(t0,X0,'o',t0(1)+k, xhat0,'-+')

图论模型与算法

图的构成：顶点集V(G)，边集E(G)，关联函数 $\phi_G$ ，其为从边映射到节点的函数。环：端点重合为一点。连杆：端点不重合的边，重边：具有两个相同端点的边。
图分为有向图和无向图
对于有向图，我们把边改名为弧。
子图：原图去掉边或者节点得到的图。
完全图：图上任意两点都是有边相连的。
连通图：从图上的一点可以到任意一点。
图的表示方法：邻接矩阵、关联矩阵，如对于右向关联矩阵: $m_{va}=\{1,-1,0\}$ 表示边a与顶点v的关系，1表示头，-1表示尾，0表示无关系。邻接矩阵 $A=(a_{uv})$ 表示是否存在定点u到v的弧。
度：图G中与v关联的边数。对于有向图，度=入度+出度。
matlab的图论工具箱：

graphallshortestpaths 求图中所有顶点之间的最短距离
graphconnredcomp 找无(有)向图的(强/弱)连通分支
graphmaxflow 计算有向图中的最大流
graphshortestpaths 求指定两点间的最短路径
。。。

计算最短路的代码示例:

[a,b,c,d,e,f]=deal(1,2,3,4,5,6);
w=[0 2 3 0 0 0
2 0 6 5 3 0
3 6 0 0 1 0
0 5 0 0 1 2
0 3 1 1 0 4
0 0 0 2 4 0];
W=sparse(w);
[dist,path,pred]=graphshortestpath(W,a,f)

此外还有网络分析的工具箱，求一些关于度的算法。

排队论

排队论的基本构成：
输入过程：顾客总体(有限or无限)，到达的类型(单个or成批)，到达时间间隔。
排队规则：顾客按怎样的规定次序接受服务。有等待制、损失制、闭合制。
服务机构：服务台的数量、服务时间服从的分布。

排队论的数量指标：
队长：系统中的平均顾客数
等待对长：系统中处于等待的顾客数量
等待时间：系统中包括顾客的平均逗留时间。
忙期：连续保持服务的时常

关于排队论中的符号表示：
A/B/C/n
A输入过程 B服务时间 C服务台数 n系统容量

一个实例
$M/M/S/\infty$
输入过程服从Poisson流
服务时间服从负指数分布
系统有S个服务台平行服务
系统容量为无穷大的等待排队系统
当S=1，即单服务台
输入：
$P\{X(t)=k\}=\frac{(\lambda t)^ke^{-\lambda t}}{k!}$ ，[0,t]时间内平均到达顾客数为 $\lambda t$
服务时间：
每个顾客接受服务的时间1服从参数为 $\mu$ 的负指数分布。 $f(t)=\mu e^{-\mu t}(t>0)$ ，每个顾客的平均服务时间为 $\frac{1}{\mu}$
系统服务强度： $\rho=\frac{\lambda}{\mu}$
无顾客的概率：
$p_0=1-\rho=1-\frac{\lambda}{\mu}$
有n个顾客的概率：
$p_n=(1-\rho)\rho^n$
平均队长： $L_s=\frac{\lambda}{\mu-\lambda}$
平均逗留时间： $W_s=\frac{1}{\mu-\lambda}$
平均等待时间： $W_q=\frac{1}{\mu-\lambda}-\frac{1}{\mu}$
对于S>0的情况：
服务能力和强度分别为 $S\mu$ 与 $\rho =\frac{\lambda}{S\mu}$ ，其他的公式直接查找参考文献使用即可。
对于单服务台：
服务时刻(i)=max{到达时刻(i)，离开时刻(i-1)}
离开时刻(i)=服务时刻(i)+服务时长(i)
等待时长(i)=离开时刻(i)-到达时刻(i)
对于多个服务台也有近似情况。