基于内点法的最优潮流计算

1、摘要内点法是一种能在可行域内部寻优的方法，即从初始内点出发，沿着中心路径方向在可行域内部直接走向最优解的方法。其中路径跟踪法是目前最具有发展潜力的一类内点算法，该方法鲁棒性强，对初值的选择不敏感，在目前电力系统优化问题中得到了广泛的应用。本文采用路径跟踪法进行最优求解，首先介绍了路径跟踪法的基本模型，并且结合具体算例，用编写的Matlab程序进行仿真分析，验证了该方法在最优潮流计算中的优越性能。关键词：最优潮流、内点法 、路径跟踪法、仿真目 次0、引言11、路径跟踪法的基本数学模型22、路径跟踪法的最优潮流求解思路33、具体算例及程序实现流程63.1、算例描述63.2、程序具体实现流程74、运

2、行结果及分析1141 运行结果1142结果分析165、结论176、编程中遇到的问题17参考文献19附录200、引言电力系统最优潮流，简称OPF（Optimal Power Flow）。OPF问题是一个复杂的非线性规划问题，要求满足待定的电力系统运行和安全约束条件下，通过调整系统中可利用控制手段实现预定目标最优的系统稳定运行状态。针对不同的应用，OPF模型课以选择不同的控制变量、状态变量集合，不同的目标函数，以及不同的约束条件，其数学模型可描述为确定一组最优控制变量u，以使目标函数取极小值，并且满足如下等式和不等式。minufx,uS.t.hx,u=0gx,u0 （0-1）其中minufx,u为

3、优化的目标函数，可以表示系统运行成本最小、或者系统运行网损最小。S.t.hx,u=0为等式约束，表示满足系统稳定运行的功率平衡。gx,u0为不等式约束，表示电源有功出力的上下界约束、节点电压上下线约束、线路传输功率上下线约束等等。电力系统最优潮流算法大致可以分为两类：经典算法和智能算法。其中经典算法主要是指以简化梯度法、牛顿法、内点法和解耦法为代表的基于线性规划和非线性规划以及解耦原则的算法，是研究最多的最优潮流算法, 这类算法的特点是以一阶或二阶梯度作为寻找最优解的主要信息。智能算法主要是指遗传算法和模拟退火发等，这类算法的特点是不以梯度作为寻优信息，属于非导数的优化方法。因此经典算法的优点

4、是能按目标函数的导数信息确定搜索方向，计算速度快，算法比较成熟，结果可信度高。缺点是对目标函数及约束条件有一定的限制，可能出现局部极小时难以收敛。而智能算法的优点是计算与导数无关，灵活性高，随机性强，缺点是算法不稳定，结果不可信，并且控制参数需凭经验给出。通过对这些常见算法的简单比较，内点法具有其优越的性能，特别是路径跟踪法，其算法收敛迅速，鲁棒性强，对初值的选择不敏感，其迭代次数与系统规模或控制变量的数目关系不大，因此本文采用该方法进行最优计算。1、 路径跟踪法的基本数学模型内点法最初的基本思路是希望通过寻优迭代过程始终在可行域内进行，因此，初始点应在可行域内，并在可行域的边界设置障碍使迭代

5、点接近边界时其目标函数迅速增大，从而保证迭代点均在可行域的内点。但是对于大规模实际问题而言，寻找初始点往往十分困难。为此许多学者长期以来致力于内点算法初始内点条件的改进。以下介绍的路径跟踪法只要求在寻优过程中松弛变量和拉格朗日乘子满足简单的大于0或者小于0的条件，即可代替原来必须在可行域内求解的要求，使得计算过程大为简化。一般可以将最优潮流模型简化为如下的非线性优化模型。Obj. minufx,u （1-1） s.t. S.t.hx,u=0 （1-2） g-gx,ug- （1-3）其中minufx,u为优化的目标函数， S.t.hx,u=0为等式约束， gx,u为不等式约束，路径跟踪内点法的基

