fun_ccv.m

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

function [f, g] = fun_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)
%FUN_CCV  Evaluates objective function & constraints for OPF.
%   [f, g] = fun_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 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;

%% 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 objective function -----
%% put Pg & Qg back in gen
gen(on, PG) = Pg * baseMVA;					%% active generation in MW
gen(on, QG) = Qg * baseMVA;					%% reactive generation in MVAR

%% compute objective value
if ncv > ng				%% no free VARs
	f = sum([Cp; Cq]);
else					%% free VARs
	f = sum(Cp);
end

%%----- evaluate 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);
	
	%% rebuild Sbus
	Sbus = makeSbus(baseMVA, bus, gen);		%% net injected power in p.u.
	
	%% evaluate power flow equations
	mis = V .* conj(Ybus * V) - Sbus;
	
	%% compute branch power flows
	Sf = V(branch(:, F_BUS)) .* conj(Yf * V);	%% complex power at "from" bus (p.u.)
	St = V(branch(:, T_BUS)) .* conj(Yt * V);	%% complex power injected at "to" bus (p.u.)
	
	%% compute generator cost constraints ( costfcn @ Pg - Cp , etc.)
	Qcc = [];
	nsegs = pcost(:, N) - 1;			%% number of cost constraints for each gen
	ncc = sum(nsegs);					%% total number of cost constraints
	Pcc = zeros(ncc, 1);
	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));
		b = yy(i1) - m .* xx(i1);
		Pcc(sum(nsegs(1:(i-1))) + [1:nsegs(i)]) = ...
							m .* gen(on(i), PG) + b - Cp(i);
	end
	if ncv > ng				%% no free VARs
		nsegs = qcost(:, N) - 1;			%% number of cost constraints for each gen
		ncc = sum(nsegs);					%% total number of cost constraints
		Qcc = zeros(ncc, 1);
		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));
			b = yy(i1) - m .* xx(i1);
			Qcc(sum(nsegs(1:(i-1))) + [1:nsegs(i)]) = ...
								m .* gen(on(i), QG) + b - Cq(i);
		end
	end
	
	%% compute constraint function values
	g = [
		%% equality constraints
		real(mis);							%% active power mismatch for all buses
		imag(mis);							%% reactive power mismatch for all buses

		%% inequality constraints (variable limits, voltage & generation)
		bus(:, VMIN) - x(j5:j6);				%% lower voltage limit for var V
		x(j5:j6) - bus(:, VMAX);				%% upper voltage limit for var V
		gen(on, PMIN) / baseMVA - x(j7:j8);		%% lower generator P limit
		x(j7:j8) - gen(on, PMAX) / baseMVA;		%% upper generator P limit
		gen(on, QMIN) / baseMVA - x(j9:j10);	%% lower generator Q limit
		x(j9:j10) - gen(on, QMAX) / baseMVA;	%% upper generator Q limit

		%% inequality constraints (line flow limits)
		abs(Sf) - branch(:, RATE_A) / baseMVA;	%% branch apparent power limits (from bus)
		abs(St) - branch(:, RATE_A) / baseMVA;	%% branch apparent power limits (to bus)
		
		%% inequality cost constraints
		Pcc;
		Qcc;
	];
end

return;