數(shù)學(xué)建模穩(wěn)定狀態(tài)模型及數(shù)學(xué)建模算法大全現(xiàn)代優(yōu)化算法簡介

第十四章穩(wěn)定狀態(tài)模型雖然動(dòng)態(tài)過程的變化規(guī)律一般要用微分方程建立的動(dòng)態(tài)模型來描述，但是對于某些實(shí)際問題，建模的主要目的并不是要尋求動(dòng)態(tài)過程每個(gè)瞬時(shí)的性態(tài)，而是研究某種意義下穩(wěn)定狀態(tài)的特征，特別是當(dāng)時(shí)間充分長以后動(dòng)態(tài)過程的變化趨勢。譬如在什么情況下描述過程的變量會(huì)越來越接近某些確定的數(shù)值，在什么情況下又會(huì)越來越遠(yuǎn)離這些數(shù)值而導(dǎo)致過程不穩(wěn)定。為了分析這種穩(wěn)定與不穩(wěn)定的規(guī)律常常不需要求解微分方程，而可以利用微分方程穩(wěn)定性理論，直接研究平衡狀態(tài)的穩(wěn)定性就行了。本章先介紹平衡狀態(tài)與穩(wěn)定性的概念，然后列舉幾個(gè)這方面的建模例子?！?微分方程穩(wěn)定性理論簡介定義1稱一個(gè)常微分方程（組）是自治的，如果方程（組）（1）中的，即在中不含時(shí)間變量。事實(shí)上，如果增補(bǔ)一個(gè)方程，一個(gè)非自治系統(tǒng)可以轉(zhuǎn)化自治系統(tǒng)，就是說，如果定義，且引入另一個(gè)變量，則方程（1）與下述方程是等價(jià)的。這就是說自治系統(tǒng)的概念是相對的。下面僅考慮自治系統(tǒng)，這樣的系統(tǒng)也稱為動(dòng)力系統(tǒng)。定義2系統(tǒng)（2）的相空間是以為坐標(biāo)的空間，特別，當(dāng)時(shí)，稱相空間為相平面?？臻g中的點(diǎn)集稱為系統(tǒng)（2）的軌線，所有軌線在相空間中的分布圖稱為相圖。定義3相空間中滿足的點(diǎn)稱為系統(tǒng)（2）的奇點(diǎn)（或平衡點(diǎn)）。奇點(diǎn)可以是孤立的，也可以是連續(xù)的點(diǎn)集。例如，系統(tǒng)（3）當(dāng)時(shí)，有一個(gè)連續(xù)的奇點(diǎn)的集合。當(dāng)時(shí)，是這個(gè)系統(tǒng)的唯一的奇點(diǎn)。下面僅考慮孤立奇點(diǎn)。為了知道何時(shí)有孤立奇點(diǎn)，給出下述定理：定理1設(shè)是實(shí)解析函數(shù)，且系統(tǒng)（2）的奇點(diǎn)。若在點(diǎn)處的Jacobian矩陣是非奇異的，則是該系統(tǒng)的孤立奇點(diǎn)。定義4設(shè)是（2）的奇點(diǎn)，稱（=1\*romani）是穩(wěn)定的，如果對于任意給定的，存在一個(gè)，使得如果，則對所有的都成立。（=2\*romanii）是漸近穩(wěn)定的，如果它是穩(wěn)定的，且。這樣，如果當(dāng)系統(tǒng)的初始狀態(tài)靠近于奇點(diǎn)，其軌線對所有的時(shí)間仍然接近它，于是說是穩(wěn)定的。另一方面，如果當(dāng)時(shí)這些軌線趨于，則是漸近穩(wěn)定的。定義5一個(gè)奇點(diǎn)不是穩(wěn)定的，則稱這個(gè)奇點(diǎn)是不穩(wěn)定的。對于常系數(shù)齊次線性系統(tǒng)（3）有下述定理。定理2設(shè)是系統(tǒng)（3）的通解。則（=1\*romani）如果系統(tǒng)（3）的系數(shù)矩陣的一切特征根的實(shí)部都是負(fù)的，則系統(tǒng)（3）的零解是漸近穩(wěn)定的。（=2\*romanii）如果的特征根中至少有一個(gè)根的實(shí)部是正的，則系統(tǒng)（3）的零解是不穩(wěn)定的。（=3\*romaniii）如果的一切特征根的實(shí)部都不是正的，但有零實(shí)部，則系統(tǒng)（3）的零解可能是穩(wěn)定的，也可能是不穩(wěn)定的，但總不會(huì)是漸近穩(wěn)定的。定理2告訴我們：系統(tǒng)（3）的零解漸近穩(wěn)定的充分必要條件是的一切特征根的實(shí)部都是負(fù)的。對于非線性系統(tǒng)，一般不可能找出其積分曲線或軌跡，也就不可能直接導(dǎo)出奇點(diǎn)的穩(wěn)定性。為克服這一困難，在奇點(diǎn)附近用一個(gè)線性系統(tǒng)來近似這個(gè)非線性系統(tǒng)，用這個(gè)近似系統(tǒng)的解來給出這個(gè)奇點(diǎn)的穩(wěn)定解。定義6設(shè)是系統(tǒng)（2）的一個(gè)孤立奇點(diǎn)。稱系統(tǒng)在點(diǎn)幾乎是線性的，如果在的Jacobian矩陣是非奇異的，即。設(shè)在的某鄰域內(nèi)連續(xù)，并有直到二階連續(xù)偏導(dǎo)數(shù)，則由多元函數(shù)的Taylor公式，可將展開成，其中是一個(gè)常數(shù)矩陣，這樣得到的線性系統(tǒng)（4）稱為系統(tǒng)（2）的線性近似。一開始，人們以為總可以用線性近似系統(tǒng)來代替所研究的原系統(tǒng)。但后來人們發(fā)現(xiàn)，這種看法是不對的，或至少說是不全面的，非線性系統(tǒng)中的許多性質(zhì)，在它的線性近似中不再保留。即使象零解穩(wěn)定性這樣一個(gè)問題，也要在一定條件下，才可用它的線性近似系統(tǒng)代替原系統(tǒng)來研究。關(guān)于這個(gè)問題，我們有下述定理：定理3如果系統(tǒng)（4）的零解是漸近穩(wěn)定的，或不穩(wěn)定的，則原系統(tǒng)的零解也是漸近穩(wěn)定的或不穩(wěn)定的。然而，如果系統(tǒng)（4）的零解是穩(wěn)定的，則原系統(tǒng)的零解是不定的，即此時(shí)不能從線性化的系統(tǒng)來導(dǎo)出原系統(tǒng)的穩(wěn)定性。系統(tǒng)（3）在其系數(shù)矩陣的行列式的條件下，可知是系統(tǒng)（3）的唯一的平衡點(diǎn)，它的穩(wěn)定性由特征方程：的根（特征根）決定。定理4設(shè)線性系統(tǒng)（3）所對應(yīng)的特征方程是其中，。設(shè)和是它的根，則當(dāng)時(shí)關(guān)于奇點(diǎn)有下述結(jié)論：（=1\*romani），是穩(wěn)定結(jié)點(diǎn)；（=2\*romanii），是穩(wěn)定退化結(jié)點(diǎn)；（=3\*romaniii），是不穩(wěn)定結(jié)點(diǎn)；（=4\*romaniv），是不穩(wěn)定退化結(jié)點(diǎn)；（=5\*romanv），是不穩(wěn)定鞍點(diǎn)；（=6\*romanvi），是穩(wěn)定焦點(diǎn)；（=7\*romanvii），是不穩(wěn)定焦點(diǎn)；（=8\*romanviii），是不穩(wěn)定中心。定理5設(shè)非線性系統(tǒng)（5）中的和滿足條件：（=1\*romani）在點(diǎn)的某鄰域內(nèi)存在連續(xù)的一階偏導(dǎo)數(shù)。（=2\*romanii）存在常數(shù)，使得，（）又設(shè)系統(tǒng)（5）的一次近似系統(tǒng)（3）的特征方程的根沒有零實(shí)部，則（5）式與（3）式的奇點(diǎn)的類型相同，并有相同的穩(wěn)定性或不穩(wěn)定性?！?再生資源的管理和開發(fā)漁業(yè)資源是一種再生資源，再生資源要注意適度開發(fā)，不能為了一時(shí)的高產(chǎn)“竭澤而漁”，應(yīng)該在持續(xù)穩(wěn)產(chǎn)的前提下追求最高產(chǎn)量或最優(yōu)的經(jīng)濟(jì)效益。這是一類可再生資源管理與開發(fā)的模型，這類模型的建立一般先考慮在沒有收獲的情況下資源自然增長模型，然后再考慮收獲策略對資源增長情況的影響。2.1資源增長模型考慮某種魚的種群的動(dòng)態(tài)。在建立模型之前，做如下的基本假設(shè)：（=1\*romani）魚群的數(shù)量本身是離散變量，談不到可微性。但是，由于突然增加或減少的只是單一個(gè)體或少數(shù)幾個(gè)個(gè)體，與全體數(shù)量相比，這種增長率是微小的。所以，可以近似地假設(shè)魚群的數(shù)量隨時(shí)間連續(xù)地，甚至是可微地變化。（=2\*romanii）假設(shè)魚群生活在一個(gè)穩(wěn)定的環(huán)境中，即其增長率與時(shí)間無關(guān)。（=3\*romaniii）種群的增長是種群個(gè)體死亡與繁殖共同作用的結(jié)果。（=4\*romaniv）資源有限的生存環(huán)境對種群的繁衍，生長有抑制作用，而且這一作用與魚群的數(shù)量是成正比的。記時(shí)刻漁場中魚量為，我們可以得到所滿足的Logistic模型：（6）其中是固有增長率，是環(huán)境容許的最大魚量。由分離變量法求得方程（6）解為（6）式有兩個(gè)平衡點(diǎn)，即，，其中是不穩(wěn)定的，在正半軸內(nèi)全局穩(wěn)定。2.2資源開發(fā)模型建立一個(gè)在捕撈情況下漁場魚量遵從的方程，分析魚量穩(wěn)定的條件，并且在穩(wěn)定的前提下，討論如何控制捕撈使持續(xù)產(chǎn)量或經(jīng)濟(jì)效益達(dá)到最大。設(shè)單位時(shí)間的捕撈量與漁場魚量成正比，比例系數(shù)表示單位時(shí)間捕撈率，可以進(jìn)一步分解分解為，稱為捕撈強(qiáng)度，用可以控制的參數(shù)如出海漁船數(shù)來度量；稱為捕撈系數(shù)，表示單位強(qiáng)度下的捕撈率。為方便取，于是單位時(shí)間的捕撈量為。常數(shù)，表示一個(gè)特定的捕撈策略，即要求捕魚者每天只能捕撈一定的數(shù)量。這樣，捕撈情況下漁場魚量滿足方程（7）這是一個(gè)一階非線性方程，且是黎卡提型的。也稱為Scheafer模型。希望知道漁場的穩(wěn)定魚量和保持穩(wěn)定的條件，即時(shí)間足夠長以后漁場魚量的趨向，并且由此確定最大持續(xù)產(chǎn)量。在平衡點(diǎn)處有，方程（7）有兩個(gè)平衡點(diǎn)，。顯然，它們均是方程的解。在的情況下，是一正平衡點(diǎn)。（7）式可改寫為（8）易知，當(dāng)時(shí)，；時(shí)，，即平衡解是不穩(wěn)定的，而是穩(wěn)定平衡解。即在捕撈強(qiáng)度的情況下，漁場魚量將穩(wěn)定在的水平，因此產(chǎn)量（單位時(shí)間的捕撈量）也將穩(wěn)定在的水平，即此時(shí)可獲得持續(xù)收獲量。當(dāng)然，當(dāng)時(shí)，，漁場魚量將逐漸減少至，這時(shí)的捕撈其實(shí)是“竭澤而漁”，當(dāng)然談不上獲得持續(xù)產(chǎn)量了。如何才能做到漁資源在持續(xù)捕撈的條件下為我們提供最大的收益？從數(shù)學(xué)上說，就是在或的條件下極大化所期望的“收益”，這里的“收益”可理解為產(chǎn)量，則問題就可以數(shù)學(xué)地?cái)⑹鰹橄率鰞?yōu)化問題：約束條件為。這里它可以歸結(jié)為的二次函數(shù)的最大值問題。簡單的推導(dǎo)不難得到最大持續(xù)捕撈強(qiáng)度為，最大持續(xù)產(chǎn)量為。捕撈強(qiáng)度是得到最大持續(xù)捕魚量的策略。2.3經(jīng)濟(jì)效益模型當(dāng)今，對魚類資源的開發(fā)和利用已經(jīng)成為人類經(jīng)濟(jì)活動(dòng)的一部分。其目的不是追求最大的漁產(chǎn)量而是最大的經(jīng)濟(jì)收益。因而一個(gè)自然的想法就是進(jìn)一步分析經(jīng)濟(jì)學(xué)行為對魚類資源開發(fā)利用的影響。如果經(jīng)濟(jì)效益用從捕撈所得的收入中扣除開支后的利潤來衡量，并且簡單地設(shè)魚的銷售單價(jià)為常數(shù)，單位捕撈強(qiáng)度（如每條出海漁船）的費(fèi)用為常數(shù)，那么單位時(shí)間的收入和支出分別為，單位時(shí)間的利潤為利潤是漁民所關(guān)注的焦點(diǎn)。因此在制定管理策略時(shí)所期望極大化的“收益”，這時(shí)就應(yīng)理解為經(jīng)濟(jì)利潤或凈收入而不是魚的產(chǎn)量。因而所討論的問題就變成了在使魚量穩(wěn)定在的約束條件下的。即求的最大值。容易求出使達(dá)到最大的捕撈強(qiáng)度為最大利潤下的漁場穩(wěn)定魚量最大利潤下漁場單位時(shí)間的持續(xù)產(chǎn)量為最大可持續(xù)凈收益與前一模型相比較可以看出，在最大效益原則下捕撈強(qiáng)度和持續(xù)產(chǎn)量均有減少，而漁場的魚量有所增加。并且，減少或增加的比例隨著捕撈成本的增長而變大，隨著銷售價(jià)格的增長而變小，這顯然是符合實(shí)際情況的。2.4種群的相互競爭模型有甲乙兩個(gè)種群，當(dāng)它們獨(dú)自在一個(gè)自然環(huán)境中生存時(shí)，數(shù)量的演變均遵從Logistic規(guī)律。記是兩個(gè)種群的數(shù)量，是它們的固有增長率，是它們的最大容量。于是，對于種群甲有其中，因子反映由于甲對有限資源的消耗導(dǎo)致的對它本身增長的阻滯作用，可解釋為相對而言單位數(shù)量的甲消耗的食物量（設(shè)食物總量為1）。當(dāng)兩個(gè)種群在同一自然環(huán)境中生存時(shí)，考察由于乙消耗同一種有限資源對甲的增長產(chǎn)生的影響，可以合理地在因子中再減去一項(xiàng)，該項(xiàng)與種群乙的數(shù)量（相對于而言）成正比，于是，種群甲增長的方程為（9）這里的意義是，單位數(shù)量乙（相對而言）消耗的供養(yǎng)甲的食物量為單位數(shù)量甲（相對）消耗的供養(yǎng)乙的食物量的倍，類似地，甲的存在也影響了乙的增長，種群乙的方程應(yīng)該是（10）對可作相應(yīng)的解釋。在兩個(gè)種群的相互競爭中，是兩個(gè)關(guān)鍵的指標(biāo)。從上面對它們的解釋可知，表示在消耗供養(yǎng)甲的資源中，乙的消耗多于甲，對可作相應(yīng)的理解。一般來說，之間沒有確定的關(guān)系，在此我們僅討論相互獨(dú)立的情形。目的是研究兩個(gè)種群相互競爭的結(jié)局，即時(shí)，的趨向，不必要解方程組（9）和（10），只需對它的平衡點(diǎn)進(jìn)行穩(wěn)定性分析。為此我們解代數(shù)方程(11)得到四個(gè)平衡點(diǎn)分別為，，，。為分析這些點(diǎn)的穩(wěn)定性，需使用相空間的技巧。首先找出在平面上使或的區(qū)域。注意到，當(dāng)時(shí)，但要使和，當(dāng)且僅當(dāng)類似地，當(dāng)且僅當(dāng)這樣我們得到在平面中,直線(12)把平面劃分為和兩個(gè)區(qū)域。類似地，對進(jìn)行分析得到（=1\*romani），當(dāng)且僅當(dāng)（=2\*romanii），當(dāng)且僅當(dāng)（=3\*romaniii）直線（13）將平面劃分為和兩個(gè)區(qū)域。兩直線（12）和（13）之間的位置關(guān)系可以由下圖的四種情況來說明。每種可能性對應(yīng)于平衡點(diǎn)的穩(wěn)定性說明如下：（=1\*romani），由圖（b）知，兩直線將平面劃分為三個(gè)區(qū)域：（14）（15）（16）可以說明不論軌線從哪個(gè)區(qū)域出發(fā)，時(shí)都將趨向。若軌線從出發(fā)，由（14）式可知隨著的增加軌線向右上方運(yùn)動(dòng)，必然進(jìn)入；若軌線從出發(fā)，由（15）式可知軌線向下方運(yùn)動(dòng)，那么它或者趨向點(diǎn)，或者進(jìn)入，但進(jìn)入是不可能的。因?yàn)?，如果設(shè)軌線在某時(shí)刻經(jīng)直線（12）式進(jìn)入，則，由式（9）、（10）不難看出。由式（15）、（16）知，故，表明在達(dá)到極小值，而這是不可能的，因?yàn)樵谥?，即一直是增加的。若軌線從出發(fā)，由（16）式可知軌線向左下方運(yùn)動(dòng)，那么它或者趨向點(diǎn)，或者進(jìn)入，而進(jìn)入后，根據(jù)上面的分析最終也將趨向。綜合上述分析可以畫出軌線示意圖。因?yàn)橹本€（12）式上，所以在（12）式上軌線方向垂直于軸；在（13）式上，軌線方向平行于軸。（=2\*romanii），類似的分析可知穩(wěn)定。（=3\*romaniii），穩(wěn)定。（=4\*romaniv），不穩(wěn)定（鞍點(diǎn)）。因?yàn)檐壘€的初始位置不同，其走向也不同或趨于或趨于。根據(jù)建模過程中的含義，可以說明點(diǎn)穩(wěn)定在生態(tài)學(xué)上的意義：①：意味著在對供養(yǎng)甲的資源的競爭中乙弱于甲，意味著在對供養(yǎng)乙的資源的競爭中甲強(qiáng)于乙，于是種群乙終將滅絕，種群甲趨向最大容量，即趨向平衡點(diǎn)。②：情況與=1\*GB3①正好相反。=3\*GB3③：因?yàn)樵诟偁幖椎馁Y源中乙較弱，而在競爭乙的資源中甲較弱，于是可以達(dá)到一個(gè)雙方共存的穩(wěn)定的平衡狀態(tài)，這是種群競爭中很少出現(xiàn)的情況。=4\*GB3④：留作習(xí)題。§3Volterra模型意大利生物學(xué)家D'Ancona曾致力于魚類種群相互制約關(guān)系的研究，在研究過程中他無意中發(fā)現(xiàn)了第一次世界大戰(zhàn)期間地中海各港口捕獲的幾種魚類占捕獲總量百分比的資料，從這些資料中他發(fā)現(xiàn)各種軟骨掠肉魚，如鯊魚、鰱魚等我們稱之為捕食者的一些不是很理想的魚占總漁獲量的百分比，在1914～1923年期間，意大利阜姆港收購的捕食者所占的比例有明顯的增加：年代1914 1915191619171918百分比11.921.422.121.236.4年代1919192019211922 1923百分比27.316.015.914.810.7他知道，捕獲的各種魚的比例基本上代表了地中海漁場中各種魚類的比例。戰(zhàn)爭中捕獲量大幅度下降，當(dāng)然使?jié)O場中食用魚（食餌）增加，以此為生的鯊魚也隨之增加。但是捕獲量的下降為什么會(huì)使鯊魚的比例增加，即對捕食者而不是對食餌有利呢？他無法解釋這個(gè)現(xiàn)象，于是求助于著名的意大利數(shù)學(xué)家V.Volterra，希望建立一個(gè)食餌—捕食者系統(tǒng)的數(shù)學(xué)模型，定量地回答這個(gè)問題。3.l形成模型為建立這樣的模型，我們分別用和記食餌和捕食者在時(shí)刻的數(shù)量。因?yàn)榇蠛Ｖ恤~類的資源豐富，可以假設(shè)如果食餌獨(dú)立生存則食餌將以增長率按指數(shù)規(guī)律增長，即有。捕食者的存在使食餌的增長率降低，設(shè)降低的程度與捕食者數(shù)量成正比，于是滿足方程（17）比例系數(shù)反映捕食者掠取食餌的能力。捕食者離開食餌無法生存，若設(shè)它獨(dú)自存在時(shí)死亡率為，即，而食餌為它提供食物的作用相當(dāng)于使死亡率降低，或使之增長。設(shè)這個(gè)作用與食餌數(shù)量成正比，于是滿足（18）比例系數(shù)反映食餌對捕食者的供養(yǎng)能力。方程（17）和（18）是在沒有人工捕獲情況下自然環(huán)境中食餌與捕食者之間的制約關(guān)系，是Volterra提出的最簡單的模型。這個(gè)模型沒有引入競爭項(xiàng)。3.2模型分析這是一個(gè)非線性模型，不能求出其解析解，所以我們還是通過平衡點(diǎn)的穩(wěn)定性分析，研究的變化規(guī)律。容易得到方程（17）和（18）的平衡點(diǎn)為,(19)當(dāng)然，平衡解對我們來說是沒有意義的。這個(gè)方程組還有一族解，和，。因此，軸和軸都是方程組（17），（18）的軌線。這意味著：方程（17）、（18）在由第一象限出發(fā)的每一個(gè)解在以后一切時(shí)間都保持在第一象限內(nèi)。當(dāng)時(shí)，方程（17）、（18）的軌線是一階方程的解曲線。用分離變量方法解得(20)是任意常數(shù)。因此，方程（17），（18）的軌線是由式（20）定義的曲線族，我們來證明這些曲線是封閉的。引理1當(dāng)時(shí)，方程（20）定義了一組封閉曲線。證明我們首先來確定當(dāng)時(shí)函數(shù)和的性狀。利用微積分方法可以作出和的圖形。如下圖所示。若它們的極大值分別記作和，則不難確定滿足，（21），（22）顯然，僅當(dāng)（20）式右端常數(shù)時(shí)相軌線才有定義。當(dāng)時(shí)，，，將式（21）和（22）與（19）式比較可知正是平衡點(diǎn)，所以是相軌線的退化點(diǎn)。為了考察時(shí)軌線的形狀，我們只需考慮的情況，其中。首先注意到：方程具有一個(gè)解和另一個(gè)解。因此，當(dāng)或時(shí)，方程沒有解。當(dāng)或時(shí)，這個(gè)方程具有唯一的解，而對于，則具有兩個(gè)解和。較小的解總是小于，較大的解總是大于。當(dāng)趨于或時(shí)，和都趨向于。因此，當(dāng)和是正數(shù)時(shí)，由（20）所定義的曲線都是封閉的。而且，這些封閉曲線中的每一條（除和以外），都不含（17）和（18）的任何平衡點(diǎn)。所以（17）和（18）的具有初始條件，的所有的解，都是時(shí)間的周期函數(shù)。也就是說，（17）和（18）的具有初始條件，的每一個(gè)解，都具有這樣的性質(zhì)：，，其中是某一正數(shù)。D'Ancona所用的數(shù)據(jù)實(shí)際上是捕食者的百分比在每一年中的平均值。因此，為了把這些數(shù)據(jù)同方程組（17）和（18）的結(jié)果進(jìn)行比較，對于（17）和（18）的任何解，，我們必須算出和的“平均值”。值得注意的是，即使還沒有準(zhǔn)確地求得和，我們?nèi)匀荒軌蛩愠鲞@些平均值。引理2設(shè)，是（17）和（18）的周期解，其周期，和的平均值定義為，這時(shí)，，。換句話說，和的平均值是平衡解。證明把（17）的兩端除以，得到，于是由于因此，，于是，。類似地，把（18）的兩端除以，由0到積分，我們得到。下面，我們考慮漁業(yè)對于上述數(shù)學(xué)模型的影響。注意到漁業(yè)使得食餌總數(shù)以速率減少，而使得捕食者的總數(shù)以速率減少。常數(shù)反映漁業(yè)的水平；即反映了海上的漁船數(shù)和下水的網(wǎng)數(shù)。因此，真實(shí)的狀態(tài)由下列修正的微分方程組來描述：（23）這個(gè)方程組同（17），（18）完全一樣（當(dāng)時(shí)），只是其中換成，而換成。因此，現(xiàn)在和的平均值是，平均說來，中等捕魚量實(shí)際上會(huì)增加食餌的數(shù)目，而減少捕食者的數(shù)目。相反，捕魚量的降低，平均說來，會(huì)增加捕食者的數(shù)目，而減少食餌的數(shù)目。這個(gè)值得注意的結(jié)果稱為Volterra原理，它解釋了D'Ancona的數(shù)據(jù)。值得注意的是Volterra模型是非常粗糙的，有興趣的讀者可以作進(jìn)一步的探討。習(xí)題十四1.單棵樹木的商品價(jià)值是由這棵樹能夠生產(chǎn)的木材體積和質(zhì)量所決定的。顯然依賴于樹木的年齡。假設(shè)曲線已知，為樹木砍伐成本。試給出砍伐樹木（更確切地說砍伐相同年齡的樹木）的最優(yōu)年齡。如果考慮到森林輪種問題，即一旦樹木從某一處砍掉，這塊土地便可以種植新樹，假定各輪種周期具有相等的長度，試建模討論最優(yōu)砍伐輪種的森林管理策略的問題。2.如果兩個(gè)種群都能獨(dú)立生存，共處時(shí)又能相互提供食物，試建立種群依存模型并討論平衡點(diǎn)及穩(wěn)定性，解釋穩(wěn)定的意義。3.如果兩個(gè)種群都不能獨(dú)立生存，但共處時(shí)可以相互提供食物，試建模以討論共處的可能性。4.如果在食餌—捕食者系統(tǒng)中，捕食者掠奪的對象只是成年的食餌，而未成年的食餌因體積太小免遭捕獲。在適當(dāng)?shù)募僭O(shè)下建立這三者之間關(guān)系的模型，求平衡點(diǎn)。第十五章常微分方程的解法建立微分方程只是解決問題的第一步，通常需要求出方程的解來說明實(shí)際現(xiàn)象，并加以檢驗(yàn)。如果能得到解析形式的解固然是便于分析和應(yīng)用的，但是我們知道，只有線性常系數(shù)微分方程，并且自由項(xiàng)是某些特殊類型的函數(shù)時(shí)，才可以肯定得到這樣的解，而絕大多數(shù)變系數(shù)方程、非線性方程都是所謂“解不出來”的，即使看起來非常簡單的方程如，于是對于用微分方程解決實(shí)際問題來說，數(shù)值解法就是一個(gè)十分重要的手段。§1常微分方程的離散化下面主要討論一階常微分方程的初值問題，其一般形式是(1)在下面的討論中，我們總假定函數(shù)連續(xù)，且關(guān)于滿足李普希茲(Lipschitz)條件，即存在常數(shù)，使得這樣，由常微分方程理論知，初值問題(1)的解必定存在唯一。所謂數(shù)值解法，就是求問題（1）的解在若干點(diǎn)處的近似值的方法，稱為問題（1）的數(shù)值解，稱為由到的步長。今后如無特別說明，我們總?cè)〔介L為常量。建立數(shù)值解法，首先要將微分方程離散化，一般采用以下幾種方法：（=1\*romani）用差商近似導(dǎo)數(shù)若用向前差商代替代入（1）中的微分方程，則得化簡得如果用的近似值代入上式右端，所得結(jié)果作為的近似值，記為，則有（2）這樣，問題（1）的近似解可通過求解下述問題（3）得到，按式（3）由初值可逐次算出。式（3）是個(gè)離散化的問題，稱為差分方程初值問題。需要說明的是，用不同的差商近似導(dǎo)數(shù)，將得到不同的計(jì)算公式。（=2\*romanii）用數(shù)值積分方法將問題（1）的解表成積分形式，用數(shù)值積分方法離散化。例如，對微分方程兩端積分，得(4)右邊的積分用矩形公式或梯形公式計(jì)算。（=3\*romaniii）Taylor多項(xiàng)式近似將函數(shù)在處展開，取一次Taylor多項(xiàng)式近似，則得再將的近似值代入上式右端，所得結(jié)果作為的近似值，得到離散化的計(jì)算公式以上三種方法都是將微分方程離散化的常用方法，每一類方法又可導(dǎo)出不同形式的計(jì)算公式。其中的Taylor展開法，不僅可以得到求數(shù)值解的公式，而且容易估計(jì)截?cái)嗾`差?！?歐拉（Euler）方法Euler方法Euler方法就是用差分方程初值問題（3）的解來近似微分方程初值問題（1）的解，即由公式（3）依次算出的近似值。這組公式求問題（1）的數(shù)值解稱為向前Euler公式。如果在微分方程離散化時(shí)，用向后差商代替導(dǎo)數(shù)，即，則得計(jì)算公式（5）用這組公式求問題（1）的數(shù)值解稱為向后Euler公式。向后Euler法與Euler法形式上相似，但實(shí)際計(jì)算時(shí)卻復(fù)雜得多。向前Euler公式是顯式的，可直接求解。向后Euler公式的右端含有，因此是隱式公式，一般要用迭代法求解，迭代公式通常為（6）2.2Euler方法的誤差估計(jì)對于向前Euler公式（3）我們看到，當(dāng)時(shí)公式右端的都是近似的，所以用它計(jì)算的會(huì)有累積誤差，分析累積誤差比較復(fù)雜，這里先討論比較簡單的所謂局部截?cái)嗾`差。假定用（3）式時(shí)右端的沒有誤差，即，那么由此算出（7）局部截?cái)嗾`差指的是，按（7）式計(jì)算由到這一步的計(jì)算值與精確值之差。為了估計(jì)它，由Taylor展開得到的精確值是（8）（7）、（8）兩式相減（注意到）得（9）即局部截?cái)嗾`差是階的，而數(shù)值算法的精度定義為：若一種算法的局部截?cái)嗾`差為，則稱該算法具有階精度。顯然越大，方法的精度越高。式（9）說明，向前Euler方法是一階方法，因此它的精度不高?！?改進(jìn)的Euler方法3.1梯形公式利用數(shù)值積分方法將微分方程離散化時(shí)，若用梯形公式計(jì)算式(4)中之右端積分，即并用代替，則得計(jì)算公式這就是求解初值問題（1）的梯形公式。直觀上容易看出，用梯形公式計(jì)算數(shù)值積分要比矩形公式好。梯形公式為二階方法。梯形公式也是隱式格式，一般需用迭代法求解，迭代公式為(10)由于函數(shù)關(guān)于滿足Lipschitz條件，容易看出其中為Lipschitz常數(shù)。因此，當(dāng)時(shí)，迭代收斂。但這樣做計(jì)算量較大。如果實(shí)際計(jì)算時(shí)精度要求不太高，用公式（10）求解時(shí)，每步可以只迭代一次，由此導(dǎo)出一種新的方法—改進(jìn)Euler法。3.2改進(jìn)Euler法按式（5）計(jì)算問題（1）的數(shù)值解時(shí)，如果每步只迭代一次，相當(dāng)于將Euler公式與梯形公式結(jié)合使用：先用Euler公式求的一個(gè)初步近似值，稱為預(yù)測值，然后用梯形公式校正求得近似值，即（11）式（11）稱為由Euler公式和梯形公式得到的預(yù)測—校正系統(tǒng)，也叫改進(jìn)Euler法。為便于編制程序上機(jī)，式（11）常改寫成（12）改進(jìn)Euler法是二階方法?！?龍格—庫塔（Runge—Kutta）方法回到Euler方法的基本思想—用差商代替導(dǎo)數(shù)—上來。實(shí)際上，按照微分中值定理應(yīng)有注意到方程就有（13）不妨記，稱為區(qū)間上的平均斜率?？梢娊o出一種斜率，（13）式就對應(yīng)地導(dǎo)出一種算法。向前Euler公式簡單地取為，精度自然很低。改進(jìn)的Euler公式可理解為取，的平均值，其中，這種處理提高了精度。如上分析啟示我們，在區(qū)間內(nèi)多取幾個(gè)點(diǎn)，將它們的斜率加權(quán)平均作為，就有可能構(gòu)造出精度更高的計(jì)算公式。這就是龍格—庫塔方法的基本思想。4.1首先不妨在區(qū)間內(nèi)仍取2個(gè)點(diǎn)，仿照（13）式用以下形式試一下（14）其中為待定系數(shù)，看看如何確定它們使（14）式的精度盡量高。為此我們分析局部截?cái)嗾`差，因?yàn)?，所以?4）可以化為（15）其中在點(diǎn)作了Taylor展開。（15）式又可表為注意到中，，可見為使誤差，只須令，，（16）待定系數(shù)滿足（16）的（15）式稱為2階龍格—庫塔公式。由于（16）式有4個(gè)未知數(shù)而只有3個(gè)方程，所以解不唯一。不難發(fā)現(xiàn)，若令，，即為改進(jìn)的Euler公式。可以證明，在內(nèi)只取2點(diǎn)的龍格—庫塔公式精度最高為2階。4.24階龍格—庫塔公式要進(jìn)一步提高精度，必須取更多的點(diǎn)，如取4點(diǎn)構(gòu)造如下形式的公式：（17）其中待定系數(shù)共13個(gè)，經(jīng)過與推導(dǎo)2階龍格—庫塔公式類似、但更復(fù)雜的計(jì)算，得到使局部誤差的11個(gè)方程。取既滿足這些方程、又較簡單的一組，可得（18）這就是常用的4階龍格—庫塔方法（簡稱RK方法）?！?線性多步法以上所介紹的各種數(shù)值解法都是單步法，這是因?yàn)樗鼈冊谟?jì)算時(shí)，都只用到前一步的值，單步法的一般形式是（19）其中稱為增量函數(shù)，例如Euler方法的增量函數(shù)為，改進(jìn)Euler法的增量函數(shù)為如何通過較多地利用前面的已知信息，如，來構(gòu)造高精度的算法計(jì)算，這就是多步法的基本思想。經(jīng)常使用的是線性多步法。讓我們試驗(yàn)一下，即利用計(jì)算的情況。從用數(shù)值積分方法離散化方程的（4）式出發(fā)，記，，式中被積函數(shù)用二節(jié)點(diǎn)，的插值公式得到（因，所以是外插。（20）此式在區(qū)間上積分可得于是得到（21）注意到插值公式（20）的誤差項(xiàng)含因子，在區(qū)間上積分后得出，故公式（21）的局部截?cái)嗾`差為，精度比向前Euler公式提高1階。若取可以用類似的方法推導(dǎo)公式，如對于有（22）其局部截?cái)嗾`差為。如果將上面代替被積函數(shù)用的插值公式由外插改為內(nèi)插，可進(jìn)一步減小誤差。內(nèi)插法用的是，取時(shí)得到的是梯形公式，取時(shí)可得（23）與（22）式相比，雖然其局部截?cái)嗾`差仍為，但因它的各項(xiàng)系數(shù)（絕對值）大為減小，誤差還是小了。當(dāng)然，（23）式右端的未知，需要如同向后Euler公式一樣，用迭代或校正的辦法處理?！?一階微分方程組與高階微分方程的數(shù)值解法6.1一階微分方程組的數(shù)值解法設(shè)有一階微分方程組的初值問題（24）若記，，，則初值問題（24）可寫成如下向量形式（25）如果向量函數(shù)在區(qū)域：連續(xù)，且關(guān)于滿足Lipschitz條件，即存在，使得對，，都有那么問題（25）在上存在唯一解。問題（25）與（1）形式上完全相同，故對初值問題（1）所建立的各種數(shù)值解法可全部用于求解問題（25）。6.2高階微分方程的數(shù)值解法高階微分方程的初值問題可以通過變量代換化為一階微分方程組初值問題。設(shè)有階常微分方程初值問題（26）引入新變量，問題（26）就化為一階微分方程初值問題（27）然后用6.1中的數(shù)值方法求解問題（27），就可以得到問題（26）的數(shù)值解。最后需要指出的是，在化學(xué)工程及自動(dòng)控制等領(lǐng)域中，所涉及的常微分方程組初值問題常常是所謂的“剛性”問題。具體地說，對一階線性微分方程組（28）其中為階方陣。若矩陣的特征值滿足關(guān)系則稱方程組（28）為剛性方程組或Stiff方程組，稱數(shù)為剛性比。對剛性方程組，用前面所介紹的方法求解，都會(huì)遇到本質(zhì)上的困難，這是由數(shù)值方法本身的穩(wěn)定性限制所決定的。理論上的分析表明，求解剛性問題所選用的數(shù)值方法最好是對步長不作任何限制。§7Matlab解法7.1Matlab數(shù)值解7.1.1非剛性常微分方程的解法Matlab的工具箱提供了幾個(gè)解非剛性常微分方程的功能函數(shù)，如ode45，ode23，ode113，其中ode45采用四五階RK方法，是解非剛性常微分方程的首選方法，ode23采用二三階RK方法，ode113采用的是多步法，效率一般比ode45高。Matlab的工具箱中沒有Euler方法的功能函數(shù)。（=1\*ROMANI）對簡單的一階方程的初值問題改進(jìn)的Euler公式為我們自己編寫改進(jìn)的Euler方法函數(shù)eulerpro.m如下：function[x,y]=eulerpro(fun,x0,xfinal,y0,n);ifnargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;fori=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end例1用改進(jìn)的Euler方法求解，，解編寫函數(shù)文件doty.m如下：functionf=doty(x,y);f=-2*y+2*x^2+2*x;在Matlab命令窗口輸入：[x,y]=eulerpro('doty',0,0.5,1,10)即可求得數(shù)值解。（=2\*ROMANII）ode23，ode45，ode113的使用Matlab的函數(shù)形式是[t,y]=solver('F',tspan，y0)這里solver為ode45，ode23，ode113，輸入?yún)?shù)F是用M文件定義的微分方程右端的函數(shù)。tspan=[t0,tfinal]是求解區(qū)間，y0是初值。例2用RK方法求解，，解同樣地編寫函數(shù)文件doty.m如下：functionf=doty(x,y);f=-2*y+2*x^2+2*x;在Matlab命令窗口輸入：[x,y]=ode45('doty',0,0.5,1)即可求得數(shù)值解。7.1.2剛性常微分方程的解法Matlab的工具箱提供了幾個(gè)解剛性常微分方程的功能函數(shù)，如ode15s，ode23s，ode23t，ode23tb，這些函數(shù)的使用同上述非剛性微分方程的功能函數(shù)。7.1.3高階微分方程解法例3考慮初值問題解（=1\*romani）如果設(shè)，那么初值問題可以寫成的形式，其中。（=2\*romanii）把一階方程組寫成接受兩個(gè)參數(shù)和，返回一個(gè)列向量的M文件F.m：functiondy=F(t,y);dy=[y(2);y(3);3*y(3)+y(2)*y(1)];注意：盡管不一定用到參數(shù)和，M—文件必須接受此兩參數(shù)。這里向量必須是列向量。（=3\*romaniii）用Matlab解決此問題的函數(shù)形式為[T,Y]=solver('F',tspan,y0)這里solver為ode45、ode23、ode113，輸入?yún)?shù)F是用M文件定義的常微分方程組，tspan=[t0tfinal]是求解區(qū)間，y0是初值列向量。在Matlab命令窗口輸入[T,Y]=ode45('F',[01],[0;1;-1])就得到上述常微分方程的數(shù)值解。這里Y和時(shí)刻T是一一對應(yīng)的，Y(:,1)是初值問題的解，Y(:,2)是解的導(dǎo)數(shù)，Y(:,3)是解的二階導(dǎo)數(shù)。例4求vanderPol方程的數(shù)值解，這里是一參數(shù)。解（=1\*romani）化成常微分方程組。設(shè)，則有（=2\*romanii）書寫M文件（對于）vdp1.m:functiondy=vdp1(t,y);dy=[y(2);(1-y(1)^2)*y(2)-y(1)];（=3\*romaniii）調(diào)用Matlab函數(shù)。對于初值，解為[T,Y]=ode45('vdp1',[020],[2;0]);（=4\*romaniv）觀察結(jié)果。利用圖形輸出解的結(jié)果：plot(T,Y(:,1),'-',T,Y(:,2),'--')title('SolutionofvanderPolEquation,mu=1');xlabel('timet');ylabel('solutiony');legend('y1','y2');例5vanderPol方程，（剛性）解(=1\*romani)書寫M文件vdp1000.m:functiondy=vdp1000(t,y);dy=[y(2);1000*(1-y(1)^2)*y(2)-y(1)];（=2\*romanii）觀察結(jié)果[t,y]=ode15s('vdp1000',[03000],[2;0]);plot(t,y(:,1),'o')title('SolutionofvanderPolEquation,mu=1000');xlabel('timet');ylabel('solutiony(:,1)');7.2常微分方程的解析解在Matlab中，符號運(yùn)算工具箱提供了功能強(qiáng)大的求解常微分方程的符號運(yùn)算命令dsolve。常微分方程在Matlab中按如下規(guī)定重新表達(dá)：符號D表示對變量的求導(dǎo)。Dy表示對變量y求一階導(dǎo)數(shù)，當(dāng)需要求變量的n階導(dǎo)數(shù)時(shí)，用Dn表示，D4y表示對變量y求4階導(dǎo)數(shù)。由此，常微分方程在Matlab中，將寫成'D2y+2*Dy=y'。7.2.1求解常微分方程的通解無初邊值條件的常微分方程的解就是該方程的通解。其使用格式為：dsolve('diff_equation')dsolve('diff_equation'，'var')式中diff_equation為待解的常微分方程，第1種格式將以變量t為自變量進(jìn)行求解，第2種格式則需定義自變量var。例6試解常微分方程解編寫程序如下：symsxydiff_equ='x^2+y+(x-2*y)*Dy=0';dsolve(diff_equ,'x')7.2.2求解常微分方程的初邊值問題求解帶有初邊值條件的常微分方程的使用格式為：dsolve('diff_equation'，'condition1，condition2，…'，'var')其中condition1，condition2，…即為微分方程的初邊值條件。例7試求微分方程，的解。解編寫程序如下：y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x')7.2.3求解常微分方程組求解常微分方程組的命令格式為：dsolve('diff_equ1，diff_equ2，…'，'var')dsolve('diff_equ1，diff_equ2，…'，'condition1，condition2，…'，'var')第1種格式用于求解方程組的通解，第2種格式可以加上初邊值條件，用于具體求解。例8試求常微分方程組：的通解和在初邊值條件為的解。解編寫程序如下：clc,clearequ1='D2f+3*g=sin(x)';equ2='Dg+Df=cos(x)';[general_f,general_g]=dsolve(equ1,equ2,'x')[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')7.2.4求解線性常微分方程組（=1\*romani）一階齊次線性微分方程組，，這里的’表示對求導(dǎo)數(shù)。是它的基解矩陣。，的解為。例9試解初值問題，解編寫程序如下：symsta=[2,1,3;0,2,-1;0,0,2];x0=[1;2;1];x=expm(a*t)*x0（=2\*romanii）非齊次線性方程組由參數(shù)變易法可求得初值問題，的解為.例10試解初值問題，。解編寫程序如下：clc,clearsymstsa=[1,0,0;2,1,-2;3,2,1];ft=[0;0;exp(t)*cos(2*t)];x0=[0;1;1];x=expm(a*t)*x0+int(expm(a*(t-s))*subs(ft,s),s,0,t);x=simple(x)習(xí)題十五1.用歐拉方法和龍格—庫塔方法求微分方程數(shù)值解，畫出解的圖形，對結(jié)果進(jìn)行分析比較。（=1\*romani）（Bessel方程，令，精確解。（=2\*romanii），冪級數(shù)解2.一只小船渡過寬為的河流，目標(biāo)是起點(diǎn)正對著的另一岸點(diǎn)。已知河水流速與船在靜水中的速度之比為。（=1\*romani）建立小船航線的方程，求其解析解。（=2\*romanii）設(shè)m，m/s，m/s，用數(shù)值解法求渡河所需時(shí)間、任意時(shí)刻小船的位置及航行曲線，作圖，并與解析解比較。第二十三章現(xiàn)代優(yōu)化算法簡介§1現(xiàn)代優(yōu)化算法簡介現(xiàn)代優(yōu)化算法是80年代初興起的啟發(fā)式算法。這些算法包括禁忌搜索（tabusearch），模擬退火（simulatedannealing），遺傳算法（geneticalgorithms），人工神經(jīng)網(wǎng)絡(luò)（neuralnetworks）。它們主要用于解決大量的實(shí)際應(yīng)用問題。目前，這些算法在理論和實(shí)際應(yīng)用方面得到了較大的發(fā)展。無論這些算法是怎樣產(chǎn)生的，它們有一個(gè)共同的目標(biāo)－求NP-hard組合優(yōu)化問題的全局最優(yōu)解。雖然有這些目標(biāo)，但NP-hard理論限制它們只能以啟發(fā)式的算法去求解實(shí)際問題。啟發(fā)式算法包含的算法很多，例如解決復(fù)雜優(yōu)化問題的蟻群算法（AntColonyAlgorithms）。有些啟發(fā)式算法是根據(jù)實(shí)際問題而產(chǎn)生的，如解空間分解、解空間的限制等；另一類算法是集成算法，這些算法是諸多啟發(fā)式算法的合成?，F(xiàn)代優(yōu)化算法解決組合優(yōu)化問題，如TSP（TravelingSalesmanProblem）問題，QAP（QuadraticAssignmentProblem）問題，JSP（Job-shopSchedulingProblem）問題等效果很好。本章我們只介紹模擬退火算法，初步介紹一下蟻群算法，其它優(yōu)化算法可以參看相關(guān)的參考資料。§2模擬退火算法2.1算法簡介模擬退火算法得益于材料的統(tǒng)計(jì)力學(xué)的研究成果。統(tǒng)計(jì)力學(xué)表明材料中粒子的不同結(jié)構(gòu)對應(yīng)于粒子的不同能量水平。在高溫條件下，粒子的能量較高，可以自由運(yùn)動(dòng)和重新排列。在低溫條件下，粒子能量較低。如果從高溫開始，非常緩慢地降溫（這個(gè)過程被稱為退火），粒子就可以在每個(gè)溫度下達(dá)到熱平衡。當(dāng)系統(tǒng)完全被冷卻時(shí)，最終形成處于低能狀態(tài)的晶體。如果用粒子的能量定義材料的狀態(tài)，Metropolis算法用一個(gè)簡單的數(shù)學(xué)模型描述了退火過程。假設(shè)材料在狀態(tài)之下的能量為，那么材料在溫度時(shí)從狀態(tài)進(jìn)入狀態(tài)就遵循如下規(guī)律：（1）如果，接受該狀態(tài)被轉(zhuǎn)換。（2）如果，則狀態(tài)轉(zhuǎn)換以如下概率被接受：其中是物理學(xué)中的波爾茲曼常數(shù)，是材料溫度。在某一個(gè)特定溫度下，進(jìn)行了充分的轉(zhuǎn)換之后，材料將達(dá)到熱平衡。這時(shí)材料處于狀態(tài)的概率滿足波爾茲曼分布：其中表示材料當(dāng)前狀態(tài)的隨機(jī)變量，表示狀態(tài)空間集合。顯然其中表示集合中狀態(tài)的數(shù)量。這表明所有狀態(tài)在高溫下具有相同的概率。而當(dāng)溫度下降時(shí)，其中且。上式表明當(dāng)溫度降至很低時(shí)，材料會(huì)以很大概率進(jìn)入最小能量狀態(tài)。假定我們要解決的問題是一個(gè)尋找最小值的優(yōu)化問題。將物理學(xué)中模擬退火的思想應(yīng)用于優(yōu)化問題就可以得到模擬退火尋優(yōu)方法?？紤]這樣一個(gè)組合優(yōu)化問題：優(yōu)化函數(shù)為，其中，它表示優(yōu)化問題的一個(gè)可行解，，表示函數(shù)的定義域。表示的一個(gè)鄰域集合。首先給定一個(gè)初始溫度和該優(yōu)化問題的一個(gè)初始解，并由生成下一個(gè)解，是否接受作為一個(gè)新解依賴于下面概率：換句話說，如果生成的解的函數(shù)值比前一個(gè)解的函數(shù)值更小，則接受作為一個(gè)新解。否則以概率接受作為一個(gè)新解。泛泛地說，對于某一個(gè)溫度和該優(yōu)化問題的一個(gè)解，可以生成。接受作為下一個(gè)新解的概率為：（1）在溫度下，經(jīng)過很多次的轉(zhuǎn)移之后，降低溫度，得到。在下重復(fù)上述過程。因此整個(gè)優(yōu)化過程就是不斷尋找新解和緩慢降溫的交替過程。最終的解是對該問題尋優(yōu)的結(jié)果。我們注意到，在每個(gè)下，所得到的一個(gè)新狀態(tài)完全依賴于前一個(gè)狀態(tài)，可以和前面的狀態(tài)無關(guān)，因此這是一個(gè)馬爾可夫過程。使用馬爾可夫過程對上述模擬退火的步驟進(jìn)行分析，結(jié)果表明：從任何一個(gè)狀態(tài)生成的概率，在中是均勻分布的，且新狀態(tài)被接受的概率滿足式（1），那么經(jīng)過有限次的轉(zhuǎn)換，在溫度下的平衡態(tài)的分布由下式給出：（2）當(dāng)溫度降為0時(shí)，的分布為：并且這說明如果溫度下降十分緩慢，而在每個(gè)溫度都有足夠多次的狀態(tài)轉(zhuǎn)移，使之在每一個(gè)溫度下達(dá)到熱平衡，則全局最優(yōu)解將以概率1被找到。因此可以說模擬退火算法可以找到全局最優(yōu)解。在模擬退火算法中應(yīng)注意以下問題：（1）理論上，降溫過程要足夠緩慢，要使得在每一溫度下達(dá)到熱平衡。但在計(jì)算機(jī)實(shí)現(xiàn)中，如果降溫速度過緩，所得到的解的性能會(huì)較為令人滿意，但是算法會(huì)太慢，相對于簡單的搜索算法不具有明顯優(yōu)勢。如果降溫速度過快，很可能最終得不到全局最優(yōu)解。因此使用時(shí)要綜合考慮解的性能和算法速度，在兩者之間采取一種折衷。（2）要確定在每一溫度下狀態(tài)轉(zhuǎn)換的結(jié)束準(zhǔn)則。實(shí)際操作可以考慮當(dāng)連續(xù)次的轉(zhuǎn)換過程沒有使?fàn)顟B(tài)發(fā)生變化時(shí)結(jié)束該溫度下的狀態(tài)轉(zhuǎn)換。最終溫度的確定可以提前定為一個(gè)較小的值，或連續(xù)幾個(gè)溫度下轉(zhuǎn)換過程沒有使?fàn)顟B(tài)發(fā)生變化算法就結(jié)束。（3）選擇初始溫度和確定某個(gè)可行解的鄰域的方法也要恰當(dāng)。2.2應(yīng)用舉例例已知敵方100個(gè)目標(biāo)的經(jīng)度、緯度如下：經(jīng)度緯度經(jīng)度緯度經(jīng)度緯度經(jīng)度緯度53.712115.304651.17580.032246.325328.275330.33136.934856.543221.418810.819816.252922.789123.104510.158412.481920.105015.45621.94510.205726.495122.122131.48478.964026.241818.176044.035613.540128.983625.987938.472220.173128.269429.001132.19105.869936.486329.72840.971828.14778.958624.663516.561823.614310.559715.117850.211110.29448.15199.532522.107518.55690.121518.872648.207716.888931.949917.63090.77320.465647.413423.778341.86713.566743.54743.906153.352426.725630.816513.459527.71335.070623.92227.630651.961222.851112.793815.73074.95688.366921.505124.090915.254827.21116.20705.144249.243016.704417.116820.035434.168822.75719.44023.920011.581214.567752.11810.40889.555911.421924.45096.563426.721328.566737.584816.847435.66199.933324.46543.16440.77756.957614.470313.636819.866015.12243.16164.242818.524514.359858.684927.148539.516816.937156.508913.709052.521115.795738.43008.464851.818123.01598.998323.644050.115623.781613.79091.951034.057423.396023.06248.431919.98575.790240.880114.297858.828914.522918.66356.743652.842327.288039.949429.511447.509924.066410.112127.266228.781227.66598.083127.67059.155614.130453.79890.219933.64900.39801.349616.835949.98166.082819.363517.662236.954523.026515.732019.569711.511817.388444.039816.263539.713928.42036.990923.180438.339219.995024.654319.605736.998024.39924.15913.185340.140020.303023.98769.403041.108427.7149我方有一個(gè)基地，經(jīng)度和緯度為（70,40）。假設(shè)我方飛機(jī)的速度為1000公里/小時(shí)。我方派一架飛機(jī)從基地出發(fā)，偵察完敵方所有目標(biāo)，再返回原來的基地。在敵方每一目標(biāo)點(diǎn)的偵察時(shí)間不計(jì)，求該架飛機(jī)所花費(fèi)的時(shí)間（假設(shè)我方飛機(jī)巡航時(shí)間可以充分長）。這是一個(gè)旅行商問題。我們依次給基地編號為1，敵方目標(biāo)依次編號為2，3，…，101，最后我方基地再重復(fù)編號為102（這樣便于程序中計(jì)算）。距離矩陣，其中表示表示兩點(diǎn)的距離，，這里為實(shí)對稱矩陣。則問題是求一個(gè)從點(diǎn)1出發(fā)，走遍所有中間點(diǎn)，到達(dá)點(diǎn)102的一個(gè)最短路徑。上面問題中給定的是地理坐標(biāo)（經(jīng)度和緯度），我們必須求兩點(diǎn)間的實(shí)際距離。設(shè)兩點(diǎn)的地理坐標(biāo)分別為，，過兩點(diǎn)的大圓的劣弧長即為兩點(diǎn)的實(shí)際距離。以地心為坐標(biāo)原點(diǎn)，以赤道平面為平面，以0度經(jīng)線圈所在的平面為平面建立三維直角坐標(biāo)系。則兩點(diǎn)的直角坐標(biāo)分別為：其中為地球半徑。兩點(diǎn)的實(shí)際距離，化簡得。求解的模擬退火算法描述如下：（1）解空間解空間可表為{}的所有固定起點(diǎn)和終點(diǎn)的循環(huán)排列集合，即其中每一個(gè)循環(huán)排列表示偵察100個(gè)目標(biāo)的一個(gè)回路，表示在第次偵察點(diǎn)，初始解可選為，本文中我們使用MonteCarlo方法求得一個(gè)較好的初始解。（2）目標(biāo)函數(shù)此時(shí)的目標(biāo)函數(shù)為偵察所有目標(biāo)的路徑長度或稱代價(jià)函數(shù)。我們要求而一次迭代由下列三步構(gòu)成：（3）新解的產(chǎn)生=1\*GB3①2變換法任選序號（）交換與之間的順序，此時(shí)的新路徑為：②3變換法任選序號和，將和之間的路徑插到之后，對應(yīng)的新路徑為（設(shè)）（4）代價(jià)函數(shù)差對于2變換法，路徑差可表示為（5）接受準(zhǔn)則如果，則接受新的路徑。否則，以概率接受新的路徑，即若大于0到1之間的隨機(jī)數(shù)則接受。（6）降溫利用選定的降溫系數(shù)進(jìn)行降溫即：，得到新的溫度，這里我們?nèi)?。?）結(jié)束條件用選定的終止溫度，判斷退火過程是否結(jié)束。若，算法結(jié)束，輸出當(dāng)前狀態(tài)。我們編寫如下的matlab程序如下：clc,clearloadsj.txt%加載敵方100個(gè)目標(biāo)的數(shù)據(jù)，數(shù)據(jù)按照表格中的位置保存在純文本文件sj.txt中x=sj(:,1:2:8);x=x(:);y=sj(:,2:2:8);y=y(:);sj=[xy];d1=[70,40];sj=[d1;sj;d1];sj=sj*pi/180;%距離矩陣dd=zeros(102);fori=1:101forj=i+1:102temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));d(i,j)=6370*acos(temp);endendd=d+d';S0=[];Sum=inf;rand('state',sum(clock));forj=1:1000S=[11+randperm(100),102];temp=0;fori=1:101temp=temp+d(S(i),S(i+1));endiftemp<SumS0=S;Sum=temp;endende=0.1^30;L=20000;at=0.999;T=1;%退火過程fork=1:L%產(chǎn)生新解c=2+floor(100*rand(1,2));c=sort(c);c1=c(1);c2=c(2);%計(jì)算代價(jià)函數(shù)值df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1));%接受準(zhǔn)則ifdf<0S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;elseifexp(-df/T)>rand(1)S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;endT=T*at;ifT<ebreak;endend%輸出巡航路徑及路徑長度S0,Sum計(jì)算結(jié)果為44小時(shí)左右。其中的一個(gè)巡航路徑如下圖所示：§3蟻群算法3.1蟻群算法簡介蟻群是自然界中常見的一種生物，人們對螞蟻的關(guān)注大都是因?yàn)椤跋伻喊峒遥煲掠辍敝惖拿裰V。然而隨著近代仿生學(xué)的發(fā)展，這種似乎微不足道的小東西越來越多地受到學(xué)者們地關(guān)注。1991年意大利學(xué)者M(jìn).Dorigo等人首先提出了蟻群算法，人們開始了對蟻群的研究：相對弱小，功能并不強(qiáng)大的個(gè)體是如何完成復(fù)雜的工作的（如尋找到食物的最佳路徑并返回等）。在此基礎(chǔ)上一種很好的優(yōu)化算法逐步發(fā)展起來。蟻群算法的特點(diǎn)是模擬自然界中螞蟻的群體行為?？茖W(xué)家發(fā)現(xiàn)，蟻群總是能夠發(fā)現(xiàn)從蟻巢到食物源的最短路徑。經(jīng)研究發(fā)現(xiàn)，螞蟻在行走過的路上留下一種揮發(fā)性的激素，螞蟻就是通過這種激素進(jìn)行信息交流。螞蟻趨向于走激素積累較多的路徑。找到最短路徑的螞蟻總是最早返回巢穴，從而在路上留下了較多的激素。由于最短路徑上積累了較多的激素，選擇這條路徑的螞蟻就會(huì)越來越多，到最后所有的螞蟻都會(huì)趨向于選擇這條最短路徑。基于螞蟻這種行為而提出的蟻群算法具有群體合作，正反饋選擇，并行計(jì)算等三大特點(diǎn)，并且可以根據(jù)需要為人工蟻加入前瞻、回溯等自然蟻所沒有的特點(diǎn)。在使用蟻群算法求解現(xiàn)實(shí)問題時(shí)，先生成具有一定數(shù)量螞蟻的蟻群，讓每一只螞蟻建立一個(gè)解或解的一部分，每只人工蟻從問題的初始狀態(tài)出發(fā)，根據(jù)“激素”濃度來選擇下一個(gè)要轉(zhuǎn)移到的狀態(tài)，直到建立起一個(gè)解，每只螞蟻根據(jù)所找到的解的好壞程度在所經(jīng)過的狀態(tài)上釋放與解的質(zhì)量成正比例的“激素”。之后，每只螞蟻又開始新的求解過程，直到尋找到滿意解。為避免停滯現(xiàn)象，引入了激素更新機(jī)制。3.2解決TSP問題的蟻群算法描述現(xiàn)以TSP問題的求解為例說明蟻群系統(tǒng)模型。首先引進(jìn)如下記號：為城市的個(gè)數(shù)；為蟻群中螞蟻的數(shù)量；為兩城市和之間距離；為時(shí)刻位于城市的螞蟻的個(gè)數(shù)，；為時(shí)刻邊弧的軌跡強(qiáng)度（即連線上殘留的信息量），且設(shè)（為常數(shù)），，；為時(shí)刻邊弧的能見度，反映由城市轉(zhuǎn)移到城市的期望程度。根據(jù)上述原理，螞蟻在運(yùn)動(dòng)過程中根據(jù)各條路徑上的信息量決定轉(zhuǎn)移方向。與真實(shí)蟻群系統(tǒng)不同，人工蟻群系統(tǒng)具有一定的記憶功能。隨著時(shí)間的推移，以前留下的信息逐漸消逝，經(jīng)個(gè)時(shí)刻，螞蟻完成一次循環(huán)，各路徑上信息量要作調(diào)整。由此得到下述的人工蟻群系統(tǒng)模型：1）設(shè)人工蟻群在并行地搜索TSP的解，并通過一種信息素做媒介相互通信，在每個(gè)結(jié)點(diǎn)上且和該結(jié)點(diǎn)相連的邊上以信息素量做搜索下一結(jié)點(diǎn)的試探依據(jù)，直到找到一個(gè)TSP問題的可行解。2）在時(shí)刻人工蟻由位置轉(zhuǎn)移至位置的轉(zhuǎn)移概率為（3）其中參數(shù)為軌跡的相對重要性（）；為能見度的相對重要性（）；為可行點(diǎn)集，即螞蟻下一步允許選擇的城市。分別反映了螞蟻在運(yùn)動(dòng)過程中所積累的信息及啟發(fā)式因子在螞蟻選擇路徑中所起的不同作用。3）當(dāng)個(gè)人工蟻按（3）式找到了可行解，則將各邊的信息量用下式修改。即調(diào)整信息量的軌跡強(qiáng)度更新方程為，（4）其中為第只螞蟻在本次循環(huán)中留在路徑上的信息量；為本次循環(huán)中路徑上的信息量的增量；參數(shù)為軌跡的持久性；為軌跡衰減度，表示信息消逝程度。對上述系統(tǒng)模型，采用人工蟻群方法求解的算法步驟可歸結(jié)為：step1：（為迭代步數(shù)或搜索次數(shù)）；各和的初始化；將個(gè)螞蟻置于個(gè)頂點(diǎn)上。step2：將各螞蟻的初始出發(fā)點(diǎn)置于當(dāng)前解集中；對每個(gè)螞蟻（）按概率轉(zhuǎn)移至下一頂點(diǎn)；將頂點(diǎn)置于當(dāng)前解集。step3：計(jì)算各螞蟻的目標(biāo)函數(shù)值，記錄當(dāng)前的最好解。step4：按更新方程修改軌跡強(qiáng)度。step5：，若預(yù)定的迭代次數(shù)且無退化行為（即找到的都是相同解），則轉(zhuǎn)step2。若為了簡化計(jì)算，增加處理較大規(guī)模的TSP問題的能力，則可將（4）式修改為：

