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

1、 Company Document number：WTUT-WT88Y-W8BBGB-BWYTT-19998基于內(nèi)點(diǎn)法的最優(yōu)潮流計(jì)算摘要內(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)，對初值的選擇不敏感，在目前電力系統(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)法 、路徑跟蹤法、仿真目 次0、引言電力系統(tǒng)最優(yōu)

2、潮流，簡稱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)。針對不同的應(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)行的功率平衡

3、。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)解的主要信息。智能算法主要是指遺傳算法和模擬退火發(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)是對目標(biāo)函數(shù)及約束條件有一定的限制，可能出現(xiàn)

4、局部極小時(shí)難以收斂。而智能算法的優(yōu)點(diǎn)是計(jì)算與導(dǎo)數(shù)無關(guān)，靈活性高，隨機(jī)性強(qiáng)，缺點(diǎn)是算法不穩(wěn)定，結(jié)果不可信，并且控制參數(shù)需憑經(jīng)驗(yàn)給出。通過對這些常見算法的簡單比較，內(nèi)點(diǎn)法具有其優(yōu)越的性能，特別是路徑跟蹤法，其算法收斂迅速，魯棒性強(qiáng)，對初值的選擇不敏感，其迭代次數(shù)與系統(tǒng)規(guī)模或控制變量的數(shù)目關(guān)系不大，因此本文采用該方法進(jìn)行最優(yōu)計(jì)算。1、 路徑跟蹤法的基本數(shù)學(xué)模型內(nèi)點(diǎn)法最初的基本思路是希望通過尋優(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)。但是對于大規(guī)模實(shí)際問題而言，尋找初始點(diǎn)往往十分困難。為此許多學(xué)

5、者長期以來致力于內(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.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=l

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

7、蹤法的最優(yōu)潮流求解思路路徑跟蹤法的最優(yōu)潮流求解過程就是對拉格朗日目標(biāo)函數(shù)求極小值問題：式中：,和均為拉格朗日乘子。該問題極小值存在的必要條件是拉格朗日函數(shù)對所有變量及乘子的偏導(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，u0時(shí)，產(chǎn)生的序列x(u)收斂至x*。建議采用：。式中稱為中心參數(shù)，一般取，在大多數(shù)場合可獲得較好的收斂效果。通過偏導(dǎo)數(shù)為0的表達(dá)式可以可得內(nèi)點(diǎn)法的修正方程為： （2-4）求解方程可得到第k次迭代的修正量，于是最優(yōu)解的一個(gè)新的近似解為： （

8、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ì)算和更新原始變量及拉格朗日乘子kKmaxGap求解修正方程，求出，輸出“計(jì)算不收斂”初始化是是否否圖1. 內(nèi)點(diǎn)法潮流計(jì)算流程圖3、 具體算例及程序?qū)崿F(xiàn)流程這部分主要有算例描述以及程序的實(shí)現(xiàn)流程兩部分，其中算例描述主要是對系統(tǒng)參數(shù)以及優(yōu)化問題進(jìn)行說明。而程序的實(shí)現(xiàn)流程主要描述的是

9、最優(yōu)潮流計(jì)算中所涉及的矩陣方程的描述。、算例描述該算例為王錫凡編寫的現(xiàn)代電力系統(tǒng)分析中的3-1的例題，是以系統(tǒng)燃料最省為最優(yōu)潮流的目標(biāo)函數(shù)。選擇該題目作為算例分析的原因是，該題目有比較詳細(xì)的解題思路以及列寫出了比較詳細(xì)的迭代結(jié)果，方便對編寫程序的運(yùn)行結(jié)果進(jìn)行比對。求如下圖所示簡化系統(tǒng)的系統(tǒng)燃料最省的最優(yōu)潮流計(jì)算：12243251:1+2+j1+除了由上圖所提供的系統(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)值為10

