matlab代码参考，潮流计算（不全对）

xiaoxiaoyang_

已于 2022-04-03 22:23:16 修改

阅读量1.1k

点赞数

文章标签： matlab

于 2022-04-03 22:22:00 首次发布

本文链接：https://blog.csdn.net/xiaoxiaoyang_/article/details/123946014

版权

2022.4.3 找到去年编写的潮流计算代码，不完全对，当作matlab笔记用。
% Y为导纳矩阵
% ykb为雅可比矩阵
% F为不平衡量
% X2存的每次的修正量
clc,clear
disp(‘电力系统极坐标下的牛顿-拉夫逊潮流计算：’)

load Data.mat
Jiedian_Data1=Jiedian_Data;
Zhilu_Data1=Zhilu_Data;
%% 说明
%支路变比不为1时，p1是1次侧节点编号；p2是k次测节点编号；
%使用了两个结构体数组。节点 jd(i),支路zl(i).

% %%
% n=input(‘请输入节点数： n=’);
% pq=input(‘请输入PV节点数：pv=’);
% pv=input(‘请输入PQ节点数：pq=’);
% esp=(‘请输入精确度：esp’);
n=9;pq=6;pv=2;esp=1e-5; %这里直接输入了。
%% 输入节点数据
A=(Jiedian_Data1); %把节点数据赋给A
[m1,n1]=size(Jiedian_Data1); %节点数据的行列数
jd=struct(‘num’,0,‘ty’,0,‘U’,0,‘P’,0,‘Q’,0,‘zkj’,0);
for i=1:m1
jd(i).num=A(i,1); %节点编号
jd(i).ty=A(i,2); %节点类型号 1为PQ; 2为PV ；3为平衡
jd(i).U=A(i,3); %节点电压
jd(i).P=A(i,4); %节点有功功率
jd(i).Q=A(i,5); %无功功率
jd(i).zkj=A(i,6); %阻抗角
end
%%输出节点数据
for i=1:m1
jd(i)
% jd(i).num
% jd(i).ty
% jd(i).U
% jd(i).P;
% jd(i).Q
% jd(i).zkj
end
% 输入支路数据
B=Zhilu_Data1;
[m2,n2]=size(Zhilu_Data1);
z1=struct(‘num’,0,‘p1’,0,‘p2’,0,‘k’,0,‘Z’,0,‘y’,0);
for i=1:m2
zl(i).num=B(i,1);
zl(i).p1=B(i,2);
zl(i).p2=B(i,3);
zl(i).k=B(i,4);
zl(i).Z=B(i,5); %支路阻抗
zl(i).y=B(i,6);%支路对地导纳
end
%% 生成节点导纳矩阵 Y
Y=zeros(m1,m1); %m1为节点数
%变压器支路线路上是 yt/k; 即 1/(zl(i).Zzl(i).k)
%1次测对地导纳为 (k-1)/kyt
%k测对地导纳为(1-k)/(KK)yt
for i=1:m1
if(zl(i).k==1)
Y(zl(i).p1,zl(i).p1) = Y(zl(i).p1,zl(i).p1)+1/zl(i).Z+zl(i).y;
Y(zl(i).p2,zl(i).p2) = Y(zl(i).p2,zl(i).p2)+1/zl(i).Z+zl(i).y;
Y(zl(i).p1,zl(i).p2) = Y(zl(i).p1,zl(i).p2) - 1/zl(i).Z;%互导
Y(zl(i).p2,zl(i).p1) = Y(zl(i).p1,zl(i).p2);
else
if (zl(i).k~=1) %处理变压器支路
yt=1/zl(i).Z;
%yt/zl(i).k^2;
Y(zl(i).p1,zl(i).p1) = Y(zl(i).p1,zl(i).p1)+yt; %1次侧
Y(zl(i).p2,zl(i).p2) = Y(zl(i).p2,zl(i).p2)+yt/(zl(i).k^2); %k测
Y(zl(i).p1,zl(i).p2) = Y(zl(i).p1,zl(i).p2) - yt/zl(i).k;%互导
Y(zl(i).p2,zl(i).p1) = Y(zl(i).p1,zl(i).p2);
end
end
end
Y %节点导纳矩阵
%% 对电压、角度赋初值
%初值在数据里已经有
%角度赋0；电压赋1
%% 计算雅可比矩阵
n=9;
Qi=0;
Pi=0;
for i=1:n
for j=1:n
%求H
if((i<=8)&&(j<=8))
if(i~=j)
H(i,j) = -jd(i).Ujd(j).U(real(Y(i,j))sin(jd(i).zkj-jd(j).zkj)…
-imag(Y(i,j))cos(jd(i).zkj-jd(j).zkj));
else if(i==j)
for k=1:n
Qi=Qi+jd(k).U(real(Y(i,k))sin(jd(i).zkj-jd(k).zkj)-imag(Y(i,k))cos(jd(i).zkj-jd(k).zkj));
end
Qi=jd(i).UQi;
H(i,i)=Qi+imag(Y(i,i))jd(i).U^2;
end
end
end
%求N 86
if((i<=8)&&(j<=6))
if(i~=j)
N(i,j) = -jd(i).Ujd(j).U(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)…
+imag(Y(i,j))sin(jd(i).zkj-jd(j).zkj));
else if(i==j)
for k=1:n
Pi=Pi+jd(k).U(real(Y(i,k))*cos(jd(i).zkj-jd(k).zkj)+imag(Y(i,k))sin(jd(i).zkj-jd(k).zkj));
end
Pi=jd(i).UPi;
N(i,i)=-Pi-real(Y(i,i))*jd(i).U^2;
end
end
end

  %求K 6*8
  if((i<=6)&&(j<=8))
      if(i~=j)
            K(i,j) = jd(i).U*jd(j).U*(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)...
         +imag(Y(i,j))*sin(jd(i).zkj-jd(j).zkj));
      else if(i==j)
             for k=1:n
                 Pi=Pi+jd(k).U*(real(Y(i,k))*cos(jd(i).zkj-jd(k).zkj)+imag(Y(i,k))*sin(jd(i).zkj-jd(k).zkj));
             end
             Pi=jd(i).U*Pi;
             K(i,i)=-Pi+real(Y(i,i))*jd(i).U^2;
          end
      end
  end
  
  %L 6*6
  if((i<=6)&&(j<=6)) 
     if(i~=j)
         L(i,j) = -jd(i).U*jd(j).U*(real(Y(i,j))*sin(jd(i).zkj-jd(j).zkj)...
         -imag(Y(i,j))*cos(jd(i).zkj-jd(j).zkj));
     else if(i==j)
             for k=1:n
                 Qi=Qi+jd(k).U*(real(Y(i,k))*sin(jd(i).zkj-jd(k).zkj)-imag(Y(i,k))*cos(jd(i).zkj-jd(k).zkj));
             end
             Qi=jd(i).U*Qi;
             L(i,i)=-Qi+imag(Y(i,i))*jd(i).U^2;
         end
     end
  end 
  
        
