基于改進(jìn)隨機(jī)梯度下降的反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì)方法

唐立森, 陳偉鋒

基于改進(jìn)隨機(jī)梯度下降的反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì)方法

唐立森, 陳偉鋒

(浙江工業(yè)大學(xué) 信息工程學(xué)院, 浙江 杭州 310023)

針對(duì)傳統(tǒng)優(yōu)化方法利用所有采樣數(shù)據(jù)進(jìn)行參數(shù)估計(jì)存在的求解困難問題，在聯(lián)立求解的框架下，通過引入隨機(jī)優(yōu)化和擴(kuò)展目標(biāo)函數(shù)，提出基于改進(jìn)隨機(jī)梯度下降的反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì)方法。該方法對(duì)多數(shù)據(jù)集的大規(guī)模系統(tǒng)進(jìn)行機(jī)理建模，基于靈敏度微分方程法獲得靈敏度矩陣，同時(shí)利用模型標(biāo)度化技術(shù)處理多狀態(tài)變量對(duì)多參數(shù)估計(jì)的同步收斂性問題。為了減小迭代過程中噪聲方差的影響，在現(xiàn)有的隨機(jī)平均梯度下降方法的基礎(chǔ)上，利用隨機(jī)擴(kuò)展目標(biāo)函數(shù)增加目標(biāo)函數(shù)中計(jì)算梯度的信息量，并給出該方法收斂的理論性分析。數(shù)值仿真結(jié)果驗(yàn)證了該方法的有效性和可行性。

參數(shù)估計(jì)；隨機(jī)優(yōu)化；擴(kuò)展目標(biāo)；靈敏度矩陣

1 前言

由于開環(huán)式的控制策略遵循相似的軌跡模式，隨著間歇反應(yīng)過程[1]多批次的運(yùn)行，雖然數(shù)據(jù)量會(huì)不斷增加，但有效信息量并不隨批次的增加呈線性增長(zhǎng)，因此間歇反應(yīng)過程往往是在有限信息量的基礎(chǔ)上進(jìn)行建模和參數(shù)估計(jì)[2]。針對(duì)多批次數(shù)據(jù)，傳統(tǒng)參數(shù)估計(jì)方法如極大似然法(maximum likelihood estimation，MLE)[3-4]、最大期望法(expectation maximization，EM)[5]等聯(lián)立所有數(shù)據(jù)進(jìn)行參數(shù)估計(jì)會(huì)導(dǎo)致問題的規(guī)模十分龐大，并且隨著批次數(shù)據(jù)逐漸增加，采用傳統(tǒng)優(yōu)化方法求解基于多批次數(shù)據(jù)的參數(shù)估計(jì)問題會(huì)在求解能力和計(jì)算效率上存在一定的問題?；跀U(kuò)展卡爾曼濾波(extended kalman filter，EKF)[6]的估計(jì)方法是有效的，但在非線性很強(qiáng)的情況下參數(shù)估計(jì)是有偏的；馬爾可夫鏈蒙特卡爾方法(markov chain monte carlo，MCMC)[7]通過抽取樣本數(shù)值逼近概率密度函數(shù)進(jìn)行參數(shù)估計(jì)，但對(duì)于具有多狀態(tài)多參數(shù)模型，它的計(jì)算代價(jià)是非常昂貴的。針對(duì)此類問題通過引入隨機(jī)梯度下降算法(stochastic gradient descent，SGD)[8]對(duì)模型中的反應(yīng)動(dòng)力學(xué)參數(shù)進(jìn)行估計(jì)。SGD源于1951年Robbins和Monro[9]提出的隨機(jī)逼近，主要用于求解大規(guī)模系統(tǒng)優(yōu)化問題以及處理機(jī)器學(xué)習(xí)任務(wù)[10-14]。隨著對(duì)隨機(jī)優(yōu)化算法的深入研究，衍生了不同版本的變體算法[15-18]，為了充分利用歷史梯度信息，隨機(jī)平均梯度算法(stochastic average gradient，SAG)[19]和加速隨機(jī)平均梯度(stochastic average gradient accelerated，SAGA)等[20]通過新梯度替代舊梯度的方式，充分考慮歷史梯度信息的同時(shí)減少了梯度計(jì)算量。但以上隨機(jī)梯度下降算法大部分以黑箱優(yōu)化器的形式使用，針對(duì)多批次數(shù)據(jù)集的反應(yīng)機(jī)理模型，Bae等[21]提出廣義化的拉普拉斯近似極大似然估計(jì)(generalization of laplace approximation maximum likelihood estimation，gLAMLE)算法，基于擴(kuò)展目標(biāo)函數(shù)，利用多批次數(shù)據(jù)進(jìn)行迭代估計(jì)，并在迭代更新中引入學(xué)習(xí)率，減小由于參數(shù)變化過快導(dǎo)致的數(shù)值震蕩，但結(jié)果仍然無法保證被估參數(shù)的收斂性，甚至可能會(huì)減緩收斂速度。本研究引用擴(kuò)展目標(biāo)函數(shù)的思想，提出基于擴(kuò)展目標(biāo)函數(shù)的隨機(jī)梯度下降算法，在聯(lián)立求解的框架下進(jìn)行參數(shù)估計(jì)，減小了單次估計(jì)的計(jì)算量，提升單次估計(jì)速度的同時(shí)保證了被估參數(shù)的收斂性。

