grad_ccv.m

关于电力系统的计算程序!最新的改进法计算潮流程序！用matlab实现的！ 适用于电力行业的相关从业人员！

function [df, dg] = grad_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)
%GRAD_CCV  Evaluates gradients of objective function & constraints for OPF.
%   [df, dg] = grad_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)

%   MATPOWER Version 2.0
%   by Ray Zimmerman, PSERC Cornell    12/12/97
%   Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/ for more info.

%%----- initialize -----
%% define named indices into data, branch matrices
[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;

%% options
alg = mpopt(11);

%% constant
j = sqrt(-1);

%% generator info
on = find(gen(:, GEN_STATUS));				%% which generators are on?
gbus = gen(on, GEN_BUS);					%% what buses are they at?

%% sizes of things
nb = size(bus, 1);
nl = size(branch, 1);
npv	= length(pv);
npq	= length(pq);
ng = npv + 1;			%% number of generators that are turned on

%% check for costs for Qg
[pcost, qcost] = pqcost(gencost, size(gen, 1), on);
if isempty(qcost)		%% set number of cost variables
	ncv = ng;			%% only Cp
else
	ncv = 2 * ng;		%% Cp and Cq
end

%% set up indexing for x
j1 = 1;			j2	= npv;				%% j1:j2	- V angle of pv buses
j3 = j2 + 1;	j4	= j2 + npq;			%% j3:j4	- V angle of pq buses
j5 = j4 + 1;	j6	= j4 + nb;			%% j5:j6	- V mag of all buses
j7 = j6 + 1;	j8	= j6 + ng;			%% j7:j8	- P of generators
j9 = j8 + 1;	j10	= j8 + ng;			%% j9:j10	- Q of generators
j11 = j10 + 1;	j12	= j10 + npv + 1;	%% j11:j12	- Cp, cost of Pg

%% grab Pg & Qg and their costs
Pg = x(j7:j8);								%% active generation in p.u.
Qg = x(j9:j10);								%% reactive generation in p.u.
Cp = x(j11:j12);							%% active costs in $/hr
if ncv > ng				%% no free VARs
	j13 = j12 + 1;	j14	= j12 + npv + 1;	%% j13:j14	- Cq, cost of Qg
	Cq = x(j13:j14);						%% reactive costs in $/hr
end

%%----- evaluate partials of objective function -----
%% compute values of objective function partials
df = [	zeros(j10, 1);				%% partial w.r.t. Va, Vm, Pg, Qg
		ones(ncv, 1)	];			%% partial w.r.t. Cp (and Cq)

%%----- evaluate partials of constraints -----
if nargout > 1
	%% reconstruct V
	Va = zeros(nb, 1);
	Va([ref; pv; pq]) = [angle(V(ref)); x(j1:j2); x(j3:j4)];
	Vm = x(j5:j6);
	V = Vm .* exp(j * Va);
	
	%% compute partials of injected bus powers
	[dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);		%% w.r.t. V
	dSbus_dPg = sparse(gbus, 1:ng, -1, nb, ng);		%% w.r.t. Pg
	dSbus_dQg = sparse(gbus, 1:ng, -j, nb, ng);		%% w.r.t. Qg
	
	%% compute partials of line flows w.r.t. V
	[dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);

	%% line limits are w.r.t apparent power, so compute partials of apparent power
	[dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
						dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St);	

	%% cost constraints w.r.t everything ( d(costfcn @ Pg - Cp) , etc.)
	dQcc_dQg = zeros(0, ng);
	dQcc_dCq = [];
	nsegs = pcost(:, N) - 1;			%% number of cost constraints for each gen
	nPcc = sum(nsegs);					%% total number of cost constraints
	dPcc_dPg = sparse([], [], [], nPcc, ng, nPcc);	%% nPcc x ng
	dPcc_dCp = sparse([], [], [], nPcc, ng, nPcc);	%% nPcc x ng
	for i = 1:ng
		xx = pcost(i,		COST:2:( COST + 2*(nsegs(i))	))';
		yy = pcost(i,	(COST+1):2:( COST + 2*(nsegs(i)) + 1))';
		i1 = 1:nsegs(i);
		i2 = 2:(nsegs(i) + 1);
		m = (yy(i2) - yy(i1)) ./ (xx(i2) - xx(i1));
		ii = sum(nsegs(1:(i-1))) + [1:nsegs(i)];
		dPcc_dPg(ii, i) = m * baseMVA;
		dPcc_dCp(ii, i) = -1 * ones(nsegs(i), 1);
	end
	nQcc = 0;
	if ncv > ng				%% no free VARs
		nsegs = qcost(:, N) - 1;			%% number of cost constraints for each gen
		nQcc = sum(nsegs);					%% total number of cost constraints
		dQcc_dQg = sparse([], [], [], nQcc, ng, nQcc);	%% nQcc x ng
		dQcc_dCq = sparse([], [], [], nQcc, ng, nQcc);	%% nQcc x ng
		for i = 1:ng
			xx = qcost(i,		COST:2:( COST + 2*(nsegs(i))	))';
			yy = qcost(i,	(COST+1):2:( COST + 2*(nsegs(i)) + 1))';
			i1 = 1:nsegs(i);
			i2 = 2:(nsegs(i) + 1);
			m = (yy(i2) - yy(i1)) ./ (xx(i2) - xx(i1));
			ii = sum(nsegs(1:(i-1))) + [1:nsegs(i)];
			dQcc_dQg(ii, i) = m * baseMVA;
			dQcc_dCq(ii, i) = -1 * ones(nsegs(i), 1);
		end
	end
	%%     [  dcc_dV      dcc_dPg   dcc_dQg        dcc_dCp         dcc_dCq ]
	dPcc = [sparse(nPcc,j6), dPcc_dPg, sparse(nPcc,ng), dPcc_dCp, sparse(nPcc,ncv-ng)];
	dQcc = [sparse(nQcc,j8), dQcc_dQg, sparse(nQcc,ng), dQcc_dCq];

	%% evaluate partials of constraints
	dg = [
		%% equality constraints
		real(dSbus_dVa(:,[pv;pq])), real(dSbus_dVm), ...
			real(dSbus_dPg), real(dSbus_dQg), sparse(nb,ncv); 	%% P mismatch
		imag(dSbus_dVa(:,[pv;pq])), imag(dSbus_dVm), ...
			imag(dSbus_dPg), imag(dSbus_dQg), sparse(nb,ncv); 	%% Q mismatch
	
		%% inequality constraints (variable limits, voltage & real generation)
		sparse(nb,j4), -speye(nb,nb), sparse(nb,2*ng+ncv);		%% Vmin for var V
		sparse(nb,j4),  speye(nb,nb), sparse(nb,2*ng+ncv);		%% Vmax for var V
		sparse(ng,j6), -speye(ng,ng), sparse(ng,ng+ncv);		%% Pmin for generators
		sparse(ng,j6),  speye(ng,ng), sparse(ng,ng+ncv);		%% Pmax for generators
		sparse(ng,j8), -speye(ng,ng), sparse(ng,ncv);			%% Qmin for generators
		sparse(ng,j8),  speye(ng,ng), sparse(ng,ncv);			%% Qmax for generators
	
		%% inequality constraints (reactive generation & line flow limits)
		dAf_dVa(:,[pv;pq]), dAf_dVm, sparse(nl,2*ng+ncv);		%% |Sf| limit
		dAt_dVa(:,[pv;pq]), dAt_dVm, sparse(nl,2*ng+ncv);		%% |St| limit

		%% inequality constraints on generator cost
		dPcc;
		dQcc;
	]';
	
	%% make full so optimization toolbox doesn't go wacky
	dg = full(dg);
end

return;