孫燕 張弛 路興龍 王靖戈 付俊

動(dòng)態(tài)優(yōu)化是約束中含有微分或差分方程的一類數(shù)學(xué)規(guī)劃問(wèn)題[1?3],其廣泛存在于電力系統(tǒng)[4?5]、石油化工[6?7]、生物工程[8?9]、清潔能源[10?11]等. 目前在化工過(guò)程中基于常微分方程模型的動(dòng)態(tài)優(yōu)化問(wèn)題不僅涉及微分方程,還包括代數(shù)方程約束,這樣的系統(tǒng)統(tǒng)稱為微分代數(shù)系統(tǒng).此外,實(shí)際的流程工業(yè)過(guò)程中,需要某些狀態(tài)變量或控制變量的函數(shù)在其全部運(yùn)行過(guò)程中或部分運(yùn)行時(shí)間內(nèi)不能超過(guò)約束的設(shè)定值,來(lái)保證工業(yè)過(guò)程的質(zhì)量和安全,這就需要考慮具有不等式路徑約束的動(dòng)態(tài)優(yōu)化問(wèn)題.

動(dòng)態(tài)優(yōu)化問(wèn)題(Dynamic optimization problem,DOP)的數(shù)值求解算法通常分為直接法[12?13]和間接法[14?15].間接法通過(guò)求解原問(wèn)題的最優(yōu)性條件(即必要條件),間接地獲得原問(wèn)題的最優(yōu)解.但是對(duì)于相對(duì)復(fù)雜的問(wèn)題,比如包含不等式約束的最優(yōu)控制問(wèn)題,由于狀態(tài)變量或者控制變量在約束中存在,求解極為困難[16].此外,復(fù)雜的非線性系統(tǒng)最優(yōu)性條件很難被定義,且兩點(diǎn)邊值問(wèn)題的求解過(guò)程中收斂域可能很小[17].直接法是通過(guò)對(duì)控制變量離散化,即控制向量參數(shù)化[18](Control vector parameterization,CVP),或控制與狀態(tài)變量同時(shí)離散化,即正交配置(Orthogonal collocation,OC),將無(wú)限維的最優(yōu)控制問(wèn)題轉(zhuǎn)化為有限維的非線性最優(yōu)化問(wèn)題.由于OC方法的計(jì)算復(fù)雜度較高,本文采用CVP方法.

針對(duì)具有不等式路徑約束的微分代數(shù)方程(Differential-algebraic equations,DAE)系統(tǒng)的動(dòng)態(tài)優(yōu)化問(wèn)題,采用直接或間接法對(duì)其進(jìn)行轉(zhuǎn)化后,仍需要對(duì)等式路徑約束和不等式路徑約束進(jìn)行處理.胡云卿等將等式路徑約束進(jìn)行微分[19],將原問(wèn)題轉(zhuǎn)化為具有不等式的常微分方程(Ordinary differential equations,ODE)動(dòng)態(tài)優(yōu)化問(wèn)題,再使用基于ODE的動(dòng)態(tài)優(yōu)化方法進(jìn)行處理,但是DAE到ODE的轉(zhuǎn)化過(guò)程中需要考慮初值條件的相容性,并且無(wú)法在理論上證明有關(guān)解的存在性和唯一性[20],此外,當(dāng)約束中變量的冪較高或者變量間耦合程度較高的情況下是不可行的.Pontryagin等將等式路徑約束轉(zhuǎn)化為兩點(diǎn)邊值問(wèn)題[21]進(jìn)行求解,該方法需要增加約束條件,增加了求解難度.Fabien將等式約束轉(zhuǎn)化為終點(diǎn)約束[22],這樣會(huì)導(dǎo)致問(wèn)題的維數(shù)變得很大,需要對(duì)轉(zhuǎn)化后的問(wèn)題進(jìn)行降維處理.所以間接地處理DAE中的等式路徑約束會(huì)引起一系列的問(wèn)題,因此本文采用后向差分法(Backward differentiation formula,BDF)[23]直接求解DAE.

