基于內(nèi)點法的最優(yōu)潮流計算(共33頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上摘要內(nèi)點法是一種能在可行域內(nèi)部尋優(yōu)的方法，即從初始內(nèi)點出發(fā)，沿著中心路徑方向在可行域內(nèi)部直接走向最優(yōu)解的方法。其中路徑跟蹤法是目前最具有發(fā)展?jié)摿Φ囊活悆?nèi)點算法，該方法魯棒性強，對初值的選擇不敏感，在目前電力系統(tǒng)優(yōu)化問題中得到了廣泛的應(yīng)用。本文采用路徑跟蹤法進(jìn)行最優(yōu)求解，首先介紹了路徑跟蹤法的基本模型，并且結(jié)合具體算例，用編寫的Matlab程序進(jìn)行仿真分析，驗證了該方法在最優(yōu)潮流計算中的優(yōu)越性能。關(guān)鍵詞：最優(yōu)潮流、內(nèi)點法 、路徑跟蹤法、仿真專心-專注-專業(yè)目 次0、引言電力系統(tǒng)最優(yōu)潮流，簡稱OPF（Optimal Power Flow）。OPF問題是一個復(fù)雜的非線性規(guī)劃

2、問題，要求滿足待定的電力系統(tǒng)運行和安全約束條件下，通過調(diào)整系統(tǒng)中可利用控制手段實現(xiàn)預(yù)定目標(biāo)最優(yōu)的系統(tǒng)穩(wěn)定運行狀態(tài)。針對不同的應(yīng)用，OPF模型課以選擇不同的控制變量、狀態(tài)變量集合，不同的目標(biāo)函數(shù)，以及不同的約束條件，其數(shù)學(xué)模型可描述為確定一組最優(yōu)控制變量u，以使目標(biāo)函數(shù)取極小值，并且滿足如下等式和不等式。minufx,uS.t.hx,u=0gx,u0 （0-1）其中minufx,u為優(yōu)化的目標(biāo)函數(shù)，可以表示系統(tǒng)運行成本最小、或者系統(tǒng)運行網(wǎng)損最小。S.t.hx,u=0為等式約束，表示滿足系統(tǒng)穩(wěn)定運行的功率平衡。gx,u0為不等式約束，表示電源有功出力的上下界約束、節(jié)點電壓上下線約束、線路傳輸功率上

3、下線約束等等。電力系統(tǒng)最優(yōu)潮流算法大致可以分為兩類：經(jīng)典算法和智能算法。其中經(jīng)典算法主要是指以簡化梯度法、牛頓法、內(nèi)點法和解耦法為代表的基于線性規(guī)劃和非線性規(guī)劃以及解耦原則的算法，是研究最多的最優(yōu)潮流算法, 這類算法的特點是以一階或二階梯度作為尋找最優(yōu)解的主要信息。智能算法主要是指遺傳算法和模擬退火發(fā)等，這類算法的特點是不以梯度作為尋優(yōu)信息，屬于非導(dǎo)數(shù)的優(yōu)化方法。因此經(jīng)典算法的優(yōu)點是能按目標(biāo)函數(shù)的導(dǎo)數(shù)信息確定搜索方向，計算速度快，算法比較成熟，結(jié)果可信度高。缺點是對目標(biāo)函數(shù)及約束條件有一定的限制，可能出現(xiàn)局部極小時難以收斂。而智能算法的優(yōu)點是計算與導(dǎo)數(shù)無關(guān)，靈活性高，隨機(jī)性強，缺點是算法不穩(wěn)定

4、，結(jié)果不可信，并且控制參數(shù)需憑經(jīng)驗給出。通過對這些常見算法的簡單比較，內(nèi)點法具有其優(yōu)越的性能，特別是路徑跟蹤法，其算法收斂迅速，魯棒性強，對初值的選擇不敏感，其迭代次數(shù)與系統(tǒng)規(guī)?；蚩刂谱兞康臄?shù)目關(guān)系不大，因此本文采用該方法進(jìn)行最優(yōu)計算。1、 路徑跟蹤法的基本數(shù)學(xué)模型內(nèi)點法最初的基本思路是希望通過尋優(yōu)迭代過程始終在可行域內(nèi)進(jìn)行，因此，初始點應(yīng)在可行域內(nèi)，并在可行域的邊界設(shè)置障礙使迭代點接近邊界時其目標(biāo)函數(shù)迅速增大，從而保證迭代點均在可行域的內(nèi)點。但是對于大規(guī)模實際問題而言，尋找初始點往往十分困難。為此許多學(xué)者長期以來致力于內(nèi)點算法初始內(nèi)點條件的改進(jìn)。以下介紹的路徑跟蹤法只要求在尋優(yōu)過程中松弛變量

