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

1、摘要內(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)化問題中得到了廣泛的應(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)法 、路徑跟蹤法、仿真 電力系統(tǒng)穩(wěn)態(tài)分析課程文論 第 頁 共 頁目 次0、引言11、路徑跟蹤法的基本數(shù)學(xué)模型22、路徑跟蹤法的最優(yōu)潮流求解思路33、具體算例及程序?qū)崿F(xiàn)流程63.1、

2、算例描述63.2、程序具體實(shí)現(xiàn)流程74、運(yùn)行結(jié)果及分析1141 運(yùn)行結(jié)果1142結(jié)果分析165、結(jié)論176、編程中遇到的問題17參考文獻(xiàn)19附錄20 電力系統(tǒng)穩(wěn)態(tài)分析課程論文 第 31 頁 共 30頁0、引言電力系統(tǒng)最優(yōu)潮流，簡稱OPF（Optimal Power Flow）。OPF問題是一個(gè)復(fù)雜的非線性規(guī)劃問題，要求滿足待定的電力系統(tǒng)運(yùn)行和安全約束條件下，通過調(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ù)取極小值，并且滿足如下

3、等式和不等式。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)電壓上下線約束、線路傳輸功率上下線約束等等。電力系統(tǒng)最優(yōu)潮流算法大致可以分為兩類：經(jīng)典算法和智能算法。其中經(jīng)典算法主要是指以簡化梯度法、牛頓法、內(nèi)點(diǎn)法和解耦法為代表的基于線性規(guī)劃和非線性規(guī)劃以及解耦原則的算法，是研究最多的最優(yōu)潮流算法, 這類算法的特點(diǎn)是以一階或二階梯度作為尋找最優(yōu)解的主要信息。智能算法主要是指遺傳算法

4、和模擬退火發(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ù)無關(guān)，靈活性高，隨機(jī)性強(qiáng)，缺點(diǎn)是算法不穩(wěn)定，結(jié)果不可信，并且控制參數(shù)需憑經(jīng)驗(yàn)給出。通過對(duì)這些常見算法的簡單比較，內(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)法最初的基本思路是希

5、望通過尋優(yōu)迭代過程始終在可行域內(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í)際問題而言，尋找初始點(diǎn)往往十分困難。為此許多學(xué)者長期以來致力于內(nèi)點(diǎn)算法初始內(nèi)點(diǎn)條件的改進(jìn)。以下介紹的路徑跟蹤法只要求在尋優(yōu)過程中松弛變量和拉格朗日乘子滿足簡單的大于0或者小于0的條件，即可代替原來必須在可行域內(nèi)求解的要求，使得計(jì)算過程大為簡化。一般可以將最優(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

6、為優(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)滿足u0，l0這樣原問題就轉(zhuǎn)化為問題A：Obj. minufx,u S.t. 然后，把目標(biāo)函數(shù)改造成障礙函數(shù)，該函數(shù)在可行域內(nèi)應(yīng)接近于原函數(shù)f(x)，而在邊界時(shí)變得很大。一次可得帶優(yōu)化問題B：obj. s.t. 其中擾動(dòng)因子或者障礙因子u0。當(dāng)l或u接近邊界時(shí)，以上函數(shù)將趨于無窮大，因此滿足以上障礙目標(biāo)函數(shù)的極小解不可能在邊界上找到。這樣

7、就通過目標(biāo)函數(shù)的變化把含不等式限制的優(yōu)化問題A變成只含等式限制優(yōu)化的問題B了，因此可以直接用拉格朗日乘子法來求解。優(yōu)化問題B的拉格朗日目標(biāo)函數(shù)為： 式中：,和均為拉格朗日乘子。因此最后簡化的求解問題就是求取上述表達(dá)式的極小解。2、 路徑跟蹤法的最優(yōu)潮流求解思路路徑跟蹤法的最優(yōu)潮流求解過程就是對(duì)拉格朗日目標(biāo)函數(shù)求極小值問題：式中：,和均為拉格朗日乘子。該問題極小值存在的必要條件是拉格朗日函數(shù)對(duì)所有變量及乘子的偏導(dǎo)數(shù)為0。即： （2-1）通過上述表達(dá)式可以解出： （2-2）定義：，稱為互補(bǔ)間隙?？傻茫?（2-3）如果x*是優(yōu)化問題A的最優(yōu)解，當(dāng)u固定時(shí)，x(u)是優(yōu)化問題B的解，那么當(dāng)Gap0，u

8、0時(shí)，產(chǎn)生的序列x(u)收斂至x*。建議采用：。式中稱為中心參數(shù)，一般取0.1，在大多數(shù)場(chǎng)合可獲得較好的收斂效果。通過偏導(dǎo)數(shù)為0的表達(dá)式可以可得內(nèi)點(diǎn)法的修正方程為： （2-4）求解方程可得到第k次迭代的修正量，于是最優(yōu)解的一個(gè)新的近似解為： （2-5）式中，和為步長： （2-6）其潮流計(jì)算的流程圖如下圖1所示，其中初始化部分包括：(1)、設(shè)置松弛變量l和u，保證l，uT0(2)、設(shè)置拉格朗日乘子w、y、z，滿足w0，Y！=0T (3)、設(shè)置優(yōu)化問題的初值。(4)、取中心參數(shù)，給定計(jì)算精度，迭代次數(shù)初值K=0。輸出最優(yōu)解，停止計(jì)算計(jì)算互補(bǔ)間隙Gap計(jì)算擾動(dòng)因子計(jì)算和更新原始變量及拉格朗日乘子kK

9、maxGap求解修正方程，求出，輸出“計(jì)算不收斂”初始化是是否否圖1. 內(nèi)點(diǎn)法潮流計(jì)算流程圖3、 具體算例及程序?qū)崿F(xiàn)流程這部分主要有算例描述以及程序的實(shí)現(xiàn)流程兩部分，其中算例描述主要是對(duì)系統(tǒng)參數(shù)以及優(yōu)化問題進(jìn)行說明。而程序的實(shí)現(xiàn)流程主要描述的是最優(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ì)。求如下圖所示簡化系統(tǒng)的系統(tǒng)燃料最省的最優(yōu)潮流計(jì)算：12243251:1.051.05:

