非線性回歸模型參數(shù)估計(jì)的輪回選擇算法

陳庭木 方兆偉 王寶祥 劉艷 邢運(yùn)高 徐波 徐大勇

摘要：生物學(xué)領(lǐng)域規(guī)律更多呈現(xiàn)不可線性化的純非線性方程關(guān)系，對(duì)此類方程的參數(shù)準(zhǔn)確估計(jì)能促進(jìn)生物學(xué)領(lǐng)域的應(yīng)用數(shù)學(xué)研究。但純非線性方程最優(yōu)擬合一直是應(yīng)用數(shù)學(xué)沒(méi)有完全解決的課題;受水稻輪回選擇育種啟發(fā)，以育種進(jìn)化方法構(gòu)建新的進(jìn)化算法，求解復(fù)雜非線性方程的擬合問(wèn)題。結(jié)果表明，通過(guò)Richard方程、房室模型中的二室模型試算，前者與NIST結(jié)果相同，后者明顯優(yōu)于原文結(jié)果;通過(guò)奶牛泌乳曲線中的Dijkstra、Wood方程各14組數(shù)據(jù)最優(yōu)擬合計(jì)算，表明本算法準(zhǔn)確有效，且計(jì)算結(jié)果優(yōu)于統(tǒng)計(jì)軟件SAS。本算法為生物學(xué)研究中更復(fù)雜的非線性方程擬合及其應(yīng)用提供了可能，為生物數(shù)學(xué)的深入發(fā)展提供了有力計(jì)算工具。

關(guān)鍵詞：最優(yōu)擬合;輪回選擇;非線性回歸

中圖分類號(hào)： S11+5文獻(xiàn)標(biāo)志碼： A

文章編號(hào)：1002-1302（2019）18-0253-07

收稿日期：2018-05-04

基金項(xiàng)目：現(xiàn)代農(nóng)業(yè)技術(shù)體系建設(shè)專項(xiàng)（編號(hào)：CARS-01-61）;江蘇省科技計(jì)劃重點(diǎn)項(xiàng)目（編號(hào)：SBE2017310472）。

作者簡(jiǎn)介：陳庭木（1977—），男，安徽蕪湖人，副研究員，主要從事水稻品質(zhì)遺傳育種與生物統(tǒng)計(jì)研究。E-mail：chentingmu@139.com。

通信作者：徐大勇，博士，研究員，主要從事水稻遺傳育種研究。E-mail：xudayong3030@sina.com。

一般線性回歸關(guān)系有成熟算法，能應(yīng)用最小二乘法準(zhǔn)確無(wú)偏估計(jì)回歸參數(shù)及其標(biāo)準(zhǔn)差，在農(nóng)業(yè)科學(xué)研究領(lǐng)域，變量間更多呈現(xiàn)非線性關(guān)系，非線性關(guān)系只有經(jīng)線性代換化為線性方程后才能應(yīng)用成熟的線性算法，且參數(shù)估計(jì)只是在線性化條件下最優(yōu)，以原非線性方程考量不是最優(yōu)[1-2]。對(duì)于不能化為線性方程的非線性方程，稱為純非線性方程。純非線性方程在生物學(xué)研究中更為普遍，如生物生長(zhǎng)衰老方程、病害發(fā)展動(dòng)態(tài)模型[3-7]、混合正態(tài)分布方程、復(fù)雜房室模型[8-11]、植物光合光響應(yīng)方程[12-14]、奶牛泌乳曲線機(jī)理方程[15]、土壤水分模型[16-18]、農(nóng)藥殘留降解模型[19-21]等。對(duì)于純非線性方程最優(yōu)擬合難度極大，傳統(tǒng)方法有最速下降法、牛頓法、共軛梯度法，近代的有麥夸法、模擬退火法、遺傳算法、縮張算法、人工神經(jīng)網(wǎng)絡(luò)、蟻群算法等。傳統(tǒng)最優(yōu)擬合方法存在對(duì)初值依賴性強(qiáng)、易陷于局優(yōu)陷阱[1-2]的缺點(diǎn)，麥夸法對(duì)簡(jiǎn)單非線性模型初值要求較低[22]，但對(duì)于復(fù)雜非線性方程（如Richard方程）極難擬合，對(duì)于房室模型中一室模型相對(duì)容易，對(duì)二室模型則會(huì)得出不恰當(dāng)?shù)慕鈁23]。其他算法多是模擬生物自然演化過(guò)程的進(jìn)化類算法。進(jìn)化算法屬于計(jì)算智能，現(xiàn)有進(jìn)化算法以遺傳算法為代表，但遺傳算法存在早熟收斂的問(wèn)題，至今未能很好克服[24]?？s張算法是一種優(yōu)秀的最優(yōu)擬合算法，不依賴導(dǎo)數(shù)支持，在已報(bào)道實(shí)例中均能達(dá)到全局最優(yōu)，但對(duì)于參數(shù)過(guò)多，如參數(shù)20個(gè)以上，則按最少試算點(diǎn)計(jì)算，每次試算量達(dá)320個(gè)試算點(diǎn)，計(jì)算難以實(shí)現(xiàn)。故應(yīng)當(dāng)進(jìn)一步研究開(kāi)發(fā)新的算法，以克服遺傳算法與縮張算法的固有缺陷。

1?材料與方法

1.1?數(shù)據(jù)來(lái)源

1.1.1?Richard方程擬合