5、和拉格朗日乘子滿足簡單的大于0或者小于0的條件，即可代替原來必須在可行域內(nèi)求解的要求，使得計算過程大為簡化。一般可以將最優(yōu)潮流模型簡化為如下的非線性優(yōu)化模型。Obj. minufx,u （1-1） s.t. S.t.hx,u=0 （1-2） g-gx,ug- （1-3）其中minufx,u為優(yōu)化的目標(biāo)函數(shù)， S.t.hx,u=0為等式約束， gx,u為不等式約束，路徑跟蹤內(nèi)點法的基本思路是：首先將式（1-3）的不等約束變成等式約束：gx,u+u=g- （1-4）gx,u-l=g- （1-5）其中松弛變量 l=l1,lrT， u=u1,urT，應(yīng)滿足u>0，l>0這樣原問題就轉(zhuǎn)化為問

6、題A：Obj. minufx,u S.t. 然后，把目標(biāo)函數(shù)改造成障礙函數(shù)，該函數(shù)在可行域內(nèi)應(yīng)接近于原函數(shù)f(x)，而在邊界時變得很大。一次可得帶優(yōu)化問題B：obj. s.t. 其中擾動因子或者障礙因子u>0。當(dāng)l或u接近邊界時，以上函數(shù)將趨于無窮大，因此滿足以上障礙目標(biāo)函數(shù)的極小解不可能在邊界上找到。這樣就通過目標(biāo)函數(shù)的變化把含不等式限制的優(yōu)化問題A變成只含等式限制優(yōu)化的問題B了，因此可以直接用拉格朗日乘子法來求解。優(yōu)化問題B的拉格朗日目標(biāo)函數(shù)為： 式中：,和均為拉格朗日乘子。因此最后簡化的求解問題就是求取上述表達(dá)式的極小解。2、 路徑跟蹤法的最優(yōu)潮流求解思路路徑跟蹤法的最優(yōu)潮流求解過

7、程就是對拉格朗日目標(biāo)函數(shù)求極小值問題：式中：,和均為拉格朗日乘子。該問題極小值存在的必要條件是拉格朗日函數(shù)對所有變量及乘子的偏導(dǎo)數(shù)為0。即： （2-1）通過上述表達(dá)式可以解出： （2-2）定義：，稱為互補間隙?？傻茫?（2-3）如果x*是優(yōu)化問題A的最優(yōu)解，當(dāng)u固定時，x(u)是優(yōu)化問題B的解，那么當(dāng)Gap0，u0時，產(chǎn)生的序列x(u)收斂至x*。建議采用：。式中稱為中心參數(shù)，一般取0.1，在大多數(shù)場合可獲得較好的收斂效果。通過偏導(dǎo)數(shù)為0的表達(dá)式可以可得內(nèi)點法的修正方程為： （2-4）求解方程可得到第k次迭代的修正量，于是最優(yōu)解的一個新的近似解為： （2-5）式中，和為步長： （2-6）其潮流

8、計算的流程圖如下圖1所示，其中初始化部分包括：(1)、設(shè)置松弛變量l和u，保證l，uT>0(2)、設(shè)置拉格朗日乘子w、y、z，滿足w<0,z>0，Y！=0T (3)、設(shè)置優(yōu)化問題的初值。(4)、取中心參數(shù)，給定計算精度，迭代次數(shù)初值K=0。輸出最優(yōu)解，停止計算計算互補間隙Gap計算擾動因子計算和更新原始變量及拉格朗日乘子k<KmaxGap<求解修正方程，求出，輸出“計算不收斂”初始化是是否否圖1. 內(nèi)點法潮流計算流程圖3、 具體算例及程序?qū)崿F(xiàn)流程這部分主要有算例描述以及程序的實現(xiàn)流程兩部分，其中算例描述主要是對系統(tǒng)參數(shù)以及優(yōu)化問題進(jìn)行說明。而程序的實現(xiàn)流程主要描述

