lhpfcomp.m

电力系统电压稳定研究的图形化软件

        % VST_LFSetParam;

			% set LF control control errors
         set(AbsErrorDisp,'String',num2str(AbsError));
         %==============================================
        	if isstr(RelError)
            set(RelErrorDisp,'String',RelError);
        	else
            set(RelErrorDisp,'String',num2str(RelError));
         end
         %==============================================
         

        set(NumIterations,'String',num2str(j*ReportCycle));
        set(IterationTime,'String',num2str(etime(clock,t0)/ReportCycle));
        %=====================================
        if (AbsError<=LFAbsTol) ...
            & ((~isstr(RelError)) ...
            & (RelError<=LFRelTol) ...
            | isstr(RelError))
         	ConvergenceFlag=1;
         
            % check the zero eigenvalue
            %=========================
            if k==NRS_Steps+1
               LF_jacob=J(2:n+1,1:n);
               eig_LF_Jacob=eig(LF_jacob);
               eig_Dyg_tip=eig(J(sub_strt+1:fn+1,sub_strt:fn));
            end
            %======================
            % Evaluate right and left eigenvector 
            % at the point of collpase
            %=====================================
            if k==NRS_Steps+1
               vpoc=v;
               wpoc=-null(J(2:n+1,1:n)');
               [mp,np]=size(XX);
            end
            %=======================================  
            
            break;
         end
         %===============================================================================
      end
      %====================================================================================
	%===============================
   if ConvergenceFlag==0
      'NRS Failed to Converge'
       
        break;
     end
     %=============================
	%=======================
    if alpha>=alphamax
        return;
     end
     %=======================

XX=[XX x];
AA=[AA alpha];
PP=[PP param];
lambda=lambda+deltalambda;
Ksys=J(2:no_gen,1:no_gen-1)-J(2:no_gen,no_gen:n)...
   *(J(no_gen+1:n+1,no_gen:n)\J(no_gen+1:n+1,1:no_gen-1));
Asys=[zeros(size(diag(gen_inertia(2:no_gen)))) ...
        diag(gen_inertia(2:no_gen))
     -Ksys, -diag(gen_damp(2:no_gen))/diag(gen_inertia(2:no_gen))];
%Asys_1=[zeros(size(diag(gen_inertia(2:no_gen)))) eye(no_gen-1)
 %    -diag(gen_inertia(2:no_gen))\Ksys -diag(gen_inertia(2:no_gen))\diag(gen_damp(2:no_gen))];


dd=eig(Asys);
[r_maxdd,s]=max(real(dd));
sys_eig=[sys_eig dd];

%if k<=NRS_Steps
 %  OPTIONS.disp=0;
  % [v_small dd_small]=eigs(Asys_1,1,'LR',OPTIONS);
   %sys_small=[sys_small dd_small];
   %sys_eig_small1=[sys_eig_small1 (v_small)];
   %v_small_gen=abs(v_small(1:no_gen-1));
	%v_small_gen=real(v_small(1:no_gen-1));
   %v_small_gen=v_small_gen/norm(v_small_gen);
   %v_small_gen=abs(v_small_gen);
	%sys_eig_small=[sys_eig_small v_small_gen];


%end



%ee=eigs(J(sub_strt+1:fn+1,sub_strt:fn),OPTIONS);
%[ind1]=find((imag(ee)==0)&(real(ee)<0));
%reg_index1=[reg_index1 length(ind1)];
%[ind2]=find((imag(ee)==0)&(real(ee)>0));
%reg_index2=[reg_index2 length(ind2)];

%Dyg_eig=[Dyg_eig ee];
%ff=eigs(Ksys,OPTIONS);											
%Ksys_eig=[Ksys_eig ff];	
							
	%=================================
    if r_maxdd<=100*eps  
    	if sign(imag(dd(s)))~=0       
        Stab=[Stab 1];
        else
        Stab=[Stab 2];
        end
     elseif r_maxdd>100*eps  
    	if sign(imag(dd(s)))~=0      
        Stab=[Stab 3];
        else
        Stab=[Stab 4];
        end   
    
    else
    	Stab=[Stab 5];	
    	
    end
    %===============================
    
 end
 %===================================================================

 'back to the NR algorithm'

% Back to NR algorithm

[n_rows,n_cols]=size(AA);
%***********************************************
%the following is used for IEEE14 bus system
del_alpha=(AA(n_cols)/(NR_steps));  
%***********************************************
%del_alpha=(alphamax/(NR_steps));

%alpha=alpha-del_alpha;

for k=1:NR_steps+(0.6)*NR_steps
	ConvergenceFlag=0;
	for j=1:round(MaxIterations/ReportCycle)
        t0=clock;
        for i=1:ReportCycle
            x0=x;
            [f,J]=eval([CurrentSystem,'(data,x,[0;param],v)']);
            J=full(J);
            delta=-sparse(J(2:n+1,1:n))\f(2:n+1);
            x=x0+delta;
         end
         
        AbsError=max(abs(x-x0));
        if x0==0
            RelError='NA';
        else
            RelError=AbsError/max(abs(x0));
        end
          
        %%set LF control control errors
        
       set(AbsErrorDisp,'String',num2str(AbsError));
        if isstr(RelError)
            set(RelErrorDisp,'String',RelError);
        else
           set(RelErrorDisp,'String',num2str(RelError));
        end

       set(NumIterations,'String',num2str(j*ReportCycle));
       set(IterationTime,'String',num2str(etime(clock,t0)/ReportCycle))

        if (AbsError<=LFAbsTol*0.001) ...
           & ((~isstr(RelError)) ...
            & (RelError<=LFRelTol*0.01) ...
             | isstr(RelError))
            ConvergenceFlag=1;
            break;
        end
    end

if ConvergenceFlag==0
   'NR fails to converge'
   
   break;
    end

    XX=[XX x];
    AA=[AA alpha];
    PP=[PP param];
   
    Ksys=J(2:no_gen,1:no_gen-1)-J(2:no_gen,no_gen:n)...
        *(J(no_gen+1:n+1,no_gen:n)\J(no_gen+1:n+1,1:no_gen-1));
    Asys=[zeros(size(diag(gen_inertia(2:no_gen))))...
        diag(gen_inertia(2:no_gen))
     -Ksys, -diag(gen_damp(2:no_gen))/diag(gen_inertia(2:no_gen))];
%  Asys_1=[zeros(size(diag(gen_inertia(2:no_gen)))) eye(no_gen-1)
 %    -diag(gen_inertia(2:no_gen))\Ksys -diag(gen_inertia(2:no_gen))\diag(gen_damp(2:no_gen))];


dd=eig(Asys);

    sys_eig=[sys_eig dd];
    %ee=eigs(J(sub_strt+1:fn+1,sub_strt:fn),OPTIONS);
    %[ind1]=find((imag(ee)==0)&(real(ee)<0));
    %reg_index1=[reg_index1 length(ind1)];
    %[ind2]=find((imag(ee)==0)&(real(ee)>0));
    %reg_index2=[reg_index2 length(ind2)];

	% Dyg_eig=[Dyg_eig ee];
    %ff=eigs(Ksys,OPTIONS);										
    %Ksys_eig=[Ksys_eig ff];
    [r_maxdd,s]=max(real(dd));

    if r_maxdd<=100*eps  
    	if sign(imag(dd(s)))~=0       
        Stab=[Stab 1];
        else
        Stab=[Stab 2];
        end
     elseif r_maxdd>100*eps  
    	if sign(imag(dd(s)))~=0       
        Stab=[Stab 3];
        else
        Stab=[Stab 4];
        end   
    
    else
    	Stab=[Stab 5];	
    	
    end
    alpha=alpha-del_alpha;
    if alpha<=0
        return;
     end
     %*****************************************************************************
     % The following investigates the sensitivity around SIB if it happens....
     % I want to identify the stability exchange from unstable to stable 
     sib=length(Stab);				% obtain length of the Stab
     sib_1=Stab(sib);				% Check stability feature
     sib_2=Stab(sib-1);				%Check stability feature of previous load value
     if sib_2~=sib_1					%check stability exchanges
        'hello:Singularity-Induced Bifurcation has occured';
        sib
        options.disp=0;
        Dyg=J(sub_strt+1:fn+1,sub_strt:fn);			%Dyg is the Jacobian of algebraic part
       [dyg_v1,dygeig1,flag1]=eigs(Dyg,1,'SM',options);		
       [sib_v1,sib_eig1,flag2]=eigs(Asys,1,'LM',options);
    end
    
  % ===================================================================  
  param=param-p*del_alpha;
  
end