2 基于多批次數(shù)據(jù)的參數(shù)估計(jì)問題

2.1 問題描述

假定間歇反應(yīng)過程中每一次反應(yīng)操作都是相同的，本研究考慮實(shí)際工況中進(jìn)料時(shí)存在隨機(jī)擾動(dòng)的影響，導(dǎo)致構(gòu)建系統(tǒng)模型時(shí)狀態(tài)初值會(huì)受隨機(jī)變量的影響而發(fā)生一定的變化，由于狀態(tài)初值發(fā)生隨機(jī)性改變，導(dǎo)致每次工況產(chǎn)生的批次數(shù)據(jù)存在一定的差異性。隨著反應(yīng)操作的不斷進(jìn)行，會(huì)生成大規(guī)模多批次的數(shù)據(jù)量，基于多批次數(shù)據(jù)進(jìn)行反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì)會(huì)使得求解優(yōu)化問題的規(guī)模十分龐大。針對(duì)此類問題考慮如下微分代數(shù)模型：

式中：(n)()屬于Rnz，表示第批次數(shù)據(jù)集對(duì)應(yīng)的微分代數(shù)模型中的狀態(tài)向量，Rnz表示nz維實(shí)數(shù)集；p(n)()和m(n)()同時(shí)屬于Rny，分別為第批次數(shù)據(jù)集對(duì)應(yīng)的輸出預(yù)測(cè)向量和輸出測(cè)量向量，Rny表示ny維實(shí)數(shù)集；屬于Rny，為輸出測(cè)量噪聲向量，其中噪聲向量中每個(gè)元素服從均值為0、方差為2的正態(tài)分布；屬于Rnz，表示進(jìn)料向量上的隨機(jī)擾動(dòng)，其中擾動(dòng)向量中每個(gè)元素服從均值為0、方差為2的正態(tài)分布；屬于Rnp，表示模型中的參數(shù)向量，Rnp表示np維實(shí)數(shù)集，且: Rnz+np?Rnz，: Rnz+np?Rny屬于可微函數(shù)；表示時(shí)間，0表示初始時(shí)刻，0表示0時(shí)刻對(duì)應(yīng)的狀態(tài)初值，上標(biāo)=1,…,，為數(shù)據(jù)總批次大小。

2.2 聯(lián)立配置點(diǎn)法

針對(duì)公式(1)中的模型，采用聯(lián)立配置點(diǎn)法進(jìn)行求解，微分狀態(tài)變量使用有限元的多項(xiàng)式來近似，其中有限元[-1,]，=1,…,nfe，滿足0<1<…

式中：p(n,s)表示第批次數(shù)據(jù)中t時(shí)刻的預(yù)測(cè)向量。第批次數(shù)據(jù)中t時(shí)刻的測(cè)量向量可以表示為m(n,s)，其中=1,…,，假設(shè)有批次數(shù)據(jù)并且每批次含有個(gè)時(shí)間采樣點(diǎn)，利用多批次測(cè)量數(shù)據(jù)，可以得到以下優(yōu)化問題用于估計(jì)模型參數(shù)：

3 改進(jìn)隨機(jī)優(yōu)化算法

