9節(jié)點(diǎn)系統(tǒng)小干擾穩(wěn)定分析

MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\r1\hSEQMTChap\r2\h9節(jié)點(diǎn)系統(tǒng)小干擾穩(wěn)定分析摘要：關(guān)鍵詞：引言電力系統(tǒng)在運(yùn)行過(guò)程中無(wú)時(shí)不遭受到一些小的干擾，例如負(fù)荷的隨機(jī)變化及隨后的發(fā)電機(jī)組調(diào)節(jié)；因風(fēng)吹引起架空線路線間距離變化從而導(dǎo)致線路等值電抗的變化，等等。這些現(xiàn)象隨時(shí)都在發(fā)生?；ヂ?lián)電力系統(tǒng)可以提高發(fā)輸電的經(jīng)濟(jì)性和可靠性，但是多個(gè)地區(qū)之間電力系統(tǒng)的多重互聯(lián)又誘發(fā)出許多新的動(dòng)態(tài)問(wèn)題，使系統(tǒng)失穩(wěn)的可能性增大。小干擾穩(wěn)定性和暫態(tài)穩(wěn)定性、電壓穩(wěn)定性、頻率穩(wěn)定性一起作為電力系統(tǒng)保持安全運(yùn)行的基本要求，其特殊重要的地位逐漸為人們所認(rèn)識(shí)。自從上世紀(jì)七十年代以來(lái)，世界上多個(gè)國(guó)家的電力系統(tǒng)都曾發(fā)生過(guò)因失去小干擾穩(wěn)定性而導(dǎo)致系統(tǒng)振蕩或電壓崩潰的嚴(yán)重事故，造成了巨大的經(jīng)濟(jì)損失，所以對(duì)小干擾穩(wěn)定性問(wèn)題的研究逐漸為人們所關(guān)注。如今，小干擾穩(wěn)定性分析已經(jīng)成為在電力系統(tǒng)的規(guī)劃和運(yùn)行中必須進(jìn)行的重要內(nèi)容之一。小干擾的發(fā)生一般不會(huì)引起系統(tǒng)結(jié)構(gòu)的變化。電力系統(tǒng)小干擾穩(wěn)定分析研究遭受小干擾后電力系統(tǒng)的穩(wěn)定性。小干擾穩(wěn)定性是指正常運(yùn)行的電力系統(tǒng)受到微小的、瞬時(shí)出現(xiàn)但又立即消失的擾動(dòng)后，恢復(fù)到它原來(lái)運(yùn)行狀況的能力，或者這種擾動(dòng)雖不消失，但可用原有的運(yùn)行狀態(tài)近似地表示可能的新運(yùn)行狀況。系統(tǒng)喪失小干擾穩(wěn)定性表現(xiàn)為以下幾種形式：(1)次同步振蕩：電力系統(tǒng)中一種特殊的運(yùn)行狀態(tài)，此時(shí)系統(tǒng)的輸電網(wǎng)絡(luò)和發(fā)電機(jī)間以一個(gè)或多個(gè)略低于系統(tǒng)同步頻率的自然頻率進(jìn)行著能量的交換。(2)低頻振蕩：由于系統(tǒng)中發(fā)電機(jī)阻尼轉(zhuǎn)矩不足而引起的增幅轉(zhuǎn)子振蕩。電力系統(tǒng)中，低頻振蕩的頻率一般在0.2~3Hz范圍內(nèi)。(3)單調(diào)不穩(wěn)定性：在特定的穩(wěn)態(tài)運(yùn)行條件下，由于缺乏同步轉(zhuǎn)矩而造成的系統(tǒng)不穩(wěn)定。在實(shí)際的大型互聯(lián)系統(tǒng)中，小干擾穩(wěn)定性問(wèn)題通常是由于阻尼不足而引起的低頻振蕩問(wèn)題，而低頻振蕩問(wèn)題又分為局部和全局兩種：(1)局部問(wèn)題可以分為以下幾種：①只與一臺(tái)發(fā)電機(jī)或一個(gè)單獨(dú)的電廠相對(duì)于系統(tǒng)剩余部分的轉(zhuǎn)子角振蕩有關(guān)，這種振蕩稱(chēng)之為局部電廠模式振蕩。②與幾臺(tái)臨近的發(fā)電機(jī)轉(zhuǎn)子之間的振蕩有關(guān)，這種振蕩稱(chēng)之為機(jī)間或站間模式振蕩。③其他可能的局部問(wèn)題包括與裝置的控制相關(guān)的模式的不確定。(2)全局問(wèn)題：是由大量發(fā)電機(jī)組之間的相互影響而造成的，因而有廣泛的影響，它們表現(xiàn)為一個(gè)區(qū)域里的一組發(fā)電機(jī)相對(duì)另一區(qū)域的一組發(fā)電機(jī)發(fā)生擺動(dòng)的振蕩，這種振蕩稱(chēng)之為區(qū)域模式振蕩。對(duì)電力系統(tǒng)小干擾穩(wěn)定性進(jìn)行分析的方法大致可以分為以下幾種：數(shù)值仿真方法、建立在線性模型基礎(chǔ)上的分析方法、小干擾穩(wěn)定域分析方法、非線性理論分析方法和計(jì)及模型不確定性的分析方法。對(duì)小干擾穩(wěn)定性問(wèn)題可以采用線性模型進(jìn)行研究，這種線性模型是將描述系統(tǒng)動(dòng)態(tài)行為的微分方程和代數(shù)方程在穩(wěn)態(tài)運(yùn)行點(diǎn)處線性化后得到的。目前，建立在線性模型基礎(chǔ)上的電力系統(tǒng)小干擾穩(wěn)定性分析方法主要有兩種：以狀態(tài)空間模型描述為基礎(chǔ)的特征值分析法和以傳遞函數(shù)矩陣為基礎(chǔ)的頻域分析法。本文采用全部特征值分析法（QR法）對(duì)9節(jié)點(diǎn)系統(tǒng)在穩(wěn)態(tài)處進(jìn)行小干擾穩(wěn)定分析。小干擾穩(wěn)定分析步驟對(duì)含有FACTS的交直流系統(tǒng)進(jìn)行小干擾穩(wěn)定分析時(shí)，其基本流程可總結(jié)如下：(1)對(duì)給定的西永穩(wěn)定運(yùn)行情況進(jìn)行潮流計(jì)算，求出系統(tǒng)各節(jié)點(diǎn)電壓、電流和功率；(2)將網(wǎng)絡(luò)方程寫(xiě)成分塊矩陣形式，并寫(xiě)出在坐標(biāo)下節(jié)點(diǎn)注入電流偏差與節(jié)點(diǎn)電壓偏差之間的線性化方程，形成網(wǎng)絡(luò)導(dǎo)納矩陣；(3)已知各負(fù)荷的功率及負(fù)荷節(jié)點(diǎn)電壓的穩(wěn)態(tài)值為、、、。根據(jù)負(fù)荷的恒阻抗模型求出負(fù)荷模型的導(dǎo)納矩陣元素、、、，用它們修改導(dǎo)納矩陣中對(duì)應(yīng)于各負(fù)荷節(jié)點(diǎn)的對(duì)角子塊；(4)首先計(jì)算出各發(fā)電機(jī)組中所有變量的初值，然后分別形成矩陣、、、、及、、。然后求出矩陣、、、，從而得到各發(fā)電機(jī)組的線性化方程。對(duì)于其他動(dòng)態(tài)元件，用同樣的方法可得到其線性化方程中的系數(shù)矩陣。直至得出系統(tǒng)中所有動(dòng)態(tài)元件的線性化方程；(5)形成矩陣、、、，進(jìn)一步計(jì)算出系統(tǒng)的狀態(tài)矩陣；(6)應(yīng)用QR法計(jì)算矩陣A的全部特征值，從而判斷系統(tǒng)在所給定的穩(wěn)態(tài)運(yùn)行情況下的小干擾穩(wěn)定性?，F(xiàn)對(duì)上述步驟進(jìn)行詳細(xì)描述。特征值分析法MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\r1\hSEQMTChap\r2\h特征值分析法是以李雅普諾夫第一定律為基礎(chǔ)，將電力系統(tǒng)模型線性化，以狀態(tài)空間法將其描述為一般的線性系統(tǒng)，把描述系統(tǒng)動(dòng)態(tài)的一階微分方程及描述網(wǎng)絡(luò)的代數(shù)方程線性化得到GOTOBUTTONZEqnNum294441REFZEqnNum294441\*Charformat\!(2.1)、GOTOBUTTONZEqnNum125948REFZEqnNum125948\*Charformat\!(2.2)的形式，根據(jù)GOTOBUTTONZEqnNum788133REFZEqnNum788133\*Charformat\!(2.3)求解狀態(tài)矩陣的特征值，然后根據(jù)電力系統(tǒng)固有模式和特征值的對(duì)應(yīng)關(guān)系，從特征值得到模式的阻尼和頻率，從特征向量得到模式和系統(tǒng)中各狀態(tài)變量的關(guān)系，從而給出表示系統(tǒng)穩(wěn)定性的定性和定量的信息。 上式中，為系統(tǒng)受擾動(dòng)后的運(yùn)行狀態(tài)相對(duì)于原平衡狀態(tài)的偏離量；為運(yùn)行變量相對(duì)于原平衡狀態(tài)的偏離量。特征值分析法主要分為全部特征值分析法和部分特征值分析法，其優(yōu)點(diǎn)包括：(1)由于特征值可以準(zhǔn)確反映出振蕩模式的頻率和阻尼以及非振蕩模式的衰減率，因此分析人員可以對(duì)系統(tǒng)的小擾動(dòng)情況有一個(gè)比較全面的了解，并能夠及時(shí)發(fā)現(xiàn)系統(tǒng)中可能存在的全部弱阻尼模式或不穩(wěn)定模式。(2)特征值分析法可以提供有關(guān)特征值與系統(tǒng)參數(shù)間的靈敏度，這將有助于揭示出系統(tǒng)小干擾穩(wěn)定性的實(shí)質(zhì)，進(jìn)而尋求改善系統(tǒng)穩(wěn)定性的對(duì)策。(3)特征值分析法不僅可以用于各種電力系統(tǒng)小干擾穩(wěn)定性問(wèn)題的分析，而且還可以用于電力系統(tǒng)中各種控制器參數(shù)的調(diào)整。利用QR法分析電力系統(tǒng)小干擾穩(wěn)定性的方法已經(jīng)發(fā)展的十分成熟，這種方法數(shù)值解穩(wěn)定，適應(yīng)性強(qiáng)，并且收斂快，是計(jì)算矩陣全部特征值的一種最有效的方法，具有魯棒性強(qiáng)、收斂速度快等特點(diǎn)。然而，QR法也有一定的局限性，一般只適合于中小規(guī)模電力系統(tǒng)的小干擾穩(wěn)定性分析，這是因?yàn)槭褂肣R法用于計(jì)算大型電力系統(tǒng)特征值時(shí)會(huì)遇到以下3個(gè)難以克服的困難：(1)QR法不能夠利用矩陣的稀疏特性，它需要將矩陣滿陣存儲(chǔ)，所需的存儲(chǔ)容量龐大；(2)使用QR法計(jì)算大型矩陣特征值時(shí)非常費(fèi)時(shí)；(3)當(dāng)矩陣的維數(shù)過(guò)高時(shí)還會(huì)導(dǎo)致QR法失效而無(wú)法得出正確的結(jié)果。對(duì)于給定的和正交陣，考慮如下迭代： 其中，每個(gè)都是正交陣，每個(gè)都是上三角陣。由歸納法可得： 這樣，每個(gè)都與相似。由于有復(fù)特征值，因此不會(huì)收斂到嚴(yán)格的、“特征值暴露”的三角陣，而只能滿足于計(jì)算稱(chēng)為實(shí)Schur分解的另一種分解。 對(duì)角均為塊或塊的分塊上三角陣稱(chēng)為擬上三角陣。實(shí)Schur分解相當(dāng)于將矩陣實(shí)歸為一個(gè)擬上三角陣。若，則存在一個(gè)正交陣使得： 其中每個(gè)是矩陣或矩陣。若是的，其元素就是的特征值；若是的，的特征值是的一對(duì)共軛的復(fù)特征值。線性化網(wǎng)絡(luò)方程為了敘述方便，將網(wǎng)絡(luò)方程寫(xiě)成分塊矩陣形式，并注意到網(wǎng)絡(luò)方程本身是線性的，因而可以直接寫(xiě)出在坐標(biāo)系下節(jié)點(diǎn)注入電流偏差與節(jié)點(diǎn)電壓偏差之間的線性化方程： 式中： 對(duì)各負(fù)荷節(jié)點(diǎn)，把負(fù)荷的恒阻抗模型給出的注入電流偏差與節(jié)點(diǎn)電壓偏差關(guān)系代入上式，即可消去負(fù)荷節(jié)點(diǎn)的電流偏差。設(shè)負(fù)荷接在節(jié)點(diǎn)，則消去該負(fù)荷后的網(wǎng)絡(luò)方程僅是對(duì)原網(wǎng)絡(luò)方程GOTOBUTTONZEqnNum267812REFZEqnNum267812\*Charformat\!(2.7)的簡(jiǎn)單修正：節(jié)點(diǎn)的電流偏差變?yōu)?，導(dǎo)納矩陣中的第個(gè)對(duì)角塊變?yōu)?，而其他?nèi)容不變。全系統(tǒng)線性化微分方程各發(fā)電機(jī)組微分方程的線性化方程可寫(xiě)成如下矩陣形式： 而定子電壓方程式的線性化方程表示為： 上式中的和為各發(fā)電機(jī)本身、軸電壓和電流分量的偏差，因此必須把它們轉(zhuǎn)換成統(tǒng)一的同步旋轉(zhuǎn)坐標(biāo)參考軸下的相應(yīng)分量，以便將它們和電力網(wǎng)絡(luò)聯(lián)系起來(lái)。對(duì)于發(fā)電機(jī)電壓、電流可推導(dǎo)如下關(guān)系： 其中： 將GOTOBUTTONZEqnNum265872REFZEqnNum265872\*Charformat\!(2.11)、GOTOBUTTONZEqnNum800204REFZEqnNum800204\*Charformat\!(2.12)帶入GOTOBUTTONZEqnNum123506REFZEqnNum123506\*Charformat\!(2.9)、GOTOBUTTONZEqnNum681572REFZEqnNum681572\*Charformat\!(2.10)消去電流偏差量可得： 其中： 對(duì)于各發(fā)電機(jī)狀態(tài)變量的個(gè)數(shù)，取決于發(fā)電機(jī)選用模型的復(fù)雜程度，根據(jù)需要確定其勵(lì)磁系統(tǒng)、調(diào)速器等模型。9節(jié)點(diǎn)系統(tǒng)描述問(wèn)題描述9節(jié)點(diǎn)電力系統(tǒng)的單線圖、支路數(shù)據(jù)、發(fā)電機(jī)參數(shù)、正常運(yùn)行情況下的系統(tǒng)潮流分別如REF_Ref471506621\h圖31、REF_Ref471506631\h表31、REF_Ref471506633\h表32、REF_Ref471506634\h表33所示。系統(tǒng)頻率為60HZ。各負(fù)荷均用恒定阻抗模擬。發(fā)電機(jī)1采用經(jīng)典模型，發(fā)電機(jī)2和3采用雙軸模型。發(fā)電機(jī)2和3均裝有自并勵(lì)靜止勵(lì)磁系統(tǒng)，其參數(shù)如下：另外，各發(fā)電機(jī)的阻尼系數(shù)均取為1.0.圖STYLEREF1\s3SEQ圖\*ARABIC\s119節(jié)點(diǎn)系統(tǒng)單線結(jié)構(gòu)圖表STYLEREF1\s3SEQ表\*ARABIC\s11支路數(shù)據(jù)首端母線名末端母線名電阻(標(biāo)幺值)電抗(標(biāo)幺值)容納之半(標(biāo)幺值)變壓器非標(biāo)準(zhǔn)變比450.0100.0850.088460.0170.0920.079570.0320.1610.153690.0390.1700.179780.00850.0720.0745890.01190.10080.1045140.00.05761.0270.00.06251.0390.00.05861.0表STYLEREF1\s3SEQ表\*ARABIC\s12發(fā)電機(jī)數(shù)據(jù)發(fā)電機(jī)母線名1147.280.00.14600.06080.06080.06088.961.02212.800.00.89580.11980.86450.19696.000.5351.0336.020.01.31250.18131.25780.25005.890.6001.0表STYLEREF1\s3SEQ表\*ARABIC\s13正常運(yùn)行情況下的系統(tǒng)潮流母線名電壓發(fā)電機(jī)負(fù)荷幅值相角/°有功功率無(wú)功功率有功功率無(wú)功功率11.0400.00000.71640.270521.02509.28001.63000.066531.02504.66480.8500-0.108641.0258-2.216850.9956-3.98881.25000.500061.0127-3.68740.90000.300071.02583.719781.01590.72751.00000.350091.03241.9667小干擾分析利用潮流結(jié)果，計(jì)算各發(fā)電機(jī)變量的初值，如REF_Ref471507017\h表34所示：表STYLEREF1\s3SEQ表\*ARABIC\s14各發(fā)電機(jī)變量的初值發(fā)電機(jī)12.271651.039180.041220.678010.287161.05664261.098440.633610.805710.931991.290151.789320.788170.62220354.136620.666070.779090.619410.561471.402900.767860.62424對(duì)于發(fā)電機(jī)1采用經(jīng)典模型，因而微分方程只有轉(zhuǎn)子運(yùn)動(dòng)方程的兩階，勵(lì)磁系統(tǒng)假定勵(lì)磁電壓保持恒定；對(duì)于發(fā)電機(jī)2和3采用雙軸模型，且均裝有自并勵(lì)靜止勵(lì)磁系統(tǒng)，因而其微分方程均為7階。自并勵(lì)靜止勵(lì)磁系統(tǒng)傳遞函數(shù)框圖如下所示：圖STYLEREF1\s3SEQ圖\*ARABIC\s12自并勵(lì)靜止勵(lì)磁系統(tǒng)傳遞函數(shù)框圖 發(fā)電機(jī)1各矩陣計(jì)算結(jié)果如下所示： 發(fā)電機(jī)2各矩陣計(jì)算結(jié)果如下所示： 發(fā)電機(jī)2各矩陣計(jì)算結(jié)果如下所示：進(jìn)一步，系統(tǒng)的線性化方程GOTOBUTTONZEqnNum294441REFZEqnNum294441\*Charformat\!(2.1)中的矩陣為： 計(jì)算得、、、如下圖所示：圖STYLEREF1\s3SEQ圖\*ARABIC\s13矩陣A圖STYLEREF1\s3SEQ圖\*ARABIC\s14矩陣B圖STYLEREF1\s3SEQ圖\*ARABIC\s15矩陣C圖STYLEREF1\s3SEQ圖\*ARABIC\s16矩陣D 據(jù)此計(jì)算狀態(tài)矩陣如下：圖STYLEREF1\s3SEQ圖\*ARABIC\s17狀態(tài)矩陣 運(yùn)用QR法計(jì)算得矩陣如下：圖STYLEREF1\s3SEQ圖\*ARABIC\s18特征根矩陣 從REF_Ref471513401\h圖38中的對(duì)角元即可得出9節(jié)點(diǎn)系統(tǒng)在穩(wěn)態(tài)下小干擾穩(wěn)定特征根。小干擾穩(wěn)定特征值分析上述特征根矩陣中每個(gè)對(duì)角元為，是矩陣或矩陣。若是的，其元素就是的特征值；若是的，的特征值是的一對(duì)共軛的復(fù)特征值。因此可求得系統(tǒng)特征根為： 其中可近似為0，即我們已知的零特征值，系統(tǒng)其他所有特征值都具有負(fù)實(shí)部，因此系統(tǒng)在給定的運(yùn)行方式下是小干擾穩(wěn)定的。心得體會(huì)附錄Small_Signal_Analysis.mclear;clc;Data1=[1400.05760;2700.06250;3900.05860;450.010.0850.088;460.0170.0920.079;570.0320.1610.153;690.0390.170.179;780.00850.0720.0745;890.01190.10080.1045];%支路數(shù)據(jù)首端母線名末端母線名電阻電抗容納之半Data2=[11.0400.71640.2705;21.0259.281.630.0665;31.0254.66480.85-0.1086;41.0258-2.216800;50.9956-3.98881.250.5;61.0127-3.68740.90.3;71.02583.719700;81.01590.727510.35;91.03241.966700];%%正常運(yùn)行情況下系統(tǒng)潮流母線名電壓幅值電壓相角/°有功無(wú)功Data3=[47.2800.1460.06080.06080.06088.9601;12.800.89580.11980.86450.196960.5351;6.0201.31250.18131.25780.255.890.61];%發(fā)電機(jī)數(shù)據(jù)TjRaXdXd'XqXq'Td0'Tq0'DData4=zeros(3,8);fort=1:3V=Data2(t,2)*cos(Data2(t,3)*pi/180)+Data2(t,2)*sin(Data2(t,3)*pi/180)*j;I=(Data2(t,4)-Data2(t,5)*j)/conj(V);EQ=V+I*(Data3(t,2)+Data3(t,5)*j);delta=atan(imag(EQ)/real(EQ));Data4(t,1)=delta*180/pi;Data4(t,2)=real(V)*cos(delta)+imag(V)*sin(delta);Data4(t,3)=real(V)*sin(delta)-imag(V)*cos(delta);Data4(t,4)=real(I)*cos(delta)+imag(I)*sin(delta);Data4(t,5)=real(I)*sin(delta)-imag(I)*cos(delta);Data4(t,6)=Data4(t,2)+Data3(t,2)*Data4(t,4)+Data3(t,3)*Data4(t,5);Data4(t,7)=Data4(t,2)+Data3(t,2)*Data4(t,4)+Data3(t,4)*Data4(t,5);Data4(t,8)=Data4(t,3)+Data3(t,2)*Data4(t,5)-Data3(t,6)*Data4(t,4);end%各發(fā)電機(jī)變量初值δVq0Vd0Iq0Id0Efq0Eq0'Ed0'[M,N]=size(Data1);n=9;nl=[568];Y=zeros(n,n);fort=1:MY(Data1(t,1),Data1(t,2))=-1/(Data1(t,3)+Data1(t,4)*j);endY=Y+Y.';fort=1:nY(t,t)=-sum(Y(t,:));fors=1:Mift==Data1(s,1)||t==Data1(s,2)Y(t,t)=Y(t,t)+Data1(s,5)*j;endend%初始網(wǎng)絡(luò)導(dǎo)納矩陣fors=1:length(nl)ift==nl(s)Y(t,t)=Y(t,t)+(Data2(t,4)-Data2(t,5)*j)/(Data2(t,2))^2;break;%恒阻抗負(fù)荷模型修正后的網(wǎng)絡(luò)導(dǎo)納矩陣endendendGG=real(Y);BB=imag(Y);YY=zeros(2*n,2*n);fort=1:nfors=1:nYY(2*t-1,2*s-1)=GG(t,s);YY(2*t-1,2*s)=-BB(t,s);YY(2*t,2*s-1)=BB(t,s);YY(2*t,2*s)=GG(t,s);endendAg=cell(1,3);BIg=cell(1,3);BVg=cell(1,3);Pg=cell(1,3);Zg=cell(1,3);Zgg=cell(1,3);RVg=cell(1,3);RIg=cell(1,3);Tg=cell(1,3);AG=cell(1,3);BG=cell(1,3);CG=cell(1,3);DG=cell(1,3);wb=120*pi;TB=10;TC=1;TA=0.02;KA=200;TR=0.03;XC=0;fort=1:3ift==1Ag{t}=[0wb;0-Data3(t,9)/Data3(t,1)];BIg{t}=[00;-(Data3(t,5)-Data3(t,4))*Data4(t,4)/Data3(t,1)-(Data4(t,7)+(Data3(t,5)-Data3(t,4))*Data4(t,5))/Data3(t,1)];BVg{t}=[00;00];Pg{t}=[00;00];Zg{t}=[-Data3(t,2)Data3(t,5);-Data3(t,4)-Data3(t,2)];Zgg{t}=pinv(Zg{t});RVg{t}=[Data4(t,2)0;-Data4(t,3)0];RIg{t}=[Data4(t,4)0;-Data4(t,5)0];Tg{t}=[sin(Data4(t,1)*pi/180)-cos(Data4(t,1)*pi/180);cos(Data4(t,1)*pi/180)sin(Data4(t,1)*pi/180)];elseAg{t}=[0wb00000;0-Data3(t,9)/Data3(t,1)-Data4(t,4)/Data3(t,1)-Data4(t,5)/Data3(t,1)000;00-1/Data3(t,7)01/Data3(t,7)00;000-1/Data3(t,8)000;0000-1/TB-(TC/TA-1)/TB-KA*TC/TA/TB;00000-1/TA-KA/TA;000000-1/TR];BIg{t}=[00;-(Data4(t,8)+(Data3(t,6)-Data3(t,4))*Data4(t,4))/Data3(t,1)-(Data4(t,7)+(Data3(t,6)-Data3(t,4))*Data4(t,5))/Data3(t,1);-(Data3(t,3)-Data3(t,4))/Data3(t,7)0;0(Data3(t,5)-Data3(t,6))/Data3(t,8);00;00;00];BVg{t}=[00;00;00;00;00;00;Data4(t,3)/Data2(t,2)/TRData4(t,2)/Data2(t,2)/TR];Pg{t}=[0001000;0010000];Zg{t}=[-Data3(t,2)Data3(t,6);