美國(guó)國(guó)家標(biāo)準(zhǔn)網(wǎng)站（NIST）提供了Richard方程的最優(yōu)擬合計(jì)算案例，殘差平方和（http：//www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml）由15對(duì)數(shù)據(jù)組成，最小離回歸平方和Qe=8 786.404 908 0。

yi=A（1+Be-kti）1/N。（1）

式中：yi為第i個(gè)觀察值;A為生長(zhǎng)極限值;B為初始生長(zhǎng)量參數(shù);k為生長(zhǎng)速度參數(shù);ti為第i個(gè)時(shí)刻;N為曲線形狀參數(shù)。

Qe=∑n-1i=0（yi-yei）2。（2）

式中：Qe為殘差平方和;yi為第i個(gè)觀察值;yei為第i個(gè)預(yù)測(cè)值。

1.1.2?二室模型

二室模型共5個(gè)參數(shù)，且滿足Km>K2>K4。摘自袁志發(fā)《多元統(tǒng)計(jì)分析》中二室模型擬合實(shí)例演算[22]，文獻(xiàn)[25-26]以麥夸法求解，初值以退層法求得。

二室模型：A2（e-k2t-e-kmt）+A4（e-k4t-e-kmt）。（3）

1.1.3?奶牛泌乳方程Dijkstra與Wood方程[25-26]

分別用下面2個(gè)公式表示。

y=Aeb（1-e-ct）/c-dt;（4）

y=Atbe-ct。（5）

式中：各字母含義見(jiàn)文獻(xiàn)[25-26]。

以2種數(shù)學(xué)模型分別對(duì)羅清堯、能本海等《中國(guó)荷斯坦奶牛第二泌乳期泌乳曲線模型的研究》《中國(guó)荷斯坦奶牛第三泌乳期泌乳曲線模型的研究》文中第2、第3胎的各7組數(shù)據(jù)重算。

1.2?輪回選擇算法

輪回選擇是經(jīng)典的育種方法，工作重點(diǎn)是選種、鑒定與自由互交。輪回選擇在有限的遺傳背景中，能最大程度地淘汰不良基因，聚合優(yōu)良基因，以產(chǎn)生更高產(chǎn)量、抗性與品質(zhì)的新種質(zhì)。在輪回選擇時(shí)，每輪選擇大量選種材料，鑒定出10%的優(yōu)良選種材料，作為自由互交親本，親本間自由雜交產(chǎn)生下輪選種材料。如此周而復(fù)始，不斷提高親本遺傳表現(xiàn)水平，多代選擇后必會(huì)獲得一群高水平的個(gè)體。

筆者受輪回選擇育種程序啟發(fā)，將殘差平方和作為目標(biāo)函數(shù)值，看作遺傳表現(xiàn)值，約束條件不滿足看作致死基因（只選擇存活個(gè)體），各自變量看作染色體，染色體以實(shí)數(shù)編碼，先按各參數(shù)給定區(qū)間隨機(jī)生成染色體組，構(gòu)成若干親本，親本間進(jìn)行自由雜交、自由突變產(chǎn)生一定量存活選種群體，對(duì)存活選種群體依目標(biāo)函數(shù)值排序，選擇若干優(yōu)良選種材料作下輪親本，親本間再自由雜交、自由突變產(chǎn)生下代選種群體，如此周而復(fù)始，直到親本間沒(méi)有遺傳變異與或目標(biāo)函數(shù)值收斂到迭代精度要求為止。經(jīng)歷多個(gè)連續(xù)世代進(jìn)化，收斂達(dá)到精度要求，稱為1輪輪回選擇進(jìn)化，此進(jìn)化算法稱為輪回選擇算法。

當(dāng)對(duì)參數(shù)取值區(qū)間完全不知情，可以設(shè)定1個(gè)很寬的區(qū)間，首輪（第1輪）進(jìn)化結(jié)束，以最優(yōu)個(gè)體自變量（染色體）取值為中心，以其絕對(duì)值為鄰域半徑，進(jìn)行第2輪輪回選擇進(jìn)化計(jì)算，之后每輪輪回選擇進(jìn)化計(jì)算，均以最優(yōu)個(gè)體自變量（染色體值）取值為中心，將區(qū)間縮小到0.618倍，重新生成隨機(jī)親本，并加入上輪最優(yōu)個(gè)體，重新作輪回選擇進(jìn)化計(jì)算，一般進(jìn)行10輪計(jì)算方能達(dá)到全局最優(yōu)解。如上組織多輪輪回選擇進(jìn)化計(jì)算，放寬初始區(qū)間要求及提高計(jì)算精度的算法稱為多重輪回選擇算法。

1.2.1?基本輪回選擇算法

輪回選擇算法基于親本自由雜交及自由突變產(chǎn)生選種群體，再?gòu)倪x種群體中優(yōu)選存活個(gè)體組成新一輪親本群體，每一個(gè)由初始親本群體形成下輪親本群體的進(jìn)化過(guò)程稱為1個(gè)進(jìn)化世代。不同進(jìn)化世代間，只要上代親本選種群體存在差異且有優(yōu)劣區(qū)別，進(jìn)化就會(huì)發(fā)生，反之進(jìn)化則會(huì)停止。每輪進(jìn)化由下列算子構(gòu)成：初始親本生成算子、雜交算子、突變算子、選擇算子、迭代終止算子。

不斷地應(yīng)用雜交算子、突變算子、選擇算子、迭代終止算子進(jìn)行世代更替，當(dāng)世代間進(jìn)化達(dá)到收斂要求，則完成了一輪輪回選擇計(jì)算。計(jì)算實(shí)踐表明，1輪輪回選擇計(jì)算結(jié)果一般不差于遺傳算法，一般表現(xiàn)更優(yōu)。

