clear all;clc;
% clf;
y=[0-66.66667i 0 0+63.49206i 0 0;
0 0-33.33333i 0 0+31.74603i 0;
0+63.49206i 0 1.45390-66.98082i -0.82987+3.11203i -0.62402+3.90015i;
0 0+31.74603i -0.82987+3.11203i 1.58459-35.73786i -0.75471+2.64150i;
0 0 -0.62402+3.90015i -0.75471+2.64150i 1.37874-6.29166i] ;%导纳矩阵
g=real(y);
b=imag(y);
u(:,1)=[1.05;1.05;1;1;1]; %各节点初始电压值
ct(:,1)=[0;0;0;0;0]; %初始角度
p0=[1;0;-0.4;-0.74;-0.32]; %各节点的输入有功功率
q0=[0;0;-0.2;-0.26;-0.16]; %各节点的输入无功功率
x(:,1)=[0;1;0;1;0;1;0]; %控制变量的初始值
e=0.0001; %迭代精度
zu(:,1)=[1;1;1;1.05]; %最优控制变量解对应为P1、U1、U2、T
for k=1:1000
for i=1:5
a=0;
for j=1:5
a=u(j,k)*(g(i,j)*cos(ct(i,k)-ct(j,k))+b(i,j)*sin(ct(i,k)-ct(j,k)))+a;
end
dp(i,k)=p0(i)-u(i,k)*a;
end
for i=3:5
z=0;
for j=1:5
z=u(j,k)*(g(i,j)*sin(ct(i,k)-ct(j,k))-b(i,j)*cos(ct(i,k)-ct(j,k)))+z;
end
dq(i,k)=q0(i)-u(i,k)*z;
end
for i=1:5
for j=1:5
if i~=j
hh(i,j)=u(i,k)*u(j,k)*(g(i,j)*sin(ct(i,k)-ct(j,k))-b(i,j)*cos(ct(i,k)-ct(j,k)));
jj(i,j)=-u(i,k)*u(j,k)*(g(i,j)*cos(ct(i,k)-ct(j,k))+b(i,j)*sin(ct(i,k)-ct(j,k)));
nn(i,j)=u(i,k)*u(j,k)*(g(i,j)*cos(ct(i,k)-ct(j,k))+b(i,j)*sin(ct(i,k)-ct(j,k)));
ll(i,j)=u(i,k)*u(j,k)*(g(i,j)*sin(ct(i,k)-ct(j,k))-b(i,j)*cos(ct(i,k)-ct(j,k)));
else c=0;
d=0;
for m=1:5
if m~=i
c=u(m,k)*(g(i,m)*sin(ct(i,k)-ct(m,k))-b(i,m)*cos(ct(i,k)-ct(m,k)))+c;
d=u(m,k)*(g(i,m)*cos(ct(i,k)-ct(m,k))+b(i,m)*sin(ct(i,k)-ct(m,k)))+d;
end
end
hh(i,i)=-u(i,k)*c;
jj(i,i)=u(i,k)*d;
nn(i,i)=u(i,k)*d+2*u(i,k)^2*g(i,i);
ll(i,i)=u(i,k)*c-2*u(i,k)^2*b(i,i);
end
end
end
jkb=[hh(3,3) nn(3,3) hh(3,4) nn(3,4) hh(3,5) nn(3,5) hh(3,1);
jj(3,3) ll(3,3) jj(3,4) ll(3,4) jj(3,5) ll(3,5) jj(3,1);
hh(4,3) nn(4,3) hh(4,4) nn(4,4) hh(4,5) nn(4,5) hh(4,1);
jj(4,3) ll(4,3) jj(4,4) ll(4,4) jj(4,5) ll(4,5) jj(4,1);
hh(5,3) nn(5,3) hh(5,4) nn(5,4) hh(5,5) nn(5,5) hh(5,1);
jj(5,3) ll(5,3) jj(5,4) ll(5,4) jj(5,5) ll(5,5) jj(5,1);
hh(1,3) nn(1,3) hh(1,4) nn(1,4) hh(1,5) nn(1,5) hh(1,1)]; %牛拉法的雅各比矩阵
dd=[dp(3,k);dq(3,k);dp(4,k);dq(4,k);dp(5,k);dq(5,k);dp(1,k)];
dx(:,k)=inv(jkb)*dd;
xzh=[1 0 0 0 0 0 0 ;0 u(3,k) 0 0 0 0 0;0 0 1 0 0 0 0;0 0 0 u(4,k) 0 0 0;0 0 0 0 1 0 0;0 0 0 0 0 u(5,k) 0;0 0 0 0 0 0 1];
x(:,k+1)=x(:,k)+xzh*dx(:,k); %潮流计算解
u(:,k+1)=[u(1,k);u(2,k);x(2,k+1);x(4,k+1);x(6,k+1)];
ct(:,k+1)=[x(7,k+1);0;x(1,k+1);x(3,k+1);x(5,k+1)];
f1=0;
for i=1:5
f1=u(j,k+1)*(g(1,j)*cos(ct(1,k+1)-ct(j,k+1))+b(1,j)*sin(ct(1,k+1)-ct(j,k+1)))+f1;
end
pg1=u(1,k+1)*f1;
f2=0;
for i=1:5
f2=u(j,k+1)*(g(2,j)*cos(ct(2,k+1)-ct(j,k+1))+b(2,j)*sin(ct(2,k+1)-ct(j,k+1)))+f2;
end
pg2=u(2,k+1)*f2;
fdx=[(702*pg1+100)*u(1,k+1)*u(3,k+1)*(g(1,3)*sin(ct(1,k+1)-ct(3,k+1))-b(1,3)*cos(ct(1,k+1)-ct(3,k+1)))+(778*pg2+100)*u(2,k+1)*u(3,k+1)*(g(2,3)*sin(ct(2,k+1)-ct(3,k+1))-b(2,3)*cos(ct(2,k+1)-ct(3,k+1)));
(702*pg1+100)*u(1,k+1)*(g(1,3)*cos(ct(1,k+1)-ct(3,k+1))+b(1,3)*sin(ct(1,k+1)-ct(3,k+1)))+(778*pg2+100)*u(2,k+1)*(g(2,3)*cos(ct(2,k+1)-ct(3,k+1))+b(2,3)*sin(ct(2,k+1)-ct(3,k+1)));
(702*pg1+100)*u(1,k+1)*u(4,k+1)*(g(1,4)*sin(ct(1,k+1)-ct(4,k+1))-b(1,4)*cos(ct(1,k+1)-ct(4,k+1)))+(778*pg2+100)*u(2,k+1)*u(4,k+1)*(g(2,4)*sin(ct(2,k+1)-ct(4,k+1))-b(2,4)*cos(ct(2,k+1)-ct(4,k+1)));
(702*pg1+100)*u(1,k+1)*(g(1,4)*cos(ct(1,k+1)-ct(4,k+1))+b(1,4)*sin(ct(1,k+1)-ct(4,k+1)))+(778*pg2+100)*u(2,k+1)*(g(2,4)*cos(ct(2,k+1)-ct(4,k+1))+b(2,4)*sin(ct(2,k+1)-ct(4,k+1)));
(702*pg1+100)*u(1,k+1)*u(5,k+1)*(g(1,5)*sin(ct(1,k+1)-ct(5,k+1))-b(1,5)*cos(ct(1,k+1)-ct(5,k+1)))+(778*pg2+100)*u(2,k+1)*u(5,k+1)*(g(2,5)*sin(ct(2,k+1)-ct(5,k+1))-b(2,5)*cos(ct(2,k+1)-ct(5,k+1)));
(702*pg1+100)*u(1,k+1)*(g(1,5)*cos(ct(1,k+1)-ct(5,k+1))+b(1,5)*sin(ct(1,k+1)-ct(5,k+1)))+(778*pg2+100)*u(2,k+1)*(g(2,5)*cos(ct(2,k+1)-ct(5,k+1))+b(2,5)*sin(ct(2,k+1)-ct(5,k+1)));
(702*pg1+100)*u(1,k+1)*u(1,k+1)*(g(1,1)*sin(ct(1,k+1)-ct(1,k+1))-b(1,1)*cos(ct(1,k+1)-ct(1,k+1)))+(778*pg2+100)*u(2,k+1)*u(1,k+1)*(g(2,1)*sin(ct(2,k+1)-ct(1,k+1))-b(2,1)*cos(ct(2,k+1)-ct(1,k+1)))];
lmd=-inv((jkb)')*fdx;
aa=g(2,1)*cos(ct(2,k+1)-ct(1,k+1))+b(2,1)*sin(ct(2,k+1)-ct(1,k+1));
tdf=[702*pg1+100;878*pg2*u(2,k+1)*aa;878*(u(1,k+1)*aa+2*u(2,k+1)*g(2,2));0];
if norm(tdf)<=e %精度判别
pg1=zu(1,k)
u1=zu(2,k)
u2=zu(3,k)
T=zu(4,k)
break;
end
bch=1.2; %步长因子
zu(:,k+1)=zu(:,k)-bch*tdf;
if zu(1,k+1)>=1.2 %不等式约束条件
zu(1,k+1)=1.2;
elseif zu(1,k+1)<=0.3
zu(1,k+1)=0.3;
end
zu(1,k+1)=1;
if zu(2,k+1)>=1.1
zu(2,k+1)=1.1;
elseif zu(2,k+1)<=1.0
end
zu(2,k+1)=1.02;
if zu(3,k+1)>=1.1
zu(3,k+1)=1.1;
elseif zu(3,k+1)<=1.0
zu(3,k+1)=1.0;
end
zu(3,k+1)=1.02;
u(1,k+1)=zu(2,k+1);
u(2,k+1)=zu(3,k+1);
p0(1)=zu(1,k+1);
end
u=u(:,k+1) %输出结果
ct=ct(:,k+1)
x=x(:,k+1)
zu=zu(:,k)
用简化梯度法进行电力系统最优潮流计算_matlab源码
版权申诉
5星 · 超过95%的资源 173 浏览量
2022-03-19
21:48:06
上传
评论 1
收藏 2KB ZIP 举报
阿里matlab建模师
- 粉丝: 3285
- 资源: 2783
最新资源
- IMG_0694.GIF
- 基于图像的三维模型重建C++源代码+文档说明(高分课程设计)
- 基于聚焦法的工件立体测量方案,根据数据进行三维重建 使用HALCON处理图像,MATLAB拟合数据+源代码+数据集+效果图
- 锄战三国村 修改:货币使用不减 v1.10(2) 原创 (中文).apk
- 基于python实现的单目双目视觉三维重建+源代码+图像图片(高分课程设计)
- 基于C+++OPENCV的全景图像拼接源码(课程设计)
- 基于Python+OpenCV对多张图片进行全景图像拼接,消除鬼影,消除裂缝+源代码+文档说明+界面截图(高分课程设计)
- 基于C++实现的全景图像拼接源码(课程设计)
- 基于SIFT特征点提取和RASIC算法实现全景图像拼接python源码+文档说明+界面截图+详细注释(95分以上课程大作业)
- 基于matlab实现眼部判别的疲劳检测系统+源代码+全部数据+文档说明+详细注释+使用说明+截图(高分课程设计)
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
- 1
- 2
前往页