3.1 隨機(jī)擴(kuò)展目標(biāo)函數(shù)

采用多批次數(shù)據(jù)進(jìn)行參數(shù)估計(jì)時(shí)，若在噪聲擾動(dòng)下每次只基于單批次數(shù)據(jù)進(jìn)行參數(shù)估計(jì)，則參數(shù)估計(jì)的精度不夠理想。Bae等[21]利用gLAMLE算法進(jìn)行參數(shù)估計(jì)時(shí)，通過設(shè)定采樣率來決定m大小，將整個(gè)數(shù)據(jù)集放入緩沖區(qū)，遍歷每一批次數(shù)據(jù)時(shí)，再從余下所有批次數(shù)據(jù)中隨機(jī)選取m批次數(shù)據(jù)集構(gòu)成擴(kuò)展目標(biāo)函數(shù)達(dá)到抵消噪聲擾動(dòng)的作用，同時(shí)基于m+1批次數(shù)據(jù)的全部樣本對(duì)參數(shù)進(jìn)行優(yōu)化更新。本研究借助gLAMLE算法中擴(kuò)展目標(biāo)函數(shù)的思想，結(jié)合隨機(jī)梯度下降設(shè)計(jì)了近似隨機(jī)梯度下降算法。

對(duì)于每一批次數(shù)據(jù)設(shè)計(jì)如下目標(biāo)函數(shù)：

類似于gLAMLE算法，采用如下隨機(jī)擴(kuò)展目標(biāo)函數(shù)：

3.2 隨機(jī)優(yōu)化算法

隨機(jī)梯度下降算法是機(jī)器學(xué)習(xí)、深度學(xué)習(xí)領(lǐng)域中比較流行的優(yōu)化算法，在隨機(jī)平均梯度下降算法的基礎(chǔ)上，通過引入3.1中隨機(jī)擴(kuò)展目標(biāo)函數(shù)的概念，抵消部分由于隨機(jī)干擾引起的估計(jì)誤差，命名為隨機(jī)擴(kuò)展目標(biāo)平均梯度下降法(stochastic extended objective average gradient descent，SEOAG)。

假定模型預(yù)測(cè)函數(shù)為

隨機(jī)擴(kuò)展目標(biāo)函數(shù)為

SEOAG算法流程如下所示：

1) 初始化參數(shù)，修正梯度向量設(shè)為v，且令初始時(shí)刻的v=0，迭代停止精度e，第0次梯度信息設(shè)為0=0，參數(shù)初值0=(1,…,np)，迭代次數(shù)=0，迭代終止次數(shù)Number，學(xué)習(xí)率為；隨機(jī)擴(kuò)展批次為m；

2) 根據(jù)初值0得到的初始靈敏度信息和預(yù)測(cè)數(shù)據(jù)，將靈敏度數(shù)據(jù)和預(yù)測(cè)數(shù)據(jù)映射到[-1, 1]，得到標(biāo)度化矩陣、，利用和對(duì)模型進(jìn)行標(biāo)度化處理；

11)輸出參數(shù)。

SEOAG算法在每一次迭代更新項(xiàng)中不僅包含新梯度信息，并且包含舊梯度信息。每次迭代更新前將樣本數(shù)據(jù)進(jìn)行隨機(jī)打亂處理，進(jìn)行順序遍歷確保充分利用了所有數(shù)據(jù)，然后通過引入隨機(jī)擴(kuò)展目標(biāo)函數(shù)的概念，增加隨機(jī)數(shù)據(jù)量達(dá)到抵消隨機(jī)擾動(dòng)的效果，此算法在保留SAG算法優(yōu)點(diǎn)的同時(shí)，抵消了一部分噪聲擾動(dòng)的影響，減小了參數(shù)變化引起的數(shù)值震蕩。

3.3 靈敏度計(jì)算

3.2節(jié)中梯度信息的計(jì)算涉及狀態(tài)變量對(duì)參數(shù)的靈敏度信息，而式(1)中約束模型整體屬于常微分方程組，因此可以采用微分方程法求解待估計(jì)參數(shù)的靈敏度信息矩陣：

對(duì)式(1a)中等式兩邊的參數(shù)進(jìn)行求導(dǎo)得到

運(yùn)用聯(lián)立配置法將式(2)、(3)、(4)代入式(11)，配置點(diǎn)方程可以寫成：

