内点法最优潮流MATLAB算法.docx
- 文档编号:29262312
- 上传时间:2023-07-21
- 格式:DOCX
- 页数:21
- 大小:17.88KB
内点法最优潮流MATLAB算法.docx
《内点法最优潮流MATLAB算法.docx》由会员分享,可在线阅读,更多相关《内点法最优潮流MATLAB算法.docx(21页珍藏版)》请在冰豆网上搜索。
内点法最优潮流MATLAB算法
内点法最优潮流MATLAB算法
clear;
%clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%数据加载
n=input('请输入要计算的节点系统(5):
')
loadNode5.txt;%节点数据
loadBranch5.txt;%支路数据
loadGenerator5.txt;%发电机数据
Node=Node5;
Branch=Branch5;
Generator=Generator5;
%节点数据处理
N=Node(:
1);%节点号
Type=Node(:
2);%节点类型
Uamp=Node(:
3);%节点电压幅值
Dlta=Node(:
4);%节点电压相角
Pd=Node(:
5);%节点负荷有功
Qd=Node(:
6);%节点负荷无功
Pg=Node(:
7);%节点出力有功
Qg=Node(:
8);%节点出力无功
Umax=Node(:
9);%节点电压幅值上限Umin=Node(:
10);%节点电压幅值下限Bc=Node(:
11);%节点补偿电容电纳值%支路数据处理
Nbr=Branch(:
1);%支路号
Nl=Branch(:
2);%支路首节点
Nr=Branch(:
3);%支路末节点
R=Branch(:
4);%支路电阻
X=Branch(:
5);%支路电抗
Z=R+1i*X;%支路阻抗=支路电阻+支路电抗Bn=Branch(:
6);%支路对地电纳
K=Branch(:
7);%支路变压器变比,0表示无变压器Ptmax=Branch(:
8);%线路传输功率上限
%发电机数据处理
Ng=Generator(:
1);%发电机序号
Nbus=Generator(:
2);%所在母线号
Pumax=Generator(:
3);%发电机有功出力上界Qumax=Generator(:
4);%发电机无功出力上界Pumin=Generator(:
5);%发电机有功出力下界Qumin=Generator(:
6);%发电机无功出力下界
a2=Generator(:
7);%燃料耗费曲线二次系数
a1=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+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)=Pumin;
gmin(ng+1:
2*ng)=Qumin;
gmin(2*ng+1:
2*ng+n)=Umin;
gmin(2*ng+n+1:
2*ng+n+nbr)=-Ptmax;%gmax
gmax=zeros(2*ng+n+nbr,1);
gmax(1:
ng)=Pumax;
gmax(ng+1:
2*ng)=Qumax;
gmax(2*ng+1:
2*ng+n)=Umax;
gmax(2*ng+n+1:
2*ng+n+nbr)=Ptmax;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成导纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y=zeros(n,n);
%%%%%%%%%%%%%%%%%%%%计算非对角元素%%%%%%%%%%%%%%%%%%%%%forii=1:
nbr
ifK(ii)==0%非变压器支路
Y(Nl(ii),Nr(ii))=-1/Z(ii);
Y(Nr(ii),Nl(ii))=Y(Nl(ii),Nr(ii));
else%变压器支路
Y(Nl(ii),Nr(ii))=-1/Z(ii)/K(ii);
Y(Nr(ii),Nl(ii))=Y(Nl(ii),Nr(ii));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%计算对角元素%%%%%%%%%%%%%%%%%%%%%%
forii=1:
n%将支路导纳加入到对角元素中
forjj=1:
nbr
ifK(jj)==0&&(Nl(jj)==ii||Nr(jj)==ii)%非变压器支路
Y(ii,ii)=Y(ii,ii)+1/Z(jj);
elseifK(jj)~=0&&(Nl(jj)==ii||Nr(jj)==ii)%变压器支路
Y(ii,ii)=Y(ii,ii)+1/Z(jj)/K(jj);
end
end
end
end
forii=1:
nbr%将对地电纳加入到对角元素中
ifK(ii)==0%非变压器支路
Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+1i*Bn(ii);
Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+1i*Bn(ii);
else%变压器支路
Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+(K(ii)-1)/K(ii)/Z(ii);
Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+(1-K(ii))/K(ii)/K(ii)/Z(ii);
end
end
forii=1:
n
Y(ii,ii)=Y(ii,ii)+i*Bc(ii);
end
G=real(Y);%电导
B=imag(Y);%电纳
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=0;%迭代次数
Kmax=150;%最大迭代次数
iteration=1e-4;%误差精度
delta=0.08;
Gap=(l'*z-u'*w)*ones(Kmax,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%whilek<50
%计算互补间隙Gap
Gap(k+1)=l'*z-u'*w;
ifGap>iteration
miu=delta*Gap(k+1)/(2*(2*ng+n+nbr));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成系数矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%相角差计算%%%%%%%%%%%%%%%%%%%%%%
theta=zeros(n,n);
forii=1:
n
forjj=1:
n
theta(ii,jj)=Dlta(ii)-Dlta(jj);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%1、等式约束雅克比矩阵%%%%%%%%%%%%%%%%
pxh=zeros(2*(ng+n),2*n);
%%%%%%%%%%%%%%%%%%%%%%%ah/aP%%%%%%%%%%%%%%%%%%%%%%%
forii=1:
ng
pxh(Ng(ii),2*Nbus(ii)-1)=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%ah/aQ%%%%%%%%%%%%%%%%%%%%%%%
forii=1:
ng
pxh(Ng(ii)+ng,2*Nbus(ii))=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%ah/ax%%%%%%%%%%%%%%%%%%%%%%%
HH=zeros(n,n);
JJ=zeros(n,n);
NN=zeros(n,n);
LL=zeros(n,n);
forii=1:
n
forjj=1:
n
ifii~=jj%i!
=j时的情况
%非对角元素
HH(ii,jj)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));
JJ(ii,jj)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));
NN(ii,jj)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));
LL(ii,jj)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));
%对角元素
HH(ii,ii)=HH(ii,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));
JJ(ii,ii)=JJ(ii,ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));
NN(ii,ii)=NN(ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));
LL(ii,ii)=LL(ii,ii)-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));
end
end
NN(ii,ii)=NN(ii,ii)-2*Uamp(ii)*G(ii,ii);
LL(ii,ii)=LL(ii,ii)+2*Uamp(ii)*B(ii,ii);
end
pxh(1+2*ng:
2:
2*(n+ng)-1,1:
2:
2*n-1)=HH';
pxh(1+2*ng:
2:
2*(n+ng)-1,2:
2:
2*n)=JJ';
pxh(2+2*ng:
2:
2*(n+ng),1:
2:
2*n-1)=NN';
pxh(2+2*ng:
2:
2*(n+ng),2:
2:
2*n)=LL';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%2、不等式约束的雅克比矩阵%%%%%%%%%%%%%%%%%%%%
%g1:
电源有功出力上下限约束
ag1aP=eye(ng,ng);
ag1aQ=zeros(ng,ng);
ag1ax=zeros(2*n,ng);
%g2:
电源无功出力上下限约束
ag2aP=zeros(ng,ng);
ag2aQ=eye(ng,ng);
ag2ax=zeros(2*n,ng);
%g3:
节点电压幅值上下限约束
ag3aP=zeros(ng,n);
ag3aQ=zeros(ng,n);
ag3ax=zeros(2*n,n);
forii=1:
n
ag3ax(2*ii,ii)=1;
end
%g4:
线路潮流上下限约束
ag4aP=zeros(ng,nbr);
ag4aQ=zeros(ng,nbr);
ag4ax=zeros(2*n,nbr);
forii=1:
n
forjj=1:
nbr
ifNl(jj)==ii
ag4ax(2*ii-1,jj)=-Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),N
r(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
ifNr(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=[ag1aPag2aPag3aPag4aP;
ag1aQag2aQag3aQag4aQ;
ag1axag2axag3axag4ax];%此即为不等式约束的雅克比矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%3、对角矩阵%%%%%%%%%%%%%%%%%%%%%%%%
L_1Z=zeros(2*ng+n+nbr,2*ng+n+nbr);
U_1W=zeros(2*ng+n+nbr,2*ng+n+nbr);
forii=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);
forii=1:
n
forjj=1:
n%元素位置为:
12
ifii~=jj%34
%对角线上与ii对应的元素
%Ap
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)));%对角线处1号元素
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)));%对角线处2号元素%%3号元素与之相等
%Aq
Aqb(2*ii-1,2*ii-1,ii)=Aqb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处1号元素
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)));%对角线处2号元素%%3号元素与之相等
%对角线上与jj对应的元素
%Ap
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)));%对角线处1号元素
Apb(2*jj-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,
jj)));%对角线处2号元素
Apb(2*jj,2*jj-1,ii)=Apb(2*jj-1,2*jj,ii);
%3号元素与2号元素相等
%Aq
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)));%对角线处1号元素
Aqb(2*jj-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%对角线处2号元素
Aqb(2*jj,2*jj-1,ii)=Aqb(2*jj-1,2*jj,ii);
%3号元素与2号元素相等
%4号元素为0
%非对角线行元素
%Ap
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)));%非对角线行处1号元素
Apb(2*ii-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处2号元素
Apb(2*ii,2*jj-1,ii)=-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处3号元素
Apb(2*ii,2*jj,ii)=-(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对
角线行处4号元素
%Aq
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)));%非对角线行处1号元素
Aqb(2*ii-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处2号元素
Aqb(2*ii,2*jj-1,ii)=Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处3号元素
Aqb(2*ii,2*jj,ii)=-(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对
角线行处4号元素
%非对角线列元素
%Ap
Apb(2*jj-1,2*ii-1,ii)=Apb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素
Apb(2*jj-1,2*ii,ii)=Apb(2*ii,2*jj-1,ii);%非对角线列处2号元素
Apb(2*jj,2*ii-1,ii)=Apb(2*ii-1,2*jj,ii);%非对角线列处3号元素
Apb(2*jj,2*ii,ii)=Apb(2*ii,2*jj,ii);%%非对角线列处4号元素
%Aq
Aqb(2*jj-1,2*ii-1,ii)=Aqb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素
Aqb(2*jj-1,2*ii,ii)=Aqb(2*ii,2*jj-1,ii);%非对角线列处2号元素
Aqb(2*jj,2*ii-1,ii)=Aqb(2*ii-1,2*jj,ii);%非对角线列处3号元素
Aqb(2*jj,2*ii,ii)=Aqb(2*ii,2*jj,ii);%%非对角线列处4号元素
end
end
%对角线上与ii对应的元素
%Ap
Apb(2*ii,2*ii-1,ii)=Apb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等
Apb(2*ii,2*ii,ii)=-2*G(ii,ii);%对角线处4号元素
%Aq
Aqb(2*ii,2*ii-1,ii)=Aqb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等
Aqb(2*ii,2*ii,ii)=2*B(ii,ii);%对角线处4号元素
end
forii=1:
n
A=A+Apb(:
:
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);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 内点法 最优 潮流 MATLAB 算法
![提示](https://static.bdocx.com/images/bang_tan.gif)