黃祥玉，齊先軍

（合肥工業(yè)大學(xué) 安徽省新能源利用與節(jié)能重點(diǎn)實(shí)驗(yàn)室，安徽 合肥 230009）

0 引言

隨著電網(wǎng)規(guī)模日益擴(kuò)大、互聯(lián)程度不斷提高，電力部門越來越多地運(yùn)用可靠性評(píng)估技術(shù)對(duì)電網(wǎng)規(guī)劃、設(shè)計(jì)和運(yùn)行的各個(gè)階段進(jìn)行定量的可靠性評(píng)估。近年來，電力系統(tǒng)可靠性評(píng)估在模型算法和工程應(yīng)用上均取得了較大進(jìn)展［1］。節(jié)點(diǎn)負(fù)荷具有較強(qiáng)的隨機(jī)性、分散性和多樣性［2］，其模型的準(zhǔn)確與否對(duì)可靠性評(píng)估結(jié)果有較大的影響，如何考慮節(jié)點(diǎn)負(fù)荷不確定性［3］和相關(guān)性對(duì)可靠性評(píng)估結(jié)果的影響，已成為一個(gè)亟待深入研究的重要課題。

在電力系統(tǒng)可靠性評(píng)估中，各節(jié)點(diǎn)負(fù)荷通常是按照完全相關(guān)或完全獨(dú)立［4］來處理的。在這2種情況下都可以通過核密度估計(jì)方法得到節(jié)點(diǎn)負(fù)荷的邊緣概率分布函數(shù)，對(duì)此函數(shù)進(jìn)行抽樣獲得系統(tǒng)負(fù)荷樣本。文獻(xiàn)［5］采用系統(tǒng)年峰荷表達(dá)系統(tǒng)負(fù)荷模型，各節(jié)點(diǎn)負(fù)荷等于系統(tǒng)年峰荷的一個(gè)固定比例，其計(jì)算結(jié)果只能反映系統(tǒng)資源最緊張情況下的抗風(fēng)險(xiǎn)能力；文獻(xiàn)［6］利用多維正態(tài)分布函數(shù)進(jìn)行負(fù)荷相關(guān)和不確定性分析，但是如果采取的假設(shè)與實(shí)際情況不相符，則依據(jù)上述負(fù)荷模型得出的計(jì)算結(jié)果和相關(guān)推論可能導(dǎo)致較大的誤差甚至錯(cuò)誤的結(jié)論。文獻(xiàn)［1］提出了一種基于多元核密度估計(jì)的相關(guān)性負(fù)荷模型，這種方法是一種基于原始數(shù)據(jù)驅(qū)動(dòng)、避免主觀經(jīng)驗(yàn)影響的先進(jìn)數(shù)學(xué)方法；然而多元核密度估計(jì)是一種抽樣效率低下的方法，為了得到比較精確的計(jì)算結(jié)果，必須以犧牲計(jì)算成本為代價(jià)。

本文采用copula函數(shù)建立相關(guān)性負(fù)荷模型，其中以非參數(shù)核密度估計(jì)方法估計(jì)copula函數(shù)參數(shù)。這種方法具有如下優(yōu)點(diǎn)：

（1）基于原始數(shù)據(jù)驅(qū)動(dòng)，可以完全避免主觀經(jīng)驗(yàn)的影響。

（2）不必假設(shè)各節(jié)點(diǎn)負(fù)荷邊緣分布函數(shù)類型和它們的聯(lián)合概率分布函數(shù)類型相同，且各節(jié)點(diǎn)負(fù)荷的邊緣分布函數(shù)類型也不必相同，這使得copula函數(shù)的應(yīng)用范圍較廣。

本文主要介紹了copula函數(shù)參數(shù)選擇、最優(yōu)copula函數(shù)選擇和基于相關(guān)性負(fù)荷模型的節(jié)點(diǎn)負(fù)荷狀態(tài)抽樣，并運(yùn)用非序貫蒙特卡洛方法編程計(jì)算，得到系統(tǒng)可靠性指標(biāo)。

1 copula函數(shù)簡(jiǎn)介