，其中此處為本次最優(yōu)路線上的邊集。

3.3人工蟻群算法性能的討論人工蟻群算法是一種基于種群的進(jìn)化算法。作為一個(gè)新興的研究領(lǐng)域，雖它還遠(yuǎn)未像GA、SA等算法那樣形成系統(tǒng)的分析方法和堅(jiān)實(shí)的數(shù)學(xué)基礎(chǔ)，但目前已有一些基本結(jié)果。在M.Dorigo三種不同的模型中，循環(huán)路徑上信息量的增量不同。1）Ant-quantitysystem模型中，2）在Ant-densitysystem模型中，3）在Ant-cyclesystem模型中，其中是反映螞蟻所留軌跡數(shù)量的常數(shù)，表示第只螞蟻在本次循環(huán)中所走路徑的長度；且時(shí)，，。算法中模型1）、2）利用的是局部信息，模型3）利用的是整體信息。人工蟻群算法中，等參數(shù)對算法性能也有很大的影響。值的大小表明留在每個(gè)結(jié)點(diǎn)上的信息量受重視的程度，值越大，螞蟻選擇以前選過的點(diǎn)的可能性越大，但過大會(huì)使搜索過早陷于局部極小點(diǎn)；的大小表明啟發(fā)式信息受重視的程度；值會(huì)影響算法的收斂速度，過大會(huì)使算法收斂于局部極小值，過小又會(huì)影響算法的收斂速度，隨問題規(guī)模的增大的值也需要隨之變化；螞蟻的數(shù)目越多，算法的全局搜索能力越強(qiáng)，但數(shù)目加大將使算法的收斂速度減慢。第二十四章時(shí)間序列模型時(shí)間序列是按時(shí)間順序排列的、隨時(shí)間變化且相互關(guān)聯(lián)的數(shù)據(jù)序列。分析時(shí)間序列的方法構(gòu)成數(shù)據(jù)分析的一個(gè)重要領(lǐng)域，即時(shí)間序列分析。時(shí)間序列根據(jù)所研究的依據(jù)不同，可有不同的分類。1．按所研究的對象的多少分，有一元時(shí)間序列和多元時(shí)間序列。2．按時(shí)間的連續(xù)性可將時(shí)間序列分為離散時(shí)間序列和連續(xù)時(shí)間序列兩種。3．按序列的統(tǒng)計(jì)特性分，有平穩(wěn)時(shí)間序列和非平穩(wěn)時(shí)間序列。如果一個(gè)時(shí)間序列的概率分布與時(shí)間無關(guān)，則稱該序列為嚴(yán)格的（狹義的）平穩(wěn)時(shí)間序列。如果序列的一、二階矩存在，而且對任意時(shí)刻滿足：（1）均值為常數(shù)（2）協(xié)方差為時(shí)間間隔的函數(shù)。則稱該序列為寬平穩(wěn)時(shí)間序列，也叫廣義平穩(wěn)時(shí)間序列。我們以后所研究的時(shí)間序列主要是寬平穩(wěn)時(shí)間序列。4．按時(shí)間序列的分布規(guī)律來分，有高斯型時(shí)間序列和非高斯型時(shí)間序列?！?確定性時(shí)間序列分析方法概述時(shí)間序列預(yù)測技術(shù)就是通過對預(yù)測目標(biāo)自身時(shí)間序列的處理，來研究其變化趨勢的。一個(gè)時(shí)間序列往往是以下幾類變化形式的疊加或耦合。（1）長期趨勢變動(dòng)。它是指時(shí)間序列朝著一定的方向持續(xù)上升或下降，或停留在某一水平上的傾向，它反映了客觀事物的主要變化趨勢。（2）季節(jié)變動(dòng)。（3）循環(huán)變動(dòng)。通常是指周期為一年以上，由非季節(jié)因素引起的漲落起伏波形相似的波動(dòng)。（4）不規(guī)則變動(dòng)。通常它分為突然變動(dòng)和隨機(jī)變動(dòng)。通常用表示長期趨勢項(xiàng)，表示季節(jié)變動(dòng)趨勢項(xiàng)，表示循環(huán)變動(dòng)趨勢項(xiàng)，表示隨機(jī)干擾項(xiàng)。常見的確定性時(shí)間序列模型有以下幾種類型：（1）加法模型（2）乘法模型（3）混合模型其中是觀測目標(biāo)的觀測記錄，，。如果在預(yù)測時(shí)間范圍以內(nèi)，無突然變動(dòng)且隨機(jī)變動(dòng)的方差較小，并且有理由認(rèn)為過去和現(xiàn)在的演變趨勢將繼續(xù)發(fā)展到未來時(shí)，可用一些經(jīng)驗(yàn)方法進(jìn)行預(yù)測，具體方法如下：1.1移動(dòng)平均法設(shè)觀測序列為，取移動(dòng)平均的項(xiàng)數(shù)。一次移動(dòng)平均值計(jì)算公式為：二次移動(dòng)平均值計(jì)算公式為：當(dāng)預(yù)測目標(biāo)的基本趨勢是在某一水平上下波動(dòng)時(shí)，可用一次移動(dòng)平均方法建立預(yù)測模型：，，其預(yù)測標(biāo)準(zhǔn)誤差為：，最近期序列值的平均值作為未來各期的預(yù)測結(jié)果。一般取值范圍：。當(dāng)歷史序列的基本趨勢變化不大且序列中隨機(jī)變動(dòng)成分較多時(shí)，的取值應(yīng)較大一些。否則的取值應(yīng)小一些。在有確定的季節(jié)變動(dòng)周期的資料中，移動(dòng)平均的項(xiàng)數(shù)應(yīng)取周期長度。選擇最佳值的一個(gè)有效方法是，比較若干模型的預(yù)測誤差。均方預(yù)測誤差最小者為好。當(dāng)預(yù)測目標(biāo)的基本趨勢與某一線性模型相吻合時(shí)，常用二次移動(dòng)平均法，但序列同時(shí)存在線性趨勢與周期波動(dòng)時(shí)，可用趨勢移動(dòng)平均法建立預(yù)測模型：，其中，。例1某企業(yè)1月～11月份的銷售收入時(shí)間序列如下表所示。取，試用簡單一次滑動(dòng)平均法預(yù)測第12月份的銷售收入，并計(jì)算預(yù)測的標(biāo)準(zhǔn)誤差。月份123456銷售收入533.8574.6606.9649.8705.1772.0月份789101112銷售收入816.4892.7963.91015.11102.7解：首先計(jì)算出，由于，，則，預(yù)測的標(biāo)準(zhǔn)誤差為計(jì)算的Matlab程序如下：y=[533.8574.6606.9649.8705.1772.0...816.4892.7963.91015.11102.7];temp=cumsum(y);mt=(temp(4:11)-[0temp(1:7)])/4y12=mt(end)ythat=mt(1:end-1);fangcha=mean((y(5:11)-ythat).^2);sigma=sqrt(fangcha)1.2指數(shù)平滑法一次移動(dòng)平均實(shí)際上認(rèn)為最近期數(shù)據(jù)對未來值影響相同，都加權(quán)；而期以前的數(shù)據(jù)對未來值沒有影響，加權(quán)為0。但是，二次及更高次移動(dòng)平均數(shù)的權(quán)數(shù)卻不是，且次數(shù)越高，權(quán)數(shù)的結(jié)構(gòu)越復(fù)雜，但永遠(yuǎn)保持對稱的權(quán)數(shù)，即兩端項(xiàng)權(quán)數(shù)小，中間項(xiàng)權(quán)數(shù)大，不符合一般系統(tǒng)的動(dòng)態(tài)性。一般說來歷史數(shù)據(jù)對未來值的影響是隨時(shí)間間隔的增長而遞減的。所以，更切合實(shí)際的方法應(yīng)是對各期觀測值依時(shí)間順序進(jìn)行加權(quán)平均作為預(yù)測值。指數(shù)平滑法可滿足這一要求，而且具有簡單的遞推形式。設(shè)觀測序列為，為加權(quán)系數(shù)，，一次指數(shù)平滑公式為：（1）假定歷史序列無限長，則有（2）（2）式表明是全部歷史數(shù)據(jù)的加權(quán)平均，加權(quán)系數(shù)分別為；顯然有由于加權(quán)系數(shù)序列呈指數(shù)函數(shù)衰減，加權(quán)平均又能消除或減弱隨機(jī)干擾的影響，所以（2）稱為一次指數(shù)平滑，類似地，二次指數(shù)平滑公式為：（3）同理，三次指數(shù)平滑公式為：（4）一般次指數(shù)平滑公式為：