end

end
H, N,K,L
ykb=[H,N;K,L];
ykb1=ykb^-1;
rank(ykb1)
%% 求解不平衡量，修正方程f(x);△P △Q
F=zeros(2*pq+pv,1); %14行，1列
for i=1:n
if(i<=pq) %求pq节点的△P △Q
%求pq节点的6个 △P F(x)的前8行

          Pi=0; %对每一个节点开始Pi赋值为0
       for j=1:n  %循环所有节点求和，得到  Pi/Ui
           Pi=Pi+jd(j).U*(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)+imag(Y(i,j))*sin(jd(i).zkj-jd(j).zkj));
       end
       Pi=Pi*jd(i).U;
       F(i)=jd(i).P/100-Pi;  % P/100是标幺值
       
       
       %求pq节点的6个 △Q F(x)的第 9-14行
       Qi=0;
       for j=1:n
           Qi=Qi+jd(j).U*(real(Y(i,j))*sin(jd(i).zkj-jd(j).zkj)-imag(Y(i,j))*cos(jd(i).zkj-jd(j).zkj));
       end
       Qi=Qi*jd(i).U;
       F(i+pq+pv)=jd(i).Q/100-Qi;   %第9-14行 Q/100 是进行标幺
end
if ( (pq<i)&&(i<=pq+pv) ) %PV节点
       Pi=0; %对每一个节点开始Pi赋值为0
       for j=1:n  %循环所有节点求和，得到  Pi/Ui
           Pi=Pi+jd(j).U*(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)+imag(Y(i,j))*sin(jd(i).zkj-jd(j).zkj));
       end
       disp('Pi Fi=:');
       Pi=Pi*jd(i).U;
       F(i)=jd(i).P/100-Pi   %第7-8行   
