一種暫態(tài)功角穩(wěn)定故障損失估計方法

劉 哲，賈宏杰

(天津大學(xué)智能電網(wǎng)教育部重點實驗室，天津 300072)

一種暫態(tài)功角穩(wěn)定故障損失估計方法

劉 哲，賈宏杰

(天津大學(xué)智能電網(wǎng)教育部重點實驗室，天津 300072)

在進行電力系統(tǒng)暫態(tài)穩(wěn)定風(fēng)險評估時，科學(xué)估算故障引起的系統(tǒng)損失是其核心問題．為此，將故障引起的損失看作將系統(tǒng)從不穩(wěn)定狀態(tài)轉(zhuǎn)化為穩(wěn)定狀態(tài)所需要的最小控制代價，借助軌跡靈敏度技術(shù)，將上述問題轉(zhuǎn)化為一個考慮暫態(tài)穩(wěn)定約束的最優(yōu)潮流問題實現(xiàn)求解．推導(dǎo)了發(fā)電機相對功角同發(fā)電機及負荷功率之間的靈敏度系數(shù)；利用靈敏度系數(shù)和發(fā)電機相對功角形成系統(tǒng)的暫態(tài)穩(wěn)定約束，將之引入最優(yōu)潮流實現(xiàn)對故障最小控制代價的求解；并且利用IEEE-14、New England-39節(jié)點等系統(tǒng)驗證了所提方法的正確性和有效性．

暫態(tài)穩(wěn)定；風(fēng)險評估；軌跡靈敏度；最優(yōu)潮流

隨著我國電力系統(tǒng)的發(fā)展，聯(lián)網(wǎng)規(guī)模的不斷擴大，系統(tǒng)中的不確定性因素不斷增加，通過風(fēng)險分析，合理評估不確定性因素對系統(tǒng)穩(wěn)定運行的影響意義重大[1]．風(fēng)險分析的難點是如何準(zhǔn)確估算故障引發(fā)系統(tǒng)失穩(wěn)后所造成的損失，利用最優(yōu)穩(wěn)定控制代價來近似估計系統(tǒng)損失，為該問題的求解提供了一種可行的途徑[1]．文獻[2]采取預(yù)先設(shè)定低壓減載策略對系統(tǒng)進行緊急控制，以獲得控制代價；文獻[3]采取發(fā)電機維修、啟動和出力調(diào)整費用以及負荷緊急控制費用來代替系統(tǒng)失穩(wěn)造成的損失；文獻[4]則利用動態(tài)安全域?qū)ο到y(tǒng)進行預(yù)防控制和緊急控制，從而得到系統(tǒng)動態(tài)不安全的風(fēng)險代價．上述方法都需針對故障預(yù)設(shè)控制策略，在實際應(yīng)用時會存在一定局限性．

軌跡靈敏度能夠在暫態(tài)過程中的每一個時刻提供系統(tǒng)豐富的運行信息，且不受元件模型限制[5]，近年來在考慮暫態(tài)穩(wěn)定約束的最優(yōu)潮流(OPF with transient stability constraints，OTS)[6-7]和暫態(tài)穩(wěn)定控制優(yōu)化(包括預(yù)防控制[8-11]和緊急控制[12])等方面得到成功應(yīng)用．筆者將故障引起系統(tǒng)失穩(wěn)后的損失，視為將系統(tǒng)從不穩(wěn)定狀態(tài)轉(zhuǎn)化為穩(wěn)定狀態(tài)所需最小控制代價，將之轉(zhuǎn)化為一個OTS問題加以求解．故障中的暫態(tài)穩(wěn)定約束通過發(fā)電機相對功角加以考慮，后者則利用軌跡靈敏度近似加以估計[13]；發(fā)電機出力調(diào)整和負荷裁減作為控制手段，不同的控制手段具有不同的控制代價(成本)．當(dāng)發(fā)電機備用容量充足時，優(yōu)先調(diào)整發(fā)電機出力(控制成本較低)來維持系統(tǒng)穩(wěn)定；而當(dāng)發(fā)電機備用容量不足時，則通過控制成本更高的負荷裁減措施來保證系統(tǒng)穩(wěn)定，其優(yōu)化過程更貼合實際．

1 軌跡靈敏度

1.1 基本原理