聯(lián)立求解方程組(6b～6d)和(12a～12c),即可得到待估計(jì)參數(shù)的靈敏度信息矩陣。

3.4 標(biāo)度化處理

機(jī)器學(xué)習(xí)中不同樣本的特征數(shù)據(jù)存在量綱的差異，由于數(shù)據(jù)間的差別可能很大，會(huì)直接影響學(xué)習(xí)的結(jié)果，因此必須對(duì)輸入和輸出數(shù)據(jù)按照一定比例進(jìn)行縮放，使之落在特定的區(qū)域內(nèi)，便于進(jìn)行結(jié)果分析，其中主要有歸一化和標(biāo)準(zhǔn)化等預(yù)處理方式。而基于隨機(jī)優(yōu)化算法對(duì)反應(yīng)機(jī)理模型進(jìn)行參數(shù)估計(jì)時(shí)，不同狀態(tài)變量對(duì)不同參數(shù)也存在量綱的差異，由于這種差異會(huì)導(dǎo)致模型對(duì)不同參數(shù)的影響程度不同，基于此類問題，本研究借鑒標(biāo)準(zhǔn)化預(yù)處理思想對(duì)模型進(jìn)行標(biāo)度化處理。整個(gè)參數(shù)估計(jì)過程可以看成[p,1(n,1)…p,ny(n,1)…p,1(n,M)…p,ny(n,M)]T=′的形式，其中輸入量為靈敏度矩陣，輸出量為[p,1(n,1)…p,ny(n,1)…p,1(n,M)…p,ny(n,M)]T，基于初始時(shí)刻0的靈敏度矩陣和第一批次的測(cè)量數(shù)據(jù)得到標(biāo)度化矩陣，利用標(biāo)度化矩陣對(duì)模型的輸入變量和輸出變量同時(shí)進(jìn)行標(biāo)度化處理后將數(shù)據(jù)映射到[-1,1]，具體標(biāo)度過程如下所示：

將轉(zhuǎn)換后的tr和ptr(n,s)分別代入式(6)和式(12)中對(duì)模型進(jìn)行標(biāo)度化處理后轉(zhuǎn)換為

4 收斂性分析

假設(shè)式(8)中任意批次的目標(biāo)函數(shù)(n,s)都是可微的，并且對(duì)應(yīng)的每一個(gè)梯度都是Lipschitz連續(xù)的，即對(duì)于參數(shù)區(qū)間內(nèi)任意的1，2(屬于Rnp)，滿足不等式

式中：L為L(zhǎng)ipschitz常數(shù)，假設(shè)任意批次的目標(biāo)函數(shù)(n,s)都滿足強(qiáng)凸性，由于SEOAG目標(biāo)函數(shù)中隨機(jī)擴(kuò)展了m個(gè)(j,s)，其中?e(n)，1，因此e(n,s)等于m個(gè)強(qiáng)凸函數(shù)之和，即e(n,s)() =(v_1,s)()+(v_2,s)()+… +(v_Nm,s)()，其中v_1,…,v_m表示從1,…,中(除以外)隨機(jī)選取的m個(gè)數(shù)，根據(jù)不等式(14)可知對(duì)于任意的1、2滿足

根據(jù)式(15)并結(jié)合向量范數(shù)三角不等式得到

由式(16)得到

則e(n,s)()也是Lipschitz連續(xù)的。

文獻(xiàn)[23]中給出了強(qiáng)凸函數(shù)的充分必要條件：對(duì)于Rnp區(qū)間內(nèi)任意的1、2滿足以下不等式關(guān)系

則由(18)式得到

根據(jù)式(19)得到

因此根據(jù)強(qiáng)凸函數(shù)的充分必要條件[23]可知e(n,s)()是強(qiáng)凸函數(shù)。進(jìn)一步地，根據(jù)文獻(xiàn)[10]可知存在常數(shù)滿足

則可以得到

5 數(shù)值實(shí)驗(yàn)分析

5.1 案例一

對(duì)于細(xì)胞反應(yīng)操作過程，進(jìn)料操作的主要目的是最大化細(xì)胞生長(zhǎng)和產(chǎn)物形成的速率，從而使產(chǎn)物形成的總速率(生產(chǎn)率)或產(chǎn)物收獲率(選擇性)最大化。通過調(diào)節(jié)限制底物、誘導(dǎo)劑、前體或中間體的投料速率和選擇適當(dāng)?shù)某跏紬l件來實(shí)現(xiàn)。給出了線性變化進(jìn)料速率的基本細(xì)胞反應(yīng)模型：

