eign6161.m

MATLAB潮流计算

a(6*gn+1:7*gn,1:gn)=-inv(te)*ke*k5;
a(6*gn+1:7*gn,3*gn+1:4*gn)=-inv(te)*ke*k6;
a(6*gn+1:7*gn,5*gn+1:6*gn)=-inv(te)*ke*k7;
a(6*gn+1:7*gn,6*gn+1:7*gn)=-inv(te);
%--------------------------------
%根据发电机模型类型修正A 阵
for n=1:gn,
   if mt(1,n)==1,                         %无g绕组,
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0
   elseif mt(1,n)==2,                     %无D+g绕组，
      a(3*gn+n,:)=0;                      %eqpp 行置0，
      a(:,2*gn+n)=a(:,2*gn+n)+a(:,3*gn+n);% eqpp 列+ eqp列
      a(:,3*gn+n)=0;                      %eqpp列置0
      
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0

   elseif mt(1,n)==3                      %无g+Q绕组
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0

      a(5*gn+n,:)=0;                      %edpp 行置0，
      a(:,5*gn+n)=0;                      %edpp列置0
   elseif mt(1,n)==4,                     %无D+g+Q绕组，
      a(3*gn+n,:)=0;                      %eqpp 行置0，
      a(:,2*gn+n)=a(:,2*gn+n)+a(:,3*gn+n);% eqpp 列+ eqp列
      a(:,3*gn+n)=0;                      %eqpp列置0
      
      a(5*gn+n,:)=a(5*gn+n,:)-a(4*gn+n,:);%edpp行-edp行
      a(4*gn+n,:)=0;                      %edp行置0
      a(:,4*gn+n)=0;                      %edp列置0

      a(5*gn+n,:)=0;                      %edpp 行置0，
      a(:,5*gn+n)=0;                      %edpp列置0
   end
end
%----------------------------------
%根据发电机模型类型压缩A 阵
dp=0;qp=0;gp=0;
for n=1:gn,
   k=mt(1,n);
   if k==1,                                      % 无g绕组
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 删edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1
   end   
   if k==2,                                      % 无D+g绕组
      a(3*gn+n-dp,:)=[];    a(:,3*gn+n-dp)=[];   % 删eqpp 行和列
      a(7*gn,7*gn)=1;       dp=dp+1;             %末元素置1
      
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 删edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1
   end
   if k==3,                                      % 无g+Q绕组
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 删edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1

      a(5*gn+n-qp-dp-gp,:)=[]; a(:,5*gn+n-qp-dp-gp)=[];%删edpp 行和列
      a(7*gn,7*gn)=1;       qp=qp+1;             %末元素置1
   end
   if k==4,      % 无D+g+Q绕组,
      tqp=0;tqpp=0;tdpp=0;
      a(3*gn+n-dp,:)=[];    a(:,3*gn+n-dp)=[];   % 删eqpp 行和列
      a(7*gn,7*gn)=1;       dp=dp+1;             %末元素置1
      
      a(4*gn+n-gp-dp,:)=[];a(:,4*gn+n-gp-dp)=[]; % 删edp 行和列
      a(7*gn,7*gn)=1;       gp=gp+1;             %末元素置1

      a(5*gn+n-qp-dp-gp,:)=[]; a(:,5*gn+n-qp-dp-gp)=[];%删edpp 行和列
      a(7*gn,7*gn)=1;       qp=qp+1;             %末元素置1
   end
