數學建模

1.啟發(fā)式算法

它主要包括禁忌搜索,模擬退火,遺傳算法,神經網絡,蟻群算法

模擬退火算法

流程圖

Metropolis準則:對于目標函數E,若E_{i+1}<E_i,則接收E_{i+1},否則依e^{-(E_{i+1}-E_i)/KT}概率接受。
模擬退火與物理退火的對應關系:

解--粒子狀態(tài)
目標函數--能量
最優(yōu)解--能量最低態(tài)
設定初溫--加熱過程
擾動--熱漲落
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//返回一個最優(yōu)秀的個體
如何對解進行編碼

編碼方法有二進制編碼,格雷編碼,實數編碼,符號編碼
編碼的方法取決于變異雜交等等。

適應度函數

一般由目標函數直接或間接改造得到,且一般是非負的

選擇

輪盤賭(比例選擇):p_i=f_i/\sum_{j=1}^Nf_j
兩兩競爭:
從父代中隨機選擇兩個個體,比較適應值,保留優(yōu)秀個體,淘汰較差個體。
排序選擇:
根據個體適應度大小進行排序,然后根據序號進行選擇。

雜交

有單點交叉和兩點交叉,根據交叉點互換部分基因。

變異

分為單點變異和換位變異

TSP問題的遺傳算法求解

對問題進行編碼:
[1,...,i...j...,n,1]
構造適應度函數:
由于距離越短越好,所以構造適應度函數為目標函數的倒數。
1/\sum_{i=1}^n dist(S_{(i)},S_{(i+1)})
選擇運算:
兩兩競爭、輪盤賭


2.元胞自動機

元胞分布于一維線性網格上,且元胞僅具有占和空兩種狀態(tài),一個元胞的狀態(tài)由兩個鄰居決定(某種規(guī)則,總共8種狀態(tài))。
對于二維的元胞,則其相當于生命游戲,元胞狀態(tài)由周圍八個鄰居決定。
將元胞自動機抽象為數學:
A=(L,d,S,N,f)

L:元胞網格空間
d:元胞空間的維數
S:有限離散的狀態(tài)集合
N:某鄰域內所有元胞的集合
f:局部映射或局部規(guī)劃

常用的二維網格:正方形網格,六邊形網格,三角形網格。
對于邊界條件:鄰居周圍是邊界,一般有定值型、周期型、反射型、吸收型。
規(guī)則:實際上是一種狀態(tài)轉移函數,類型為總和型與合法型。

交通問題

總路程: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,所以方程變?yōu)椋?img class="math-inline" src="https://math.jianshu.com/math?formula=(%5Crho%20v)_x%2B%5Crho_t%3D0" alt="(\rho v)_x+\rho_t=0" mathimg="1">
速度函數:
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)時間序列

時間序列:預測對象按照時間順序排列而成的序列。
時間預測:根據時序過去的變化規(guī)律,推測今后趨勢。
時間序列變化形式:

長期變動趨勢 T_t
季節(jié)變動 S_t
循環(huán)變動 C_t
不規(guī)則變動 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),關聯(lián)函數\phi_G,其為從邊映射到節(jié)點的函數。環(huán):端點重合為一點。連桿:端點不重合的邊,重邊:具有兩個相同端點的邊。
圖分為有向圖和無向圖
對于有向圖,我們把邊改名為弧。
子圖:原圖去掉邊或者節(jié)點得到的圖。
完全圖:圖上任意兩點都是有邊相連的。
連通圖:從圖上的一點可以到任意一點。
圖的表示方法:鄰接矩陣、關聯(lián)矩陣,如對于右向關聯(lián)矩陣:m_{va}=\{1,-1,0\}表示邊a與頂點v的關系,1表示頭,-1表示尾,0表示無關系。鄰接矩陣A=(a_{uv})表示是否存在定點u到v的弧。
度:圖G中與v關聯(lián)的邊數。對于有向圖,度=入度+出度。
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成批),到達時間間隔。
排隊規(guī)則:顧客按怎樣的規(guī)定次序接受服務。有等待制、損失制、閉合制。
服務機構:服務臺的數量、服務時間服從的分布。

排隊論的數量指標:
隊長:系統(tǒng)中的平均顧客數
等待對長:系統(tǒng)中處于等待的顧客數量
等待時間:系統(tǒng)中包括顧客的平均逗留時間。
忙期:連續(xù)保持服務的時常

關于排隊論中的符號表示:
A/B/C/n
A輸入過程 B服務時間 C服務臺數 n系統(tǒng)容量

一個實例
M/M/S/\infty
輸入過程服從Poisson流
服務時間服從負指數分布
系統(tǒng)有S個服務臺平行服務
系統(tǒng)容量為無窮大的等待排隊系統(tǒng)
當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}
系統(tǒng)服務強度:\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)
對于多個服務臺也有近似情況。

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容