如圖1～ 3所示為*=1.0′10-6以及Number為12 000時(shí)，選取不同學(xué)習(xí)率時(shí)不同算法對(duì)反應(yīng)動(dòng)力學(xué)參數(shù)=(xs，ps)的估計(jì)情況。參數(shù)的初值0=(8.0, 4.0)，圖1中，當(dāng)學(xué)習(xí)率取0.000 1時(shí)，由于學(xué)習(xí)率偏小導(dǎo)致算法迭代估計(jì)速度較慢，其中SEOAG算法明顯受學(xué)習(xí)率影響較小，參數(shù)估計(jì)值收斂至真值附近，而其余的算法均未收斂；各個(gè)算法單次參數(shù)估計(jì)的時(shí)間損耗分別為SGD：0.093 8 s，MBGD：0.103 2 s，SAG：0.101 5 s，SEOAG：0.115 4 s，表明相比其他算法，SEOAG在目標(biāo)函數(shù)中擴(kuò)展了m批次的數(shù)據(jù)量，導(dǎo)致SEOAG單次估計(jì)速度會(huì)偏慢一些；圖2中當(dāng)學(xué)習(xí)率取0.001時(shí)，由于學(xué)習(xí)率選取適中，各個(gè)算法參數(shù)均能迭代估計(jì)到真值附近，而SEOAG算法比其他算法參數(shù)估計(jì)的收斂速度更快；增大學(xué)習(xí)率會(huì)增加收斂速度，但同時(shí)會(huì)在參數(shù)迭代估計(jì)中產(chǎn)生數(shù)值震蕩，圖3中當(dāng)學(xué)習(xí)率取0.01時(shí)，學(xué)習(xí)率過大導(dǎo)致不同算法對(duì)參數(shù)的估計(jì)值表現(xiàn)出明顯的數(shù)值震蕩，而SEOAG算法相比其他算法，有明顯降低震蕩的效果；具體參數(shù)估計(jì)值如表1～ 3所示。

圖1 不同算法對(duì)Yxs，Yps迭代估計(jì)過程的比較(a=0.000 1)

圖2 不同算法對(duì)Yxs，Yps迭代估計(jì)過程的比較(a=0.001)

圖3 不同算法對(duì)Yxs，Yps迭代估計(jì)過程的比較(a=0.01)

為了防止參數(shù)估計(jì)值受隨機(jī)性的影響，表1～ 3中分別列出了選取不同學(xué)習(xí)率時(shí)，不同算法分別運(yùn)行10次，取10次終止時(shí)刻參數(shù)估計(jì)的平均值。表1中由于學(xué)習(xí)率取值過小導(dǎo)致SGD、MBGD、SAG算法在終止時(shí)刻的估計(jì)值偏差較大、平均相對(duì)誤差增大，而SEOAG算法的參數(shù)估計(jì)值偏差最小、估計(jì)精度最好；從表2和3中的數(shù)據(jù)可知，增大學(xué)習(xí)率時(shí)，相比其他3種算法，SEOAG算法的參數(shù)估計(jì)精度也屬于中上水平；從圖1～ 3以及表1～3中的數(shù)據(jù)可知，SEOAG算法可以減小算法對(duì)學(xué)習(xí)率選取的依賴性，當(dāng)學(xué)習(xí)率取值較大時(shí)，SEOAG在加快收斂速度的同時(shí)可以減緩數(shù)值震蕩、減小估計(jì)誤差；當(dāng)學(xué)習(xí)率取值過小時(shí)，SEOAG在保證參數(shù)估計(jì)精度的同時(shí)可以明顯增加收斂速度。

表1 不同算法對(duì)Yxs，Yps估計(jì)值的比較(a=0.000 1)

表2 不同算法對(duì)Yxs，Yps估計(jì)值的比較(a=0.001)

表3 不同算法對(duì)Yxs，Yps估計(jì)值的比較(a=0.01)

5.2 案例二

考慮在一個(gè)體積固定的容器內(nèi)發(fā)生的兩步化學(xué)反應(yīng)過程：