針對(duì)本算法，以C++編程語(yǔ)言，以面向?qū)ο蟮木幊谭椒ㄔO(shè)計(jì)RSBase類，將常用參數(shù)設(shè)計(jì)成類的屬性，各算子設(shè)計(jì)成類的私有方法，參數(shù)輸入與輸出方法設(shè)計(jì)為公有方法，以利于使用者改寫(xiě)與調(diào)用。類中將Fx（目標(biāo)函數(shù)定義）、Gx（約束函數(shù)定義）設(shè)定為純虛函數(shù)，每次調(diào)用必須根據(jù)使用者特有的統(tǒng)計(jì)模型改寫(xiě)，同樣，參數(shù)輸出方法Para也設(shè)計(jì)成了虛函數(shù)，方便改寫(xiě)。另外設(shè)計(jì)了函數(shù)求極值及積分的公有方法，利于對(duì)使用者定義的復(fù)雜數(shù)學(xué)模型求函數(shù)極值點(diǎn)、極值及區(qū)間積分。軟件著作權(quán)《連農(nóng)最優(yōu)生長(zhǎng)擬合類庫(kù)軟件V1.0》（登記號(hào)：2147408）包含了本算法全部C++代碼。算法流程詳見(jiàn)圖1。

1.2.2?多重輪回選擇算法

一輪輪回選擇計(jì)算多由于區(qū)間不恰當(dāng)或區(qū)間過(guò)大，一般不能達(dá)到全局最優(yōu)解，如要求更高計(jì)算精度，可在上輪計(jì)算結(jié)果基礎(chǔ)上精算。每輪初始親本生成邊界值由空間收斂算子給出，算子算法描述如下：

tmp=0.618repeat

repeatrepeat+1

if（tmp<0.001）tmp0.001

xa=optX-tmp×|optX|

xb=optX+tmp×|optX|。

對(duì)每個(gè)參數(shù)獨(dú)立計(jì)算取值區(qū)間，由xa、xb計(jì)算產(chǎn)生初始親本群體及表征突變程度上限。收斂倍率0.618是筆者經(jīng)過(guò)數(shù)百次試算由經(jīng)驗(yàn)給出，一般取值區(qū)間（0.5，1）。repeat初始值為0，optX為參數(shù)上輪最優(yōu)估計(jì)值，tmp最小值暫定為0.001，最佳值有待深入研究確定。

每輪最優(yōu)親本由上輪最優(yōu)選種個(gè)體遺傳，能保障輪次間收斂與持續(xù)進(jìn)化，參數(shù)取值區(qū)間的縮小，能更大程度地減少早熟收斂的可能性。多重輪回選擇算法流程見(jiàn)圖2。

1.2.3?輪回選擇算法屬性構(gòu)造與各算子設(shè)計(jì)

1.2.3.1?屬性構(gòu)造

本算法以C++語(yǔ)言設(shè)計(jì)，程序執(zhí)行前，要求定義一些基本的參數(shù)，以利程序安排必要的內(nèi)存空間用于計(jì)算。重要參數(shù)如下：

Nx為參數(shù)變量個(gè)數(shù);Repeat為多重輪回選擇計(jì)算輪數(shù)計(jì)數(shù)器;nP、nS為親本與選種群體規(guī)模;iG、iS、maxG為進(jìn)化世代計(jì)數(shù)、選種個(gè)體計(jì)數(shù)及最大進(jìn)化代數(shù);isint為整型變量標(biāo)志指針，規(guī)定變量取實(shí)數(shù)還是整數(shù);xa、xb參數(shù)取值區(qū)間指針，初始群體由之計(jì)算產(chǎn)生，突變程度也由之表征;optX、optValue最優(yōu)個(gè)體染色體取值指針及其目標(biāo)函數(shù)值，作為算法輸出;Parent、valueP親本群體染色體取值及對(duì)應(yīng)目標(biāo)函數(shù)值指針;Select、valueS選種群體染色體取值及對(duì)應(yīng)目標(biāo)函數(shù)值指針。

1.2.3.2?初始群體構(gòu)造算子

第i個(gè)變量依據(jù)（xa，i，xb，i）區(qū)間隨機(jī)生成親本參數(shù)，如整型標(biāo)志變量isinti=0將相應(yīng)變量取整，每個(gè)變量的區(qū)間、整型約束獨(dú)立處理。初始群體由存活個(gè)體構(gòu)成。初始群體的構(gòu)造，不像其他算法那樣過(guò)度依賴初始值，給研究者一個(gè)寬松的條件，之后在第1輪最優(yōu)結(jié)果獲得后，再據(jù)其修正取值區(qū)間，區(qū)間長(zhǎng)為最優(yōu)個(gè)體相應(yīng)自變量取值絕對(duì)值的2倍，以防初始區(qū)間不恰當(dāng)造成的隨機(jī)生成樣本分布不均勻，之后每輪區(qū)間長(zhǎng)度縮小到上輪0.618倍，當(dāng)區(qū)間長(zhǎng)達(dá)到原區(qū)間長(zhǎng)的0.001倍時(shí)不再縮小。每輪區(qū)間收斂能保證最優(yōu)解區(qū)間的收斂，減少進(jìn)化難度，區(qū)間愈小，進(jìn)化愈快。

1.2.3.3?雜交算子

