jiance.m

粒子群优化算法(PSO)是一种进化计算技术(evolutionary computation)

function jiance(x)
%
%   testsopf(casename, which_flow)
%
%   Runs a full Newton's method power flow where casename is the name of
%   the m-file (without the .m extension) containing the power flow data,
%  
%   which_flow ==> 1: opf; 0: load flow
%
%   You can change options vector (see 'help mpoption' for details).
%   Uses default 'case' if 1st parameter is not given. 

%% --- define names for columns to data matrices

PQ=1; PV=2; REF=3; NONE=4; BUS_I=1; BUS_TYPE=2; PD=3; QD=4; GS=5; BS=6; BUS_AREA=7;
VM=8; VA=9; BASE_KV=10;	ZONE=11; VMAX=12; VMIN=13; LAM_P=14; LAM_Q=15; MU_VMAX=16; MU_VMIN=17; 

GEN_BUS=1; PG=2; QG=3; QMAX=4; QMIN=5; VG=6; MBASE=7; GEN_STATUS=8; PMAX=9; PMIN=10;
MU_PMAX=11; MU_PMIN=12; MU_QMAX=13; MU_QMIN=14; % MU: KT multiplier

F_BUS=1;  T_BUS=2;  BR_R=3;   BR_X=4; BR_B=5; % standard notation (in input)
RATE_A=6; RATE_B=7; RATE_C=8; TAP=9;  SHIFT=10; BR_STATUS=11; % standard notation (in input)
PF=12; QF=13; PT=14; QT=15; MU_SF=16; MU_ST=17; % MU_SF: idx of MU on MVA lim at f bus (in opf)

% --- read data 

[baseMVA, bus, gen, branch, area, gencost] = feval('case30',x); 

  
% --- convert bus numbering to internal bus numbering

i2e	= bus(:, BUS_I);
e2i = zeros(max(i2e), 1);
e2i(i2e) = [1:size(bus, 1)]';
bus(:, BUS_I)	= e2i( bus(:, BUS_I)	);
gen(:, GEN_BUS)	= e2i( gen(:, GEN_BUS)	);
branch(:, F_BUS)= e2i( branch(:, F_BUS)	);
branch(:, T_BUS)= e2i( branch(:, T_BUS)	);

% --- preliminary calculations

%[ref, pv, pq] = bustypes(bus, gen);   % sort out bus types

ref	= find(bus(:, BUS_TYPE) == REF);	%% reference bus index
pv	= find(bus(:, BUS_TYPE) == PV );	%% PV bus indices
pq	= find(bus(:, BUS_TYPE) == PQ );	%% PQ bus indices
if isempty(ref); fprinft('\n -- OOPS, no reference bus, check data! -- \n'); return; end


% -------- get constants, initial states ----------

nb = size(bus, 1);        nl = size(branch, 1);        npv     = length(pv);     npq     = length(pq);    ng = npv + 1;             
ref	= find(bus(:, BUS_TYPE) == REF);      pvpq = [1:ref-1, ref+1:size(bus, 1)];   gbus = gen(:, GEN_BUS);  

[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);   % build admitance matrix

% --- call newtonpf.m
 
   %%gen(:, PG)=x;    %%---………………………………………………………………………％

   verbose = 0;

   % compute complext bus injections
   gbus = gen(:, GEN_BUS);				% at what buses are they at?
   Sbus = -(bus(:, PD) + sqrt(-1) * bus(:, QD));	%% power injected by loads
   Sbus(gbus) = Sbus(gbus) + gen(:, PG) + sqrt(-1) * gen(:, QG)  ; %% plus generation
   Sbus = Sbus /baseMVA; % convert to p.u system

   V0	= bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
   V0(gen(:, GEN_BUS)) = gen(:, VG) ./ abs(V0(gen(:, GEN_BUS))).* V0(gen(:, GEN_BUS));

   [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, verbose);

   [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);
   
   %%gen(:,PG)
   
   %%branch(:, PF)
   Pg = gen(:, PG)/baseMVA
   Qg = gen(:, QG)/baseMVA;        Va = angle(V(pvpq));      Vm = abs(V);
  
   %%------evalution the objective function-----------%%
   
   Sbus = -(bus(:, PD) + sqrt(-1) * bus(:, QD));	                               %% power injected by loads
   Sbus(gbus) = Sbus(gbus) + gen(:, PG) + sqrt(-1) * gen(:, QG);   %% plus generation
   Sbus = Sbus /baseMVA;                                                          % convert to p.u system
   mis = V .* conj(Ybus * V) - Sbus; 		               %% evaluate power flow equations
   Sf = V(branch(:, F_BUS)) .* conj(Yf * V);	                               %% power injected at "from" bus
   St = V(branch(:, T_BUS)) .* conj(Yt * V);	                               %% power injected at "to" bus

   g1= [  real(mis);                                                  
          imag(mis) ];
   g2= [  bus(:, VMIN)         - Vm;         %% Vmim
          Vm                   - bus(:, VMAX)];     
   g3= [  gen(:, PMIN)/baseMVA - Pg;		          %% Pgmin
          Pg                   - gen(:, PMAX)/baseMVA] ;   
   g4= [  gen(:, QMIN)/baseMVA - Qg;		          %% Qgmin
          Qg                   - gen(:, QMAX)/baseMVA ] ;   
   g5= [  abs(Sf)              - branch(:, RATE_A)/baseMVA; %% brch apparent pwr lim
          abs(St)              - branch(:, RATE_A)/baseMVA ]; %% brch apparent pwr lim  
      
   max(0,g2)
   
   max(0,g3)
   
   vv=sum(max(0,g2))
      
   vg=sum(max(0,g3))
      
   vq=sum(max(0,g4))
      
   vs=sum(max(0,g5))
   
   %%------Frist  the objective---total  generator  cost ----**
   F1=0.0;
   for j = 1:ng;   
   totalcost(j) = gencost(j,5)*(Pg(j,:)*Pg(j,:))+gencost(j,6)*Pg(j,:)+gencost(j,7); 
   end
   F1=sum(totalcost)
  %% F11=F1+kv*vv+kg*vg+kq*vq+ks*vs
   
   %%-----Sencond the objective ---emission----%%
   F2=0.0;
   for j = 1:ng;     
   %%emission(j) = polyval( gencost( j, 8:10) , Pg(j,:) )+gencost( j,11)*exp(gencost( j, 12)*Pg(j,:));
   end
   %%F2=sum(emission);
   
   %%-----third the objective ---reactive power optimal----%%
   %%sum(gen(:,PG))-sum(bus(:,PD))
   %%sum(abs(Vm-1))
   w=1.5;
   %%F3=sum(gen(:,PG))-sum(bus(:,PD));  %%+w*sum(abs(Vm-1))
   
 
return;