quadratic_bid_ninebus.m

可以进行多节点多平衡节点的电力系统潮流计算

function quadratic_results = quadratic_pro(casename)

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
	VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
	GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;

[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
	RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch;

[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, N, COST] = idx_cost;


%% read data & convert to internal bus numbering
[baseMVA, bus, gen, branch, area, gen_bid, load_bid] = feval(casename);
[i2e, bus, gen, branch] = ext2int(bus, gen, branch);


%% get bus index lists of each type of bus
[ref, pv, pq] = bustypes(bus, gen);

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

alg = 1;
%%%发电机成本系数
%if gencost(1, MODEL) == 2
%    c2 = 5;
%    c1 = 6;
%    c0 = 7;
%end


%% start timer
t0 = clock;


node = bus(:, BUS_I);
node_n = length(node);

gen_node = gen(:, GEN_BUS);     %发电机节点号
gen_n= length(gen_node);        %发电机节台数

branch_linef = branch(:,F_BUS);
branch_n = length(branch_linef);   %支路数



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Bp = makeB(baseMVA, bus, branch, alg);
Bp = Bp(2:node_n,2:node_n);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  min     1/2 x'H x + f'x
%%  Sub to  A . x <= b
%%          Aeq . x  = beq
%%          lb <= x <= ub
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%temp = Bp(:, [pv; pq])';
%Bp = temp(:, [pv; pq])';                    %去掉V,theta节点
%%－－－－－－－－－－准备参数－－－－－－－－－－－－－－－－－－－－－－－－－－

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%准备 H ,f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n_zero = zeros(node_n-1,1);     %比节点数小1的列向量
H_n = gen_n + node_n + node_n - 1;
H_zero = zeros(H_n, H_n);
gencost_c2 = gen_bid(:, 2);
loadcost_c2 = load_bid(:,3);
for(i=1:1:gen_n)
    H_zero(i,i) = gencost_c2(i);
end
for(i=1:1:node_n)
    H_zero(gen_n+i,gen_n+i) = -loadcost_c2(i);
end
H = (baseMVA^2)*H_zero;         
f = baseMVA*[gen_bid(:, 1);-load_bid(:,2);n_zero];                %成本第二项，c1为列向量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%准备Aeq, beq%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pg_n = ones(1,gen_n);                       %gen_n:发电机台数,作为发电与负荷等式Pg系数矩阵
pd_n = ones(1,node_n);

pg_co = zeros(node_n,gen_n);                %先形成一个pg系数的0矩阵
for(i=1:1:gen_n)                            %发电机对应节点不符
    pg_co(gen_node(i), i) = -1;             %添加系数
end
%pg_co(:, 1) = 0;
%pg_co(:, ref) = 0;                          %把关于平衡节点的列系数变为0
%%%%%%%%%%%%%%%%%%%%%%%%%Ａ区，或整个区域％％％％％％
pg_co = pg_co(2:node_n, :);                 %去掉关于平衡节点的行，直接取去掉第一行的矩阵，默认为节点1是平衡节点



pd_co = eye(node_n);
pd_co = pd_co(2:node_n, :);


net_Aeq = [-pg_n, pd_n, n_zero'];           %网络平衡系数
node_Aeq = [pg_co,pd_co,Bp];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
net_b = 0;                               %平衡方程等式右边为零，忽略网损original
node_b = n_zero;                         %%比节点数小1的零列向量original



Aeq = [net_Aeq; node_Aeq];
beq = [net_b; node_b];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%准备A, b %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A_theta = zeros(2*(branch_n),node_n);%theta矩阵形式
A_pg = zeros(2*(branch_n),(gen_n+node_n));
branch_fbus = branch(:,F_BUS);%（支路起始节点号）
branch_tbus = branch(:,T_BUS);%（支路末节点号）
%赋值
for ii = 1:1:branch_n
    A_theta(ii,branch_fbus(ii)) = 1;
    A_theta(ii,branch_tbus(ii)) = -1;
end
for jj = 1:1:branch_n
    A_theta(branch_n+jj,branch_fbus(jj)) = -1;
    A_theta(branch_n+jj,branch_tbus(jj)) = 1;
end

%去掉平衡节点的theta
A_theta = A_theta(:,2:node_n);
b_once = (branch(:,RATE_A)/baseMVA).*branch(:,BR_X);

A = [A_pg,A_theta];
b = [b_once;b_once];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%准备 lb,ub %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta_max = ones(node_n-1,1); %用有名值算，theta扩大baseMVA搜索

lb = [gen_bid(:,3)/baseMVA; load_bid(:,4)/baseMVA; -theta_max];            %下限
ub = [gen_bid(:,4)/baseMVA; load_bid(:,5)/baseMVA; theta_max];             %上限
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
%x0 = [gen(:, PG);bus(:, PD);n_zero]/baseMVA;
[x,fval,exitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub);%有约束
%[x,fval,exitflag,output,lambda] = quadprog(H,f,[],[],Aeq,beq,lb,ub);%无约束



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%不计网损结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = length(x);
socialcost = fval
gen_MW = x(1:gen_n)*baseMVA
load_MW = x(gen_n+1:gen_n+node_n)*baseMVA
theta = x(gen_n+node_n+1:n);


theta = [0;theta]
branch_p = (theta(branch(:,F_BUS)) - theta(branch(:,T_BUS))) ./ branch(:,BR_X);
%branch_p = 100*branch_p;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%计及网损%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T = A.P 形成A%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Bp_inverse = inv(Bp);
%B_inverse = [n_zero';Bp_inverse];
%A_L = B_inverse(branch(:,F_BUS),:) - B_inverse(branch(:,T_BUS),:);
%for i = 1:1:branch_n
%    A_L(i,:) = A_L(i,:)/branch(i,BR_X);
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 网损微增率%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%R_T = 2 * (branch(:,BR_R) .* branch_p);
%for i=1:1:29
%    loss_increment(i) = sum(R_T.*A_L(:,i));
%end
%loss_increment = [0,loss_increment];


%loss1 = sum(branch(:,BR_R).*(branch_p.*branch_p));


%theta_ij = theta(branch(:,F_BUS)) - theta(branch(:,T_BUS));
%ploss = (theta_ij.*theta_ij).*( branch(:,BR_R) ./ (branch(:,BR_R).*branch(:,BR_R)+branch(:,BR_X).*branch(:,BR_X)));
%loss2 = sum(ploss);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


et = etime(clock, t0);

computingtime = et

branch_p_a = baseMVA*branch_p
branch_p = [branch(:,RATE_A),baseMVA*branch_p]%支路额定功率－支路实际功率

node_price(1) = lambda.eqlin(1)/baseMVA;    %参考节点电价
for(i=2:1:node_n)
    node_price(i) = (lambda.eqlin(1) + lambda.eqlin(i))/baseMVA;
end

multiplier = [bus(:,BUS_I),lambda.eqlin]
nodeprice = [bus(:,BUS_I),node_price']

%[x,fval,exitflag,output,lambda] = quadprog(H,f,A,b,[],[],lb,ub,x0)
return