end
dq=dp+qp+gp;                                     %删去的行列总数
a(7*gn-dq+1:7*gn,:)=[];                          %删去最后dq行和列
a(:,7*gn-dq+1:7*gn)=[];
%====================================
%计算特征根,自然振荡频率,阻尼比eig0,wn,k
[v1,d]=eig(a);                  %左特征向量，特征根对角阵
eig0=diag(d);                   %特征根
%====================================
%计算右特征向量w和参与矩阵pwv
w=inv(v1)';                     %右特征向量
pwv=abs(w.*v1);                 %参与矩阵
%====================================
n=7*gn-dq;                     %以七阶模型为基础
pwvm=0.;
for k=1:n
   pwvm(1,k)=max(pwv(:,k));
   for kk=1:n
      if pwvm(1,k)==pwv(kk,k);
         pwvm(2,k)=kk;
      end
   end
   %pwv(:,k)=pwv((:,k)/pwvm(1,k);
end
pwvm=pwvm';
for k=1:n
   pwvm(k,3)=eig0(k);
end
%找出和特征值次相关的变量位置及其相关度
%注意此时变量是以变量为一簇的
pwvm(:,5)=(min(pwv))';
for i=1:n
   for l=1:n
      wox=pwv(l,i);
      if wox<pwvm(i,1)&wox>pwvm(i,5);
         pwvm(i,5)=wox;
         pwvm(i,4)=l;
      end
   end
end
%=====================================
%机电模式相关比pmm和电气模式相关比pee
n=7*gn-dq;                     %以七阶模型为基础
pmm=(sum(pwv(1:2*gn,:))./sum(pwv(2*gn+1:n,:)))';
format long,
pee=1./pmm;
%====================================
%pmm(n-2:n,1)=pmm(n-2:n,1)*0;
%=====================================
ewk=0;%（1特征值编号、2特征值、3阻尼比、4频率、5强相关度、
          %6强相关变量位置、7次相关度、8次相关变量位置）
for n=1:7*gn-dq                     %以七阶模型为基础
   ewk(n,1)=n;                      %特征根序号
end
ewk(:,2)=eig0;                      %特征根
ewk(:,3)=-real(eig0)./abs(eig0);    %--特征根实部/特征根的模值
ewk(:,4)=abs(imag(eig0)/(2*3.1416));%频率
ewk(:,5:9)=pwvm;
%====================-------------------------------======================
for k=1:7*gn-dq% 参与矩阵规格化
   pwv(:,k)=pwv(:,k)./pwvm(k,1);
end
%==========================================
format short
%分离发电机的参数
%-----------------------------------------------------
ef=[1,1,1,1,1,1,1,1,1,1];   %各发电机励磁系统的阶数
%统计发电机的各个参数的个数及各发电机的阶数fact
eqpp=zeros(1,2);
ed=zeros(1,2);%edp
edpp=zeros(1,2);
fact=zeros(1,gn);
for i=1:gn           
   if mt(i)==1       %D+Q
      eqpp(1,1)=eqpp(1,1)+1;
      edpp(1,1)=edpp(1,1)+1;
      fact(i)=5+ef(i);
   elseif mt(i)==2   %Q
      edpp(1,1)=edpp(1,1)+1;
      fact(i)=4+ef(i);
   elseif mt(i)==3   %D
      eqpp(1,1)=eqpp(1,1)+1;
      fact(i)=4+ef(i);
   elseif mt(i)==4
      fact(i)=3+ef(i);
   elseif mt(i)==5 %D+Q+g
      eqpp(1,1)=eqpp(1,1)+1;
      edpp(1,1)=edpp(1,1)+1;
      ed(1,1)=ed(1,1)+1;
      fact(i)=6+ef(i);
   end
end
%---------------------------------------
%----------------------------
%确定发电机的各个变量的起始位置
eqpp(1,2)=3*gn+1;
ed(1,2)=eqpp(1,1)+eqpp(1,2);
edpp(1,2)=ed(1,1)+ed(1,2);
%励磁系统的阶数及起始位置
begin=edpp(1,1)+edpp(1,2);
%-----------------------------   
%行置换
num=0;
for i=1:gn
   b(num+1,:)=a(i,:);
   b(num+2,:)=a(gn+i,:);
   b(num+3,:)=a(2*gn+i,:);
   if mt(i)==1
      b(num+4,:)=a(eqpp(1,2),:);
      eqpp(1,2)=eqpp(1,2)+1;
      b(num+5,:)=a(edpp(1,2),:);
      edpp(1,2)=edpp(1,2)+1;
      %加入励磁系统变量,多阶的话就在此添加循环.就是num+5+l(励磁的循环次数)
      b(num+6,:)=a(begin,:);
      begin=begin+ef(i);
      num=num+fact(i);
   elseif mt(i)==2
      b(num+4,:)=a(edpp(1,2),:);
      edpp(1,2)=edpp(1,2)+1;
       %加入励磁系统变量,多阶的话就在此添加循环
      b(num+5,:)=a(begin,:);
      begin=begin+ef(i);
      num=num+fact(i);
   elseif mt(i)==3
      b(num+4,:)=a(eqpp(1,2),:);
      eqpp(1,2)=eqpp(1,2)+1;
       %加入励磁系统变量,多阶的话就在此添加循环
      b(num+5,:)=a(begin,:);
      begin=begin+ef(i);
      num=num+fact(i);
   elseif mt(i)==4
      %加入励磁系统变量,多阶的话就在此添加循环
      b(num+4,:)=a(begin,:);
      begin=begin+ef(i);
      num=num+fact(i);
   elseif mt(i)==5
      b(num+4,:)=a(eqpp(1,2),:);
      eqpp(1,2)=eqpp(1,2)+1;
      b(num+5,:)=a(ed(1,2),:);
      ed(1,2)=ed(1,2)+1;
      b(num+6,:)=a(edpp(1,2),:);
      edpp(1,2)=edpp(1,2)+1;
       %加入励磁系统变量,多阶的话就在此添加循环
      b(num+7,:)=a(begin,:);
      begin=begin+ef(i);
      num=num+fact(i);
   end
end
%=====----------------
%列置换
global af;
af=zeros(size(a));%分离后的系统矩阵

%确定发电机的各个变量的起始位置
eqpp(1,2)=3*gn+1;
ed(1,2)=eqpp(1,1)+eqpp(1,2);
edpp(1,2)=ed(1,1)+ed(1,2);
%励磁系统的阶数及起始位置
%ef=[1,1,1];
begin=edpp(1,1)+edpp(1,2);
num=0;
for i=1:gn
   af(:,num+1)=b(:,i);
   af(:,num+2)=b(:,gn+i);
   af(:,num+3)=b(:,2*gn+i);
   if mt(i)==1
      af(:,num+4)=b(:,eqpp(1,2));
      eqpp(1,2)=eqpp(1,2)+1;
      af(:,num+5)=b(:,edpp(1,2));
      edpp(1,2)=edpp(1,2)+1;
      %加入励磁系统变量,多阶的话就在此添加循环
      af(:,num+6)=b(:,begin);
      begin=begin+ef(i);
      num=num+5+ef(i);
   elseif mt(i)==2
      af(:,num+4)=b(:,edpp(1,2));
      edpp(1,2)=edpp(1,2)+1;
       %加入励磁系统变量,多阶的话就在此添加循环
      af(:,num+5)=b(:,begin);
      begin=begin+ef(i);
      num=num+4+ef(i);
   elseif mt(i)==3
      af(:,num+4)=b(:,eqpp(1,2));
      eqpp(1,2)=eqpp(1,2)+1;
       %加入励磁系统变量,多阶的话就在此添加循环
      af(:,num+5)=b(:,begin);
      begin=begin+ef(i);
      num=num+4+ef(i);
   elseif mt(i)==4
      %加入励磁系统变量,多阶的话就在此添加循环
      af(:,num+4)=b(:,begin);
      begin=begin+ef(i);
      num=num+3+ef(i);
   elseif mt(i)==5
      af(:,num+4)=b(:,eqpp(1,2));
      eqpp(1,2)=eqpp(1,2)+1;
      af(:,num+5)=b(:,ed(1,2));
      ed(1,2)=ed(1,2)+1;
      af(:,num+6)=b(:,edpp(1,2));
      edpp(1,2)=edpp(1,2)+1;
       %加入励磁系统变量,多阶的话就在此添加循环
      af(:,num+7)=b(:,begin);
      begin=begin+ef(i);
      num=num+6+ef(i);
   end
end
%===========
%--------------------------利用模糊的方法对系数矩阵Ki进行聚类分析------------------
k100=abs(k1);
s=size(k1);
a1=s(1);
a2=s(2);
%模糊数学最大最小法(a1行数,a2列数)
%几何平均最小法
for i=1:a1
   for c=1:a1
      s1=0;
      s2=0;
      for k=1:a2
         if k100(i,k)<k100(c,k)
            s1=s1+k100(i,k);
         else
            s1=s1+k100(c,k);
         end
         s2=s2+(k100(c,k)*k100(i,k))^0.5;
      end
   wwww(i,c)=s1/s2;   
   end
end

%传递闭包
zs=ceil(log2(-1));
for k=1:zs
   for i=1:a1
      for c=1:a1
         s2=0;
         for h=1:a1
            if wwww(i,h)<wwww(h,c)
               s1=wwww(i,h);
            else
               s1=wwww(h,c);
            end
            if s1>s2
               s2=s1;
            end
         end
         wwww(i,c)=s2;
      end
   end
end
wwww1=wwww;
for i=1:a1
   for c=1:a1
      if wwww1(i,c)<0.6
         wwww1(i,c)=0;
      else
         wwww1(i,c)=1;
      end
   end
end
%==============================================================
global order gnum;
gnum=[5,1];
order=4;
bn=0;
[eig1,pwvm1,pmm1,n1]=jie_ou(1,bn);
bn=bn+n1;
[eig2,pwvm2,pmm2,n2]=jie_ou(2,bn);

%#############SMA法###############
%从计算结果读入初值(初值应考虑到是不是机电模式--虚部是不是在1.26--15.7之间）
value_first1=[eig1;eig2];
[num_all,z]=size(value_first1);
num_yes=0;
for k=1:num_all
   if real(value_first1(k))>0|...%是不是有不稳定的值
         (abs(imag(value_first1(k)))>1.25&abs(imag(value_first1(k)))<15.8)
      %初值考虑到是不是机电模式--虚部是不是在1.26--15.7之间
      num_yes=num_yes+1;
      value_first(num_yes,1)=value_first1(k);
   end
end
%是机电模式的特征值初值
%倪以信SMA法
no=zeros(gn,gn);
one=eye(gn,gn);
%收敛限制
ersp=1e-6;
%机电模式的个数
[n,n1]=size(value_first);
mm=zeros(n,1);%计算收敛前计算的次数
for k=1:n
value=value_first(k);
l=1;
while l==1
%各控制器的传递函数
Ge=ke/(one+te*value);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%计算简化二阶的系统阵
L1s=(one+value*tdp)+k3+Ge*(k6+Gpss*k2);
L2s=-k4-Ge*(k5-Gpss*k1);
Ks=k1+k2*inv(L1s)*L2s;
Ds=D+Ggov+k2*inv(L1s)*Ge*Gpss;
Ar=[no,one*100*pi;-inv(tj)*Ks,-inv(tj)*Ds];
[vot,val]=eig(Ar);
value0=diag(val);
[min_erro,h_s]=min(abs(value-value0));
if min_erro<ersp
      l=0;
end
votor=vot(:,h_s);
value=value0(h_s);
mm(k,1)=mm(k,1)+1;
end %while
%取出真实的计算结果
value_last(k,1)=value;
votor_last(:,k)=votor;
end %for