設(shè)隨機(jī)向量X＝（X1，…，XN）的概率分布函數(shù)為F（x1，…，xN），由聯(lián)合分布函數(shù)得到變量的邊緣分布函數(shù)，在聯(lián)合分布中去除邊緣分布的信息，只剩下相關(guān)結(jié)構(gòu)的信息。如果存在函數(shù)C，滿足（1）式，則稱C是聯(lián)合分布函數(shù)F的copula，也稱C為隨機(jī)向量（X1，…，XN）的copula［7］。定義c為對(duì)應(yīng)copula函數(shù)的密度函數(shù)。（1）式為：

其中，θ為參數(shù)向量；Fi（·）（i＝1，…，N）為邊緣分布函數(shù)。

copula函數(shù)實(shí)質(zhì)上是一種將邊緣分布函數(shù)和聯(lián)合分布函數(shù)連接起來的紐帶函數(shù)，又稱連接函數(shù)。應(yīng)用copula函數(shù)實(shí)現(xiàn)相關(guān)性建模一般有如下幾個(gè)步驟：① 預(yù)選幾個(gè)常用的copula函數(shù)；②根據(jù)原始數(shù)據(jù)估計(jì)各常用copula函數(shù)的參數(shù)；③ 從所有預(yù)選的copula函數(shù)中選擇最優(yōu)copula函數(shù)。

2 負(fù)荷相關(guān)性建模

2.1 常用的copula函數(shù)

相關(guān)性分析中常用的copula函數(shù)主要有橢球copula類和阿基米德copula類。橢球copula類主要包括Gaussian copula和t copula（本文不考慮t copula，只做簡(jiǎn)單介紹）；阿基米德copula主要包括Clayton copula、Frank copula和Gumbel copula。

（1）Gaussian copula函數(shù)。表達(dá)式為：

其中，ui＝Fi（xi）（i＝1，2，…，N）；φρ為相關(guān)系數(shù)矩陣ρ的多元標(biāo)準(zhǔn)正態(tài)分布函數(shù)；φ－1為標(biāo)準(zhǔn)正態(tài)分布函數(shù)的逆函數(shù)。

（2）t copula函數(shù)。表達(dá)式為：

其中，ui＝Fi（xi）（i＝1，2，…，N）；Tρ，v為相關(guān)系數(shù)矩陣ρ、自由度v的多元t分布函數(shù)；為自由度是vi的標(biāo)準(zhǔn)t分布函數(shù)的逆函數(shù)。

（3）阿基米德copula函數(shù)。表達(dá)式為：

其中，ui＝Fi（xi）（i＝1，2，…，N）；φ為對(duì)應(yīng)阿基米德copula的生成元；φ－1為生成元的逆函數(shù)。3種常見阿基米德copula的生成元見表1所列，表中α為相依系數(shù)，表征變量間的相依程度。

表1 阿基米德copula函數(shù)生成元

2.2 負(fù)荷相關(guān)性建模中copula函數(shù)參數(shù)選擇

常用的copula函數(shù)模型的參數(shù)估計(jì)方法包括精確極大似然法（Exact Maximum Likelihood method，EML）、邊緣分布推導(dǎo)法（Inference Functions for Margins method，IFM）、規(guī)范極大似然法（Canonical Maximum Likelihood method、CML）、基于核密度的極大似然法（Maximum Likelihood based on Kernel density method，MLK）和Genest and Revest法［8］，為了編程計(jì)算的簡(jiǎn)便和精確，本文采用MLK法計(jì)算copula函數(shù)參數(shù)。

假設(shè)系統(tǒng)中有N個(gè)負(fù)荷節(jié)點(diǎn)，1a中每個(gè)小時(shí)各負(fù)荷節(jié)點(diǎn)的負(fù)荷數(shù)據(jù)已知。第i個(gè)負(fù)荷節(jié)點(diǎn)在第t個(gè)小時(shí)的負(fù)荷為xit（i＝1，2，…，N；t＝1，2，…，T），T為1a小時(shí)時(shí)段數(shù)。

（1）利用核密度估計(jì)函數(shù)將負(fù)荷變量xit進(jìn)行概率變換，得到服從均勻分布的變量uit為：

其中，hi為第i個(gè)負(fù)荷節(jié)點(diǎn)對(duì)應(yīng)窗寬；K（·）為核函數(shù)。因?yàn)樵谧顑?yōu)窗寬的情況下核函數(shù)的選擇對(duì)密度估計(jì)結(jié)果影響不大，所以通常選取核函數(shù)為高斯核。此時(shí)，窗寬可以由hi＝1.059σi求得，其中，σi為第i個(gè)節(jié)點(diǎn)負(fù)荷數(shù)據(jù)標(biāo)準(zhǔn)差，則（5）式可變?yōu)椋?/p>