電力系統(tǒng)軌跡靈敏度方法最早由Pai等[5]提出，它認為暫態(tài)過程中的系統(tǒng)狀態(tài)變量同系統(tǒng)參數(shù)和狀態(tài)變量初值近似滿足線性關(guān)系，而線性系數(shù)即為系統(tǒng)的軌跡靈敏度．

電力系統(tǒng)動態(tài)可用微分-代數(shù)方程表示為

式中：x為系統(tǒng)狀態(tài)變量；y為系統(tǒng)代數(shù)變量；λ為系統(tǒng)控制變量；x0、y0分別表示狀態(tài)變量和代數(shù)變量的初值．將式(1)對λ 求偏導(dǎo)可得

式中：xλ、yλ分別為系統(tǒng)軌跡對控制變量λ的靈敏度矩陣，分別為微分方程組對狀態(tài)變量、代數(shù)變量和控制變量的偏導(dǎo)；分別為代數(shù)方程組對狀態(tài)變量、代數(shù)變量和控制變量的偏導(dǎo)數(shù)．根據(jù)系統(tǒng)初始運行點可求得軌跡靈敏度的初試條件，將式(1)和式(2)聯(lián)立求解即可得到系統(tǒng)各時刻的軌跡靈敏度．

1.2 發(fā)電機功角同其有功出力間的靈敏度

當(dāng)發(fā)電機采取經(jīng)典模型時，電力系統(tǒng)動態(tài)方程為

式中：δi、ωi、Mi、Pg,i、Ei′分別為發(fā)電機轉(zhuǎn)子角、轉(zhuǎn)速、轉(zhuǎn)動慣量、機械功率和內(nèi)電勢；Gii、Gij為僅保留發(fā)電機內(nèi)節(jié)點后，系統(tǒng)導(dǎo)納矩陣元素的實部；分別為為系統(tǒng)導(dǎo)納矩陣元素的虛部；m為發(fā)電機節(jié)點數(shù)．

對式(3)兩端分別求發(fā)電機k有功出力Pg,k的偏導(dǎo)數(shù)可得

由發(fā)電機機端功率方程可得

式中：iI˙為發(fā)電機注入網(wǎng)絡(luò)的電流；iS、iα分別為發(fā)電機輸出的視在功率和功率因數(shù)角；Vi、θi分別為發(fā)電機機端電壓幅值和相角．

發(fā)電機內(nèi)電勢iE′˙可表示為

因此，有

式(9)兩端分別對Pg,k求偏導(dǎo)數(shù)，得

式(7)兩端分別對Pg,k求偏導(dǎo)數(shù)，得

發(fā)電機節(jié)點無功注入方程為

式中n為系統(tǒng)節(jié)點數(shù)．

式(12)兩端分別對Pg,k求偏導(dǎo)數(shù)得

又知穩(wěn)態(tài)時的發(fā)電機機械功率為

式(15)兩端分別對Pg,k求偏導(dǎo)數(shù)，可得

1.3 發(fā)電機功角同負荷有功間的靈敏度

假設(shè)負荷采取恒阻抗模型，式(3)兩端分別對負荷節(jié)點k的有功負荷PL,k求偏導(dǎo)數(shù)，得

式(18)兩端分別對PL,k求偏導(dǎo)數(shù)，得

對于恒阻抗負荷而言，等值導(dǎo)納可表示為

當(dāng)發(fā)電機采用較復(fù)雜模型時，依然可以推導(dǎo)出發(fā)電機功角同發(fā)電機有功出力、負荷有功和無功間的靈敏度關(guān)系[14]．但對于復(fù)雜模型，為便于計算，可利用數(shù)值方法進行軌跡靈敏度的求解[13]．

當(dāng)系統(tǒng)控制參數(shù)發(fā)生變化時，利用已求得的發(fā)電機絕對功角的軌跡靈敏度，可近似估計發(fā)電機ij之間相對功角的變化量，即

2 基于靈敏度的電力系統(tǒng)暫態(tài)穩(wěn)定最小控制成本模型

2.1 暫態(tài)穩(wěn)定最小控制成本模型

將軌跡靈敏度和最優(yōu)潮流相結(jié)合，求解非線性優(yōu)化問題，即可得暫態(tài)穩(wěn)定最小控制成本，其模型為

約束條件為

