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

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

2、問(wèn)題，要求滿足待定的電力系統(tǒng)運(yùn)行和安全約束條件下，通過(guò)調(diào)整系統(tǒng)中可利用控制手段實(shí)現(xiàn)預(yù)定目標(biāo)最優(yōu)的系統(tǒng)穩(wěn)定運(yùn)行狀態(tài)。針對(duì)不同的應(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)運(yùn)行成本最小、或者系統(tǒng)運(yùn)行網(wǎng)損最小。S.t.hx,u=0為等式約束，表示滿足系統(tǒng)穩(wěn)定運(yùn)行的功率平衡。gx,u0為不等式約束，表示電源有功出力的上下界約束、節(jié)點(diǎn)電壓上下線約束、線路傳輸功率上

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

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

5、和拉格朗日乘子滿足簡(jiǎn)單的大于0或者小于0的條件，即可代替原來(lái)必須在可行域內(nèi)求解的要求，使得計(jì)算過(guò)程大為簡(jiǎn)化。一般可以將最優(yōu)潮流模型簡(jiǎn)化為如下的非線性優(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)點(diǎn)法的基本思路是：首先將式（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這樣原問(wèn)題就轉(zhuǎn)化為問(wèn)

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

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

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

9、的是最優(yōu)潮流計(jì)算中所涉及的矩陣方程的描述。3.1、算例描述該算例為王錫凡編寫的現(xiàn)代電力系統(tǒng)分析中的3-1的例題，是以系統(tǒng)燃料最省為最優(yōu)潮流的目標(biāo)函數(shù)。選擇該題目作為算例分析的原因是，該題目有比較詳細(xì)的解題思路以及列寫出了比較詳細(xì)的迭代結(jié)果，方便對(duì)編寫程序的運(yùn)行結(jié)果進(jìn)行比對(duì)。求如下圖所示簡(jiǎn)化系統(tǒng)的系統(tǒng)燃料最省的最優(yōu)潮流計(jì)算：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ī)有功無(wú)功出力上下界和燃料耗費(fèi)曲線 參數(shù)（表3-2）。若不作特殊說(shuō)明，所有數(shù)據(jù)都是以標(biāo)幺值形式給出，功率基準(zhǔn)值為100MVA，母線電壓上下界分別為1.1和0.9。表3-1線路傳輸功率邊界支路號(hào)首末端母線號(hào)線路傳輸功率邊界11-2221-30.6532-3242-4653-55表3-2 發(fā)電機(jī)數(shù)據(jù)發(fā)電機(jī)序號(hào)母線號(hào)出力上界出力下界燃料耗費(fèi)曲線參數(shù)有功無(wú)功有功無(wú)功二次系數(shù)一次系數(shù)常數(shù)14831-350.4395200.43551200.6425851-2.1200.55500.7451857.203.2、程序具體實(shí)現(xiàn)過(guò)程針對(duì)上述系統(tǒng)，首先我

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

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

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

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

20、于要將1節(jié)點(diǎn)電壓抬高至其下界以滿足不等式約束的要求而產(chǎn)生的副作用。但是網(wǎng)損的增加并不影響目標(biāo)函數(shù)的優(yōu)化。整個(gè)系統(tǒng)的燃料費(fèi)用與不優(yōu)化的潮流計(jì)算相比仍然減少了243.76$。表4-4 各有功源有功和無(wú)功源無(wú)功出力發(fā)電機(jī)序號(hào)母線序號(hào)有功出力無(wú)功出力燃料費(fèi)用/$OPFPFOPFPFOPFPF145.50565.001.77801.883113833.063463.80252.15682.57942.61942.29943870.134483.15總計(jì)7.66247.57944.39744.13057703.207946.95表4-5 各節(jié)點(diǎn)電壓向量母線序號(hào)電壓幅值電壓相角/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ì)的，具有算法收斂迅速，魯棒性強(qiáng)，對(duì)初值的選擇不敏感，其迭代次數(shù)與系統(tǒng)規(guī)?；蚩刂谱兞康臄?shù)目關(guān)系不大等特點(diǎn)。即在該簡(jiǎn)化模型中迭代次數(shù)為23次，但是由于其迭代次數(shù)是與系統(tǒng)規(guī)模關(guān)系不大，對(duì)IEEE30、IEEE118節(jié)點(diǎn)系統(tǒng)的計(jì)算結(jié)果其迭代次數(shù)始終保持在21到27次之間。另外將水火電

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

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

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

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

26、9); %燃料耗費(fèi)曲線常數(shù)項(xiàng)%計(jì)算參數(shù)初始化%n=length(N); %節(jié)點(diǎn)個(gè)數(shù)ng=length(Ng); %發(fā)電機(jī)臺(tái)數(shù)nbr=length(Nbr); %之路個(gè)數(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; %拉格朗日乘子%計(jì)算不等式約束的上下限%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; %計(jì)算互補(bǔ)間隙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的情況 %非對(duì)角元素 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); %對(duì)角元素 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、對(duì)角矩陣% 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)*(