end

end
F
%% 算 △θ △U F=ykbX
% 存一个Ui矩阵 [U];
for i= 1:pq %pq 节点的U
Upq(i)=jd(i).U;
end
%Upq=Upq’;
% 存完。
X=zeros(2pq+pv,1);
disp(‘X(0)=’);
X=(ykb^-1)F % ykb 1414 F 14*1 或者X2=ykb\F

%% 迭代
disp(‘开始迭代’);
a=0; %设置迭代次数
X1=[X(1:8,:);X(9:14).Upq] %X1为角度、电压增量
disp(‘X1=’);disp(max(X1))
disp(‘进入循环’);
while(max(abs(X1))>esp)
%X=ykb^-1F; %求△Xi
a=a+1
disp(‘迭代次数a=’);disp(a);
%X1(:,a+1)
X2(:,a)=[X(1:8,:);X(9:14).*Upq]; % 存△θi Vi

for i=1:n
    if(i<=pq)
    jd(i).zkj=jd(i).zkj-X2(i,a);
    jd(i).U=jd(i).U-X2(i+pq+pv,a);
    else if( (pq<i) && (i<pq+pv)) % PV  7 8 节点
            disp('jd(i).zkj');disp (jd(i).zkj);disp(X2(i,a));
            jd(i).zkj=jd(i).zkj-X2(i,a);
            disp('输出jd(i).zkj');
            disp(jd(i).zkj);
        end
    end
end  % 节点的角度、电压更新了 Xi+1 
%形成电压矩阵
for i= 1:pq  %pq 节点的U
Upq(i)=jd(i).U;
end

% Upq=Upq’

 %求新的ykb
              n=9;
            Qi=0;
            Pi=0;
            for i=1:n
                for j=1:n
                   %求H
                  if((i<=8)&&(j<=8)) 
                     if(i~=j)
                         H(i,j) = -jd(i).U*jd(j).U*(real(Y(i,j))*sin(jd(i).zkj-jd(j).zkj)...
                         -imag(Y(i,j))*cos(jd(i).zkj-jd(j).zkj));
                     else if(i==j)
                             for k=1:n
                                 Qi=Qi+jd(k).U*(real(Y(i,k))*sin(jd(i).zkj-jd(k).zkj)-imag(Y(i,k))*cos(jd(i).zkj-jd(k).zkj));
                             end
                             Qi=jd(i).U*Qi;
                             H(i,i)=Qi+imag(Y(i,i))*jd(i).U^2;
                         end
                     end
                  end   
                  %求N 8*6
                  if((i<=8)&&(j<=6))
                      if(i~=j)
                         N(i,j) = -jd(i).U*jd(j).U*(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)...
                         +imag(Y(i,j))*sin(jd(i).zkj-jd(j).zkj));
                      else if(i==j)
                             for k=1:n
                                 Pi=Pi+jd(k).U*(real(Y(i,k))*cos(jd(i).zkj-jd(k).zkj)+imag(Y(i,k))*sin(jd(i).zkj-jd(k).zkj));
                             end
                             Pi=jd(i).U*Pi;
                             N(i,i)=-Pi-real(Y(i,i))*jd(i).U^2;
                          end
                      end
                  end


                  %求K 6*8
                  if((i<=6)&&(j<=8))
                      if(i~=j)
                            K(i,j) = jd(i).U*jd(j).U*(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)...
                         +imag(Y(i,j))*sin(jd(i).zkj-jd(j).zkj));
                      else if(i==j)
                             for k=1:n
                                 Pi=Pi+jd(k).U*(real(Y(i,k))*cos(jd(i).zkj-jd(k).zkj)+imag(Y(i,k))*sin(jd(i).zkj-jd(k).zkj));
                             end
                             Pi=jd(i).U*Pi;
                             K(i,i)=-Pi+real(Y(i,i))*jd(i).U^2;
                          end
                      end
                  end

                  %L 6*6
                  if((i<=6)&&(j<=6)) 
                     if(i~=j)
                         L(i,j) = -jd(i).U*jd(j).U*(real(Y(i,j))*sin(jd(i).zkj-jd(j).zkj)...
                         -imag(Y(i,j))*cos(jd(i).zkj-jd(j).zkj));
                     else if(i==j)
                             for k=1:n
                                 Qi=Qi+jd(k).U*(real(Y(i,k))*sin(jd(i).zkj-jd(k).zkj)-imag(Y(i,k))*cos(jd(i).zkj-jd(k).zkj));
                             end
                             Qi=jd(i).U*Qi;
                             L(i,i)=-Qi+imag(Y(i,i))*jd(i).U^2;
                         end
                     end
                  end 


                end
            end           