式中：物料A和物料B反應(yīng)生成物料C，物料C反應(yīng)生成物料D，1、2表示反應(yīng)比率，A的初始濃度A(0)=1.5 mol×L-1，的初始濃度B(0)=1.0mol×L-1。假設(shè)該反應(yīng)中物料C和物料D的濃度是可測(cè)量的，則輸出量C=C+C；D=D+D，其中C、D都服從均值為0、方差為0.012的正態(tài)分布；假設(shè)A物料濃度初值存在隨機(jī)擾動(dòng)，服從均值為0、方差為0.012的正態(tài)分布。反應(yīng)比率真值=(1,2)=(0.5, 2)，假設(shè)該反應(yīng)在實(shí)際生產(chǎn)過程中是批量反應(yīng)的。基于多批次數(shù)據(jù)下運(yùn)用聯(lián)立求解框架重復(fù)案例一中求解的步驟，比較不同隨機(jī)優(yōu)化算法對(duì)反應(yīng)動(dòng)力學(xué)參數(shù)1、2估計(jì)結(jié)果。

圖4給出了各個(gè)物料反應(yīng)濃度B的模擬數(shù)據(jù)，反應(yīng)時(shí)間為10 s，通過在A物料濃度初始時(shí)刻增加隨機(jī)擾動(dòng)以及輸出量添加給定噪聲擾動(dòng)，模擬出=20的批次數(shù)據(jù)，*=1.0′10-6，Number=12 000，采樣率設(shè)為0.2，則m=′0.2=4；參數(shù)的初值0=(0.05, 3.5)，如圖5～7所示，當(dāng)選擇不同學(xué)習(xí)率時(shí)，不同算法對(duì)反應(yīng)動(dòng)力學(xué)參數(shù)1、2估計(jì)效果圖。為了防止參數(shù)估計(jì)值受隨機(jī)性的影響，表4～ 6列出了選擇不同學(xué)習(xí)率時(shí)，分別運(yùn)行10次，對(duì)參數(shù)1、2估計(jì)值取平均值。

圖4 k1、k2為真值時(shí)各物質(zhì)濃度的曲線

圖5 不同算法對(duì)k1，k2迭代估計(jì)過程的比較(a=0.001)

圖6 不同算法對(duì)k1，k2迭代估計(jì)過程的比較(a=0.01)

圖7 不同算法對(duì)k1，k2迭代估計(jì)過程的比較(a=0.1)

從圖5中可知學(xué)習(xí)率取0.001時(shí)，學(xué)習(xí)率取值過小導(dǎo)致各個(gè)算法對(duì)參數(shù)1、2的估計(jì)偏差較大，其中MBGD的樣本批次大小為4，SEOAG算法在相同條件下，參數(shù)估計(jì)值的收斂效果更好；圖5中各算法單次參數(shù)估計(jì)時(shí)間損耗分別為SGD：0.120 6 s，MBGD：0.128 0 s，SAG：0.134 7 s，SEOAG：0.139 4 s，表明相比其他算法，由于SEOAG在目標(biāo)函數(shù)中擴(kuò)展了m批次的數(shù)據(jù)量，導(dǎo)致SEOAG單次估計(jì)速度會(huì)偏慢一些；圖6表示學(xué)習(xí)率取0.01時(shí)，各個(gè)算法收斂速度加快，而SEOAG算法則明顯收斂速度更快；圖7表示學(xué)習(xí)率取0.1時(shí)，各個(gè)算法最后均能估計(jì)到真值附近、但由于學(xué)習(xí)率過大會(huì)導(dǎo)致明顯的數(shù)值震蕩，而SEOAG算法相比其他算法，有明顯降低震蕩的效果。表4～ 6中列出了選取不同的學(xué)習(xí)率時(shí)，各個(gè)算法分別運(yùn)行10次后，取參數(shù)估計(jì)的平均值以及估計(jì)值的平均相對(duì)誤差，表4和5中學(xué)習(xí)率取值過小，導(dǎo)致這4種算法在終止時(shí)刻參數(shù)估計(jì)值偏差都比較大，而基于表中數(shù)據(jù)可以發(fā)現(xiàn)SEOAG算法中估計(jì)值的平均相對(duì)誤差最?。划?dāng)學(xué)習(xí)率增大時(shí)，基于表6中的數(shù)據(jù)可知SEOAG算法在估計(jì)精度上優(yōu)勢(shì)不太明顯，但在取10次估計(jì)值的平均相對(duì)誤差也是最小的?；谏鲜鰣D5～7以及表4～6中的數(shù)據(jù)可知SEOAG算法可以減小算法對(duì)學(xué)習(xí)率取值的依賴性，在學(xué)習(xí)率取值過小時(shí)，SEOAG算法相比其他3種算法，收斂速度更快、更容易收斂到真值附近；在學(xué)習(xí)率取值較大時(shí)，可以減緩數(shù)值震蕩、減小估計(jì)值的平均相對(duì)誤差。