對(duì)于不等式路徑約束,文獻(xiàn)[24]采用多重打靶(Multiple shooting,MS)法,將不等式路徑約束轉(zhuǎn)化為分段的點(diǎn)約束,每次打靶均需要給出對(duì)應(yīng)打靶分段內(nèi)該點(diǎn)約束涉及的控制變量和各狀態(tài)變量合理的范圍,各變量段間的連續(xù)還需要一個(gè)合適的系數(shù)保證,不然很難得到理想的結(jié)果.文獻(xiàn)[25]對(duì)不等式路徑約束的積極與否進(jìn)行判斷,若不等式路徑約束為積極約束,則結(jié)合微分方程在DAE求解器中進(jìn)行求解,否則不考慮該不積極約束;不等式路徑約束為積極約束時(shí),DAE階次可能很高,所以該方法對(duì)DAE求解器要求較高.Jacobson等提出的松弛變量法[26],通過(guò)增加松弛變量將不等式路徑約束轉(zhuǎn)化為相應(yīng)的等式路徑約束,但是這種方法難以處理不等式路徑約束數(shù)比控制變量數(shù)多的問(wèn)題;Vassiliadis等[27]通過(guò)離散化將約束項(xiàng)在問(wèn)題時(shí)間約束段內(nèi)選擇有限個(gè)內(nèi)點(diǎn)進(jìn)行離散化,通過(guò)懲罰函數(shù)法將約束違反的積分設(shè)為0來(lái)保證約束的無(wú)違反,但是此方法可能會(huì)在NLP問(wèn)題中引進(jìn)較多的約束條件;Floudas提出的凸函數(shù)近似法[28]主要應(yīng)用于靜態(tài)系統(tǒng),Chachuat等[29]將其改進(jìn)應(yīng)用在動(dòng)態(tài)系統(tǒng)的最優(yōu)控制上,此方法求得的數(shù)值解雖然可以嚴(yán)格滿足不等式路徑約束,但算法復(fù)雜度較高;Rehbock等提出的精確罰函數(shù)法[30]能夠成功求解具有不等式路徑約束的復(fù)雜最優(yōu)控制問(wèn)題[31],該方法雖然在算法復(fù)雜度上有所改善但是不能嚴(yán)格滿足不等式約束,且約束違反程度不能確定;Liu等[32?33]提出了一種新穎的光滑化精確懲罰函數(shù)法,將不等式路徑約束轉(zhuǎn)化為一個(gè)光滑化懲罰項(xiàng)加到目標(biāo)函數(shù)中,將問(wèn)題轉(zhuǎn)變成一個(gè)無(wú)約束優(yōu)化問(wèn)題;由于目標(biāo)函數(shù)表達(dá)式會(huì)因不等式約束的違反程度隨時(shí)間而改變,所以如何獲取問(wèn)題轉(zhuǎn)化后目標(biāo)函數(shù)的梯度信息是一個(gè)難點(diǎn).2005年Chen在之前的CVP-OC混合方法基礎(chǔ)上進(jìn)一步提出了有限收斂法[34],問(wèn)題初始求解時(shí)不考慮不等式路徑約束的最優(yōu)控制問(wèn)題,將初始所得解做為初始解以此優(yōu)化原問(wèn)題,這樣可以使得優(yōu)化過(guò)程得以優(yōu)化.然而該方法受限于帶有不等式路徑約束的ODE動(dòng)態(tài)優(yōu)化問(wèn)題,對(duì)于帶有不等式路徑約束的非0階DAE動(dòng)態(tài)系統(tǒng)的適用性有待探究.

因此,本文針對(duì)具有不等式路徑約束的DAE系統(tǒng)的動(dòng)態(tài)優(yōu)化問(wèn)題,采用只對(duì)控制變量進(jìn)行參數(shù)化的CVP方法,將無(wú)限維的DOP轉(zhuǎn)化為有限維的DOP,仍然保留系統(tǒng)的動(dòng)態(tài)特性.通過(guò)分點(diǎn)離散化方法將不等式路徑約束轉(zhuǎn)化為有限的內(nèi)點(diǎn)約束進(jìn)行處理,設(shè)計(jì)了在指定的路徑約束違反容忍度下通過(guò)有限步迭代獲得KKT(Karush Kuhn Tucker)最優(yōu)點(diǎn)的算法.然后利用BDF方法求解DAE方程組,再采用序列二次規(guī)劃法(Sequential quadratic programming,SQP)獲得滿足不等式路徑約束的最優(yōu)解.最后將該算法應(yīng)用在工業(yè)化工問(wèn)題中,仿真結(jié)果表明了該算法可以獲得最優(yōu)控制軌跡并有著良好的路徑約束效果及收斂性.

