郭瑋林,鮮 勇,張大巧,凌王輝
(火箭軍工程大學(xué) 七系,西安 710025)
由于高超聲速飛行器助推段彈道設(shè)計(jì)約束條件多,動(dòng)力學(xué)過(guò)程復(fù)雜,終端入軌條件苛刻,就目前研究情況而言,實(shí)現(xiàn)助推段彈道優(yōu)化設(shè)計(jì)本身就是一個(gè)難點(diǎn)問(wèn)題[1-3],若要實(shí)現(xiàn)助推段彈道的快速高精度計(jì)算滿(mǎn)足高超聲速飛行器機(jī)動(dòng)發(fā)射要求更是一個(gè)巨大的挑戰(zhàn)。
文獻(xiàn)[4]針對(duì)高超聲速飛行器助推段彈道設(shè)計(jì)方面進(jìn)行了研究,但彈道計(jì)算速度和精度有待提高。文獻(xiàn)[5]聯(lián)合間接法和直接法可在3 min內(nèi)計(jì)算生成高超聲速飛行器助推段最優(yōu)軌跡,但彈道計(jì)算耗時(shí)較長(zhǎng)未達(dá)到快速計(jì)算的要求。文獻(xiàn)[6]主要是運(yùn)用牛頓迭代法實(shí)現(xiàn)對(duì)助推段彈道的快速計(jì)算,但牛頓法對(duì)初值敏感且易陷入局部最優(yōu)。而L-M算法作為牛頓法的改進(jìn),融合了牛頓迭代局部快速收斂和梯度下降法全局迭代迅速的優(yōu)點(diǎn),通過(guò)引入阻尼因子,能夠有效避免雅可比矩陣奇異時(shí)迭代不收斂問(wèn)題[7]。
BP神經(jīng)網(wǎng)絡(luò)是一種多層前饋網(wǎng)絡(luò),網(wǎng)絡(luò)結(jié)構(gòu)簡(jiǎn)單實(shí)時(shí)性好,只要有足夠的隱層和隱節(jié)點(diǎn),就可以逼近任意的非線(xiàn)性映射關(guān)系[8]?;贐P神經(jīng)網(wǎng)絡(luò)建立發(fā)射點(diǎn)及終端入軌點(diǎn)狀態(tài)量與彈道參數(shù)的映射關(guān)系,可以實(shí)現(xiàn)彈道諸元參數(shù)快速預(yù)測(cè),為后續(xù)L-M方法數(shù)值尋優(yōu)計(jì)算提供彈道參數(shù)初值。
因此本文采用基于BP神經(jīng)網(wǎng)絡(luò)和L-M方法的聯(lián)合算法實(shí)現(xiàn)助推段彈道快速計(jì)算。首先綜合考慮各項(xiàng)約束條件設(shè)計(jì)了助推段飛行程序和彈道優(yōu)化模型;然后利用 BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)計(jì)算諸元參數(shù)初值,并且建立了基于BP神經(jīng)網(wǎng)絡(luò)和L-M方法的聯(lián)合數(shù)值計(jì)算模型;最后通過(guò)仿真實(shí)驗(yàn)驗(yàn)證了聯(lián)合算法能夠有效地實(shí)現(xiàn)助推段彈道的快速高精度計(jì)算。
以三級(jí)固體火箭作為高超聲速飛行器助推段的運(yùn)載器,假定地球是以角速度ω自轉(zhuǎn)的橢球,且飛行器運(yùn)動(dòng)姿態(tài)能夠達(dá)到瞬時(shí)平衡。
高超聲速飛行器的助推段關(guān)機(jī)點(diǎn)高度相比彈道導(dǎo)彈較低,所以一級(jí)飛行程序采取二次攻角轉(zhuǎn)彎的設(shè)計(jì)方式,即飛行器分別在跨音速前和跨音速后進(jìn)行一次攻角轉(zhuǎn)彎,跨音速時(shí)零攻角飛行??v向俯仰程序設(shè)計(jì)如下:
式中:φcx表示飛行程序角;0~t1為垂直上升段;t1~t2為跨音速前攻角轉(zhuǎn)彎段;t2~t3為跨音速飛行段,攻角為零;t3~t4為跨音速后攻角轉(zhuǎn)彎段;t4~tk1為一二級(jí)分離前的等程序瞄準(zhǔn)段;a1、a2、a~1、a~2為攻角設(shè)計(jì)參數(shù);θ為彈道傾角,ωz為地球自轉(zhuǎn)角速度分量。
二級(jí)和三級(jí)飛行程序采用分段線(xiàn)性化方法進(jìn)行設(shè)計(jì),具體形式如公式(3),其中:為二級(jí)等程序飛行段;為二級(jí)轉(zhuǎn)彎段,k1為該段的程序轉(zhuǎn)彎角速率;為級(jí)間分離段;為三級(jí)轉(zhuǎn)彎段,為該段的程序轉(zhuǎn)彎角速率;為三級(jí)等程序飛行段,tk3為三級(jí)關(guān)機(jī)點(diǎn)。
1)助推段約束條件
由于受到彈體結(jié)構(gòu)、執(zhí)行機(jī)構(gòu)性能的限制,必須限制飛行器轉(zhuǎn)彎角速率、過(guò)載和動(dòng)壓,以滿(mǎn)足執(zhí)行機(jī)構(gòu)性能約束和穩(wěn)定飛行的要求。為減小被敵方紅外偵測(cè)系統(tǒng)識(shí)別的可能性,飛行器低彈道飛行時(shí),高度須低于 90 km;同時(shí),為順利實(shí)現(xiàn)入軌飛行,助推段關(guān)機(jī)點(diǎn)還必須滿(mǎn)足終端高度、速度和當(dāng)?shù)貜椀纼A角約束。具體約束形式如下所示:
2)優(yōu)化變量的確定
所以最終高超聲速飛行器助推段彈道參數(shù)確定為3個(gè)優(yōu)化變量,具體如下:
式中,Tcz為一級(jí)發(fā)動(dòng)機(jī)垂直上升段工作時(shí)間。
3)目標(biāo)函數(shù)
為滿(mǎn)足高超聲速飛行器終端入軌條件,采取罰函數(shù)方式將終端約束轉(zhuǎn)換為目標(biāo)函數(shù),即:
式中,a、b、c為權(quán)系數(shù)。
以某發(fā)射區(qū)域?yàn)槔M(jìn)行仿真分析(假定地球?yàn)闄E球,所以發(fā)射點(diǎn)經(jīng)度改變對(duì)終端參數(shù)無(wú)影響)。為實(shí)現(xiàn)高超聲速飛行器助推段彈道快速計(jì)算,采用 BP神經(jīng)網(wǎng)絡(luò)算法,通過(guò)訓(xùn)練樣本,建立發(fā)射點(diǎn)緯度、高程、射向及關(guān)機(jī)點(diǎn)高度、速度、彈道傾角與助推段彈道諸元參數(shù)的函數(shù)關(guān)系式,即建立以下映射關(guān)系:
根據(jù)BP神經(jīng)網(wǎng)絡(luò)基本原理,若確定了網(wǎng)絡(luò)層數(shù)、各層節(jié)點(diǎn)數(shù)、傳遞函數(shù)、初始權(quán)系數(shù)和閥值、學(xué)習(xí)方法等,BP網(wǎng)絡(luò)結(jié)構(gòu)也就得以確定。
由式(7)可知,神經(jīng)網(wǎng)絡(luò)輸入層節(jié)點(diǎn)為6個(gè),輸出層節(jié)點(diǎn)為3個(gè),所以高超聲速飛行器彈道模型為六輸入三輸出的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)。
根據(jù)Kosmogorov定理,一個(gè)具有合適參數(shù)和結(jié)構(gòu)的3層BP神經(jīng)網(wǎng)絡(luò)可以完成任意的n維到m維的映射,所以可以選取3層BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行仿真計(jì)算。隱層傳遞函數(shù)選用logsig函數(shù),輸出層傳遞函數(shù)選用purelin線(xiàn)性函數(shù),學(xué)習(xí)方法選用L-M方法。
采用網(wǎng)絡(luò)結(jié)構(gòu)增長(zhǎng)型方法進(jìn)行仿真計(jì)算,確定隱層節(jié)點(diǎn)數(shù)量。通過(guò)測(cè)試誤差的比較,選取隱層節(jié)點(diǎn)數(shù)為12個(gè)。最終確定的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示。
圖1 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 BP neural network structure
隨機(jī)選取該發(fā)射區(qū)域發(fā)射點(diǎn)緯度、高程和射向信息,利用優(yōu)化算法求解滿(mǎn)足飛行約束和入軌終端約束的高超聲速飛行器助推段彈道,從而得出5000組訓(xùn)練樣本。由于粒子群算法具有收斂速度快、易實(shí)現(xiàn)的特點(diǎn)[9-11],所以采用粒子群智能算法對(duì)彈道諸元參數(shù)進(jìn)行優(yōu)化求解,生成訓(xùn)練樣本。
利用圖1神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)對(duì)樣本進(jìn)行訓(xùn)練,得到BP神經(jīng)網(wǎng)絡(luò)輸入層至隱層的權(quán)值W1和閾值b1以及隱層至輸出層的權(quán)值W2和閾值b2。因此可以得到發(fā)射點(diǎn)位置射向、關(guān)機(jī)點(diǎn)高度、速度、彈道傾角與彈道諸元參數(shù)的函數(shù)關(guān)系式。根據(jù) BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和各層傳遞函數(shù)可得:
由式(8)和式(9)推導(dǎo)BP神經(jīng)網(wǎng)絡(luò)函數(shù)解析式的最終形式如下:
由式(10)可以快速預(yù)測(cè)計(jì)算發(fā)射區(qū)域任意發(fā)射點(diǎn)位相應(yīng)助推段彈道諸元參數(shù)。
訓(xùn)練神經(jīng)網(wǎng)絡(luò)的根本目的是確保訓(xùn)練好的網(wǎng)絡(luò)模型對(duì)非訓(xùn)練樣本具有好的泛化能力。為測(cè)試以上神經(jīng)網(wǎng)絡(luò)模型性能,再次隨機(jī)選取發(fā)射區(qū)域500組測(cè)試樣本,測(cè)試樣本發(fā)射點(diǎn)位為高超聲速飛行器助推段終端要求的標(biāo)準(zhǔn)高度、標(biāo)準(zhǔn)速度和標(biāo)準(zhǔn)彈道傾角為將代入BP神經(jīng)網(wǎng)絡(luò)解析式(10)即可計(jì)算得到彈道諸元參數(shù),將該諸元參數(shù)代入高超聲速飛行器標(biāo)準(zhǔn)彈道得到 BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的助推段終端高度、速度和彈道傾角與標(biāo)準(zhǔn)終端高度、速度和彈道傾角對(duì)比得到 BP神經(jīng)網(wǎng)絡(luò)終端高度計(jì)算偏差、速度計(jì)算偏差和彈道傾角計(jì)算偏差,如圖2~4所示。基于BP神經(jīng)網(wǎng)絡(luò)方法的樣本仿真統(tǒng)計(jì)結(jié)果如表1所示。
圖2 BP神經(jīng)網(wǎng)絡(luò)終端高度計(jì)算偏差Fig.2 Calculated deviation of terminal height based on BP neural network
圖3 BP神經(jīng)網(wǎng)絡(luò)終端速度計(jì)算偏差Fig.3 Calculated deviation of terminal velocity based on BP neural network
圖4 BP神經(jīng)網(wǎng)絡(luò)終端彈道傾角計(jì)算偏差Fig.4 Calculated deviation of terminal inclination angle based on BP neural network
表1 基于BP神經(jīng)網(wǎng)絡(luò)的樣本仿真結(jié)果Tab.1 Simulation results of samples based on BP neural network
由圖2至圖4和表1的仿真結(jié)果可知,基于BP神經(jīng)網(wǎng)絡(luò)求解的高超聲速飛行器助推段彈道具有一定的精度,但入軌精度并不高,特別是終端高度最大入軌偏差達(dá)到386.92 m,高度平均偏差也有121.56 m,所以?xún)H僅基于 BP神經(jīng)網(wǎng)絡(luò)計(jì)算飛行器助推段彈道精度是不夠的。
BP神經(jīng)網(wǎng)絡(luò)解算的助推段彈道終端參數(shù)偏差較大,不能滿(mǎn)足高精度入軌要求,但可以利用 BP神經(jīng)網(wǎng)絡(luò)解析式快速預(yù)測(cè)計(jì)算諸元參數(shù),將其作為后續(xù)數(shù)值尋優(yōu)計(jì)算的彈道參數(shù)初值,加快數(shù)值計(jì)算的收斂速度,以求得滿(mǎn)足終端參數(shù)精度要求的助推段彈道。
牛頓法具有收斂速度快、精度高等優(yōu)點(diǎn),但對(duì)初值敏感且易陷入局部最優(yōu),而L-M算法作為牛頓法的改進(jìn),通過(guò)引入阻尼因子λ,能夠有效避免雅可比矩陣奇異或病態(tài)時(shí)迭代不收斂情況發(fā)生。因此可以利用BP神經(jīng)網(wǎng)絡(luò)快速預(yù)測(cè)彈道諸元參數(shù)初值,然后采用L-M算法進(jìn)一步數(shù)值尋優(yōu)計(jì)算,以解決機(jī)動(dòng)發(fā)射條件下高超聲速飛行器彈道快速計(jì)算的問(wèn)題。
下面研究建立基于BP神經(jīng)網(wǎng)絡(luò)與L-M算法的聯(lián)合數(shù)值尋優(yōu)計(jì)算模型。
設(shè)助推段終端高度偏差Δhf、速度偏差Δvf和彈道傾角偏差ΔfΘ是諸元參數(shù)Tcz、k1、k2的三維偏差函數(shù)f(x),即
式中,x為狀態(tài)變量,將f(x)在x處進(jìn)行一階泰勒展開(kāi)可得如下公式:
式中,x0為彈道諸元參數(shù)初值,Δx為迭代步長(zhǎng),Q為助推段終端參數(shù)偏差對(duì)諸元參數(shù)的偏導(dǎo)數(shù)矩陣,可采用彈道求差法計(jì)算求解。Q矩陣形式如下:
其中,彈道諸元參數(shù)初值x0的計(jì)算采用BP神經(jīng)網(wǎng)絡(luò)解析式(10)求解,即得:
助推段終端偏差函數(shù)的最小二乘模型為
將式(12)代入式(15)得到聯(lián)合算法數(shù)值計(jì)算模型如下:
由文獻(xiàn)[12]可知式(16)的解可表示為
式中,λ為L(zhǎng)-M算法的阻尼因子,λI為阻尼項(xiàng)。
在迭代計(jì)算過(guò)程中:當(dāng)助推段終端參數(shù)偏差較大時(shí),可以選取較大的λ,使得L-M算法具有梯度下降法下降量大、迭代迅速的特點(diǎn);當(dāng)助推段終端參數(shù)偏差較小時(shí),可以選取較小的λ,使得L-M算法具有牛頓法二階收斂的特點(diǎn),從而快速收斂到最優(yōu)解。
Yamashita和Fukushima在文獻(xiàn)[13]證明了若阻尼因子時(shí),那么在一個(gè)比非奇異性要弱的局部誤差有界條件下,L-M算法將保持二次收斂。因此針對(duì)阻尼因子λ的選取,可采用如下自適應(yīng)策略:
式中,1δ、2δ為終端參數(shù)偏差量范數(shù)的閥值,且0<δ2<δ1。
基于聯(lián)合算法計(jì)算模型,根據(jù)式(17)進(jìn)行多次迭代數(shù)值尋優(yōu)計(jì)算,最終解得滿(mǎn)足助推段終端入軌精度要求的彈道諸元參數(shù)。
基于BP神經(jīng)網(wǎng)絡(luò)與L-M方法的聯(lián)合算法計(jì)算流程如圖5所示,聯(lián)合算法具體計(jì)算步驟如下:
1)將機(jī)動(dòng)發(fā)射點(diǎn)緯度、高程和射向信息與終端入軌點(diǎn)高度、速度和彈道傾角參數(shù)輸入到訓(xùn)練好的BP神經(jīng)網(wǎng)絡(luò)解析式中,計(jì)算得到彈道諸元初始值
2)再將諸元參數(shù)初值代入L-M算法數(shù)值計(jì)算模型中進(jìn)行迭代尋優(yōu)計(jì)算,生成新的彈道諸元參數(shù)
4)若助推段終端參數(shù)滿(mǎn)足飛行器終端入軌精度要求或迭代達(dá)到最大次數(shù),則退出計(jì)算,輸出結(jié)果;否則,返回步驟2)。
圖5 聯(lián)合算法計(jì)算流程圖Fig.5 Flow chart of the joint algorithm
隨機(jī)選取發(fā)射區(qū)域 500組發(fā)射點(diǎn)位作為仿真樣本,采用基于BP神經(jīng)網(wǎng)絡(luò)和L-M方法的聯(lián)合算法對(duì)高超聲速飛行器助推段彈道進(jìn)行仿真計(jì)算得到樣本助推段終端參數(shù)偏差值?;诼?lián)合算法的終端高度、速度和彈道傾角計(jì)算偏差分別如圖6~8所示,助推段彈道計(jì)算耗時(shí)Ts仿真如圖 9所示,基于聯(lián)合算法的樣本仿真統(tǒng)計(jì)結(jié)果如表2所示。本次仿真實(shí)驗(yàn)計(jì)算在普通微機(jī)上實(shí)現(xiàn)。
圖6 聯(lián)合算法終端高度計(jì)算偏差Fig.6 Calculated deviation of terminal height based on the joint algorithm
圖7 聯(lián)合算法終端速度計(jì)算偏差Fig.7 Calculated deviation of terminal velocity based on the joint algorithm
圖8 聯(lián)合算法終端彈道傾角計(jì)算偏差Fig.8 Calculated deviation of terminal inclination angle based on the joint algorithm
圖9 彈道計(jì)算耗時(shí)仿真圖Fig.9 Simulation of trajectory calculation time
由圖6至圖9和表2仿真結(jié)果可知,經(jīng)仿真計(jì)算生成的500條助推段彈道均滿(mǎn)足終端約束要求,并且終端高度偏差在2 m以?xún)?nèi),終端速度偏差在0.1 m/s以?xún)?nèi),終端彈道傾角偏差在0.01°以?xún)?nèi),保證高超聲速飛行器以較高的精度入軌。同時(shí)助推段彈道計(jì)算最大耗時(shí)為 2.84 s,平均耗時(shí)0.8211 s,彈道計(jì)算效率高、速度快。所以基于BP神經(jīng)網(wǎng)絡(luò)和L-M方法的聯(lián)合算法能夠較好地滿(mǎn)足機(jī)動(dòng)發(fā)射條件下高超聲速飛行器助推段彈道快速高精度計(jì)算的要求。
表2 基于聯(lián)合算法的樣本仿真結(jié)果Tab.2 Simulation results of samples based on joint algorithm
本文研究了機(jī)動(dòng)發(fā)射條件下的高超聲速飛行器助推段彈道快速計(jì)算問(wèn)題。基于 BP神經(jīng)網(wǎng)絡(luò)建立了發(fā)射點(diǎn)及終端入軌點(diǎn)狀態(tài)量與助推段彈道諸元參數(shù)的函數(shù)解析式,并采用 BP神經(jīng)網(wǎng)絡(luò)方法預(yù)測(cè)彈道參數(shù)初值,最終推導(dǎo)出基于BP神經(jīng)網(wǎng)絡(luò)和L-M算法的聯(lián)合數(shù)值尋優(yōu)計(jì)算模型。由仿真結(jié)果可知,采用聯(lián)合算法計(jì)算得到的助推段彈道終端高度、速度和彈道傾角偏差分別在2 m、0.1 m/s和0.01°以?xún)?nèi),彈道計(jì)算最大耗時(shí)不超過(guò)3.0 s。結(jié)果表明基于BP神經(jīng)網(wǎng)絡(luò)和L-M方法的聯(lián)合算法對(duì)于機(jī)動(dòng)發(fā)射高超聲速飛行器助推段彈道快速高精度計(jì)算具有較強(qiáng)的魯棒性和適應(yīng)性。本文通過(guò)引入神經(jīng)網(wǎng)絡(luò)算法,利用BP神經(jīng)網(wǎng)絡(luò)為L(zhǎng)-M算法數(shù)值尋優(yōu)計(jì)算提供參數(shù)初值,提高了數(shù)值算法的智能化水平,較好地解決了算法計(jì)算速度的問(wèn)題。
參考文獻(xiàn)(References):
[1]Wu G, Liu L, Wang Y, et al. Ascent trajectory optimization of hypersonic vehicle based on improved Particle Swarm algorithm[C]//IEEE Chinese Automation Congress. 2016: 115-120.
[2]Dalle D J, Torrez S M, Driscoll J F, et al. Minimum-fuel ascent of a hypersonic vehicle using surrogate optimization[J]. Journal of Aircraft, 2014, 51(6): 1973-1986.
[3]Feng L, Liu L, Wang Y. Trajectories optimization of hypersonic vehicle based on a hybrid optimization algorithm of PSO and SQP[C]//IEEE Control and Decision Conference. 2015: 4518-4522.
[4]任京濤. 助推滑翔導(dǎo)彈上升段多終端約束彈道設(shè)計(jì)及制導(dǎo)方法研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué)碩士論文, 2013.Ren J T. Research on reconfigurable control system design of heavy launch vehicle[D]. Harbin: Harbin Institute of Technology Master Degree Dissertation,2013.
[5]崔乃剛, 傅瑜, 盧寶剛. 助推-滑翔飛行器助推段最優(yōu)軌跡快速生成方法[J]. 彈道學(xué)報(bào), 2013,25(2): 33-38.Cui N G, Fu Y, Lu B G. A fast generation method of optional trajectory of boost phase for boost-glide vehicle [J]. Journal of ballistics, 2013, 25(2): 33-38.
[6]劉開(kāi)封, 陳穎, 何念念, 等. 基于正交實(shí)驗(yàn)與牛頓方法的滑翔導(dǎo)彈助推彈道設(shè)計(jì)[J]. 導(dǎo)彈與航天運(yùn)載技術(shù), 2014, 334(4): 50-54.Liu K F, Chen Y, He N N, et al. Boost trajectory design of boost-glide missile based on orthogonal experimental design and Newton's iterative method[J]. Missiles and Space Vehicles, 2014, 334(4): 50-54.
[7]Ma C F, Jiang L H. Some research on Levenberg-Marquardt method for the nonlinear equations[J].Applied Mathematics and Computation, 2007, 184:1032-1040.
[8]鮮勇, 李少朋, 李邦杰. 基于BP神經(jīng)網(wǎng)絡(luò)的固體導(dǎo)彈耗盡關(guān)機(jī)姿態(tài)調(diào)制方法的研究[J]. 兵工學(xué)報(bào),2015, 36 (4): 668-673.Xian Y, Li S P, Li B J. An approach to attitude angle modulation of solid missiles under the condition of depleted shutdown based on BP neural network[J]. Acta Armamentarii, 2015, 36(4): 668-673.
[9]Ran M P, Wang Q. Spacecraft rendezvous trajectory optimization method based on EPSO[J]. Journal of Astronautics, 2013, 34(9): 1195-1201.
[10]Pontani M, Conway B A. Optimal finite-thrust rendezvous trajectories found via particle swarm algorithm[J]. Journal of Spacecraft and Rockets,2014, 94(1): 434-445.
[11]Zhang W, Huang X, Gao X Z, et al. Multi-objective longitudinal trajectory optimization for hypersonic reentry glide vehicle based on PSO algorithm[C]//IEEE Control Conference. 2015: 2350-2356.
[12]崔乃剛, 張亮, 韋常柱, 等. 可重復(fù)使用運(yùn)載器大姿態(tài)機(jī)動(dòng)自抗擾控制[J]. 中國(guó)慣性技術(shù)學(xué)報(bào),2017, 25(3): 387-394.Cui N G, Zhang L, Wei C Z, et al. Active disturbance rejection control for reusable launch vehicle with large attitude maneuver[J]. Journal of Chinese Inertial Technology, 2017, 25(3): 387-394.
[13]Yamashita N, Fukushima M. On the rate of convergence of the Levenberg-Marquardt method[J]. Computing (Supplement 15), 2001: 239-249.