上述模型中式(24)表示針對某故障的控制成本，包括上調(diào)發(fā)電機出力的成本、下調(diào)發(fā)電機出力成本和裁減負荷成本．為簡單起見，本模型將3種控制代價均考慮為控制量的線性函數(shù)，進一步亦可采用二次、線性分段成本函數(shù)[15]，其求解過程與上述模型類似；式(25)和式(26)分別為有功和無功潮流方程；式(27)為線路潮流約束；式(28)為節(jié)點電壓約束；式(29)、式(30)和式(31)分別為發(fā)電機有功出力上調(diào)量約束、下調(diào)量約束和無功約束；式(32)表示負荷裁減量約束；式(33)表示暫態(tài)功角約束．

為簡單起見，對于負荷節(jié)點，在進行負荷有功裁減時，負荷無功裁減量計算式為即在控制過程中，負荷功率因數(shù)考慮為恒定，若負荷模型為其他形式，將其考慮成最優(yōu)控制模型的約束條件即可同理求解，式(34)中θ表示負荷的功率因數(shù)角向量．

式(24)～式(33)屬于典型的含非線性約束最優(yōu)化模型，求解過程較為簡單，求解的流程如圖1所示，利為收斂判據(jù)，當(dāng)2個條件同時滿足時，表明系統(tǒng)已無可調(diào)發(fā)電機和負荷，則程序結(jié)束．

圖1 功角穩(wěn)定控制流程Fig.1 Flow chart of angle stability control

2.2 優(yōu)化模型的幾點說明

(1) 對于系統(tǒng)中的平衡機，只有在其他發(fā)電機和負荷調(diào)整不能滿足系統(tǒng)潮流平衡時，才參與潮流調(diào)整，為達此目的，將平衡機的有功出力調(diào)整費用取為普通發(fā)電機調(diào)整費用的2倍．

(2) 模型中的式(33)為系統(tǒng)的暫態(tài)穩(wěn)定約束，當(dāng)發(fā)電機相對功角超過δlim時即認為系統(tǒng)發(fā)生暫態(tài)失穩(wěn)．不妨假設(shè)系統(tǒng)中有m個發(fā)電機和l個負荷，則有Cm2個發(fā)電機相對功角，其中有c個相對功角絕對值大于δlim，進而有

(3) 用圖2簡單解釋式(33)中tsen時刻的選取原則，該圖給出了一個3機簡單系統(tǒng)的相對功角：2號機和1、3號機的功角相互擺開，圖中只畫出了δ23和δ13，曲線δ12和δ23非常近似，但不重合．在t＞后，δ12、δ23絕對值將大于δlim，而相對功角δ13的絕對值一直小于δlim．假設(shè)在δ12、δ23曲線中，δ12最先達到δsen，此時刻記為，為了保證在時刻之前所有相對功角都位于[?δlim,δlim]內(nèi)，取tsen=．δsen和δlim是預(yù)先設(shè)定的一個相對功角的閾值，不同系統(tǒng)閾值的選取可能不同[10]，δlim越小表明穩(wěn)定區(qū)域越小，暫態(tài)穩(wěn)定最優(yōu)控制成本將越大，考慮到功角穩(wěn)定問題大多是單擺失穩(wěn)，因此δsen一般選取在2π到10π之間即可．本文暫不考慮預(yù)設(shè)閾值對算法的影響，對同一系統(tǒng)δsen和δlim取值相同．

圖2 發(fā)電機相對功角曲線示意Fig.2 Sketch of relative rotor angle curves

(4) 在式(33)中，有

式中

靈敏度矩陣的構(gòu)成式為

由于負荷按照恒功率因數(shù)進行處理，因此將發(fā)電機相對功角同無功負荷的靈敏度按照相應(yīng)的功率因數(shù)折合到有功負荷靈敏度中表示發(fā)電機i和j之間的相對功角同負荷q間的靈敏度，且有

式中qθ為負荷q的功率因數(shù)角．

3 算例分析

以IEEE-14和New England-39節(jié)點系統(tǒng)為例，對本文中方法進行驗證，計算中取δlim＝π[10]，ε=0.05MW．對于14節(jié)點系統(tǒng)，發(fā)電機采用次暫態(tài)模型，負荷采取恒阻抗負荷，δsen＝5π；對于39節(jié)點系統(tǒng)，發(fā)電機采用經(jīng)典模型，負荷采用恒阻抗模型，δsen＝10π．除平衡節(jié)點以外其他發(fā)電機控制費用：=10 K/MW，平衡節(jié)點的控制費用：=＝20 K/MW，負荷控制費用：CL＝30 K/MW，K為某種貨幣單位．