本文的主要貢獻(xiàn)如下:

1)采用后向差分法直接求解DAE,避免將DAE化為ODE的過(guò)程中會(huì)產(chǎn)生的解的相容性、維數(shù)增加等前文所述問(wèn)題.2)所采用的處理路徑約束的方法不會(huì)改變目標(biāo)函數(shù)的結(jié)構(gòu)形式.并且可以在指定路徑約束違反的容忍度的條件下,經(jīng)過(guò)有限步迭代得到最優(yōu)解.

本文結(jié)構(gòu)如下:第1節(jié)對(duì)具有不等式路徑約束的DAE系統(tǒng)的最優(yōu)控制問(wèn)題進(jìn)行描述.第2節(jié)闡述本文所采用的求解動(dòng)態(tài)優(yōu)化問(wèn)題和處理不等式路徑約束的算法.第3節(jié)將算法應(yīng)用于工業(yè)化工問(wèn)題.最后對(duì)本文所做工作進(jìn)行了總結(jié)并提出了對(duì)未來(lái)可以改進(jìn)的方面進(jìn)行了展望.

1 問(wèn)題描述

本文考慮的具有不等式路徑約束的DAE系統(tǒng)的最優(yōu)控制問(wèn)題,根據(jù)Biegler[35]給出的動(dòng)態(tài)優(yōu)化問(wèn)題描述形式添加不等式路徑約束條件,記為P1,給出數(shù)學(xué)描述如下:

s.t.

其中,x(t)∈Rnx是狀態(tài)變量,y(t)∈Rny是代數(shù)變量,v(t)∈Rnu是控制變量且上下邊界約束分別為ul和uu,這些變量都是時(shí)間標(biāo)量t∈[t0,tf]的函數(shù);動(dòng)態(tài)過(guò)程開始時(shí)刻t0和結(jié)束時(shí)刻tf都是固定的;動(dòng)態(tài)系統(tǒng)特性由常微分方程(2)組成;其初始可行狀態(tài)為xx0;動(dòng)態(tài)優(yōu)化問(wèn)題的最優(yōu)解需要滿足等式路徑約束(3)、不等式路徑約束(4);目標(biāo)函數(shù)式(1)由終值項(xiàng)φ[x(tf),y(tf)]和被積函數(shù)為‘[t,x(t),y(t),v(t)]的積分項(xiàng)組成.其中,等式路徑約束(3)和常微分方程(2)組成DAE.為不失一般性,假設(shè)這里的DAE對(duì)各變量均有連續(xù)二階偏導(dǎo)數(shù),且index為1.

這些DAE系統(tǒng)具有更大的靈活性:可以是線性的也可以是非線性的;不等式路徑約束(4)可以有多個(gè),若有多個(gè)不等式路徑約束,其處理方法同(4).

該最優(yōu)控制問(wèn)題可以簡(jiǎn)單描述為:在可行初始狀態(tài)(5)條件下,尋找最優(yōu)的控制軌跡v?(t),使得式(2)和式(3)描述的系統(tǒng)滿足約束(4),且使得目標(biāo)函數(shù)(1)的值最小.

2 主要結(jié)果

本文針對(duì)具有不等式路徑約束的DAE系統(tǒng)的最優(yōu)控制問(wèn)題,首先利用CVP技術(shù)將無(wú)限維的DOP轉(zhuǎn)化為有限維的DOP,將原問(wèn)題轉(zhuǎn)化為仍保留系統(tǒng)動(dòng)態(tài)特性的具有不等式路徑約束的DOP.對(duì)于系統(tǒng)的動(dòng)態(tài)方程,采用BDF直接求解,避免了將DAE轉(zhuǎn)化為ODE或兩點(diǎn)邊值問(wèn)題.利用逐點(diǎn)離散法處理不等式路徑約束,將不等式約束在違反處轉(zhuǎn)化為有限的點(diǎn)約束,在有限步的迭代下,得到滿足用戶指定的路徑約束違反容忍度的KKT最優(yōu)點(diǎn).算法的主要結(jié)構(gòu)如圖1所示.

