qlimitpflow.m

电力系统潮流计算的程序源代码

function [baseMVA,bus0,gen0,Ybus,bus,gen,J,tag]= Qlimitpflow(casename);
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                1 开始潮流程序的一些转换和运算                           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------—1.1 
%% 对本函数输入参数的说明（即[MVAbase, bus, gen, branch, success, et] = mypf(casename, mpopt, fname)的说明）
if nargin>2
error('输入变量太多');%casename只能输入一个文件名
end
%% 
%[baseMVA, bus, gen, branch] = feval(casename);  %%feval函数可以参考help feval
%%
lmda=0.7;
%[baseMVA, bus, gen, branch] = feval('IEEE118');
% [baseMVA, bus, gen, branch] = feval('IEEE57');
[baseMVA, bus, gen, branch] = feval('IEEE30');%feval里调用的文件名应与casename一致
%%---------------------------------------------------
bus0=bus;
gen0=gen;
bustype0=bus0(:,2);
pv0=find(bustype0==2);
pq0=find(bustype0==1);
%--------—1.1 

%-----1.2对节点类型处理的补充，发电机不在运行状态的节点作为PQ节点
%% 检查发电机状态
bus_gen_status = zeros(size(bus(:,1)));
bus_gen_status(gen(:, 1)) = gen(:,8);
%% 区分节点型式
ref	= find(bus(:,2) == 3 & bus_gen_status);	
pv	= find(bus(:,2) == 2  & bus_gen_status);	
pq	= find(bus(:, 2) == 1 | ~bus_gen_status);	
%% 选择一条新的参考节点 ( slack意外停机状态下 )
if isempty(ref)
	ref = pv(1);				%% 选择第一条PV节点
	pv = pv(2:length(pv));		%% 从节点列表中除去这个节点