9、的是最優(yōu)潮流計算中所涉及的矩陣方程的描述。3.1、算例描述該算例為王錫凡編寫的現(xiàn)代電力系統(tǒng)分析中的3-1的例題，是以系統(tǒng)燃料最省為最優(yōu)潮流的目標(biāo)函數(shù)。選擇該題目作為算例分析的原因是，該題目有比較詳細(xì)的解題思路以及列寫出了比較詳細(xì)的迭代結(jié)果，方便對編寫程序的運行結(jié)果進(jìn)行比對。求如下圖所示簡化系統(tǒng)的系統(tǒng)燃料最省的最優(yōu)潮流計算：12243251:1.051.05:1j0.015j0.03j0.25j0.250.08+j0.300.1+j0.350.04+j0.25j0.25j0.253.7+j1.32+j11.6+j0.8除了由上圖所提供的系統(tǒng)母線負(fù)荷功率數(shù)據(jù)、線路參數(shù)和變壓器之路參數(shù)數(shù)據(jù)、變壓器便

10、比數(shù)據(jù)之外，以下順序給出了線路傳輸功率邊界（表3-1），發(fā)電機(jī)有功無功出力上下界和燃料耗費曲線 參數(shù)（表3-2）。若不作特殊說明，所有數(shù)據(jù)都是以標(biāo)幺值形式給出，功率基準(zhǔn)值為100MVA，母線電壓上下界分別為1.1和0.9。表3-1線路傳輸功率邊界支路號首末端母線號線路傳輸功率邊界11-2221-30.6532-3242-4653-55表3-2 發(fā)電機(jī)數(shù)據(jù)發(fā)電機(jī)序號母線號出力上界出力下界燃料耗費曲線參數(shù)有功無功有功無功二次系數(shù)一次系數(shù)常數(shù)14831-350.4395200.43551200.6425851-2.1200.55500.7451857.203.2、程序具體實現(xiàn)過程針對上述系統(tǒng)，首先我

11、們先列寫出該算例的數(shù)學(xué)模型和有關(guān)計算公式。在該算例中，共有5個節(jié)點，相應(yīng)的狀態(tài)量為：系統(tǒng)中有2臺發(fā)電機(jī)，沒有其他無功源，因此控制變量為：應(yīng)該指出，此處發(fā)電機(jī)和無功源的編號與及誒單編號無關(guān)，是獨立編號的。這是因為系統(tǒng)中一個節(jié)點可能接有多臺發(fā)電機(jī)的緣故。因此系統(tǒng)中總變量共有14個：最優(yōu)潮流的數(shù)學(xué)模型為：目標(biāo)函數(shù)：約束條件：每個節(jié)點有兩個潮流方程，共有10個等式約束條件，對非發(fā)電機(jī)而言： （i=1,2,3）對發(fā)電機(jī)節(jié)點： （i=4,5）式中：表示第k臺發(fā)電機(jī)接在節(jié)點i上。不等式約束共有14個條件，分別是： 根據(jù)以上模型可以形成修正方程。該方程包括形成等式左邊的系數(shù)矩陣和等式右邊的常數(shù)項兩部分。1、形

12、成系數(shù)矩陣1）、等式約束的雅克比矩陣等式右端包括3個子矩陣：其中：其中：式中：i為發(fā)電機(jī)的序號；j為節(jié)點號；表示第i臺發(fā)電機(jī)是在節(jié)點j上的。（潮流計算中的雅克比矩陣）2）、不等式約束的雅克比矩陣式中：、和依次表示電源有功出力的上下界約束，無功電源出力的上下界約束，節(jié)點電壓賦值的上下界約束和線路潮流約束。，式中：第行列元素為1，其他元素均為0。3）、對角矩陣4）、海森伯矩陣這是最復(fù)雜的部分，共包含四項。有上述推導(dǎo)已經(jīng)可以得到其中的第四項為：其余三項是：目標(biāo)函數(shù)的海森伯矩陣、等式約束海森伯矩陣與拉格朗日乘子y的乘積和不等式約束海森伯矩陣與拉格朗日乘子的乘積。2、形成常數(shù)項均可根據(jù)定義直接求得?？梢?/p>