6、本思路是：首先将式（1-3）的不等约束变成等式约束：gx,u+u=g- （1-4）gx,u-l=g- （1-5）其中松弛变量 l=l1,lrT， u=u1,urT，应满足u0，l0这样原问题就转化为问题A：Obj. minufx,u S.t. 然后，把目标函数改造成障碍函数，该函数在可行域内应接近于原函数f(x)，而在边界时变得很大。一次可得带优化问题B：obj. s.t. 其中扰动因子或者障碍因子u0。当l或u接近边界时，以上函数将趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到。这样就通过目标函数的变化把含不等式限制的优化问题A变成只含等式限制优化的问题B了，因此可以直接用拉

7、格朗日乘子法来求解。优化问题B的拉格朗日目标函数为： 式中：,和均为拉格朗日乘子。因此最后简化的求解问题就是求取上述表达式的极小解。2、 路径跟踪法的最优潮流求解思路路径跟踪法的最优潮流求解过程就是对拉格朗日目标函数求极小值问题：式中：,和均为拉格朗日乘子。该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0。即： （2-1）通过上述表达式可以解出： （2-2）定义：，称为互补间隙。可得： （2-3）如果x*是优化问题A的最优解，当u固定时，x(u)是优化问题B的解，那么当Gap0，u0时，产生的序列x(u)收敛至x*。建议采用：。式中称为中心参数，一般取0.1，在大多数场合可

8、获得较好的收敛效果。通过偏导数为0的表达式可以可得内点法的修正方程为： （2-4）求解方程可得到第k次迭代的修正量，于是最优解的一个新的近似解为： （2-5）式中，和为步长： （2-6）其潮流计算的流程图如下图1所示，其中初始化部分包括：(1)、设置松弛变量l和u，保证l，uT0(2)、设置拉格朗日乘子w、y、z，满足w0，Y！=0T (3)、设置优化问题的初值。(4)、取中心参数，给定计算精度，迭代次数初值K=0。输出最优解，停止计算计算互补间隙Gap计算扰动因子计算和更新原始变量及拉格朗日乘子kKmaxGap求解修正方程，求出，输出“计算不收敛”初始化是是否否图1. 内点法潮流计算流程图3

9、、 具体算例及程序实现流程这部分主要有算例描述以及程序的实现流程两部分，其中算例描述主要是对系统参数以及优化问题进行说明。而程序的实现流程主要描述的是最优潮流计算中所涉及的矩阵方程的描述。3.1、算例描述该算例为王锡凡编写的现代电力系统分析中的3-1的例题，是以系统燃料最省为最优潮流的目标函数。选择该题目作为算例分析的原因是，该题目有比较详细的解题思路以及列写出了比较详细的迭代结果，方便对编写程序的运行结果进行比对。求如下图所示简化系统的系统燃料最省的最优潮流计算：12243251:1.051.05:1j0.015j0.03j0.25j0.250.08+j0.300.1+j0.350.04+j

10、0.25j0.25j0.253.7+j1.32+j11.6+j0.8除了由上图所提供的系统母线负荷功率数据、线路参数和变压器之路参数数据、变压器便比数据之外，以下顺序给出了线路传输功率边界（表3-1），发电机有功无功出力上下界和燃料耗费曲线 参数（表3-2）。若不作特殊说明，所有数据都是以标幺值形式给出，功率基准值为100MVA，母线电压上下界分别为1.1和0.9。表3-1线路传输功率边界支路号首末端母线号线路传输功率边界11-2221-30.6532-3242-4653-55表3-2 发电机数据发电机序号母线号出力上界出力下界燃料耗费曲线参数有功无功有功无功二次系数一次系数常数14831-3

11、50.4395200.43551200.6425851-2.1200.55500.7451857.203.2、程序具体实现过程针对上述系统，首先我们先列写出该算例的数学模型和有关计算公式。在该算例中，共有5个节点，相应的状态量为：系统中有2台发电机，没有其他无功源，因此控制变量为：应该指出，此处发电机和无功源的编号与及诶单编号无关，是独立编号的。这是因为系统中一个节点可能接有多台发电机的缘故。因此系统中总变量共有14个：最优潮流的数学模型为：目标函数：约束条件：每个节点有两个潮流方程，共有10个等式约束条件，对非发电机而言： （i=1,2,3）对发电机节点： （i=4,5）式中：表示第k台发电