end
%-----1.2对节点类型处理的补充，发电机不在运行状态的节点作为PQ节点


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%        重要部分，建立节点导纳阵          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%---------makeYbus,建立节点导纳阵
j = sqrt(-1);
nb = size(bus, 1);			%% 统计节点数目
nl = size(branch, 1);		%% 统计支路数目
%% 检查节点标号是否按照一定顺序
if any(bus(:,1) ~= [1:nb]')
	error('节点必须按照节点规定的顺序出现')
end

stat = branch(:,11);			    		%% 运行中的支路
Ys = stat ./ (branch(:, 3) + j * branch(:, 4)); %%支路的导纳	
Bc = stat .* branch(:, 5);					%% 支路的充电电纳
tap = ones(nl, 1);								%% 缺省tap置1,nl行1列
i = find(branch(:, 9));						%% 寻找非零TAP(ratio)
tap(i) = branch(i, 9);						%% 用非零TAP(ratio)替换
tap = tap .* exp(j*pi/180 * branch(:,10));	%% 考虑相位调整(若有移相器）
tap2= tap.*conj(tap);

Ybus = sparse(branch(:,1), branch(:,2), -Ys./conj(tap), nb, nb) + ...
         sparse(branch(:,2), branch(:,1), -Ys./tap, nb, nb) + ...
         sparse(branch(:,1), branch(:,1), (Ys+j*Bc/2)./tap2,nb, nb) + ...
         sparse(branch(:,2), branch(:,2), Ys+j*Bc/2,nb, nb);
Gs=bus(:,5)/ baseMVA;
Bs=bus(:,6)/ baseMVA;
Ybus = Ybus + sparse(1:nb,1:nb,Gs+j*Bs, ...
                         nb, nb);
%%————————————导纳阵完成--------------


%% 发电机信息
on = find(gen(:, 8));				%% 判断哪些机组在运行
%% 对 V 和 Pg 根据case文件的数据进行初始化
V0=ones(nb,1);
V0(gen(on, 1)) = gen(on, 6) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%              开始运行潮流               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t0 = clock;     %% 计算迭代时间（开始计时）
    %% 参数
    tol		= 1e-6;
    max_it	= 100;
    %% 初始化
    j = sqrt(-1);
    converged = 0;
    conter = 0;
    pv_change=0;
    V = V0;
    Va = angle(V);
    Vm = abs(V);
    pv_pre=pv;
    pv_post=pv;
    %--------------makeSbus,计算节点注入和输出功率
%% 发电机信息
on = find(gen(:, 8));				%% 哪个发电机在运行？
gbus = gen(on, 1);					%% 在哪个母线？
ngbus=size(gbus,1);
tag=zeros(2,ngbus);
note=0;

    %% 进行Newton迭代
   while (~converged & conter < 100)
	      %% 更新迭代指针
	   

%% --------------------形成节点注入功率矢量-------------------------


%Sbus = -(bus(:,3) + j * bus(:, 4));						%% 负荷注入的功率
%Sbus(gbus) = Sbus(gbus) + gen(on, 2) + j * gen(on, 3);	%% 备用发电机
%Sbus = Sbus / baseMVA;										%% 转换成 p.u.

 %%-------------注入功率的形成-------------------------------------------
Sbus = -(bus(:,3) + j * bus(:, 4))*(1+lmda);%% 负荷注入的功率
Pg_pump=zeros(size(gbus,1),1);
Pg=gen(:,2);
pump_location=find(gen(:,2)<0);
Pg_pump(pump_location)=Pg(pump_location);
Pg(pump_location)=0;
Sbus(gbus) = Sbus(gbus) +(Pg)*(1+lmda)+Pg_pump +j*gen(:,3);	     %% 备用发电机
Sbus = Sbus / baseMVA;										%% 转换成 p.u.   

        %% 计算 F(x)
	      mis = V .* conj(Ybus * V) - Sbus;
	      F = [	real(mis(pv));
			    real(mis(pq));
			    imag(mis(pq))];
         %检查是否收敛
         
	normF = norm(F, inf);
	if note==0&normF<tol       %note＝1为pv节点有转换，若为0则无转换
		converged = 1;
		fprintf('\nNewton迭代法经过 %d 次迭代收敛.\n', conter); break;
    end    
        
%--------------------------	Jacobian--------------
	      %% 计算 dSbus-dV,节点注入功率对对电压的变化率
	      j = sqrt(-1);
          n = length(V);
          Ibus = Ybus * V;
         conter = conter + 1;

          if issparse(Ybus)			%% 稀疏矩阵处理 (如果 Ybus 稀疏矩阵)
	         diagV		= spdiags(V, 0, n, n);
	         diagIbus	= spdiags(Ibus, 0, n, n);
	         diagVnorm	= spdiags(V./abs(V), 0, n, n);
          else						
	         diagV		= diag(V);
	         diagIbus	= diag(Ibus);
	         diagVnorm	= diag(V./abs(V));
          end

          dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm;
          dSbus_dVa = j * diagV * conj(diagIbus - Ybus * diagV);

          %%为了节省计算时间，利用临时的置换矩阵	
          % 	j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
          % 	j12 = real(dSbus_dVm([pv; pq], pq));
          % 	j21 = imag(dSbus_dVa(pq, [pv; pq]));
          % 	j22 = imag(dSbus_dVm(pq, pq));

	      temp = real(dSbus_dVa(:, [pv; pq]))';
	      j11 = temp(:, [pv; pq])';
	      temp = real(dSbus_dVm(:, pq))';
	      j12 = temp(:, [pv; pq])';
	      temp = imag(dSbus_dVa(:, [pv; pq]))';
	      j21 = temp(:, pq)';
	      temp = imag(dSbus_dVm(:, pq))';
	      j22 = temp(:, pq)';
   
	      J = [	j11 j12;
			    j21 j22;	];   %%   Jacobian 阵 (2n-r)*(2n*r)

	      %% 计算迭代步长
	      dx = -(J \ F);
          
        %% 建立修正V的索引      
    npv	= length(pv);
    npq	= length(pq);
    n1 = 1;			n2 = npv;			%% n1:n2 - PV节点电压相角n1=1;n2=npv;
    n3 = n2 + 1;	n4 = n2 + npq;		%% n3:n4 - PQ节点电压相角n3=npv+1;n4=npv+npq;
    n5 = n4 + 1;	n6 = n4 + npq;		%% n5:n6 - PQ节点电压值n5=npv+npq+1;j6=npv+npq+npq
          

	      %% 修正电压
   	      Va(pv) = Va(pv) + dx(n1:n2);
	      Va(pq) = Va(pq) + dx(n3:n4);
	      Vm(pq) = Vm(pq) + dx(n5:n6);
	      V = Vm .* exp(j * Va);
         
          %--------------------------	Jacobian----------------------------
         


 pv_pre=pv;
%--------------------------------------------------------

%%进行无功越限处理
V0m=abs(V0);
temp = Ybus.';
Sg = V(gbus) .* conj(temp(:, gbus).' * V);
Qg = zeros(size(gbus,1), 1);			
Qg_max=zeros(size(gbus,1),1);
Qg_min=zeros(size(gbus,1),1);
bustype=zeros(size(bus,1),1);
Qg=imag(Sg) * baseMVA + bus(gbus, 4)*(1+lmda);	
Qg_max=gen(:,4);
Qg_min=gen(:,5);
bustype=bus(:,2);
ngbus=size(gbus,1);
note=0;
% if abs(Qg-Qg_max)<tol
 %    Qg=Qg_max;
 %end
% if abs(Qg-Qg_min)<tol
 %    Qg=Qg_min;
 %end
for i=1:ngbus
    if bustype(gbus(i))==2
        if Qg(i)>Qg_max(i)
            gen(i,3)=Qg_max(i);
            bus(gbus(i),2)=1;
            pv	= find(bus(:,2) == 2  & bus_gen_status);
            pq	= find(bus(:, 2) == 1 | ~bus_gen_status);
            tag(1,i)=gbus(i);
            tag(2,i)=1;
            note=1;
        end
        if  Qg(i)<Qg_min(i)
            gen(i,3)=Qg_min(i);
            bus(gbus(i),2)=1;
            pv	= find(bus(:,2) == 2  & bus_gen_status);	
            pq	= find(bus(:, 2) == 1 | ~bus_gen_status);
            tag(1,i)=gbus(i);
            tag(2,i)=-1;
            note=1;
        end
    end
end

for i=1:ngbus
    if bustype(gbus(i))==1&tag(2,i)==-1
        if Vm(gbus(i))<V0m(gbus(i))
            if Qg(i)>Qg_max(i)
                gen(i,3)=Qg_max(i);
                tag(2,i)=1;
                note=1
            end
            if Qg(i)<Qg_min(i)
                note=1;
            end
            if Qg(i)<=Qg_max(i)&Qg(i)>=Qg_min(i)
                bus(gen(i,1),2)=2;
                pv	= find(bus(:,2) == 2  & bus_gen_status);	
                pq	= find(bus(:, 2) == 1 | ~bus_gen_status);
                Vm(gbus(i))=V0m(gbus(i));
                tag(2,i)=0;
                note=1;
            end
        end
    end
end

for i=1:ngbus
    if bustype(gbus(i))==1&tag(2,i)==1
        if Vm(gbus(i))>V0m(gbus(i))
            if Qg(i)<Qg_min(i)
                gen(i,3)=Qg_min(i);
                tag(2,i)=-1;
                note=1
            end
            if Qg(i)>Qg_max(i)
                note=1;
            end
            if Qg(i)<=Qg_max(i)&Qg(i)>=Qg_min(i)
                bus(gen(i,1),2)=2;
                pv	= find(bus(:,2) == 2  & bus_gen_status);	
                pq	= find(bus(:, 2) == 1 | ~bus_gen_status);
                Vm(gbus(i))=V0m(gbus(i));
                tag(2,i)=0;
                note=1;
            end
        end
    end
end

pv_post=pv;
end


if ~converged
		fprintf('\nNewton迭代法经过 %d 次迭代不收敛.\n', conter);
end
bustype=bus(:,2);
%-----------Jacobian完成--------------------

%--------------------------越限的节点--------------------------------
tag;
reachlimit=pv0(find(ismember(pv0,pv)==0));
pv_uplimit=reachlimit(find((Vm(reachlimit)-V0m(reachlimit))<0));
pv_lowlimit=reachlimit(find((Vm(reachlimit)-V0m(reachlimit))>0));


%%---------------修正节点、发电机------------------------
j = sqrt(-1);
nl = size(branch, 1);		
%% 初始化返回值
bus		= bus;
gen		= gen;
branch	= branch;

%%----- 修正节点电压 -----
bus(:, 9) = abs(V);
bus(:, 10) = angle(V) * 180 / pi;

%%----- 修正所有发电机Qg和浮动节点的Pg -----
%% 发电机信息
refgen = find(gen(on, 1) == ref);		%% 哪个是基准发电机？

%% 计算总的节点注入或输出功率
temp = Ybus.';
Sg = V(gbus) .* conj(temp(:, gbus).' * V);

%% 修正所有发电机的Qg
gen(:, 3) = zeros(size(gen, 1), 1);				%% Qg置零,行size(gen, 1),列1
gen(on, 3) = Qg;	%% 注入 Q + 基本 Qd

%% 修正所有浮动节点的Pg
gen(on(refgen), 2) = real(Sg(refgen)) * baseMVA + bus(ref, 3);	%% 注入 P + 基本 Pd

%%%%----------------------------------------------------------
et = etime(clock, t0);
%--------------------------------------


return;