13、表示為：當(dāng)知道目標(biāo)函數(shù)梯度矢量之后，再根據(jù)以上等式和不等式約束的雅克比矩陣公式就可以求得。4、 運行結(jié)果及分析41 運行結(jié)果以下對該算例的尋優(yōu)過程用數(shù)字加以說明，設(shè)4、5節(jié)點發(fā)電機(jī)均能有算法調(diào)節(jié)其出力。在初始化過程中各變量初值根據(jù)實際問題自行設(shè)置的，我們給出所用變量的處置如下：節(jié)點電壓；平衡節(jié)點；發(fā)電機(jī)出力有功取其邊界值；松弛因子，當(dāng)收斂條件時，需要迭代進(jìn)行23次（例題所給出的迭代次數(shù)為17次）。表 4-1 迭代過程中各節(jié)點電壓增量的變化情況迭代次數(shù)1-0.-0.-0.26638-0.-0.20.0.0.0.-0.32.1.3.1.3.40.-0.2.-0.2.50.0.-0.189910.-

14、0.6-0.-0.-0.31384-0.-0.7-0.-0.-0.52191-0.-0.80.50336-0.1.-0.1.90.-0.052280.-0.0.10-0.-0.-0.13409-0.-0.11-0.-0.-0.135-0.02272-0.120.0.0.-0.0.130.-0.0.0966-0.0.140.-0.0.-0.0.15-0.-0.-0.016551.24E-05-0.16-0.-0.-0.01235-0.-0.17-0.-0.-0.00473-0.-0.180.-0.002230.-0.-0.190.-0.0.-0.-0.200.-0.0.-0.-0.210.-0.

15、0.0.-0.224.42E-08-7.14E-088.44E-061.71E-05-6.67E-0823-9.43E-097.65E-091.38E-072.40E-081.58E-08表 4-2 迭代過程中各節(jié)點相角增量的變化情況迭代次數(shù)110.10.9.10.9.20.-0.0.-0.0.30.-0.0.-0.0.40.0.0.0.-0.5-28.46696-28.48845-28.5226-28.58592-28.632.35.29.147235.29.7-0.-0.-0.30025-0.-0.8-0.-0.-0.3624-0.07659-0.9-0.0.-0.271170.-0.10

16、6054.9484-6300.9746058.9143059.0885-121.11-0.0.-0.137520.-0.12-131.1933-133.1411-129.232-131.1213-131.21147130.0.0.0.0.140.-0.-0.111210.-0.150.0.03296-0.048286.43E-02-0.160.0.0.-0.0.17-0.-0.-0.18778-0.-0.18-6.68934-6.67987-6.6504-6.62907-6.19-0.06977-0.-0.00065-0.-0.20-0.00417-0.-0.00209-0.00978-0.2

17、10.973 E-060.966 E-060.97 E-060.9718 E-060.9797 E-0622-1.42 E-08-1.43 E-08-1.43 E-08-1.43 E-08-1.41 E-08231.01E-091.01E-091.1 E-081.01 E-091.01 E-09表 4-3 迭代過程中有功源有功、無功源無功增量的變化情況迭代次數(shù)有功源有功出力增量無功源無功出力增量16.0.0.2.21.0.3.-3.30.-0.-11.7288-1.4-0.-0.-15.29238.50.0.3.-1.61.0.3.0.70.0.4.-0.8-0.0.-6.524714.90.

18、0.-3.013663.100.0.1.0.110.0.1.473390.120.-0.-3.813813.13-0.0.-0.458660.140.-0.-0.472941.150.-0.0.3.58E-01160.0.0.0.17-0.0810.-0.247990.18-0.0.-0.22899-0.19-0.36020.-0.22777-0.20-0.0.0.-0.21-0.0.0.-0.22-1.41E-041.34E-042.33E-04-1.17E-0423-2.60E-062.54E-063.21E-08-1.69E-07將各次迭代過程中Gap變化情況繪制成曲線，可以顯示出路勁跟