（2）令L（θ）＝ln（c（θ）），采用極大似然估計(jì)法估計(jì)copula函數(shù)密度函數(shù)c的參數(shù)向量θ為：

2.3 負(fù)荷相關(guān)性建模中最優(yōu)copula函數(shù)的確定

最優(yōu)copula的選擇方式有很多，沒有統(tǒng)一的標(biāo)準(zhǔn)，主要可以分為圖形法和解析法，本文采用更加直觀準(zhǔn)確的解析法。在選擇最優(yōu)copula的解析法中最常用的是歐氏距離法，即基于經(jīng)驗(yàn)copula（Empirical copula，EMC）函數(shù)的最短距離法選取最優(yōu)copula函數(shù)，EMC可表示為：

其中，I（·）為指示函數(shù)，若括號(hào)內(nèi)條件滿足，則I的取值為1，否則為0；x＝｛（，…），t＝1，2，…，T｝為容量j是T的N維觀測(cè)樣本；為順序統(tǒng)計(jì)量且1≤t1，…，tN≤T。

EMC與理論 copula（Theoretical copula，THC）函數(shù)之間的歐氏距離可通過（9）式計(jì)算，即

根據(jù)上述計(jì)算結(jié)果，選擇d最小的copula函數(shù)作為最優(yōu)copula，用來描述多元隨機(jī)變量間的相依結(jié)構(gòu)。

2.4 相關(guān)性負(fù)荷的狀態(tài)抽樣

根據(jù)負(fù)荷的最優(yōu)copula函數(shù)，可以對(duì)負(fù)荷狀態(tài)抽樣。Gaussian copula函數(shù)和阿基米德copula函數(shù)所對(duì)應(yīng)的抽樣方法不同。

（1）Gaussian copula函數(shù)抽樣方法［9］。抽樣步驟如下：① 計(jì)算負(fù)荷隨機(jī)向量的相關(guān)系數(shù)矩陣ρ，對(duì)其進(jìn)行Cholesky分解，得到下三角矩陣A，如果相關(guān)系數(shù)矩陣是正定的，則存在唯一確定的N×N矩陣，使得ρ＝AAT；② 生成由N個(gè)服從標(biāo)準(zhǔn)正態(tài)分布的變量所構(gòu)成的隨機(jī)向量w＝（w1，w2，…，wN），令z＝Aw；③ 令ui＝φ（zi）（i＝1，2，…，N），得到服從N元Gaussian copula函數(shù)分布的隨機(jī)向量取值（u1，u2，…，uN）。

（2）阿基米德copula函數(shù)抽樣方法［10］。設(shè)阿基米德copula函數(shù)為C，其生成元為φ，（U1，U2，…，UN）是服從聯(lián)合分布C的隨機(jī)向量，根據(jù)Marshall和Olkin提出的算法，如果存在一個(gè)分布函數(shù)F，滿足F（0）＝0，且其拉氏變換與生成元的反函數(shù)相等，即φ－1＝L［F］，可以采用如下方法產(chǎn)生（U1，U2，…，UN）的取值（u1，u2，…，uN）：① 生成分布函數(shù)G（·）的隨機(jī)數(shù)z，其中G＝L－1［φ－1］，L－1（φ－1）為φ－1的反拉氏變換；② 生成N個(gè)服從（0，1）均勻獨(dú)立同分布的隨機(jī)數(shù)w1，w2，…，wN；③ 令ui＝G（－ln（wi）/z）（i＝1，…，N），得到服從生成元φ的N元阿基米德copula函數(shù)分布的隨機(jī)向量取值（u1，u2，…，uN）。

由上述仿真方法得到滿足給定相關(guān)結(jié)構(gòu)的隨機(jī)數(shù)向量后，通過（10）式的反變換可以得到服從聯(lián)合概率分布的多元負(fù)荷隨機(jī)數(shù)為：

其中，F(xiàn)i為節(jié)點(diǎn)i上負(fù)荷對(duì)應(yīng)的邊緣分布函數(shù)。