10、0MVA，母線電壓上下界分別為和。表3-1線路傳輸功率邊界支路號首末端母線號線路傳輸功率邊界11-2221-332-3242-4653-55表3-2 發(fā)電機(jī)數(shù)據(jù)發(fā)電機(jī)序號母線號出力上界出力下界燃料耗費(fèi)曲線參數(shù)有功無功有功無功二次系數(shù)一次系數(shù)常數(shù)14831-325851、程序具體實(shí)現(xiàn)過程針對上述系統(tǒng)，首先我們先列寫出該算例的數(shù)學(xué)模型和有關(guān)計(jì)算公式。在該算例中，共有5個(gè)節(jié)點(diǎn)，相應(yīng)的狀態(tài)量為：系統(tǒng)中有2臺發(fā)電機(jī)，沒有其他無功源，因此控制變量為：應(yīng)該指出，此處發(fā)電機(jī)和無功源的編號與及誒單編號無關(guān)，是獨(dú)立編號的。這是因?yàn)橄到y(tǒng)中一個(gè)節(jié)點(diǎn)可能接有多臺發(fā)電機(jī)的緣故。因此系統(tǒng)中總變量共有14個(gè)：最優(yōu)潮流的數(shù)學(xué)模

11、型為：目標(biāo)函數(shù)：約束條件：每個(gè)節(jié)點(diǎn)有兩個(gè)潮流方程，共有10個(gè)等式約束條件，對非發(fā)電機(jī)而言： （i=1,2,3）對發(fā)電機(jī)節(jié)點(diǎn)： （i=4,5）式中：表示第k臺發(fā)電機(jī)接在節(jié)點(diǎn)i上。不等式約束共有14個(gè)條件，分別是：根據(jù)以上模型可以形成修正方程。該方程包括形成等式左邊的系數(shù)矩陣和等式右邊的常數(shù)項(xiàng)兩部分。1、形成系數(shù)矩陣1）、等式約束的雅克比矩陣等式右端包括3個(gè)子矩陣：其中：其中：式中：i為發(fā)電機(jī)的序號；j為節(jié)點(diǎn)號；表示第i臺發(fā)電機(jī)是在節(jié)點(diǎn)j上的。（潮流計(jì)算中的雅克比矩陣）2）、不等式約束的雅克比矩陣式中：、和依次表示電源有功出力的上下界約束，無功電源出力的上下界約束，節(jié)點(diǎn)電壓賦值的上下界約束和線路潮

12、流約束。，式中：第行列元素為1，其他元素均為0。3）、對角矩陣4）、海森伯矩陣這是最復(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é)果以下對該算例的尋優(yōu)過程用數(shù)字加以說明，設(shè)4、5節(jié)點(diǎn)發(fā)電機(jī)均能有算法調(diào)節(jié)其出力。在初始化過程中各變量初值根據(jù)實(shí)際問題自行設(shè)置的，我們給出所用變量的處置如下：節(jié)點(diǎn)電壓；平衡節(jié)點(diǎn)

13、；發(fā)電機(jī)出力有功取其邊界值；松弛因子，當(dāng)收斂條件時(shí)，需要迭代進(jìn)行23次（例題所給出的迭代次數(shù)為17次）。表 4-1 迭代過程中各節(jié)點(diǎn)電壓增量的變化情況迭代次數(shù)1234567891011121314151617181920212223表 4-2 迭代過程中各節(jié)點(diǎn)相角增量的變化情況迭代次數(shù)123456789101112131415161718192021 E-06 E-06 E-06 E-06 E-0622 E-08 E-08 E-08 E-08 E-0823 E-08 E-09 E-09表 4-3 迭代過程中有功源有功、無功源無功增量的變化情況迭代次數(shù)有功源有功出力增量無功源無功出力增量1234

14、567891011121314151617181920212223將各次迭代過程中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é)果進(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。從表中可以看出

15、，由于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ì)算相比仍然減少了$。表4-4 各有功源有功和無功源無功出力發(fā)電機(jī)序號母線序號有功出力無功出力燃料費(fèi)用/$OPFPFOPFPFOPFPF1425總計(jì)表4-5 各節(jié)點(diǎn)電壓向量母線序號電壓幅值電壓相角/radOPFPFOPFPF12345005、 結(jié)論路徑跟蹤法在電力系統(tǒng)中應(yīng)用與求解線性規(guī)劃和非線性規(guī)劃模型中，還是比較有優(yōu)勢的，具有算法收斂迅

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

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

