function [df, dg, d2f] = OTSgra(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)
%OTSGRA  Evaluates gradients of objective function & constraints for OPF.
%   [df, dg] = OTSgra(x, baseMVA, bus, gen, gencost, branch, Ybus,
%                          Yf, Yt, V, ref, pv, pq, mpopt)
%   Also, if a third output argument is specified, it will compute the 2nd
%   derivative matrix for the objective function.

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

%% 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 partials of objective function -----
%% generator info
on = find(gen(:, GEN_STATUS));				%% which generators are on?
gbus = gen(on, GEN_BUS);					%% what buses are they at?

%% compute values of objective function partials
[pcost, qcost] = pqcost(gencost, size(gen, 1), on);
df_dPg = zeros(ng, 1);
df_dQg = zeros(ng, 1);
for i = 1:ng
	df_dPg(i) = polyval(polyder( pcost( i, COST:(COST+pcost(i, N)-1) ) ), Pg(i)*baseMVA) * baseMVA;	%% w.r.t p.u. Pg
end
if ~isempty(qcost)					%% Qg is not free
	for i = 1:ng
		df_dQg(i) = polyval(polyder( qcost( i, COST:(COST+qcost(i, N)-1) ) ), Qg(i)*baseMVA) * baseMVA;	%% w.r.t p.u. Qg
	end
end
df = [	zeros(j6, 1);				%% partial w.r.t. Va & Vm
		df_dPg;						%% partial w.r.t. Pg
		df_dQg	];					%% partial w.r.t. Qg

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

	%% make full so optimization toolbox doesn't go wacky
	dg = full(dg);

	%% compute 2nd derivative of cost
	if nargin > 2
		d2f_dPg2 = zeros(ng, 1);
		d2f_dQg2 = zeros(ng, 1);
		for i = 1:ng
			d2f_dPg2(i) = polyval(polyder(polyder( pcost( i, COST:(COST+pcost(i, N)-1) ) )), Pg(i)*baseMVA) * baseMVA^2;	%% w.r.t p.u. Pg
		end
		if ~isempty(qcost)					%% Qg is not free
			for i = 1:ng
				d2f_dQg2(i) = polyval(polyder(polyder( qcost( i, COST:(COST+qcost(i, N)-1) ) )), Qg(i)*baseMVA) * baseMVA^2;	%% w.r.t p.u. Qg
			end
		end
		i = [j7:j10]';
		d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2]);
	end
end

return;