每條染色體獨(dú)立雜交，每次雜交隨機(jī)選擇2個(gè)不同的親本，以[0，1]間隨機(jī)數(shù)α為因子一次雜交形成2個(gè)選種個(gè)體。有整型約束的變量變換為整型。雜交如同輪回選擇中的自由互交，個(gè)體都由自由交配形成，保障了遺傳多樣性。雜交采用實(shí)數(shù)形式表示，由式（6）、式（7）計(jì)算：

x1=αPi+（1-α）Pj;（6）

x2=（1-α）Pi+αPj。（7）

1.2.3.4?誘變算子?因?yàn)檎T變產(chǎn)生變異及雜交產(chǎn)生變異是本算法的2種變異方式，誘變概率過(guò)低會(huì)導(dǎo)致雜交變異產(chǎn)生的變異個(gè)體顯著多于誘變產(chǎn)生的變異個(gè)體，為相對(duì)均衡起見(jiàn)，不能將誘變概率定太低。αi∈[0，1]間隨機(jī)數(shù)，用于計(jì)算第i個(gè)變量突變程度，βi∈[0，1]間隨機(jī)數(shù)，表示第i個(gè)自變量是否突變，表征誘變發(fā)生的概率。當(dāng)βi<0.1時(shí)進(jìn)行突變處理，否則不作突變處理，整體突變概率很高。計(jì)算實(shí)踐表明，突變?cè)谶M(jìn)化過(guò)程中起重要作用。如果αi<0.5，采用式（8）。

x=P-αi[（xi）max-xai]。（8）

式中：（xi）max為當(dāng)前世代親本群體第i個(gè)自變量的最大值。如果αi≥0.5，采用式（9）。

x=P+αi[xbi-（xi）max]。（9）

式中：a、b表示屬性值的左、右邊界值。

1.2.3.5?選擇算子

選擇分為2個(gè)部分。一部分選擇是致死選擇，即約束函數(shù)Gx=0，為致死個(gè)體，不能參與進(jìn)化計(jì)算;另一部分選擇對(duì)存活個(gè)體選優(yōu)汰劣，由選種群體的排序功能完成，優(yōu)選nP個(gè)個(gè)體加入下輪親本群體。選擇算子優(yōu)選出的解總是可行解，這與運(yùn)籌學(xué)中的外點(diǎn)法不同，同時(shí)又避免了內(nèi)點(diǎn)法初始可行點(diǎn)難以選擇的矛盾。迭代過(guò)程中不會(huì)因?yàn)楹瘮?shù)性態(tài)差而產(chǎn)生大幅度偏差。

1.2.3.6?空間收斂算子

當(dāng)解空間過(guò)大時(shí)，會(huì)造成進(jìn)化計(jì)算獲取的采樣點(diǎn)不一定能覆蓋最優(yōu)解空間，容易造成早熟收斂，如對(duì)不包含解空間的空間進(jìn)行割舍，能集中計(jì)算能力尋找最優(yōu)解。具體方法見(jiàn)“1.1.2”節(jié)。

1.2.3.7?迭代終止算子

迭代終止算子是判斷計(jì)算能否終止的方法，當(dāng)進(jìn)化代數(shù)達(dá)到上限時(shí)，終止計(jì)算;或親本沒(méi)有遺傳變異時(shí)，終止計(jì)算;或親本群體目標(biāo)函數(shù)值變異標(biāo)準(zhǔn)差小于指定精度值，提前結(jié)束本輪計(jì)算。當(dāng)群體變異存在，但目標(biāo)函數(shù)值已經(jīng)基本沒(méi)有變化時(shí)，應(yīng)提前終止;如有整數(shù)約束時(shí)，常存在多個(gè)最優(yōu)解。

1.3?經(jīng)典算法介紹

高斯-牛頓法應(yīng)用各一階導(dǎo)數(shù)構(gòu)造一個(gè)正交矩陣，建立一個(gè)最優(yōu)搜索方向，當(dāng)目標(biāo)函數(shù)全局性態(tài)單一、有類似正定二次型的曲面形態(tài)且初始搜索點(diǎn)處對(duì)應(yīng)正交矩陣可逆，則計(jì)算過(guò)程很快，速度優(yōu)于所有現(xiàn)代算法。

麥夸法是在高斯-牛頓法基礎(chǔ)上發(fā)展起來(lái)的一種近代算法。當(dāng)初始搜索點(diǎn)處對(duì)應(yīng)正交矩陣不可逆，通過(guò)將正交矩陣主對(duì)角線元素加上一足夠小正數(shù)，可以改變矩陣性態(tài)，可以適當(dāng)放寬高斯-牛頓法的初始值條件，計(jì)算量相對(duì)其成倍提高。

牛頓法以初始點(diǎn)處黑塞矩陣計(jì)算優(yōu)化點(diǎn)，當(dāng)黑塞矩陣可逆時(shí)，迭代計(jì)算非常快，如目標(biāo)函數(shù)為正定二次型可一次迭代成功，當(dāng)黑塞矩陣不可逆時(shí)，類似于麥夸法，將黑塞矩陣主對(duì)角線元素加上一足夠小正數(shù)，此法稱為修正牛頓法。牛頓法與修正牛頓法要計(jì)算各種二階導(dǎo)數(shù)，計(jì)算復(fù)雜，不利于編程計(jì)算。二者與前述兩法相比，存在同樣的問(wèn)題，當(dāng)目標(biāo)函數(shù)性態(tài)很差時(shí)，存在相當(dāng)多的局部最優(yōu)，則計(jì)算很難達(dá)到全局最優(yōu)，有的甚至沒(méi)有可行解。

