fun_std.m

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

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

%% constant
j = sqrt(-1);

%% 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

%% 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

%% grab Pg & Qg
Pg = x(j7:j8);								%% active generation in p.u.
Qg = x(j9:j10);								%% reactive generation in p.u.

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

%% 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
[pcost, qcost] = pqcost(gencost, size(gen, 1), on);
f = sum(	[totcost(pcost, gen(on, PG)); ...	%% cost of Pg
			 totcost(qcost, gen(on, QG)) ]	);	%% cost of Qg, empty if no qcost

%%----- 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 injected at "from" bus (p.u.)
	St = V(branch(:, T_BUS)) .* conj(Yt * V);	%% complex power injected at "to" bus (p.u.)
	
	%% 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)
	];
end

return;