3 發(fā)輸電系統(tǒng)可靠性評(píng)估

發(fā)輸電系統(tǒng)可靠性評(píng)估過程主要包括4個(gè)方面，即確定元件失效模型和負(fù)荷模型、選擇系統(tǒng)狀態(tài)、識(shí)別并分析系統(tǒng)問題以及進(jìn)行可靠性指標(biāo)計(jì)算［11］。本文選擇非序貫蒙特卡洛方法計(jì)算發(fā)輸電系統(tǒng)可靠性指標(biāo)，具體流程如圖1所示。選擇系統(tǒng)狀態(tài)時(shí)，在元件方面，發(fā)電機(jī)、輸電線路和變壓器仍采用兩狀態(tài)離散分布模型進(jìn)行抽樣；在負(fù)荷方面，基于copula函數(shù)建立計(jì)及節(jié)點(diǎn)負(fù)荷相關(guān)性的負(fù)荷模型，得到節(jié)點(diǎn)負(fù)荷的聯(lián)合概率分布函數(shù)，即最優(yōu)copula函數(shù)，再按最優(yōu)copula函數(shù)對(duì)負(fù)荷進(jìn)行抽樣。

圖1 基于相關(guān)性負(fù)荷模型的發(fā)輸電系統(tǒng)可靠性評(píng)估流程

4 算例

針對(duì)IEEE－RBTS系統(tǒng)建立相關(guān)性負(fù)荷模型進(jìn)行測(cè)試。IEEE－RBTS系統(tǒng)有5個(gè)負(fù)荷節(jié)點(diǎn)，在進(jìn)行電力系統(tǒng)可靠性評(píng)估時(shí)，節(jié)點(diǎn)負(fù)荷數(shù)據(jù)通常按完全相關(guān)處理，但在實(shí)際情況中各節(jié)點(diǎn)的年負(fù)荷峰值不會(huì)在同一時(shí)刻出現(xiàn)，各節(jié)點(diǎn)負(fù)荷的變化趨勢(shì)也不完全一致。為了驗(yàn)證本文所提相關(guān)性負(fù)荷建模方法的正確性，對(duì)IEEE－RBTS系統(tǒng)原始負(fù)荷數(shù)據(jù)進(jìn)行如下修改：①2～6號(hào)節(jié)點(diǎn)分別采用黑色冶金、紡織、機(jī)械、農(nóng)業(yè)和市政生活的典型日負(fù)荷曲線［12］；② 因?yàn)镮EEE－RBTS系統(tǒng)將一年分成夏、冬、春/秋3個(gè)時(shí)間段，而文獻(xiàn)［12］只有冬、夏時(shí)段的典型日負(fù)荷曲線，故假定春/秋時(shí)段的日負(fù)荷曲線可以由冬、夏時(shí)段的取平均值得到；③ 因?yàn)槲墨I(xiàn)［12］沒有每季節(jié)周末的日負(fù)荷曲線，所以假設(shè)周末也采用工作日日負(fù)荷曲線。

因?yàn)?號(hào)節(jié)點(diǎn)基礎(chǔ)負(fù)荷最大，對(duì)系統(tǒng)EENS指標(biāo)影響最大，而2、3號(hào)負(fù)荷節(jié)點(diǎn)地理位置較近，氣候條件相似，負(fù)荷具有較大的相關(guān)性，故對(duì)2、3號(hào)節(jié)點(diǎn)負(fù)荷采用基于copula函數(shù)的相關(guān)性負(fù)荷模型，其他節(jié)點(diǎn)負(fù)荷采用獨(dú)立負(fù)荷模型。

5個(gè)節(jié)點(diǎn)負(fù)荷對(duì)應(yīng)的核密度估計(jì)函數(shù)窗寬見表2所列。

表2 各節(jié)點(diǎn)負(fù)荷數(shù)據(jù)對(duì)應(yīng)窗寬

2、3號(hào)節(jié)點(diǎn)負(fù)荷原始數(shù)據(jù)所對(duì)應(yīng)的各預(yù)選copula函數(shù)參數(shù)和歐氏距離d見表3所列。

表3 copula函數(shù)參數(shù)和歐氏距離