縮張算法是由揚(yáng)州大學(xué)顧世梁等提出的最優(yōu)化算法，通過(guò)擴(kuò)張步與收縮步組成一輪循環(huán)，初始值由隨機(jī)數(shù)生成，經(jīng)過(guò)多輪的縮張循環(huán)，可以較快地計(jì)算得到全局最優(yōu)解，在低維問(wèn)題處理上本算法目前是最優(yōu)方法，但面對(duì)高維問(wèn)題時(shí)，其試算點(diǎn)數(shù)過(guò)多，計(jì)算時(shí)間相當(dāng)長(zhǎng)，計(jì)算難以實(shí)現(xiàn)。

1.4?各算法間比較指標(biāo)

1.4.1?誤差方差

各算法比較最關(guān)鍵指標(biāo)為殘差平方和，好的算法獲得估計(jì)參數(shù)更加適合所研究模型，應(yīng)變量變異平方和中誤差平方和占比最小化。均方誤差（mean squared error，MSE）是衡量“平均誤差”的一種較方便的方法，表征實(shí)際值與預(yù)測(cè)值的差異程度，可由殘差平方和計(jì)算得到。均方根誤差（root mean squared error，RMSE）是均方誤差的算術(shù)平方根，也稱作標(biāo)準(zhǔn)誤差（standard error）。

1.4.2?AIC指標(biāo)

AIC是赤池弘次提出的衡量統(tǒng)計(jì)模型擬合優(yōu)良性的一種標(biāo)準(zhǔn)，可以權(quán)衡所估計(jì)模型的復(fù)雜度和此模型擬合數(shù)據(jù)的優(yōu)良性。

AIC=2k+nln（Qe/n）。（10）

式中：k是參數(shù)的數(shù)量;n為觀察數(shù);Qe為殘差平方和。

1.4.3?計(jì)算時(shí)間

計(jì)算時(shí)間是衡量算法效率的重要指標(biāo)，但算法有效性往往比算法效率更重要，如牛頓法是非線性最優(yōu)化問(wèn)題的最快算法，但計(jì)算往往不能得到正確結(jié)果，現(xiàn)代遺傳算法計(jì)算時(shí)間遠(yuǎn)多于牛頓法，但獲得了廣泛的應(yīng)用?？s張算法計(jì)算時(shí)間也遠(yuǎn)多于牛頓法，但它是處理低維非線性問(wèn)題的最佳算法。可見(jiàn)，計(jì)算時(shí)間在算法間比較意義相對(duì)較小，更重要的是比較計(jì)算可行性、準(zhǔn)確性與魯棒性。

1.4.4?對(duì)初值的依賴性?好的算法不僅要有精確的計(jì)算結(jié)果，還不能依賴初值的精確性，如牛頓法對(duì)初值要求較高，只有好的初值才會(huì)收斂到精確結(jié)果，麥夸法相對(duì)于牛頓法則對(duì)初值依賴性有所下降，但不當(dāng)?shù)某踔颠€會(huì)造成不收斂[22-23]。縮張算法對(duì)初值基本沒(méi)有要求，本研究算法對(duì)初值也沒(méi)有要求，對(duì)初值沒(méi)有依賴的算法稱為沒(méi)有初值依賴性。

2?結(jié)果與分析

2.1?Richard方程最優(yōu)擬合結(jié)果

應(yīng)用輪回選擇算法經(jīng)歷11輪計(jì)算，第0輪經(jīng)歷10代進(jìn)化，第1輪經(jīng)歷17代進(jìn)化，…，第10輪經(jīng)歷7代進(jìn)化（表1），殘差平方和Qe達(dá)到本算法最優(yōu)值8 786.404 913 187 29，與NIST網(wǎng)站公布的全球最優(yōu)結(jié)果8 786.404 908 0前8位有效數(shù)字完全相同，AIC分別為49.515 799 32、49.515 799 31。經(jīng)多次計(jì)算均達(dá)到相近計(jì)算結(jié)果（前9位有效數(shù)字一直相同），本次計(jì)算第6輪已達(dá)到計(jì)算精度要求，之后RMSE變化很小，且參數(shù)估計(jì)值變化也相當(dāng)小，有的輪次間不能產(chǎn)生進(jìn)化。本次嘗試隨機(jī)生成初值再作麥夸法求解，當(dāng)隨機(jī)點(diǎn)數(shù)達(dá)5 000個(gè)時(shí)，其最優(yōu)結(jié)果仍不及本算法。本算法屬性值構(gòu)造：親本群體規(guī)模100，每代雜交規(guī)模10 000，誘變規(guī)模10 000，迭代收斂精度要求1.0×10-10。A初始區(qū)間[0，1 000]，B初始區(qū)間[0，1 000]，K初始區(qū)間[0，10]，N初始區(qū)間[0，10]。建議10輪為佳，如提前進(jìn)入全局最優(yōu)值，之后輪次的進(jìn)化世代數(shù)很少，進(jìn)化很快。平均每個(gè)數(shù)據(jù)集單個(gè)模型計(jì)算用時(shí)不超過(guò)120 s。

擬合殘差總體為可接受水平，擬合決定系數(shù)0.991 8，相關(guān)系數(shù)0.995 9。殘差分析見(jiàn)表2。

2.2?房室模型中二室模型擬合