圖1 算法主要結(jié)構(gòu)圖Fig.1 Structure diagram of the algorithm

本節(jié)首先簡(jiǎn)要介紹了控制向量參數(shù)化方法及其處理后得到的有限維的DOP.然后在逐點(diǎn)離散法中,路徑約束被轉(zhuǎn)化為點(diǎn)約束的形式,同時(shí)給出了不等式路徑約束的DAE系統(tǒng)最優(yōu)控制問(wèn)題的具體求解算法.并給出了分點(diǎn)離散法收斂性的證明.

2.1 控制向量參數(shù)化

CVP是求解動(dòng)態(tài)優(yōu)化問(wèn)題的一種重要策略.CVP方法的核心思想[36]是:利用基函數(shù)將控制變量進(jìn)行參數(shù)化,將原無(wú)限維最優(yōu)控制問(wèn)題轉(zhuǎn)化為一個(gè)含有動(dòng)態(tài)過(guò)程的有限維優(yōu)化問(wèn)題后,再用傳統(tǒng)的DOP求解方法進(jìn)行求解.即將控制向量v(t)被離散為n段而狀態(tài)變量保持連續(xù),每個(gè)時(shí)間分段內(nèi)用帶參數(shù)的多項(xiàng)式(如常數(shù)、線性或二次多項(xiàng)式).關(guān)于控制變量參數(shù)化方法的收斂性,文獻(xiàn)[37?38]都進(jìn)行了相應(yīng)的證明.

首先不考慮不等式路徑約束,將動(dòng)態(tài)優(yōu)化問(wèn)題的時(shí)間域[t0,tf]劃分為N個(gè)控制階段,

在某一時(shí)間分段[tk?1,tk],k=1,2,···,N內(nèi),根據(jù)函數(shù)逼近理論,控制變量可以表示為低階多項(xiàng)式參數(shù)依賴的表達(dá)形式,

其中,j表示原問(wèn)題有 1,2,···,j,···,nu個(gè)控制變量,為多項(xiàng)式基函數(shù),M為基函數(shù)階數(shù).實(shí)際計(jì)算和應(yīng)用中控制向量的參數(shù)化可以有多種參數(shù)依賴形式,如圖2所示常用的有分段常量逼近策略(如圖2實(shí)線)、分段線性逼近策略(如圖2虛線)、分段二次多項(xiàng)式逼近策略(如圖2點(diǎn)虛線)等.由于計(jì)算機(jī)采用數(shù)字信號(hào),非常符合分段常量策略形式,同時(shí)這種策略簡(jiǎn)便易行、生成參數(shù)數(shù)目較少.因此本文采用分段常量逼近策略進(jìn)行控制變量參數(shù)化,即采用為基函數(shù).

圖2 控制變量分段示例Fig.2 Illustration of control variable pro files

從而整個(gè)時(shí)間段內(nèi)的控制變量可以表示為:

其中,X[tk?1,tk)(t) 為開關(guān)函數(shù),k=1,2,···,N,具體形式如下:

進(jìn)而每一分段內(nèi),控制變量可以表示為:

利用上述中CVP思想對(duì)原問(wèn)題(1)～(7)進(jìn)行處理,得到如下有限維動(dòng)態(tài)優(yōu)化問(wèn)題,記為P2:

s.t.

其中,u∈U是v(t)離散化后的有限維控制向量,并且U∈RN是個(gè)非空N維向量.其他變量的定義同問(wèn)題P1.

對(duì)于參數(shù)化后的DOP問(wèn)題可使用成熟的基于梯度的優(yōu)化算法進(jìn)行求解.本文首先采用BDF直接對(duì)微分代數(shù)方程組進(jìn)行求解,進(jìn)而得到帶有有限維控制向量的目標(biāo)函數(shù),再用BFGS迭代法計(jì)算拉格朗日函數(shù)的Hessian矩陣的正定擬牛頓近似值,利用SQP算法便可得到最優(yōu)的控制參數(shù)組合.