10、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ù)、變壓器便比數(shù)據(jù)之外，以下順序給出了線路傳輸功率邊界（表3-1），發(fā)電機(jī)有功無功出力上下界和燃料耗費(fèi)曲線 參數(shù)（表3-2）。若不作特殊說明，所有數(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ā)電

11、機(jī)序號(hào)母線號(hào)出力上界出力下界燃料耗費(fèi)曲線參數(shù)有功無功有功無功二次系數(shù)一次系數(shù)常數(shù)14831-350.4395200.43551200.6425851-2.1200.55500.7451857.203.2、程序具體實(shí)現(xiàn)過程針對(duì)上述系統(tǒng)，首先我們先列寫出該算例的數(shù)學(xué)模型和有關(guān)計(jì)算公式。在該算例中，共有5個(gè)節(jié)點(diǎn)，相應(yīng)的狀態(tài)量為：系統(tǒng)中有2臺(tái)發(fā)電機(jī)，沒有其他無功源，因此控制變量為：應(yīng)該指出，此處發(fā)電機(jī)和無功源的編號(hào)與及誒單編號(hào)無關(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è)等式約

12、束條件，對(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、形成系數(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）、不等式約束的雅克比矩陣式中：、和依次表示電源有功出力的上下界約束，無功電源出力的上下界約束，節(jié)點(diǎn)電壓賦值的上下界約束和線路潮流約束。，式中：第行列元素為1，其他元素均為0。3）、對(duì)角矩陣4

13、）、海森伯矩陣這是最復(fù)雜的部分，共包含四項(xiàng)。有上述推導(dǎo)已經(jīng)可以得到其中的第四項(xiàng)為：其余三項(xiàng)是：目標(biāo)函數(shù)的海森伯矩陣、等式約束海森伯矩陣與拉格朗日乘子y的乘積和不等式約束海森伯矩陣與拉格朗日乘子的乘積。2、形成常數(shù)項(xiàng)均可根據(jù)定義直接求得?？梢员硎緸椋寒?dāng)知道目標(biāo)函數(shù)梯度矢量之后，再根據(jù)以上等式和不等式約束的雅克比矩陣公式就可以求得。4、 運(yùn)行結(jié)果及分析41 運(yùn)行結(jié)果以下對(duì)該算例的尋優(yōu)過程用數(shù)字加以說明，設(shè)4、5節(jié)點(diǎn)發(fā)電機(jī)均能有算法調(diào)節(jié)其出力。在初始化過程中各變量初值根據(jù)實(shí)際問題自行設(shè)置的，我們給出所用變量的處置如下：節(jié)點(diǎn)電壓；平衡節(jié)點(diǎn)；發(fā)電機(jī)出力有功取其邊界值；松弛因子，當(dāng)收斂條件時(shí)，需要迭代進(jìn)行

14、23次（例題所給出的迭代次數(shù)為17次）。表 4-1 迭代過程中各節(jié)點(diǎn)電壓增量的變化情況迭代次數(shù)1-0.342655-0.167831-0.26638-0.229198-0.2680356620.5334240.81165980.1218280.7537347-0.0262123332.94758741.72816763.3253111.43776433.00165032240.9044296-0.2274022.314238-0.4104682.40614206250.03175650.2589346-0.189910.2779252-0.231011026-0.301923-0.117806

15、-0.31384-0.085802-0.291482167-0.396482-0.116134-0.52191-0.068149-0.5177984180.50336-0.0459651.092422-0.1343191.12869775190.2101792-0.052280.570784-0.0951050.60512496910-0.150019-0.039233-0.13409-0.021024-0.1143501811-0.141952-0.039185-0.135-0.02272-0.1228649120.35872590.00295170.638045-0.0505250.663

16、350104130.0328824-0.0071680.0966-0.0139860.101132708140.0452686-0.0075060.127766-0.01478900.037015-0.005358-0.016551.24E-05-0.0074075916-0.030654-0.005288-0.01235-0.001191-0.0056632617-0.012631-0.003095-0.00473-0.004306-0.00333363180.0016418-0.002230.002955-0.004447-0.00095998190.001716

17、6-0.0023980.001914-0.004628-0.00253085200.0005951-0.0008470.000732-0.000783-0.0008939210.000136-0.0001950.0001870.0001657-0.0002038224.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é)點(diǎn)相角增量的變化情況迭代次數(shù)110.16611310.6710059.95933310.7655089.98744180920.23178