19、蹤法最優(yōu)潮流計算的收斂特性，如圖4-1所示：圖 4-1 5節(jié)點系統(tǒng)最優(yōu)潮流內(nèi)點法收斂特性圖4-2為5節(jié)點系統(tǒng)最優(yōu)潮流計算結(jié)果截圖，其中包括迭代次數(shù)、燃料總費用、發(fā)電機(jī)有功無功出力、各節(jié)點電壓幅值與相角、以及各支路有功功率。（注：結(jié)果中的值為標(biāo)幺值，功率的基準(zhǔn)值為100MVA）42結(jié)果分析將最優(yōu)潮流計算的結(jié)果和普通潮流計算結(jié)果進(jìn)行比較，其中PF表示為普通潮流計算。普通潮流計算中，發(fā)電機(jī)不會調(diào)節(jié)其出力。即4節(jié)點為PQ節(jié)點，5節(jié)點為平衡節(jié)點。見表4-4和表4-5。從表中可以看出，由于4機(jī)組比5機(jī)組的燃料曲線系數(shù)小，因此4機(jī)組有功出力增加，5機(jī)組有功出力減少。同時系統(tǒng)的網(wǎng)損、無功功率都有所增加。這是由

20、于要將1節(jié)點電壓抬高至其下界以滿足不等式約束的要求而產(chǎn)生的副作用。但是網(wǎng)損的增加并不影響目標(biāo)函數(shù)的優(yōu)化。整個系統(tǒng)的燃料費用與不優(yōu)化的潮流計算相比仍然減少了243.76$。表4-4 各有功源有功和無功源無功出力發(fā)電機(jī)序號母線序號有功出力無功出力燃料費用/$OPFPFOPFPFOPFPF145.50565.001.77801.883113833.063463.80252.15682.57942.61942.29943870.134483.15總計7.66247.57944.39744.13057703.207946.95表4-5 各節(jié)點電壓向量母線序號電壓幅值電壓相角/radOPFPFOPFPF1

21、0.90000.8822-0.007-0.0834021.10001.07790.40490.3116031.08181.0364-0.0571-0.747341.06971.05000.478670.3116051.1001.0500005、 結(jié)論路徑跟蹤法在電力系統(tǒng)中應(yīng)用與求解線性規(guī)劃和非線性規(guī)劃模型中，還是比較有優(yōu)勢的，具有算法收斂迅速，魯棒性強，對初值的選擇不敏感，其迭代次數(shù)與系統(tǒng)規(guī)?；蚩刂谱兞康臄?shù)目關(guān)系不大等特點。即在該簡化模型中迭代次數(shù)為23次，但是由于其迭代次數(shù)是與系統(tǒng)規(guī)模關(guān)系不大，對IEEE30、IEEE118節(jié)點系統(tǒng)的計算結(jié)果其迭代次數(shù)始終保持在21到27次之間。另外將水火電

22、最優(yōu)潮流問題分解為火電最優(yōu)潮流子問題和水電子問題，提供了有效的協(xié)調(diào)算法。但是對于路徑跟蹤法的影響因素還是比較多的，比如初始點的選擇、迭代步長的選取、壁壘參數(shù)的調(diào)整、離散變量的處理等等，如果選取不當(dāng)可能會出現(xiàn)不恰當(dāng)?shù)慕Y(jié)果，因此還需要研究者們做大量的工作。6、 編程中遇到的問題上圖為例題所給出的迭代次數(shù)與仿真結(jié)果，可以看出除了迭代次數(shù)與書中不一致以外，迭代的結(jié)果基本上完全一樣。但在程序的實現(xiàn)過程中，很難按照書上的流程編寫出結(jié)果一致的程序，因為程序中有大量的矩陣計算，而且也不能將書上列寫的矩陣直接翻譯成MATLAB語言，需要進(jìn)行一些不同的處理方式，所以需要在網(wǎng)上去尋找一些算法的實現(xiàn)方式。另外在編程中