2.2 不等式路徑約束處理

2.2.1 算法步驟

本文利用文獻(xiàn)[34]中的分點(diǎn)離散法來(lái)處理不等式路徑約束,記分段點(diǎn)集合PN={t0,t1,···,tN?1,tN≡tf},ti?1

其中,tj,j=1,2,···,M是時(shí)間域內(nèi)M個(gè)違反時(shí)間段上的離散形式.給出求解具有不等式路徑約束的DAE系統(tǒng)最優(yōu)控制問(wèn)題求解具體算法如下:

初始化.設(shè)置初始控制參數(shù)u(0)、初始化不等式路徑約束違反點(diǎn)集合T(0)??T;有限且任意小的容忍度參數(shù)η>0;分段參數(shù)N∈N?;迭代次數(shù)p=0;迭代精度δ>0;

步驟1.先解決不帶有不等式路徑約束的DOP問(wèn)題P2,i.e.解式(13)～(15)、式(17)～(19)構(gòu)成的動(dòng)態(tài)優(yōu)化問(wèn)題;

步驟2.利用BDF法求解式(14),(15)和式(17)構(gòu)成的DAE初值問(wèn)題,得到xx(t)的值;

步驟3.在容忍度 下檢查h(t)≤η,t∈[t0,tf]的違反情況.將每個(gè)CVP分段中違反量最大的點(diǎn)對(duì)應(yīng)的點(diǎn)約束加到DOP的約束中,并且檢查已有違反點(diǎn)的違反情況,若不存在違反則去除該點(diǎn),更新;

步驟4.使用BFGS迭代法計(jì)算拉格朗日函數(shù)的Hessian矩陣的正定擬牛頓近似值獲取關(guān)于參數(shù)向量的梯度信息;

步驟5.通過(guò)SQP方法,求解QP子問(wèn)題,生成新的迭代步;

步驟6.若在容忍度η下h(ti)≤η,其中,且,則算法結(jié)束;否則重新求解帶現(xiàn)有點(diǎn)約束的NLP問(wèn)題,p=p+1,繼續(xù)步驟2.

若問(wèn)題P2中的不等式路徑約束存在違反,則會(huì)循環(huán)迭代求解帶不同點(diǎn)約束的DOP問(wèn)題,通過(guò)SQP方法調(diào)整控制變量u,使得點(diǎn)約束的違反程度有效降低,直至收斂至容忍度η內(nèi).不等式路徑約束到點(diǎn)約束的轉(zhuǎn)換是收斂的(詳細(xì)內(nèi)容參考第2.2.2節(jié)),因此通過(guò)有限步迭代,能找到滿足指定容忍度的KKT最優(yōu)點(diǎn).

2.2.2 收斂性證明

在原DOP問(wèn)題經(jīng)過(guò)CVP處理后的N個(gè)控制段上,對(duì)于分點(diǎn)離散法的收斂性分析,本文給出以下證明:

對(duì)于DAE問(wèn)題P1的不等式路徑約束h(t)=均有以下假設(shè):

假設(shè)1.h(t)至少是一階可導(dǎo)的.

假設(shè)2.h(t)對(duì)于t是Lipschitz連續(xù)的,即存在一個(gè)常數(shù)K,使得對(duì)所有(t)和(),在域D={t0≤t≤tf,|x|≤∞,|u|≤∞}上都有下式成立:

假設(shè)3.若h(t)對(duì)于t不是Lipschitz連續(xù)的,記不連續(xù)點(diǎn)集合為Pnon且Pnon∈PM,則h(t)滿足假設(shè)2,其中tPnon.

定理1.若問(wèn)題P1的任意一個(gè)不等式路徑約束h(t)≤0都滿足上述假設(shè),則對(duì)于該h(t)存在有限的點(diǎn)集使得h(t)≤η:η>0,t∈[t0,tf].