從表3中歐氏距離大小可知，2、3號(hào)節(jié)點(diǎn)負(fù)荷數(shù)據(jù)對(duì)應(yīng)的最優(yōu)copula函數(shù)為Gaussian copula。為了對(duì)比負(fù)荷原始數(shù)據(jù)與copula函數(shù)的分布特性，將2、3號(hào)節(jié)點(diǎn)負(fù)荷序列 ｛（x2t，x3t），t＝1，2，…，T｝按照（5）式進(jìn)行變換，得到變換序列｛（u2t，u3t），t＝1，2，…，T｝。變換序列的散點(diǎn)圖、相同相關(guān)系數(shù)下的Gaussian copula函數(shù)模擬散點(diǎn)圖、Gaussian copula密度函數(shù)圖分別如圖2～圖4所示。

圖2 變換序列散點(diǎn)圖

圖3 Gaussian copula函數(shù)模擬散點(diǎn)圖

圖4 Gaussian copula密度函數(shù)圖

由圖2可以看出，2、3號(hào)節(jié)點(diǎn)負(fù)荷數(shù)據(jù)具有較強(qiáng)對(duì)稱相關(guān)性，中間分布較為均勻，上、下尾分布較為集中，這與圖3中的模擬數(shù)據(jù)散點(diǎn)圖基本一致，圖4中的Gaussian copula密度函數(shù)也具有這些特點(diǎn)，所以選擇Gaussian copula作為最優(yōu)copula是合適的。

利用Gaussian copula函數(shù)抽樣方法產(chǎn)生2、3號(hào)節(jié)點(diǎn)負(fù)荷的抽樣樣本，抽樣200萬次得到模擬數(shù)據(jù)的統(tǒng)計(jì)量見表4所列。

表4 負(fù)荷抽樣數(shù)據(jù)統(tǒng)計(jì)量

此時(shí)，2、3號(hào)節(jié)點(diǎn)的負(fù)荷抽樣樣本的變異系數(shù)［13］分別為1.214×10－4和1.609×10－4。如果要進(jìn)一步提高計(jì)算精度，可以增加抽樣次數(shù)或采用減小方差技術(shù)。

2、3號(hào)節(jié)點(diǎn)負(fù)荷原始數(shù)據(jù)的統(tǒng)計(jì)量見表5所列。比較表4和表5，可以看出模擬數(shù)據(jù)和原始數(shù)據(jù)的統(tǒng)計(jì)量非常接近，說明了基于copula函數(shù)的相關(guān)性負(fù)荷模型和算法的有效性。

表5 負(fù)荷原始數(shù)據(jù)統(tǒng)計(jì)量

基于copula函數(shù)負(fù)荷模型分析負(fù)荷相關(guān)性對(duì)IEEE－RBTS系統(tǒng)可靠性指標(biāo)的影響。采用非序貫蒙特卡洛方法計(jì)算2種情形（計(jì)及節(jié)點(diǎn)負(fù)荷相關(guān)性和節(jié)點(diǎn)負(fù)荷完全獨(dú)立）下的發(fā)輸電系統(tǒng)可靠性指標(biāo)，見表6所列。

表6 發(fā)輸電系統(tǒng)可靠性指標(biāo)

由表6中計(jì)算結(jié)果可以看出系統(tǒng)LOLP指標(biāo)和LOLE指標(biāo)在考慮2、3號(hào)節(jié)點(diǎn)的負(fù)荷相關(guān)性和完全獨(dú)立2種情況下是一樣的，即是否考慮2、3號(hào)節(jié)點(diǎn)的負(fù)荷相關(guān)性對(duì)系統(tǒng)LOLP指標(biāo)和LOLE指標(biāo)沒有影響。而系統(tǒng)EPNS指標(biāo)和EENS指標(biāo)卻有明顯不同，負(fù)荷相關(guān)抽樣情況下系統(tǒng)EENS指標(biāo)比負(fù)荷獨(dú)立抽樣情況大4.326%，這與文獻(xiàn)［1］的結(jié)論基本一致。雖然在本文算例中系統(tǒng)EENS指標(biāo)在2種情況下差別不大，但是若計(jì)算更大規(guī)模電網(wǎng)的可靠性，考慮和不考慮負(fù)荷相關(guān)性所對(duì)應(yīng)的EENS指標(biāo)可能會(huì)相差很大，此時(shí)不考慮負(fù)荷相關(guān)性所造成的誤差不能忽略不計(jì)。而且考慮負(fù)荷相關(guān)性情況下的可靠性指標(biāo)更加符合實(shí)際情況。