23、發(fā)現(xiàn)，王錫凡的電力系統(tǒng)一書中，也有一些公式書寫有誤的現(xiàn)象，所以需要對公式進(jìn)行一定的驗證推導(dǎo)，因此很大程度上會出現(xiàn)問題。另外對于路徑跟蹤法來說，其初始點的選擇、迭代步長的選取、壁壘參數(shù)的調(diào)整等都對計算結(jié)果又一定的影響。而在本程序中，我選取的初始點以及拉格朗日因子也都與例題所提供的有些許不同，如果選擇和書中相同的參數(shù)，計算結(jié)果就會出現(xiàn)問題，這其中的原因，我的猜測是可能在用matalb語言實現(xiàn)過程中，細(xì)節(jié)方面可能與例題所展示的有所出入，由于時間關(guān)系，暫時還未找到原因。參考文獻(xiàn)1 張伯明.高等電力網(wǎng)絡(luò)分析M.清華大學(xué)出版社,2007.2 王錫凡.現(xiàn)代電力系統(tǒng)分析M. 北京科學(xué)出版社,2003 3 張江

24、紅.最優(yōu)潮流算法綜述J.華北電力，2010,074 赫玉國.一種基于KarmarKar內(nèi)點法的最優(yōu)潮流算法J.中國機(jī)電工程學(xué)報，1996,11附錄clcclear%讀取計算參數(shù)%Branch=load('Branch.txt'); %支路參數(shù)文檔Node=load('Node.txt'); %節(jié)點參數(shù)文檔Generator=load('Generator.txt'); %發(fā)電機(jī)參數(shù)文檔%支路數(shù)據(jù)提取Nbr=Branch(:,1); %支路號Nl=Branch(:,2); %支路首節(jié)點Nr=Branch(:,3); %支路末節(jié)點%節(jié)點數(shù)據(jù)提取N=N

25、ode(:,1); %節(jié)點號Type=Node(:,2); %節(jié)點類型Uamp=Node(:,3); %節(jié)點電壓幅值Dlta=Node(:,4); %節(jié)點電壓相角Pd=Node(:,5); %節(jié)點負(fù)荷有功Qd=Node(:,6); %節(jié)點負(fù)荷無功Pg=Node(:,7); %節(jié)點出力有功Qg=Node(:,8); %節(jié)點出力無功%發(fā)電機(jī)數(shù)據(jù)提取Ng=Generator(:,1); %發(fā)電機(jī)序號Nbus=Generator(:,2); %所在母線號a2=Generator(:,7); %燃料耗費曲線二次系數(shù)a1=Generator(:,8); %燃料耗費曲線一次系數(shù)a0=Generator(:,

26、9); %燃料耗費曲線常數(shù)項%計算參數(shù)初始化%n=length(N); %節(jié)點個數(shù)ng=length(Ng); %發(fā)電機(jī)臺數(shù)nbr=length(Nbr); %之路個數(shù)x=zeros(2*(ng+n),1); %控制變量+狀態(tài)量x(1:ng)=Pg(Nbus);x(ng+1:2*ng)=Qg(Nbus);x(2*ng+2):2:2*(ng+n)=Uamp;x(2*ng+1):2:2*(ng+n)-1)=Dlta;l=0.8*ones(2*ng+n+nbr,1); %松弛變量u=1.1*ones(2*ng+n+nbr,1); %松弛變量w=-1.5*ones(2*ng+n+nbr,1); %拉格