18、66-0.0563220.443469-0.1597860.48707519430.3617911-0.8085570.725769-0.8068860.25162957140.15066350.16947270.0832580.2051569-0.015875295-28.46696-28.48845-28.5226-28.58592-28.3776798632.47947435.22075729.147235.15486129.067434157-0.470667-0.412112-0.30025-0.363094-0.181197188-0.467826-0.345623-0.3624-

19、0.07659-0.68697569-0.0863110.0376297-0.271170.0365628-0.42343914106054.9484-6300.9746058.9143059.0885-121.10432511-0.1588960.0202972-0.137520.0226728-0.1302101412-131.1933-133.1411-129.232-131.1213-131.21147130.4096510.41544330.4634270.41828270.489259247140.0008763-0.001778-0.111210.0019211-0.101304

20、29150.02917830.03296-0.048286.43E-02-0.03148972160.17167570.10090020.006839-0.0628160.02257510817-0.214753-0.185676-0.18778-0.118917-0.1764744218-6.68934-6.67987-6.6504-6.62907-6.638764219-0.06977-0.045184-0.00065-0.011004-0.0887467620-0.00417-0.006299-0.00209-0.00978-0.0002933210.973 E-060.966 E-06

21、0.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.36212440.93787560.4865922.042616921.14615720.75196063.407153-3.18149530.2643491-0.817387-11.7288-1.3212614-0.263123-0.1

22、13899-15.29238.53730250.99069920.24277053.702688-1.67716861.07181020.30747043.1565150.539120170.55939840.72454994.392399-0.9614078-0.2176880.3751495-6.524714.529966490.36618140.078903-3.013663.0701114100.63010530.01895251.6192940.7513961110.14071980.41629161.473390.4541527120.2992523-0.178165-3.8138

23、13.430855713-0.2616330.3695489-0.458660.6443822140.3030544-0.117622-0.472941.1194583150.2827053-0.1053590.4451043.58E-01160.14023780.00272430.3407110.280824817-0.0810.1301595-0.247990.083329618-0.3330460.3030621-0.22899-0.09383819-0.36020.3228062-0.22777-0.11854420-0.0486180.03392570.044021-0.039712

24、21-0.0082120.00447790.009156-0.00069422-1.41E-041.34E-042.33E-04-1.17E-0423-2.60E-062.54E-063.21E-08-1.69E-07將各次迭代過程中Gap變化情況繪制成曲線，可以顯示出路勁跟蹤法最優(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ī)有功無功出力、各節(jié)點(diǎn)電壓幅值與相角、以及各支路有功功率。（注：結(jié)果中的值為標(biāo)幺值，功率的基準(zhǔn)值為100MVA）42結(jié)果分析將最優(yōu)潮流計(jì)算的結(jié)果和普通潮流計(jì)算結(jié)

25、果進(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)損、無功功率都有所增加。這是由于要將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 各有功源有功和無功源無功出力發(fā)電機(jī)序號(hào)母線序號(hào)有功出力無功出力燃料費(fèi)用/$OPFPFOPFPFOPFPF145.50565.001.77801.

26、883113833.063463.80252.15682.57942.61942.29943870.134483.15總計(jì)7.66247.57944.39744.13057703.207946.95表4-5 各節(jié)點(diǎn)電壓向量母線序號(hào)電壓幅值電壓相角/radOPFPFOPFPF10.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ì)