12、机接在节点i上。不等式约束共有14个条件，分别是： 根据以上模型可以形成修正方程。该方程包括形成等式左边的系数矩阵和等式右边的常数项两部分。1、形成系数矩阵1）、等式约束的雅克比矩阵等式右端包括3个子矩阵：其中：其中：式中：i为发电机的序号；j为节点号；表示第i台发电机是在节点j上的。（潮流计算中的雅克比矩阵）2）、不等式约束的雅克比矩阵式中：、和依次表示电源有功出力的上下界约束，无功电源出力的上下界约束，节点电压赋值的上下界约束和线路潮流约束。，式中：第行列元素为1，其他元素均为0。3）、对角矩阵4）、海森伯矩阵这是最复杂的部分，共包含四项。有上述推导已经可以得到其中的第四项为：其余三项是：

13、目标函数的海森伯矩阵、等式约束海森伯矩阵与拉格朗日乘子y的乘积和不等式约束海森伯矩阵与拉格朗日乘子的乘积。2、形成常数项均可根据定义直接求得。可以表示为：当知道目标函数梯度矢量之后，再根据以上等式和不等式约束的雅克比矩阵公式就可以求得。4、 运行结果及分析41 运行结果以下对该算例的寻优过程用数字加以说明，设4、5节点发电机均能有算法调节其出力。在初始化过程中各变量初值根据实际问题自行设置的，我们给出所用变量的处置如下：节点电压；平衡节点；发电机出力有功取其边界值；松弛因子，当收敛条件时，需要迭代进行23次（例题所给出的迭代次数为17次）。表 4-1 迭代过程中各节点电压增量的变化情况迭代次数

14、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.-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.-

15、0.-0.00473-0.-0.180.-0.002230.-0.-0.190.-0.0.-0.-0.200.-0.0.-0.-0.210.-0.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 迭代过程中各节点相角增量的变化情况迭代次数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.1472

16、35.29.7-0.-0.-0.30025-0.-0.8-0.-0.-0.3624-0.07659-0.9-0.0.-0.271170.-0.106054.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.65

17、04-6.62907-6.19-0.06977-0.-0.00065-0.-0.20-0.00417-0.-0.00209-0.00978-0.210.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 迭代过程中有功源有功、无功源无功增量的变化情况迭代次数有功源有功出力增量无功源无功出力增量16.0.0.2.21.0.3.-3.30.-0.-11

18、.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.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-

19、1.17E-0423-2.60E-062.54E-063.21E-08-1.69E-07将各次迭代过程中Gap变化情况绘制成曲线，可以显示出路劲跟踪法最优潮流计算的收敛特性，如图4-1所示：图 4-1 5节点系统最优潮流内点法收敛特性图4-2为5节点系统最优潮流计算结果截图，其中包括迭代次数、燃料总费用、发电机有功无功出力、各节点电压幅值与相角、以及各支路有功功率。（注：结果中的值为标幺值，功率的基准值为100MVA）42结果分析将最优潮流计算的结果和普通潮流计算结果进行比较，其中PF表示为普通潮流计算。普通潮流计算中，发电机不会调节其出力。即4节点为PQ节点，5节点为平衡节点。见表4-4和表

20、4-5。从表中可以看出，由于4机组比5机组的燃料曲线系数小，因此4机组有功出力增加，5机组有功出力减少。同时系统的网损、无功功率都有所增加。这是由于要将1节点电压抬高至其下界以满足不等式约束的要求而产生的副作用。但是网损的增加并不影响目标函数的优化。整个系统的燃料费用与不优化的潮流计算相比仍然减少了243.76$。表4-4 各有功源有功和无功源无功出力发电机序号母线序号有功出力无功出力燃料费用/$OPFPFOPFPFOPFPF145.50565.001.77801.883113833.063463.80252.15682.57942.61942.29943870.134483.15总计7.66

21、247.57944.39744.13057703.207946.95表4-5 各节点电压向量母线序号电压幅值电压相角/radOPFPFOPFPF10.90000.8822-0.007-0.0834021.10001.07790.40490.3116031.08181.0364-0.0571-0.747341.06971.05000.478670.3116051.1001.0500005、 结论路径跟踪法在电力系统中应用与求解线性规划和非线性规划模型中，还是比较有优势的，具有算法收敛迅速，鲁棒性强，对初值的选择不敏感，其迭代次数与系统规模或控制变量的数目关系不大等特点。即在该简化模型中迭代次数为