二室模型共5個(gè)參數(shù)，且滿足要求Km>K2>K4。對(duì)于有約束的擬合，傳統(tǒng)非線性擬合難度極大。給定初值不當(dāng)，不能產(chǎn)生符合約束的解，而初值選擇難度又很大。輪回選擇算法可將約束作為致死因子，在存活個(gè)體中優(yōu)選親本雜交、突變，產(chǎn)生新的選種群體，遴選最優(yōu)個(gè)體，求解最優(yōu)擬合參數(shù)。以式（3）模型與數(shù)據(jù)試算（表3）。本算法屬性值構(gòu)造：親本群體規(guī)模100，每代雜交規(guī)模10 000，誘變規(guī)模10 000，迭代收斂精度要求1.0×10-10。初始區(qū)間A2∈[0，5 000]、K2∈[0，10]、Km∈[0，50]、A4∈[0，50]、K4∈[0，50]，經(jīng)11輪計(jì)算，得 Qe=6 731.829 722 431 19（AIC=47.548 071 21），明顯低于文獻(xiàn)[22]中的最優(yōu)結(jié)果6 761.911 38（AIC=47.575 180 15），且完全符合約束條件。

以上擬合數(shù)據(jù)說(shuō)明，K4=0，不適合二室模型擬合，更適合一室模型擬合，但原文中K4明顯不為0，不能說(shuō)明一室模型不合適;A4與原文差異較大，其他3個(gè)參數(shù)差異較小。本算法計(jì)算更準(zhǔn)確。同Richard模型計(jì)算相似，第6輪基本不再有大的進(jìn)化。最優(yōu)擬合殘差見(jiàn)表4，決定系數(shù)0.996 4，相關(guān)系數(shù)0.998 2。

2.3?最優(yōu)數(shù)學(xué)模型選擇

奶牛泌乳方程，國(guó)內(nèi)外已報(bào)道有100多種，以“1.1.3”節(jié)模型與數(shù)據(jù)計(jì)算。

第2胎泌乳數(shù)據(jù)擬合殘差平方和（表5）作模型、組別兩因素?zé)o重復(fù)方差分析，表明模型間差異顯著（F值9.083 321 499，無(wú)差別顯著概率0.023 583 05），Dijkstra顯著優(yōu)于Wood模型。第3胎泌乳數(shù)據(jù)擬合殘差平方和（表5）作模型、組別兩因素?zé)o重復(fù)方差分析，也表明模型間差異顯著（F值11.093 4，無(wú)差別顯著概率0.015 793），Dijkstra也顯著優(yōu)于Wood模型。這與文獻(xiàn)[25-26]結(jié)論完全相反，對(duì)比表7數(shù)據(jù)發(fā)現(xiàn)，Wood模型參數(shù)文獻(xiàn)[25-26]SAS計(jì)算結(jié)果與本研究輪回選擇計(jì)算結(jié)果相似，但本算法所有殘差平方和均小于SAS結(jié)果，方差分析表明，2種方法間差異F值4.282 694 075，無(wú)差別概率0.058 986 506，接近顯著水平，2種算法殘差平方和差異越大，則參數(shù)估計(jì)值差異越大，表明輪回選擇算法在處理簡(jiǎn)單非線性模型方面優(yōu)于SAS。表6顯示，Dijkstra模型數(shù)據(jù)差異較大，少數(shù)組別相近，多數(shù)組別差異大，且擬合殘差14組數(shù)據(jù)均小于SAS結(jié)果，表明輪回選擇算法在處理復(fù)雜非線性模型方面也優(yōu)于SAS軟件系統(tǒng)的非線性分析模塊（一般采用高斯-牛頓法）。文獻(xiàn)[25-26]14組試驗(yàn)數(shù)據(jù)分別用2種模型計(jì)算，形成28組擬合數(shù)據(jù)，擬合殘差平方和Qe均表明本算法較SAS非線性擬合優(yōu)秀，但AIC結(jié)果2個(gè)模型間差別不明顯，可能AIC指標(biāo)結(jié)合了參數(shù)個(gè)數(shù)信息，實(shí)質(zhì)上泌乳數(shù)據(jù)需要更精確的預(yù)測(cè)結(jié)果，以Qe比較更為恰當(dāng)。

3?討論與結(jié)論

3.1?輪回選擇算法的優(yōu)點(diǎn)

3.1.1?魯棒性

輪回選擇算法采用親本自由雜交、自由突變產(chǎn)生選種群體，從選種群體中優(yōu)選個(gè)體組成新的親本群體，不斷地選優(yōu)汰劣，實(shí)現(xiàn)全局尋優(yōu)計(jì)算。本算法不采用導(dǎo)數(shù)計(jì)算，不要求目標(biāo)函數(shù)及約束函數(shù)連續(xù)可導(dǎo)，另可以針對(duì)整數(shù)約束求解，不受一般最優(yōu)化問(wèn)題的實(shí)數(shù)連續(xù)約束的制約。當(dāng)模型添加了區(qū)間限制（如非負(fù)），不等式約束限制（如房室模型中的二室、三室模型）等時(shí)，有特別大的靈活性。

3.1.2?適合解混合非線性模型與高難度生物學(xué)模型