證明. 原DOP問(wèn)題經(jīng)CVP處理后,時(shí)間域T∈[t0,tf],i=0,1,···,N?1會(huì)被劃分為N段相等或者不等的子時(shí)間段Ti=[ti,ti+1],i=0,1,···,N?1.在每個(gè)子時(shí)間段Ti上求解DAE初值問(wèn)題時(shí),該子時(shí)間段將被離散為M個(gè)小時(shí)間分段,k=0,1,···,M?1,其中=ti,.對(duì)于任意一個(gè)小時(shí)間分段,取由h(t)在上連續(xù),在內(nèi)可導(dǎo),那么在內(nèi)至少存在一點(diǎn),使等式成立,當(dāng)時(shí),對(duì)等式進(jìn)一步處理有:

去絕對(duì)值進(jìn)行移項(xiàng),即

綜上,在整個(gè)時(shí)間域[t0,tf]上的每個(gè)子時(shí)間段,i=0,1,···,N(η)?1上都有上式成立,因此有

其中,N(η)是η的函數(shù),η是式(22)滿足的條件.對(duì)于有限的η,N(η)也是有限的.

對(duì)于問(wèn)題P2的每一個(gè)不等式路徑約束均可此類證明.

3 仿真研究

本節(jié)采用兩個(gè)DAE系統(tǒng)的最優(yōu)控制問(wèn)題對(duì)提出的算法進(jìn)行仿真研究.仿真環(huán)境為:Intel i7-4790,3.6GHz CPU和DDMM1/1600 MHZ 8GB內(nèi)存.在仿真中求解DAE方程組求解精度設(shè)為10?6;使用SQP方法解NLP問(wèn)題,求解精度設(shè)為10?4;內(nèi)點(diǎn)約束的收斂容忍度為η=10?6.

3.1 催化劑混合問(wèn)題

由Guun等提出的催化劑混合問(wèn)題[39],一直以來(lái)被研究者研究用于動(dòng)態(tài)優(yōu)化問(wèn)題的研究[40?41].問(wèn)題可描述為:在長(zhǎng)度為lf的管式反應(yīng)器中,兩種物質(zhì)A、B在管內(nèi)發(fā)生反應(yīng).反應(yīng)式為A?B→C.其中B是副產(chǎn)物,C是目標(biāo)產(chǎn)物,反應(yīng)速率隨催化劑沿著管長(zhǎng)的混合率而變化.為了得到更多的目標(biāo)產(chǎn)物C,該動(dòng)態(tài)優(yōu)化控制問(wèn)題的求解目標(biāo)是獲得催化劑最佳配比方案.其數(shù)學(xué)模型如下式所示:

s.t.

其中,狀態(tài)變量x1(l)、x2(l)、x3(l)分別表示為A、B、C的摩爾分率,初始狀態(tài)為x1(0)=1,x2(0)=0,x3(0)=0;控制變量u(l)表示催化劑沿管長(zhǎng)的混合率,滿足的約束為0≤u(l)≤1;l在這里表示管長(zhǎng).這里取l0=0,lf=12,仿真中取u0=0.17.

根據(jù)本文給出的具有不等式路徑約束的DAE系統(tǒng)最優(yōu)控制問(wèn)題的求解算法,若不等式在[l0,lf]上存在違反,則將其轉(zhuǎn)化為內(nèi)點(diǎn)約束x2(li)?0.075≤η.圖3和圖4是控制變量取為N=24時(shí)得到的控制軌跡和狀態(tài)軌跡.圖5為x2(l)的狀態(tài)軌跡.可以看出在整個(gè)長(zhǎng)度域l∈[l0,lf]內(nèi)都滿足不等式路徑約束,說(shuō)明本文提出的方法在路徑約束處理方面有著很好的效果.

仿真中分別取N=24、29、34,得到問(wèn)題最優(yōu)目標(biāo)函數(shù)和求解時(shí)間如表1.由表1中可以看出,本文方法求解的時(shí)間與文獻(xiàn)[18]相比,省時(shí)至少20%,精確度方面會(huì)因?yàn)閰?shù)有所浮動(dòng),但是目標(biāo)函數(shù)最差在N=29時(shí),本文所得目標(biāo)函數(shù)也已經(jīng)達(dá)到文獻(xiàn)[18]的99.97%.