disp('新的雅可比矩阵： ');disp(a);
disp(ykb);

%求F
                    F=zeros(2*pq+pv,1); %14行，1列
                for i=1:n
                    if(i<=pq) %求pq节点的△P  △Q
                           %求pq节点的6个 △P F(x)的前8行
                           Pi=0; %对每一个节点开始Pi赋值为0
                           for j=1:n  %循环所有节点求和，得到  Pi/Ui
                               Pi=Pi+jd(j).U*(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)+imag(Y(i,j))*sin(jd(i).zkj-jd(j).zkj));
                           end
                           Pi=Pi*jd(i).U;
                           F(i)=jd(i).P/100-Pi;  % P/100是标幺值

                           %求pq节点的6个 △Q F(x)的第 9-14行
                           Qi=0;
                           for j=1:n
                               Qi=Qi+jd(j).U*(real(Y(i,j))*sin(jd(i).zkj-jd(j).zkj)-imag(Y(i,j))*cos(jd(i).zkj-jd(j).zkj));
                           end
                           Qi=Qi*jd(i).U;
                           F(i+pq+pv)=jd(i).Q/100-Qi;   %第9-14行 Q/100 是进行标幺
                    end
                    if ( (pq<i)&&(i<=pq+pv) ) %PV节点
                           Pi=0; %对每一个节点开始Pi赋值为0
                           for j=1:n  %循环所有节点求和，得到  Pi/Ui
                               Pi=Pi+jd(j).U*(real(Y(i,j))*cos(jd(i).zkj-jd(j).zkj)+imag(Y(i,j))*sin(jd(i).zkj-jd(j).zkj));
                           end

% disp(‘Pi Fi=:’);
Pi=Pi*jd(i).U;
F(i)=jd(i).P/100-Pi %第7-8行
end
end

    disp('新的F(x)');disp(F);
    X=(ykb^-1)*F; %求△θi+1 △ Ui+1/Ui+1
    X1=[X(1:8,:);X(9:14).*Upq];  %△θi+1  △Ui+1
    disp('新的△θi+1  △Ui+1');
    disp(a);
    disp(X1)

end

欢迎使用Markdown编辑器

你好！这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Markdown编辑器, 可以仔细阅读这篇文章，了解一下Markdown的基本语法知识。

新的改变

我们对Markdown编辑器进行了一些功能拓展与语法支持，除了标准的Markdown编辑器功能，我们增加了如下几点新功能，帮助你用它写博客：

全新的界面设计 ，将会带来全新的写作体验；
在创作中心设置你喜爱的代码高亮样式，Markdown 将代码片显示选择的高亮样式 进行展示；
增加了 图片拖拽 功能，你可以将本地的图片直接拖拽到编辑区域直接展示；
全新的 KaTeX数学公式 语法；
增加了支持甘特图的mermaid语法¹ 功能；
增加了 多屏幕编辑 Markdown文章功能；
增加了 焦点写作模式、预览模式、简洁写作模式、左右区域同步滚轮设置 等功能，功能按钮位于编辑区域与预览区域中间；
增加了 检查列表 功能。

功能快捷键

撤销：Ctrl/Command + Z
重做：Ctrl/Command + Y
加粗：Ctrl/Command + B
斜体：Ctrl/Command + I
标题：Ctrl/Command + Shift + H
无序列表：Ctrl/Command + Shift + U
有序列表：Ctrl/Command + Shift + O
检查列表：Ctrl/Command + Shift + C
插入代码：Ctrl/Command + Shift + K
插入链接：Ctrl/Command + Shift + L
插入图片：Ctrl/Command + Shift + G
查找：Ctrl/Command + F
替换：Ctrl/Command + G