如非線性模型中部分自變量存在整數(shù)約束，則稱為混合非線性模型，解此類問(wèn)題，傳統(tǒng)優(yōu)化算法無(wú)能為力。運(yùn)籌學(xué)中采用罰函數(shù)法中的外點(diǎn)法求解，初始值易找，但所得常為不可行解，內(nèi)點(diǎn)法初始點(diǎn)不易尋找，所得解均為可行解，但可行域不一定唯一，不能保證解為所有可行域中的最優(yōu)解，實(shí)踐應(yīng)用不便，如采用輪回選擇算法，易于實(shí)現(xiàn)，如將Richard方程中N要求取整數(shù)，求解中僅加整數(shù)標(biāo)志即可。

由于非線性模型擬合一般采用導(dǎo)數(shù)支持的傳統(tǒng)算法，當(dāng)模型復(fù)雜時(shí)，會(huì)在最優(yōu)解附近存在很多局部?jī)?yōu)化陷阱，擬合難度極大，如Richard方程，在美國(guó)NIST網(wǎng)站中將其歸為高難度數(shù)組。植物生長(zhǎng)過(guò)程不是一個(gè)單調(diào)生長(zhǎng)過(guò)程，還包括衰老過(guò)程，僅用生長(zhǎng)方程來(lái)擬合，不能解析衰老過(guò)程，衰老過(guò)程還干擾生長(zhǎng)過(guò)程的參數(shù)擬合。如設(shè)計(jì)生長(zhǎng)衰老方程，如式（11），將生長(zhǎng)與衰老過(guò)程相疊加，則可以解析生長(zhǎng)過(guò)程與衰老過(guò)程，還能找到二者平衡點(diǎn)，揚(yáng)州大學(xué)孫成明在博士論文中提及此想法[27]，但一直未有應(yīng)用報(bào)道，主要原因是參數(shù)擬合難度遠(yuǎn)高于Richard方程，筆者研究發(fā)現(xiàn)，輪回選擇算法可以高效準(zhǔn)確擬合此類方程，將另文報(bào)道。

W=A11+B1e-kt-A21+e-k（t-t0）;（11）

A2（e-k2t-e-kmt）;（12）

A2（e-k2t-e-kmt）N。（13）

式（12）方程用于描述生物對(duì)藥物的吸收清除過(guò)程，稱為一室模型，在應(yīng)用中一室模型擬合精度一般不高，但給方程括號(hào)項(xiàng)加一指數(shù)成分后變?yōu)樗膮?shù)方程，形如式（13），擬合精度顯著提高，也極大地提高了擬合難度，用輪回選擇算法也能準(zhǔn)確擬合，將另文報(bào)道。此類復(fù)雜的生物學(xué)方程還有很多，意義很大，難度很高，限制了在生物學(xué)科研中的應(yīng)用，輪回選擇算法的成功將拓展對(duì)復(fù)雜生物學(xué)模型的開(kāi)發(fā)與應(yīng)用，促進(jìn)生物學(xué)理論研究的深入。

3.1.3?每輪重新初始化，防止了解樣本分布不均勻、早熟收斂問(wèn)題

遺傳算法與模擬退火算法在本算法提出前早有應(yīng)用，但其早熟收斂缺陷不易克服，主要是進(jìn)化計(jì)算過(guò)程中對(duì)最優(yōu)參數(shù)空間的不自覺(jué)舍棄[22]。輪回選擇算法因采用了類似輪回選擇育種中自由互交的變異方式，單輪計(jì)算出現(xiàn)早熟收斂的概率較低，采用多輪優(yōu)化技術(shù)，每輪優(yōu)化均對(duì)最優(yōu)解空間作分割，區(qū)間長(zhǎng)縮小到0.618倍，當(dāng)早期輪次未計(jì)算出全局最優(yōu)解時(shí)，可在以后的輪次出現(xiàn)全局最優(yōu)解，在足夠的計(jì)算輪后總能找到全局最優(yōu)解。如提高計(jì)算精度（高精度計(jì)算支持）還可以突破現(xiàn)有最優(yōu)記錄。從表1和表3結(jié)果可知，隨著計(jì)算輪次的增加，進(jìn)化增益很小，一般沒(méi)必要超過(guò)10輪，完全可以處理一般科研應(yīng)用問(wèn)題的求解。

3.1.4?隱并行性

如同遺傳算法，本算法避開(kāi)了最優(yōu)方向及最優(yōu)步長(zhǎng)的搜索，也有計(jì)算的隱并行性，可實(shí)現(xiàn)多線程的并行計(jì)算，大幅提高計(jì)算時(shí)間效率。

3.2?算法優(yōu)化

輪回選擇算法計(jì)算，對(duì)不同的模型，最優(yōu)算法參數(shù)不同，采用相同的參數(shù)也完全可以計(jì)算出各種類型模型的全局最優(yōu)解，但計(jì)算效率有區(qū)別，針對(duì)不同的模型可以探討最優(yōu)參數(shù)配置問(wèn)題，如以全局最優(yōu)解解出計(jì)算時(shí)間為目標(biāo)函數(shù)值進(jìn)行復(fù)雜的非線性規(guī)劃，可以計(jì)算最優(yōu)參數(shù)配置問(wèn)題，此問(wèn)題有待深入研究。本算法有待改進(jìn)，促進(jìn)算法的完善與推廣。

參考文獻(xiàn)：

[1]顧世梁，惠大豐，莫惠棟. 非線性方程最優(yōu)擬合的縮張算法[J]. 作物學(xué)報(bào)，1998，24（5）：513-519.

[2]顧世梁，徐辰武，蒯建敏. 用縮張算法最優(yōu)擬合非線性方程的檢驗(yàn)[J]. 揚(yáng)州大學(xué)學(xué)報(bào)（自然科學(xué)版），2001，1（3）：16-19.