圖3 控制曲線Fig.3 Control pro file

圖4 路徑約束下的狀態(tài)曲線Fig.4 The status curve with path constraint

圖5 路徑約束下x2的狀態(tài)曲線Fig.5 The status curve of x2with path constraint

表1 催化劑混合問(wèn)題測(cè)試結(jié)果Table 1 Results of catalyst mixing problem

3.2 青霉素分批補(bǔ)料發(fā)酵

在文獻(xiàn)[42]中描述的一個(gè)青霉素分批補(bǔ)料發(fā)酵問(wèn)題(Fed-batch penicillin fermentation),該問(wèn)題具有很強(qiáng)的非線性和奇異性.該動(dòng)態(tài)優(yōu)化問(wèn)題描述如下:

其中,目標(biāo)函數(shù)J是總收益.X是生物量濃度(g/L),P是現(xiàn)有青霉素產(chǎn)物數(shù)量(每升活性),具體增長(zhǎng)率μ,具體產(chǎn)品形成速率θ和反應(yīng)堆的體積V(L)都是系統(tǒng)的狀態(tài)變量.底物濃度S(g/L)是控制變量.進(jìn)料流量F為1666.67L/h.該模型初始狀態(tài)為 [X(0),P(0),V(0),μ(0)]=[1,0,250000,0]. 起始時(shí)刻[t0,tf]=[0,150].控制變量S∈[0.001,0.5].

采用文獻(xiàn)[18]中處理不等式路徑約束的方法與本文方法對(duì)比,仿真中取N=12、14、16分別得到最優(yōu)目標(biāo)函數(shù)和求解時(shí)間如表2所示.由表2中可以看出,本文方法求解的時(shí)間約為文獻(xiàn)[18]的十分之一,精確度方面會(huì)因?yàn)閰?shù)有所浮動(dòng),但是目標(biāo)函數(shù)最差在N=12時(shí),本文所得目標(biāo)函數(shù)也已經(jīng)達(dá)到文獻(xiàn)[18]的99.3%.

仿真時(shí)設(shè)置u0=0.5.取控制變量N=12時(shí),圖6為得到的控制軌跡,圖7是狀態(tài)軌跡,其中狀態(tài)變量X(t)≤41.從圖中可以看出本文所提方法的有效性.

表2 青霉素分批補(bǔ)料發(fā)酵問(wèn)題測(cè)試結(jié)果Table 2 Results of fed-batch penicillin fermentation

圖6 控制曲線Fig.6 Control pro file

圖7 路徑約束下的狀態(tài)曲線Fig.7 The status curve with path constraint

綜上兩個(gè)仿真實(shí)驗(yàn),本文方法允許路徑約束違反在用戶指定范圍之內(nèi),在問(wèn)題允許存在誤差或者實(shí)際上確實(shí)存在誤差的情況下,以一定的準(zhǔn)確度換取求解速度,滿足應(yīng)用需求.

4 結(jié)論

本文針對(duì)帶有不等式路徑約束的DAE系統(tǒng)的動(dòng)態(tài)優(yōu)化問(wèn)題,提出了一種直接求解該類型問(wèn)題的框架,首先,通過(guò)CVP將無(wú)限維的動(dòng)態(tài)優(yōu)化問(wèn)題轉(zhuǎn)化為有限維的動(dòng)態(tài)優(yōu)化問(wèn)題;然后,利用分點(diǎn)離散方法,設(shè)計(jì)了一個(gè)可以在一定容忍度下滿足路徑約束的算法,并在理論上論證了該算法可以在有限步數(shù)的迭代下收斂;最后,利用催化劑混合問(wèn)題和壓力限定批量反應(yīng)堆問(wèn)題仿真研究驗(yàn)證了本文所提方法的有效性.

本文對(duì)控制變量進(jìn)行的是等分離散化,但對(duì)局部變化較大的控制信號(hào),等分離散化不能準(zhǔn)確地反應(yīng)真實(shí)的控制信號(hào).因此本文下一步的工作將研究根據(jù)控制信號(hào)曲線自適應(yīng)地進(jìn)行CVP離散化.