表4 不同算法對(duì)k1，k2估計(jì)值的比較(a=0.001)

表5 不同算法對(duì)k1，k2估計(jì)值的比較(a=0.01)

表6 不同算法對(duì)k1，k2估計(jì)值的比較(a=0.1)

6 結(jié)論

為了解決多批次數(shù)據(jù)下傳統(tǒng)優(yōu)化方法存在求解困難的情況，本研究通過對(duì)模型進(jìn)行標(biāo)度化處理解決了不同數(shù)據(jù)間的量綱差異，然后引入擴(kuò)展目標(biāo)函數(shù)提出了基于改進(jìn)隨機(jī)梯度下降的參數(shù)估計(jì)方法，并且給出了算法收斂的理論性分析。通過對(duì)主流的隨機(jī)優(yōu)化算法進(jìn)行數(shù)值實(shí)驗(yàn)對(duì)比，驗(yàn)證了所提出的SEOAG算法在估計(jì)精度、收斂速度以及受學(xué)習(xí)率影響方面的優(yōu)越性。

[1] LIM H C, HENRY C, SHIN H S. Fed-batch cultures (principles and applications of semi-batch bioreactors) [M]. New York: Cambridge University Press, 2013.

[2] MCLEAN K, MCAULEY K B. Mathematical modelling of chemical processes—obtaining the best model predictions and parameter estimates using identifiability and estimability procedures [J]. Canadian Journal of Chemical Engineering, 2012, 90(2): 351-366.

[3] KARIMI H, MCAULEY K B. A maximum-likelihood method for estimating parameters, stochastic disturbance intensities and measurement noise variances in nonlinear dynamic models with process disturbances [J]. Computers & Chemical Engineering, 2014, 67(4): 178-198.

[4] DUIJN M, GILE K J, HANDCOCK M S. A framework for the comparison of maximum pseudo-likelihood and maximum likelihood estimation of exponential family random graph models [J]. Social Networks, 2009, 31(1): 52-62.

[5] CANNARILE F, COMPARE M, ROSSI E,. A fuzzy expectation maximization based method for estimating the parameters of a multi-state degradation model from imprecise maintenance outcomes [J]. Annals of Nuclear Energy, 2017, 110(17): 739-752.

[6] HE L, HU M K, WEI Y J,. State of charge estimation by finite difference extended Kalman filter with HPPC parameters identification [J]. Science China (Technological Sciences), 2020, 63(3): 410-421.

[7] XIA W, DAI X X, FENG Y. Bayesian-MCMC-based parameter estimation of stealth aircraft RCS models [J]. Chinese Physics B, 2015, 24(12): 622-628.

[8] LI X L. Preconditioned stochastic gradient descent [J]. IEEE Transaction on Neural Networks and Learning Systems, 2018, 29(5): 1454-1466.

[9] ROBBINS H, MONRO S. A stochastic approximation method [J]. Annals of Mathematical Statistics, 1951, 22(3): 400-407.

[10] VAIDYA J, YU H, JIANG X Q. Privacy-preserving SVM classification [J]. Knowledge and Information Systems, 2008, 14(2): 161-178.

[11] SONG W J, ZHU J K, LI Y,. Image alignment by online robust PCA via stochastic gradient descent [J]. IEEE Transactions on Circuits and Systems for Video Technology, 2016, 26(7): 1241-1250.

[12] KALANTZIS V, VASSILIS K. A spectral newton-schur algorithm for the solution of symmetric generalized eigenvalue problems [J]. Electronic Transactions on Numerical Analysis, 2020, 52: 132-153.

