基于markov鏈的長(zhǎng)江水質(zhì)預(yù)測(cè)模型

基于markov鏈的長(zhǎng)江水質(zhì)預(yù)測(cè)模型

馬爾可夫預(yù)測(cè)方法是馬爾可夫鏈在預(yù)測(cè)領(lǐng)域的應(yīng)用方法。起初，這種方法在水文、天然氣、地震等方面得到了廣泛應(yīng)用，此后，經(jīng)濟(jì)學(xué)家在研究生產(chǎn)力和預(yù)測(cè)業(yè)務(wù)利潤(rùn)方面發(fā)揮了作用。在馬爾可夫預(yù)測(cè)第167個(gè)非常重要的問(wèn)題之一是一對(duì)狀態(tài)轉(zhuǎn)換的概率分布。傳統(tǒng)的估計(jì)方法是將n狀態(tài)的觀察時(shí)間和系統(tǒng)從當(dāng)前狀態(tài)移到另一個(gè)狀態(tài)的概率中進(jìn)行估計(jì)。當(dāng)系統(tǒng)從當(dāng)前狀態(tài)轉(zhuǎn)移到另一個(gè)狀態(tài)時(shí)，需要通過(guò)概率分布概率，然后計(jì)算一階躍概率矩陣。但事實(shí)上，我們不知道系統(tǒng)在n狀態(tài)下觀察到的總概率，而系統(tǒng)在不同時(shí)間內(nèi)在n狀態(tài)下的概率。1轉(zhuǎn)移m步后狀態(tài)概率向量預(yù)測(cè)模型馬爾可夫鏈?zhǔn)侵赶到y(tǒng)的未來(lái)狀態(tài),僅與現(xiàn)在的狀態(tài)有關(guān),而與以前狀態(tài)無(wú)關(guān)的隨機(jī)過(guò)程,它具有無(wú)后效性特點(diǎn).馬爾可夫鏈的預(yù)測(cè)就是根據(jù)某些變量的現(xiàn)在狀態(tài)及其變化趨向,預(yù)測(cè)其在未來(lái)某一特定期間內(nèi)可能出現(xiàn)的狀態(tài),從而為決策提供依據(jù).設(shè)一個(gè)系統(tǒng)有n個(gè)狀態(tài),在時(shí)刻t的n個(gè)狀態(tài)設(shè)為s1,s2,…,sn,在時(shí)刻t位于狀態(tài)si,而在時(shí)刻t+1轉(zhuǎn)移到狀態(tài)sj的概率為pij,i,j=1,2,…,n,稱矩陣P=(pij)n×n為一步轉(zhuǎn)移概率矩陣.根據(jù)概率轉(zhuǎn)移矩陣的定義,其中的元素滿足:pij≥0,n∑j=1pij=1,i=1,2,?,npij≥0,∑j=1npij=1,i=1,2,?,n,即矩陣P(pij)n×n中每一行的元素之和等于1.設(shè)在初始時(shí)刻t=0時(shí),系統(tǒng)的n種狀態(tài)s1,s2,…,sn的概率向量為α(0)=(p0(1),p0(2),…,p0(n)),則系統(tǒng)在時(shí)刻t=1的概率向量α(1)=(p1(1),p1(2),…,p1(n))為(p1(1),p1(2),?,p1(n))=(p0(1),p0(2),?,p0(n))[p11?p1n??pn1?pn](p1(1),p1(2),?,p1(n))=(p0(1),p0(2),?,p0(n))????p11?pn1??p1n?pn????即α(1)=α(0)P,按上式遞推,可以得到系統(tǒng)在時(shí)刻t=m的n種狀態(tài)的概率向量α(m)=(pm(1),pm(2),…,pm(n))為α(m)=α(0)Pm,此式即為轉(zhuǎn)移m步后狀態(tài)概率向量預(yù)測(cè)模型.2步狀態(tài)轉(zhuǎn)移概率矩陣計(jì)算傳統(tǒng)的估算方法178是已知系統(tǒng)存在n種狀態(tài)S={s1,s2,…,sn},假設(shè)在N次觀察中,系統(tǒng)處于第i狀態(tài)共有ni次,顯然Ν=n∑i=1niN=∑i=1nni,設(shè)系統(tǒng)當(dāng)前時(shí)刻處于第i狀態(tài)si在下一個(gè)時(shí)刻轉(zhuǎn)移到第j狀態(tài)狀態(tài)sj的次數(shù)為nij次,顯然ni=n∑j=1nijni=∑j=1nnij,用頻率估計(jì)概率的方法可以估算一步狀態(tài)轉(zhuǎn)移概率矩陣P=(pij)n×n,其中pij=Ρ{X=sj|X=si}=nijni,i,j=1,2,?,npij=P{X=sj|X=si}=nijni,i,j=1,2,?,n顯然pij∈?n∑j=1pij=1,i=1,2,?,n.pij∈?∑j=1npij=1,i=1,2,?,n.然而歷史統(tǒng)計(jì)資料有時(shí)并沒(méi)有給出系統(tǒng)處于n種狀態(tài)的總的觀察次數(shù)N,只給出系統(tǒng)處于n種狀態(tài)下的概率,此時(shí)我們無(wú)法按著上述方法來(lái)估計(jì)一步狀態(tài)轉(zhuǎn)移概率矩陣.因此有必要從另外的途徑獲取一步狀態(tài)轉(zhuǎn)移概率矩陣的估算方法.為了獲得較為精確的一步轉(zhuǎn)移概率矩陣,利用最優(yōu)化的思想,即在m個(gè)時(shí)刻中要使實(shí)際狀態(tài)的概率向量與理論計(jì)算的狀態(tài)的概率向量的誤差平方和達(dá)到最小為準(zhǔn)則,為此可建立最優(yōu)化模型.設(shè)α(t)=(pt(1),pt(2),…,pt(n))是時(shí)刻t系統(tǒng)在n個(gè)狀態(tài)下的概率向量,t=0,1,2,…,m,設(shè)一步狀態(tài)轉(zhuǎn)移概率矩陣為P=(pij)n×n.實(shí)際上由于客觀環(huán)境的變化,相鄰時(shí)刻的一步轉(zhuǎn)移概率矩陣并不完全相同,因此α(t+1)與α(t)P之間總存在誤差,t=0,1,2,…,m-1,由誤差平方和達(dá)到最小的準(zhǔn)則,建立如下最優(yōu)化模型,記為模型(1).minf(Ρ)=m-1∑t=0∥α(t+1)-α(t)Ρ∥2=m-1∑t=0(α(t+1)-α(t)Ρ)(α(t+1)-α(t)Ρ)Τs.t{n∑j=1pij=1,i=1,2,?,npij≥0,i,j=1,2,?,n其中約束條件是由一步轉(zhuǎn)移概率矩陣的定義得來(lái)的.模型(1)是一個(gè)二次規(guī)劃問(wèn)題,可以有現(xiàn)成的算法求解,也可用Matlab軟件求解.3t下面以2005年全國(guó)數(shù)學(xué)模型競(jìng)賽的A題為例,說(shuō)明本模型的應(yīng)用.已知1995—2004年枯水期長(zhǎng)江全流域Ⅰ類、Ⅱ類、Ⅲ類、Ⅳ類、Ⅴ類、劣Ⅴ類這6類水質(zhì)所占的比例,如表1所示.現(xiàn)要對(duì)長(zhǎng)江未來(lái)10年的水質(zhì)污染的發(fā)展趨勢(shì)作出預(yù)測(cè).為此可建立長(zhǎng)江水質(zhì)污染的馬爾可夫鏈趨勢(shì)預(yù)測(cè)的一步轉(zhuǎn)移概率矩陣估計(jì)的最優(yōu)化模型.設(shè)枯水期長(zhǎng)江全流域水質(zhì)在第t年屬于Ⅰ類、Ⅱ類、Ⅲ類、Ⅳ類、Ⅴ類、劣Ⅴ類這6類狀態(tài)的比例向量分別為α(t)=(pt(1),pt(2),…,pt(6)),t=0,1,2,…,9.設(shè)P=(pij)6×6為6類狀態(tài)的一步轉(zhuǎn)移概率矩陣,根據(jù)誤差平方和達(dá)到最小的準(zhǔn)則,把表一中數(shù)據(jù)代入模型(1),建立如下最優(yōu)化模型:minf(Ρ)=9∑t=0(α(t+1)-α(t)Ρ)(α(t+1)-α(t)Ρ)Τs.t{6∑j=1pij=1,i=1,2,?,6pij≥0,i,j=1,2,?,6用Matlab軟件求解,得Ρ=[0.46070.00000.53930.00000.00000.00000.03450.27780.50490.06120.06030.06140.00000.27050.33080.32300.00000.01420.00001.00000.00000.00000.00000.00000.00000.11070.17470.00000.00000.1746]由下式α(9+k)=α(9)Ρk,k=1,2,?,10可以對(duì)長(zhǎng)江未來(lái)10年的水質(zhì)污染屬于Ⅰ類、Ⅱ類、Ⅲ類、Ⅳ類、Ⅴ類、劣Ⅴ類這6類狀態(tài)的比例向量作出預(yù)測(cè),預(yù)測(cè)結(jié)果見(jiàn)表2:從預(yù)測(cè)計(jì)算結(jié)果可以看出,枯水期長(zhǎng)江全流域水質(zhì)屬于Ⅳ類、Ⅴ類、劣Ⅴ類這6類狀態(tài)的比例并沒(méi)有發(fā)生根本性的減少,