22、23次，但是由于其迭代次数是与系统规模关系不大，对IEEE30、IEEE118节点系统的计算结果其迭代次数始终保持在21到27次之间。另外将水火电最优潮流问题分解为火电最优潮流子问题和水电子问题，提供了有效的协调算法。但是对于路径跟踪法的影响因素还是比较多的，比如初始点的选择、迭代步长的选取、壁垒参数的调整、离散变量的处理等等，如果选取不当可能会出现不恰当的结果，因此还需要研究者们做大量的工作。6、 编程中遇到的问题上图为例题所给出的迭代次数与仿真结果，可以看出除了迭代次数与书中不一致以外，迭代的结果基本上完全一样。但在程序的实现过程中，很难按照书上的流程编写出结果一致的程序，因为程序中有大量

23、的矩阵计算，而且也不能将书上列写的矩阵直接翻译成MATLAB语言，需要进行一些不同的处理方式，所以需要在网上去寻找一些算法的实现方式。另外在编程中发现，王锡凡的电力系统一书中，也有一些公式书写有误的现象，所以需要对公式进行一定的验证推导，因此很大程度上会出现问题。另外对于路径跟踪法来说，其初始点的选择、迭代步长的选取、壁垒参数的调整等都对计算结果又一定的影响。而在本程序中，我选取的初始点以及拉格朗日因子也都与例题所提供的有些许不同，如果选择和书中相同的参数，计算结果就会出现问题，这其中的原因，我的猜测是可能在用matalb语言实现过程中，细节方面可能与例题所展示的有所出入，由于时间关系，暂时还

24、未找到原因。参考文献1 张伯明.高等电力网络分析M.清华大学出版社,2007.2 王锡凡.现代电力系统分析M. 北京科学出版社,2003 3 张江红.最优潮流算法综述J.华北电力，2010,074 赫玉国.一种基于KarmarKar内点法的最优潮流算法J.中国机电工程学报，1996,11附录clcclear%读取计算参数%Branch=load(Branch.txt); %支路参数文档Node=load(Node.txt); %节点参数文档Generator=load(Generator.txt); %发电机参数文档%支路数据提取Nbr=Branch(:,1); %支路号Nl=Branch(:

25、,2); %支路首节点Nr=Branch(:,3); %支路末节点%节点数据提取N=Node(:,1); %节点号Type=Node(:,2); %节点类型Uamp=Node(:,3); %节点电压幅值Dlta=Node(:,4); %节点电压相角Pd=Node(:,5); %节点负荷有功Qd=Node(:,6); %节点负荷无功Pg=Node(:,7); %节点出力有功Qg=Node(:,8); %节点出力无功%发电机数据提取Ng=Generator(:,1); %发电机序号Nbus=Generator(:,2); %所在母线号a2=Generator(:,7); %燃料耗费曲线二次系数a1

26、=Generator(:,8); %燃料耗费曲线一次系数a0=Generator(:,9); %燃料耗费曲线常数项%计算参数初始化%n=length(N); %节点个数ng=length(Ng); %发电机台数nbr=length(Nbr); %之路个数x=zeros(2*(ng+n),1); %控制变量+状态量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+

27、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; %拉格朗日乘子%计算不等式约束的上下限%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

28、)=-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); %生成导纳矩阵%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

29、,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;%迭代次数%主程序%while k1e-3 miu=0.1*Gap(k+1)/(2*(2*ng+n+nbr); %形成系数矩阵% theta=zeros(n,n); for ii=1:n for jj=1:n theta(ii,jj)=Dlta(ii)-Dlta(jj); en

30、d 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)*(G(ii,jj) *

31、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(ii,ii)=H1(ii

32、,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(ii,jj) *sin(

33、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、不等式约束的雅克比矩阵% agaP=

34、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(jj) *Uamp(Nr(

35、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)*(G(Nl(jj),N

36、r(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=zeros(2*ng+n+n

37、br,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 j

38、j=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)*Uamp(jj)*(

39、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,jj)* sin(th

40、eta(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*jj-1,2*jj,ii)

41、; 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(i