[3]劉維紅，林茂松，李紅梅，等. 人工接種測(cè)定水稻干尖線蟲(chóng)在水稻上的病害發(fā)展動(dòng)態(tài)[J]. 中國(guó)農(nóng)業(yè)科學(xué)，2007，40（12）：2734-2740.

[4]劉維紅. 水稻干尖線蟲(chóng)在“小穗頭”上的危害及人工接種測(cè)定病害發(fā)展動(dòng)態(tài)[D]. 南京：南京農(nóng)業(yè)大學(xué)，2007.

[5]彭懷俊，顧?鋼，紀(jì)成燦，等. 烤煙根系土壤中青枯病菌動(dòng)態(tài)與田間病害發(fā)生發(fā)展的關(guān)系[J]. 湖南農(nóng)業(yè)大學(xué)學(xué)報(bào)（自然科學(xué)版），2005，31（4）：384-387.

[6]曹?靜. 設(shè)施條件下馬鈴薯晚疫病流行主導(dǎo)因素及病害防治的初步研究[D]. 保定：河北農(nóng)業(yè)大學(xué)，2002.

[7]柏自琴. 中國(guó)柑橘黃龍病發(fā)生動(dòng)態(tài)及其病原菌亞洲種群分化研究[D]. 重慶：西南大學(xué)，2012.

[8]曾文藝，孫曉穎. 藥物動(dòng)力學(xué)房室模型的改進(jìn)[J]. 北京師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2012，48（1）：6-10.

[9]王曉麗，趙?晶. 傳染病房室模型的建立及求解方法[J]. 山東輕工業(yè)學(xué)院學(xué)報(bào)（自然科學(xué)版），2006，20（3）：88-90.

[10]何紹雄，曹鑒萍，黃團(tuán)華. 非線性房室模型參數(shù)計(jì)算法[J]. 數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，1986（3）：147-152.

[11]楊博文，劉?萍. 藥物動(dòng)力學(xué)二房室模型[J]. 哈爾濱師范大學(xué)（自然科學(xué)學(xué)報(bào)），2016，32（1）：13-15，66.

[12]韓?剛，趙?忠. 不同土壤水分下4種沙生灌木的光合光響應(yīng)特性[J]. 生態(tài)學(xué)報(bào)，2010，30（15）：4019-4026.

[13]閆小紅，尹建華，段世華，等. 四種水稻品種的光合光響應(yīng)曲線及其模型擬合[J]. 生態(tài)學(xué)雜志，2013，32（3）：604-610.

[14]王榮榮，夏江寶，楊吉華，等. 貝殼砂生境干旱脅迫下杠柳葉片光合光響應(yīng)模型比較[J]. 植物生態(tài)學(xué)報(bào)，2013，37（2）：111-121.

[15]賈先波，陳仕毅，王?杰，等. 中國(guó)荷斯坦牛泌乳曲線數(shù)學(xué)模型擬合及應(yīng)用[J]. 南京農(nóng)業(yè)大學(xué)學(xué)報(bào)，2016，39（1）：133-138.

[16]李洪文，高煥文. 保護(hù)性耕作土壤水分模型[J]. 中國(guó)農(nóng)業(yè)大學(xué)學(xué)報(bào)，1996，1（2）：25-30.

[17]于?浕. 土壤水分特征曲線Gardner模型參數(shù)的預(yù)報(bào)模型研究[D]. 太原：太原理工大學(xué)，2017.

[18]張?巖，朱?巖，張建軍， 等. 林地土壤水分模型SWUF在晉西黃土高原的適用性[J]. 林業(yè)科學(xué)，2012，48（5）：8-14.

[19]吳?鵬，秦智偉，周秀艷，等. 蔬菜農(nóng)藥殘留研究進(jìn)展[J]. 東北農(nóng)業(yè)大學(xué)學(xué)報(bào)，2011，42（1）：138-144.

[20]陳振德，陳雪輝，馮明祥，等. 毒死蜱在菠菜中的殘留動(dòng)態(tài)研究[J]. 農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào)，2005，24（4）：728-731.

[21]張韓杰，閆艷春. 農(nóng)藥殘留及微生物在農(nóng)藥降解中的應(yīng)用與展望[J]. 湖北植保，2004（1）：31-35.

[22]袁志發(fā)，周靜芋. 多元統(tǒng)計(jì)分析[M]. 北京：科學(xué)出版社，2002：247-248.

[23]周靜芋，宋世德，袁志發(fā)，等. 用麥夸法進(jìn)行室分析曲線擬合[J]. 西北農(nóng)業(yè)大學(xué)學(xué)報(bào)，1996，24（1）：75-78.

[24]范青武. 遺傳算法動(dòng)態(tài)演化過(guò)程的研究[D]. 北京：北京工業(yè)大學(xué)，2007.

[25]羅清堯，熊本海，馬?毅，等. 中國(guó)荷斯坦奶牛第二泌乳期泌乳曲線模型的研究[J]. 中國(guó)農(nóng)業(yè)科學(xué)，2010，43（23）：4910-4916.

[26]熊本海，馬?毅，羅清堯，等. 中國(guó)荷斯坦奶牛第三泌乳期泌乳曲線模型的研究[J]. 中國(guó)農(nóng)業(yè)科學(xué)，2011，44（2）：402-408.

[27]孫成明. FACE水稻生長(zhǎng)發(fā)育模擬模型研究[D]. 揚(yáng)州：揚(yáng)州大學(xué)，2006.