27、朗日乘子z=ones(2*ng+n+nbr,1); %拉格朗日乘子y=zeros(2*n,1); %拉格朗日乘子y(1:2:2*n-1)=1e-3;y(2:2:2*n)=-1e-3; %拉格朗日乘子%計算不等式約束的上下限%gmin的下界值gmin=zeros(2*ng+n+nbr,1);gmin(1:ng)=Generator(:,5);gmin(ng+1:2*ng)=Generator(:,6);gmin(2*ng+1:2*ng+n)=Node(:,10);gmin(2*ng+n+1:2*ng+n+nbr)=-Branch(:,8); %gmax得上界值gmax=zeros(2*ng+n+

28、nbr,1);gmax(1:ng)=Generator(:,3);gmax(ng+1:2*ng)=Generator(:,4);gmax(2*ng+1:2*ng+n)=Node(:,9);gmax(2*ng+n+1:2*ng+n+nbr)=Branch(:,8); %生成導(dǎo)納矩陣%Y=zeros(n);for k=1:nbr t1=Branch(k,2);t2=Branch(k,3);R=Branch(k,4);X=Branch(k,5);ban=Branch(k,6);K=Branch(k,7); Y(t1,t1)=Y(t1,t1)+1/(R+j*X)+j*ban; Y(t1,t2)=Y(t

29、1,t2)-1/(K*(R+j*X); Y(t2,t1)=Y(t2,t1)-1/(K*(R+j*X); Y(t2,t2)=Y(t2,t2)+1/(K*K*(R+j*X)+j*ban;endG=real(Y);B=imag(Y);k=0;%迭代次數(shù)%主程序%while k<150 Gap(k+1)=l'*z-u'*w; %計算互補間隙Gap if Gap>1e-3 miu=0.1*Gap(k+1)/(2*(2*ng+n+nbr); %形成系數(shù)矩陣% theta=zeros(n,n); for ii=1:n for jj=1:n theta(ii,jj)=Dlta(ii

30、)-Dlta(jj); end end %1、等式約束雅克比矩陣% hx=zeros(2*(ng+n),2*n); %ah/aP% for ii=1:ng hx(Ng(ii),2*Nbus(ii)-1)=1; end %ah/aQ% for ii=1:ng hx(Ng(ii)+ng,2*Nbus(ii)=1; end %ah/ax% H1=zeros(n,n); J1=zeros(n,n); N1=zeros(n,n); L1=zeros(n,n); for ii=1:n for jj=1:n if ii=jj%i!=j的情況 %非對角元素 H1(ii,jj)=-Uamp(ii)*Uamp(j

31、j)*(G(ii,jj) *sin(theta(ii,jj) -B(ii,jj)*cos(theta(ii,jj); J1(ii,jj)=Uamp(ii)*Uamp(jj)*(G(ii,jj) *cos(theta(ii,jj) +B(ii,jj)*sin(theta(ii,jj); N1(ii,jj)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj) +B(ii,jj)*sin(theta(ii,jj); L1(ii,jj)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj) -B(ii,jj)*cos(theta(ii,jj); %對角元素 H

32、1(ii,ii)=H1(ii,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj) *sin(theta(ii,jj)-B(ii,jj)*cos(theta(ii,jj); J1(ii,ii)=J1(ii,ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj) *cos(theta(ii,jj)+B(ii,jj)*sin(theta(ii,jj); N1(ii,ii)=N1(ii,ii)-Uamp(jj)*(G(ii,jj) *cos(theta(ii,jj) +B(ii,jj)*sin(theta(ii,jj); L1(ii,ii)=L1(ii,ii)-Uamp(jj)*(

33、G(ii,jj) *sin(theta(ii,jj) -B(ii,jj)*cos(theta(ii,jj); end end N1(ii,ii)=N1(ii,ii)-2*Uamp(ii)*G(ii,ii); L1(ii,ii)=L1(ii,ii)+2*Uamp(ii)*B(ii,ii); end hx(1+2*ng:2:2*(n+ng)-1,1:2:2*n-1)=H1' hx(1+2*ng:2:2*(n+ng)-1,2:2:2*n)=J1' hx(2+2*ng:2:2*(n+ng),1:2:2*n-1)=N1' hx(2+2*ng:2:2*(n+ng),2:2:2*n)

34、=L1' %2、不等式約束的雅克比矩陣% agaP=eye(ng,ng) zeros(ng,ng) zeros(ng,n) zeros(ng,nbr); agaQ=zeros(ng,ng) eye(ng,ng) zeros(ng,n) zeros(ng,nbr); ag1ax=zeros(2*n,ng); ag2ax=zeros(2*n,ng); ag3ax=zeros(2*n,n); for ii=1:n ag3ax(2*ii,ii)=1; end ag4ax=zeros(2*n,nbr); for ii=1:n for jj=1:nbr if Nl(jj)=ii ag4ax(2*i

35、i-1,jj)=-Uamp(Nl(jj) *Uamp(Nr(jj)*(G(Nl(jj),Nr(jj) *sin(theta(Nl(jj),Nr(jj)-B(Nl(jj),Nr(jj) *cos(theta(Nl(jj),Nr(jj); ag4ax(2*ii,jj)=Uamp(Nr(jj)*(G(Nl(jj),Nr(jj)* cos(theta(Nl(jj),Nr(jj)+B(Nl(jj),Nr(jj) *sin(theta(Nl(jj),Nr(jj) -2*Uamp(Nl(jj)*G(Nl(jj),Nr(jj); end if Nr(jj)=ii ag4ax(2*ii-1,jj)=Uamp(N

36、l(jj) *Uamp(Nr(jj)*(G(Nl(jj),Nr(jj) *sin(theta(Nl(jj),Nr(jj)-B(Nl(jj),Nr(jj) *cos(theta(Nl(jj),Nr(jj); ag4ax(2*ii,jj)=Uamp(Nl(jj)*(G(Nl(jj),Nr(jj) *cos(theta(Nl(jj),Nr(jj)+B(Nl(jj),Nr(jj) *sin(theta(Nl(jj),Nr(jj); end end end pxg=agaP ; agaQ;ag1ax ag2ax ag3ax ag4ax; %3、對角矩陣% L_1Z=zeros(2*ng+n+nbr,2*

37、ng+n+nbr); U_1W=zeros(2*ng+n+nbr,2*ng+n+nbr); for ii=1:2*ng+n+nbr L_1Z(ii,ii)=z(ii)/l(ii); U_1W(ii,ii)=w(ii)/u(ii); end %4、海森伯矩陣% %將海森伯矩陣分為四塊H1,H2,H3,H4 %H1% A2=diag(a2); H1=zeros(2*(ng+n),2*(ng+n); H1(1:ng,1:ng)=2*A2; %H2% H2=zeros(2*(ng+n),2*(ng+n); A=zeros(2*n,2*n); Apb=zeros(2*n,2*n,n); Aqb=zero

38、s(2*n,2*n,n); for ii=1:n for jj=1:n if ii=jj Apb(2*ii-1,2*ii-1,ii)= Apb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj)+B(ii,jj)*sin(theta(ii,jj); Apb(2*ii-1,2*ii,ii)=Apb(2*ii-1,2*ii,ii) +Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj)-B(ii,jj)*cos(theta(ii,jj); Aqb(2*ii-1,2*ii-1,ii)= Aqb(2*ii-1,

39、2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj)-B(ii,jj)*cos(theta(ii,jj); Aqb(2*ii-1,2*ii,ii)=Aqb(2*ii-1,2*ii,ii) -Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj)+B(ii,jj)*sin(theta(ii,jj); Apb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj) *cos(theta(ii,jj)+B(ii,jj) *sin(theta(ii,jj); Apb(2*jj-1,2*jj,i

40、i)=-Uamp(ii)*(G(ii,jj)* sin(theta(ii,jj)-B(ii,jj) *cos(theta(ii,jj); Apb(2*jj,2*jj-1,ii)=Apb(2*jj-1,2*jj,ii); Aqb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj) *sin(theta(ii,jj)-B(ii,jj) *cos(theta(ii,jj); Aqb(2*jj-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)* cos(theta(ii,jj)+B(ii,jj) *sin(theta(ii,jj); Aqb(2*jj,

41、2*jj-1,ii)=Aqb(2*jj-1,2*jj,ii); Apb(2*ii-1,2*jj-1,ii)=-Uamp(ii)*Uamp(jj)*(G(ii,jj) *cos(theta(ii,jj)+B(ii,jj) *sin(theta(ii,jj); Apb(2*ii-1,2*jj,ii)= Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj)-B(ii,jj)*cos(theta(ii,jj); Apb(2*ii,2*jj-1,ii)= -Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj)-B(ii,jj)*cos(theta(ii,jj); Apb(2*ii,2*jj,ii)=-(G(ii,jj)*cos(theta(ii,jj) +B(ii,jj)*sin(theta(ii,jj); Aqb(2*ii-1,2*jj-1,ii)= -Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj)-B(ii,jj) *cos(theta(ii,jj); Aqb(2*ii-1,2*jj,ii)= -Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj)+B(ii,jj)*sin(theta(ii,jj); Aqb(2*ii,2*jj-1,ii)= Uamp(jj)*(