5 結(jié)論

本文通過copula函數(shù)建立相關(guān)性負(fù)荷模型，得到基于相關(guān)性負(fù)荷模型和傳統(tǒng)不考慮負(fù)荷相關(guān)性的發(fā)輸電系統(tǒng)可靠性指標(biāo)，并將2種情況下的可靠性指標(biāo)進(jìn)行比較，驗(yàn)證了基于copula函數(shù)的相關(guān)性負(fù)荷模型的正確性。

（1）通過基于copula函數(shù)的相關(guān)性負(fù)荷模型抽樣得到的節(jié)點(diǎn)負(fù)荷數(shù)據(jù)能夠很好地保持原始數(shù)據(jù)的統(tǒng)計(jì)特性，是對(duì)負(fù)荷相關(guān)性比較適合的一種模擬方式。

（2）基于相關(guān)性負(fù)荷模型的發(fā)輸電系統(tǒng)可靠性評(píng)估能夠得到更加準(zhǔn)確的可靠性指標(biāo)。傳統(tǒng)的基于獨(dú)立負(fù)荷模型的發(fā)輸電系統(tǒng)可靠性指標(biāo)相對(duì)于本文基于相關(guān)性負(fù)荷模型的可靠性指標(biāo)會(huì)過高評(píng)估系統(tǒng)的可靠性。

在考慮多元負(fù)荷相關(guān)性的情況下，相關(guān)性分析中常用的阿基米德copula函數(shù)參數(shù)將會(huì)變得十分復(fù)雜，現(xiàn)有條件難以實(shí)現(xiàn)。故基于copula函數(shù)的多元負(fù)荷相關(guān)性情況還需要進(jìn)一步的研究。

［1］張夏菲.非參數(shù)核密度估計(jì)負(fù)荷模型在電網(wǎng)可靠性評(píng)估中的應(yīng)用［D］.重慶：重慶大學(xué)，2010.

［2］周 文，賀仁睦.電力負(fù)荷建模問題研究綜述［J］.現(xiàn)代電力，1999，16（2）：83－89.

［3］史新紅，齊先軍，王治國(guó)，等.基于UGF的發(fā)電系統(tǒng)區(qū)間可靠性評(píng)估及其仿射算法改進(jìn)［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2014，37（3）：286－291.

［4］丁 明，李生虎，洪 梅，等.電力系統(tǒng)概率分析中的K均值聚類負(fù)荷模型［J］.電力系統(tǒng)自動(dòng)化，1999，23（19）：51－54.

［5］Subcommittee P M.IEEE reliability test system［J］.IEEE Transactions on Power Apparatus and Systems，1979，6（3）：2047－2054.

［6］Li W Y，Billinton R.Effect of bus load uncertainty and cor－relation in composite system adequacy evaluation［J］.IEEE Transactions on Power Systems，1991，6（4）：1522－1529.

［7］羅俊鵬.copula理論及其在金融分析中的應(yīng)用研究［D］.天津：天津大學(xué)，2005.

［8］趙崇新.基于多元copula－SV－Var模型的開放式基金投資組合風(fēng)險(xiǎn)測(cè)度［D］.江門：五邑大學(xué)，2012.

［9］李玉敦.計(jì)及相關(guān)性的風(fēng)速模型及其在發(fā)電系統(tǒng)可靠性評(píng)估中的應(yīng)用［D］.重慶：重慶大學(xué)，2012.

［10］王麗芳.copula分布估計(jì)算法［M］.北京：機(jī)械工業(yè)出版社，2012：67－68.

［11］李文沅.電力系統(tǒng)風(fēng)險(xiǎn)評(píng)估模型、方法和應(yīng)用［M］.北京：科學(xué)出版社，2006：97－98.

［12］電力工業(yè)部電力規(guī)劃設(shè)計(jì)總院.電力系統(tǒng)設(shè)計(jì)手冊(cè)［M］.北京：中國(guó)電力出版社，1998：28－29.

［13］李生虎.電力系統(tǒng)可靠性評(píng)估的算法研究及軟件設(shè)計(jì)［D］.合肥：合肥工業(yè)大學(xué)，1999.