27、的，具有算法收斂迅速，魯棒性強(qiáng)，對(duì)初值的選擇不敏感，其迭代次數(shù)與系統(tǒng)規(guī)模或控制變量的數(shù)目關(guān)系不大等特點(diǎn)。即在該簡化模型中迭代次數(shù)為23次，但是由于其迭代次數(shù)是與系統(tǒng)規(guī)模關(guān)系不大，對(duì)IEEE30、IEEE118節(jié)點(diǎn)系統(tǒng)的計(jì)算結(jié)果其迭代次數(shù)始終保持在21到27次之間。另外將水火電最優(yōu)潮流問題分解為火電最優(yōu)潮流子問題和水電子問題，提供了有效的協(xié)調(diào)算法。但是對(duì)于路徑跟蹤法的影響因素還是比較多的，比如初始點(diǎn)的選擇、迭代步長的選取、壁壘參數(shù)的調(diào)整、離散變量的處理等等，如果選取不當(dāng)可能會(huì)出現(xiàn)不恰當(dāng)?shù)慕Y(jié)果，因此還需要研究者們做大量的工作。6、 編程中遇到的問題上圖為例題所給出的迭代次數(shù)與仿真結(jié)果，可以看出除了

28、迭代次數(shù)與書中不一致以外，迭代的結(jié)果基本上完全一樣。但在程序的實(shí)現(xiàn)過程中，很難按照書上的流程編寫出結(jié)果一致的程序，因?yàn)槌绦蛑杏写罅康木仃囉?jì)算，而且也不能將書上列寫的矩陣直接翻譯成MATLAB語言，需要進(jìn)行一些不同的處理方式，所以需要在網(wǎng)上去尋找一些算法的實(shí)現(xiàn)方式。另外在編程中發(fā)現(xiàn)，王錫凡的電力系統(tǒng)一書中，也有一些公式書寫有誤的現(xiàn)象，所以需要對(duì)公式進(jìn)行一定的驗(yàn)證推導(dǎo)，因此很大程度上會(huì)出現(xiàn)問題。另外對(duì)于路徑跟蹤法來說，其初始點(diǎn)的選擇、迭代步長的選取、壁壘參數(shù)的調(diào)整等都對(duì)計(jì)算結(jié)果又一定的影響。而在本程序中，我選取的初始點(diǎn)以及拉格朗日因子也都與例題所提供的有些許不同，如果選擇和書中相同的參數(shù)，計(jì)算結(jié)果

29、就會(huì)出現(xiàn)問題，這其中的原因，我的猜測(cè)是可能在用matalb語言實(shí)現(xiàn)過程中，細(xì)節(jié)方面可能與例題所展示的有所出入，由于時(shí)間關(guān)系，暫時(shí)還未找到原因。參考文獻(xiàn)1 張伯明.高等電力網(wǎng)絡(luò)分析M.清華大學(xué)出版社,2007.2 王錫凡.現(xiàn)代電力系統(tǒng)分析M. 北京科學(xué)出版社,2003 3 張江紅.最優(yōu)潮流算法綜述J.華北電力，2010,074 赫玉國.一種基于KarmarKar內(nèi)點(diǎn)法的最優(yōu)潮流算法J.中國機(jī)電工程學(xué)報(bào)，1996,11附錄clcclear%讀取計(jì)算參數(shù)%Branch=load(Branch.txt); %支路參數(shù)文檔Node=load(Node.txt); %節(jié)點(diǎn)參數(shù)文檔Generator=loa

30、d(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=Node(:,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ù)荷無功Pg=Node(:,7); %節(jié)點(diǎn)出力有功Qg=Node(:,8); %節(jié)點(diǎn)出力無功%發(fā)電機(jī)數(shù)據(jù)提取Ng=Generator(:,1);

31、%發(fā)電機(jī)序號(hào)Nbus=Generator(:,2); %所在母線號(hào)a2=Generator(:,7); %燃料耗費(fèi)曲線二次系數(shù)a1=Generator(:,8); %燃料耗費(fèi)曲線一次系數(shù)a0=Generator(:,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

32、*(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); %拉格朗日乘子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)=Gener

33、ator(:,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+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

34、);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(t1,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 k1e-3 miu=0.1*Gap(k+1)/(2*(2*ng+n+nbr); %形成系數(shù)矩陣% theta=zer

35、os(n,n); for ii=1:n for jj=1:n theta(ii,jj)=Dlta(ii)-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

36、 if ii=jj%i!=j的情況 %非對(duì)角元素 H1(ii,jj)=-Uamp(ii)*Uamp(jj)*(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)*s

37、in(theta(ii,jj) -B(ii,jj)*cos(theta(ii,jj); %對(duì)角元素 H1(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,j

38、j)*sin(theta(ii,jj); L1(ii,ii)=L1(ii,ii)-Uamp(jj)*(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

39、*n-1)=N1; hx(2+2*ng:2:2*(n+ng),2:2:2*n)=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

40、for jj=1:nbr if Nl(jj)=ii ag4ax(2*ii-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 i

41、f Nr(jj)=ii ag4ax(2*ii-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(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

42、; %3、對(duì)角矩陣% L_1Z=zeros(2*ng+n+nbr,2*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=zeros(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)*