3.1 IEEE-14節(jié)點系統(tǒng)

假設(shè)線路2-4靠近節(jié)點2處發(fā)生三相短路，臨界切除時間tcr＝0.29 s，設(shè)此故障由于控制設(shè)備原因，其實際切除時間為：tcl＝1.5tcr＝0.435 s，節(jié)點電壓最大值均取1.15 p.u.，節(jié)點電壓最小值均取0.85 p.u.，此時系統(tǒng)出現(xiàn)如圖3所示的失穩(wěn)情況．

圖3 發(fā)電機功角曲線(IEEE-14)Fig.3 Rotor angle curves(IEEE-14)

從圖3中可以清晰地看到發(fā)電機被分成了{1，2}和{3，4，5}2組．

情景1 當(dāng)發(fā)電機有功備用充足時，不需要調(diào)整負荷，只通過調(diào)整發(fā)電機出力即可實現(xiàn)系統(tǒng)穩(wěn)定，控制結(jié)果如表1所示，最小控制代價為1,225.10,K．調(diào)整后各發(fā)電機功角曲線如圖4所示，從圖4可以看出，經(jīng)過發(fā)電機調(diào)整系統(tǒng)可保持暫態(tài)穩(wěn)定．

表1 情景1的發(fā)電機控制結(jié)果(IEEE-14)Tab.1 Results of generator control for scenario 1(IEEE-14)

圖4 情景1調(diào)整后發(fā)電機功角曲線(IEEE-14)Fig.4 Rotor angle curves after control for scenario 1 (IEEE-14)

情景2 當(dāng)發(fā)電機備用不足時，僅靠調(diào)整發(fā)電機，將不能達到系統(tǒng)穩(wěn)定的目的，此時還應(yīng)伴隨負荷裁減控制．作為示例，不妨假設(shè)G3、G4、G5號發(fā)電機有功出力的最大值由原先的100,MW變成20,MW，通過優(yōu)化計算，發(fā)電機和負荷的調(diào)整結(jié)果如表2和表3所示．此時，發(fā)電機和負荷總的控制費用為3,294.67,K，調(diào)整后發(fā)電機的功角曲線如圖5所示，可以看出系統(tǒng)亦能保持穩(wěn)定，但由于發(fā)電機有功儲備較小，需要進行負荷裁減，情景2的控制代價要高于情景1．

3.2 New England-39節(jié)點系統(tǒng)

假設(shè)線路22-21靠近節(jié)點22處發(fā)生三相短路，臨界切除時間tcr＝0.13 s，因某種原因，系統(tǒng)實際故障切除時間變?yōu)椋簍cl＝1.5tcr＝0.2,s，節(jié)點電壓最大值取1.15 p.u.，節(jié)點電壓最小值取0.85,p.u.，此時系統(tǒng)出現(xiàn)如圖6所示的失穩(wěn)情況，從圖6中可以看出10臺發(fā)電機被分為{1，3，4，5，8，9}、{2}、{6，7}、{10}4組．

表2 情景2的發(fā)電機控制結(jié)果(IEEE-14)Tab.2 Results of generator control for scenario 2(IEEE-14)

表3 情景2的負荷控制結(jié)果(IEEE-14)Tab.3 Results of load control for scenario 2(IEEE-14)

圖5 情景2調(diào)整后發(fā)電機功角曲線(IEEE-14)Fig.5 Rotor angle curves after control for scenario 2 (IEEE-14)

圖6 發(fā)電機功角曲線(New England-39)Fig.6 Rotor angle curves(New England-39)

情景1 當(dāng)發(fā)電機備用充足時，發(fā)電機調(diào)整結(jié)果如表4所示．由于發(fā)電機均具有充足的有功調(diào)解備用，因此故障后，通過調(diào)整發(fā)電機出力即可保證系統(tǒng)恢復(fù)穩(wěn)定運行，系統(tǒng)最小控制成本為9,564.39,K，控制后系統(tǒng)功角曲線如圖7所示．