[13] ZHAO R B, TAN V. A unified convergence analysis of the multiplicative update algorithm for regularized nonnegative matrix factorization [J]. IEEE Transactions on Signal Processing, 2018, 66(1): 129-138.

[14] DUCHI J C, FENG R. Stochastic methods for composite optimization problems [J]. SIAM Journal on Optimization, 2017, 28(4): 3229-3259.

[15] QIAN N. On the momentum term in gradient descent learning algorithms [J]. Neural Networks, 1999, 12(1): 145-151.

[16] MAZHAR H, HEYN T, NEGRUT D,. Using Nesterov's method to accelerate multibody dynamics with friction and contact [J]. ACM Transactions on Graphics, 2015, 34(3): 1-14.

[17] JOHNSON R, ZHANG T. Accelerating stochastic gradient descent using predictive variance reduction [J]. News in Physiological Sciences, 2013, 1(3): 315-323.

[18] SATO H, KASAI H, MISHRA B. Riemannian stochastic variance reduced gradient algorithm with retraction and vector transport [J]. SIAM Journal on Optimization, 2019, 29(2): 1444-1472.

[19] SCHMIDT M, LE R N, BACH F. Minimizing finite sums with the stochastic average gradient [J]. Mathematical Programming, 2017, 162(1): 83-112.

[20] COURTY N, GONG X, VANDEL J,. SAGA: sparse and geometry-aware non-negative matrix factorization through non-linear local embedding [J]. Machine Learning, 2014, 97(1): 205-226.

[21] BAE J, JEONG D H, LEE J M. Ranking-based parameter subset selection for nonlinear dynamics with stochastic disturbances under limited data [J]. Industrial & Engineering Chemistry Research, 2020, 59(50): 21854-21868.

[22] CHEN W F, SHAO Z J, BIEGLER L T. A bilevel NLP sensitivity-based decomposition for dynamic optimization with moving finite elements [J]. AIChE Journal, 2014, 60(3): 966-979.

[23] BOYD S, VANDENBERGHE L. Convex optimization [M]. New York: Cambridge University Press, 2004.

[24] QIAN Q, JIN R, YI J F,. Efficient distance metric learning by adaptive sampling and mini-batch stochastic gradient descent (SGD) [J]. Machine Learning, 2015, 99(3): 353-372.

Estimation of reaction kinetic parameters based on modified stochastic gradient descent

TANG Li-sen, CHEN Wei-feng

(School of Information Engineering, Zhejiang University of Technology, Hangzhou 310023,China)

Considering the solution difficulty of conventional optimization algorithm in parameter estimation using all sampled data, a reaction kinetic parameter estimation method based on modified stochastic gradient descent was proposed by introducing stochastic optimization and extended objective function in the framework of simultaneous solution. Firstly, the mechanism of large-scale system with multiple data sets was modeled, and the sensitivity matrix was obtained based on the sensitivity differential equation method, and the model scaling technique was used to deal with the simultaneous convergence problem of multi-state variables to multi-parameter estimation. In order to reduce the influence of noise variance in the iterative process, based on the existing stochastic average gradient descent method, the stochastic extended objective function was applied to increase the amount of information for calculating the gradient in the objective function, and the theoretical convergence of the method was given. Relevant numerical simulation results have verified the effectiveness and feasibility of the proposed method.

parameter estimation; stochastic optimization; extended objective; sensitivity matrix

TQ03

A

10.3969/j.issn.1003-9015.2022.03.015

1003-9015(2022)03-0426-11

2021-06-11；

2021-08-17。

國家重點(diǎn)研發(fā)計(jì)劃(2017YFE0106700)；國家自然科學(xué)基金(61873242)。

唐立森 (1997-)，男，江西上饒人，浙江工業(yè)大學(xué)碩士生。

陳偉鋒，E-mail：wfchen@zjut.edu.cn

唐立森, 陳偉鋒. 基于改進(jìn)隨機(jī)梯度下降的反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì)方法[J]. 高?；瘜W(xué)工程學(xué)報(bào), 2022, 36(3): 426-436.

：TANG Li-sen,CHEN Wei-feng. Estimation of reaction kinetic parameters based on modified stochastic gradient descent [J]. Journal of Chemical Engineering of Chinese Universities, 2022, 36(3): 426-436.