18、中的原因，我的猜測是可能在用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(); %支路參數(shù)文檔Node=load(); %節(jié)點(diǎn)參數(shù)文檔Generator=load(); %發(fā)電機(jī)參數(shù)文檔%支路數(shù)據(jù)提取Nbr=Bra

19、nch(:,1); %支路號Nl=Branch(:,2); %支路首節(jié)點(diǎn)Nr=Branch(:,3); %支路末節(jié)點(diǎn)%節(jié)點(diǎn)數(shù)據(jù)提取N=Node(:,1); %節(jié)點(diǎn)號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); %發(fā)電機(jī)序號Nbus=Generator(:,2); %所在母線號a2=Gene

20、rator(:,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ī)臺數(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=*ones(2*ng+n+nbr,1); %

21、松弛變量u=*ones(2*ng+n+nbr,1); %松弛變量w=*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)=Generator(:,6);gmin(2*ng+1:2*ng+n)=Node(:,10);gmin(2*ng

22、+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);R=Branch(k,4);X=Branch(k,5);ban=Branch(k,6);K=Br

23、anch(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=*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)-

24、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(jj)

25、*(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); %對角元素 H1(

26、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)*(G(

27、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)=L1; %2、不等式約束的

28、雅克比矩陣% 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*ii-1,jj)=-Uamp(Nl(j

29、j) *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(Nl(jj) *Uamp(Nr(jj)

30、*(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*ng+n+nbr); U_1W=ze

31、ros(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

32、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,2*ii-1,ii)+Uamp(ii

33、)*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,ii)=-Uamp(ii)*(G(ii

34、,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,2*jj-1,ii)=Aqb(2*j

35、j-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)=-(

36、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)*(G(ii,jj)*cos(theta(ii,jj)+B(ii,jj

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

38、i); Aqb(2*jj-1,2*ii,ii)=Aqb(2*ii,2*jj-1,ii); Aqb(2*jj,2*ii-1,ii)=Aqb(2*ii-1,2*jj,ii); Aqb(2*jj,2*ii,ii)=Aqb(2*ii,2*jj,ii); end end Apb(2*ii,2*ii-1,ii)=Apb(2*ii-1,2*ii,ii); Apb(2*ii,2*ii,ii)=-2*G(ii,ii); Aqb(2*ii,2*ii-1,ii)=Aqb(2*ii-1,2*ii,ii); Aqb(2*ii,2*ii,ii)=2*B(ii,ii); end for ii=1:n A=A+Apb(:,:

39、,ii)*y(2*ii-1)+Aqb(:,:,ii)*y(2*ii); end H2(2*ng+1:2*(ng+n),2*ng+1:2*(ng+n)=A; %H3% H3=zeros(2*(ng+n),2*(ng+n); A3=zeros(2*n,2*n); Apc=zeros(2*n,2*n,nbr); for ii=1:nbr %對角線上ii Apc(2*Nl(ii)-1,2*Nl(ii)-1,ii)= -Uamp(Nl(ii)*Uamp(Nr(ii)* (G(Nl(ii),Nr(ii)*cos(theta(Nl(ii),Nr(ii)+B(Nl(ii),Nr(ii)*sin(theta(N

40、l(ii),Nr(ii); Apc(2*Nl(ii)-1,2*Nl(ii),ii)= -Uamp(Nr(ii)*(G(Nl(ii),Nr(ii) *sin(theta(Nl(ii),Nr(ii)-B(Nl(ii),Nr(ii)*cos(theta(Nl(ii),Nr(ii); Apc(2*Nl(ii),2*Nl(ii)-1,ii)= Apc(2*Nl(ii)-1,2*Nl(ii),ii); Apc(2*Nl(ii),2*Nl(ii),ii)=-2*G(Nl(ii),Nr(ii); %對角線上jj Apc(2*Nr(ii)-1,2*Nr(ii)-1,ii)=-Uamp(Nl(ii)*Uamp(Nr(ii) *(G(Nl(ii),Nr(ii)*cos(theta(Nl(ii),Nr(ii) +B(Nl(ii),Nr(ii)*sin(theta(Nl(ii),Nr(ii); Apc(2*Nr(ii)-1,2*Nr(ii),ii)= Uamp(Nl(ii)*(G(Nl(ii),N