情景2 假設(shè)10號發(fā)電機有功出力最大值由1,600,MW變?yōu)?,400,MW，系統(tǒng)將沒有足夠的發(fā)電機備用容量可調(diào)，因此同樣的故障出現(xiàn)后，單純依靠發(fā)電機調(diào)整無法保證系統(tǒng)穩(wěn)定，必須輔以負荷裁減，此時系統(tǒng)的控制方案和最小控制成本示于表5和表6．不難看出，由于發(fā)電機缺乏足夠的有功備用，系統(tǒng)的故障控制成本升高變?yōu)?1,128.62,K，控制后的系統(tǒng)功角曲線示于圖8．

圖7 情景1調(diào)整后的發(fā)電機功角曲線(New England-39)Fig.7 Rotor angle curves after control for scenario 1 (New England-39)

表4 情景1的發(fā)電機控制結(jié)果(New England-39)Tab.4 Results of generator control for scenario 1(New England-39)

表5 情景2發(fā)電機控制結(jié)果(New England-39)Tab.5 Results of generator control for scenario 2(New England-39)

表6 情景2負荷控制結(jié)果(New England-39)Tab.6 Results of load control for scenario 2(New England-39)

圖8 情景2調(diào)整后發(fā)電機功角曲線(New England-39)Fig.8 Rotor angle curves after control for scenario 2 (New England-39)

4 結(jié) 語

本文中提出了一種基于軌跡靈敏度的暫態(tài)功角穩(wěn)定故障損失估計方法，將故障損失用將系統(tǒng)從不穩(wěn)定狀態(tài)轉(zhuǎn)化為穩(wěn)定狀態(tài)所需要最小控制代價加以考慮．利用軌跡靈敏度技術(shù)近似估計系統(tǒng)參數(shù)變化后發(fā)電機相對功角的變化量，進而將暫態(tài)穩(wěn)定約束引入最優(yōu)潮流模型，以發(fā)電機調(diào)整和負荷裁減作為控制手段，實現(xiàn)最小控制代價的求解．本文尚未考慮故障隨機性的影響，但經(jīng)過大量算例驗證(包括不同的故障線路、不同的故障類型和切除時間等)，算法對所有故障具有普遍適用性，這為進一步開展電力系統(tǒng)暫態(tài)風(fēng)險評估研究奠定了基礎(chǔ)．

［1］ 薛禹勝，劉 強，董朝陽，等. 關(guān)于暫態(tài)穩(wěn)定不確定性分析的評述[J]. 電力系統(tǒng)自動化，2007，31(14)：1-6.

Xue Yusheng，Liu Qiang，Dong Zhaoyang，et al. A review of non-deterministic analysis for power system transient stability[J]. Automation of Electric Power Systems，2007，31(14)：1-6(in Chinese).

［2］ 趙珊珊，仲悟之，張東霞，等. 暫態(tài)電壓穩(wěn)定風(fēng)險評估方法及應(yīng)用[J]. 電力系統(tǒng)自動化，2009，33(19)：1-4.

Zhao Shanshan，Zhong Wuzhi，Zhang Dongxia，et al. A transient voltage instability risk assessment method and its application[J]. Automation of Electric Power Systems，2009，33(19)：1-4(in Chinese).

［3］ Vittal V，McCalley J D，Van Acker V，et al. Transient instability risk assessment[C]//Proceedings of Power Engineering Society Summer Meeting. Edmonton，Canada，1999：206-211.

［4］ 余貽鑫，王東濤. 輸電系統(tǒng)動態(tài)安全風(fēng)險評估與優(yōu)化[J]. 中國科學(xué) E輯：技術(shù)科學(xué)，2009，39(2)：286-292.

Yu Yixin，Wang Dongtao. Dynamic security risk assessment and optimization of power transmission system[J]. Science in China Series E：Technological Sci-ences，2009，39(2)：286-292(in Chinese).

［5］ Hiskens I A，Pai M A，Trajectory sensitivity analysis of hybrid systems[J]. IEEE Transactions on Circuits and Systems，2000，47(2)：204-220.

［6］ 孫景強，房大中，錘德成. 暫態(tài)穩(wěn)定約束下的最優(yōu)潮流[J]. 中國電機工程學(xué)報，2005，25(12)：12-17.

Sun Jingqiang，F(xiàn)ang Dazhong，Chui Decheng. Optimal power flow with transient stability constraints [J]. Proceedings of the CSEE，2005，25(12)：12-17(in Chinese).

［7］ 李傳棟，楊金剛，陳家榮. 暫態(tài)穩(wěn)定約束下的最優(yōu)潮流[J]. 電力系統(tǒng)及其自動化學(xué)報，2009，21(6)：45-50.

Li Chuandong，Yang Jingang，Chen Jiarong. Optimal power flow with transient stability constraints[J]. Proceedings of the CSU-EPSA，2009，21(6)：45-50(in Chinese).

［8］ Nguyen T B，Pai M A. Dynamic security-constrained rescheduling of power systems using trajectory sensitivities[J]. IEEE Transactions on Power Systems，2003，18(2)：848-854.

［9］ 孫景強，房大中，周保容. 基于軌跡靈敏度的電力系統(tǒng)動態(tài)安全預(yù)防控制算法研究[J]. 電網(wǎng)技術(shù)，2004，28(21)：26-30.

Sun Jingqiang，F(xiàn)ang Dazhong，Zhou Baorong. Study on preventive control algorithm for dynamic security of power systems based on trajectory sensitivity method[J]. Power System Technology，2004，28(21)：26-30(in Chinese).

［10］ Fang Dazhong，Yang Xiaodong，Sun Jingqiang，et al. An optimal generation rescheduling approach for transient stability enhancement[J]. IEEE Transactions on Power Systems，2007，22(1)：386-394.

［11］ 周鯤鵬，陳允平. 運用軌跡靈敏度的電力系統(tǒng)動態(tài)安全控制[J]. 電網(wǎng)技術(shù)，2003，27(12)：46-50.

Zhou Kunpeng，Chen Yunping. Power system dynamic security control based on trajectory sensitivity[J]. Power System Technology，2003，27(12)：46-50(in Chinese).

［12］ Zima M，Andersson G. Stability assessment and emergency control method using trajectory sensitivities[C]//Proceedings of IEEE Bologna Power Tech Conference. Italy，2003：189-194.

［13］ Nguyen T T，Karimishad A. Fast and accurate method for dynamic security-constrained economic dispatch via sensitivity analysis[C]// Proceedings of IEEE Power Engineer-ing Society General Meeting． Tampa，F(xiàn)L，USA，2007：1-8.

［14］ Nguyen T B，Pai M A，Hiskens I A. Direct computation of critical clearing time using trajectory sensitivities[C]//Proc of IEEE Power Engineering Society Summer Meeting. Seattle，WA，USA，2000：604-608.

［15］ Wang Hongye，Murillo-Sanchez C E，Zimmeraman R D，et al. On computational issues of market-based optimal power flow[J]. IEEE Transactions on Power Systems，2007，22(3)：1185-1193.

An Evaluation Method for Transient Stability Contingency Cost

LIU Zhe，JIA Hong-jie
(Key Laboratory of Smart Grid of Ministry of Education，Tianjin University，Tianjin 300072，China)

In power system transient stability risk assessment，it is a core issue to properly estimate the cost consequences caused by instability events. In this paper，the cost consequences of a contingency are considered as the minimum control cost to move the system from instability to stability under the contingency. Trajectory sensitivity is employed to properly consider the transient stability constraints in the optimal power flow(OPF).The minimum objective of the OPF is the determination of contingency cost. The sensitivities of the relative rotor angle with respect to the output of generator and the dispatchable load are firstly derived. The information is then used to imitate the transient stability constraints in the OPF. The correctness and effectiveness of the proposed method in this paper is finally validated by some systems，such as IEEE 14-bus system and New England 39-bus system.

transient stability；risk assessment；trajectory sensitivity；optimal power flow

TM712

A

0493-2137(2011)09-0759-09

2010-09-20；

2010-11-18．

國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)資助項目(2009CB219701，2010CB234600)；天津市科技發(fā)展計劃資助項目(09JCZDJC 25000)；“十一五”國家科技支撐計劃重大資助項目(2006BAJ03A06)；博士學(xué)科點專項科研基金資助項目.

劉 哲（1985— ），男，博士研究生，liuzhe@tju.edu.cn.

賈宏杰，hjjia@tju.edu.cn.