function [f,df,ddf] = costfmin(x, baseMVA, bus, gen, gencost, branch, areas, Ybus, Yf, Yt, mpopt, parms, ccost)
%COSTFMIN  Evaluates objective function, gradient and Hessian for OPF.
%   [f, df, ddf] = costfmin(x, baseMVA, bus, gen, gencost, branch, areas, ...
%                           Ybus, Yf, Yt, mpopt, parms, ccost)

%   MATPOWER
%   $Id: costfmin.m,v 1.2 2004/08/23 20:56:09 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   and Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ 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;

%% unpack parameters

nb = parms(1);
ng = parms(2);
nl = parms(3);
ny = parms(4);
nx = parms(5);
nvl = parms(6);
nz = parms(7);
nxyz = parms(8);
thbas = parms(9);
thend = parms(10);
vbas = parms(11);
vend = parms(12);
pgbas = parms(13);
pgend = parms(14);
qgbas = parms(15);
qgend = parms(16);
ybas = parms(17);
yend = parms(18);
zbas = parms(19);
zend = parms(20);
pmsmbas = parms(21);
pmsmend = parms(22);
qmsmbas = parms(23);
qmsmend = parms(24);
sfbas = parms(25);
sfend = parms(26);
stbas = parms(27);
stend = parms(28);

%% grab Pg & Qg
Pg = x(pgbas:pgend);                %% active generation in p.u.
Qg = x(qgbas:qgend);                %% reactive generation in p.u.

%% put Pg & Qg back in gen
gen(:, PG) = Pg * baseMVA;          %% active generation in MW
gen(:, QG) = Qg * baseMVA;          %% reactive generation in MVAr

%%----- evaluate objective function -----

%% compute objective value
%[pcost, qcost] = pqcost(gencost, size(gen, 1));
%pqcoststack = [pcost; qcost];
% use totcost only on polynomial cost; in the minimization problem
% formulation, pwl cost is the sum of the y variables.
ipol = find(gencost(:, MODEL) == POLYNOMIAL);    % poly MW and MVAr costs
ipwl = find(gencost(:, MODEL) == PW_LINEAR);     % pw_lin MW and MVAr costs
xx = [ gen(:, PG); gen(:, QG)];                  % 
if ~isempty(ipol)
  f = sum( totcost(gencost(ipol, :), xx(ipol)) );  %% cost of poly P or Q
else
  f = 0;
end
if ny+nz  > 0
  f = f + ccost * x;
end

%%----- evaluate cost gradient -----
if nargout > 1
  df_dPgQg = zeros(2*ng, 1);
  for i = ipol'
    df_dPgQg(i)= baseMVA * ...  %% w.r.t p.u. Pg
     polyval(polyder(gencost(i, COST:(COST+gencost(i, N)-1))), xx(i));
  end
  df = [  zeros(vend, 1);  % partial w.r.t. Va & Vm
          df_dPgQg;        % partial w.r.t. polynomial cost Pg and Qg
          zeros(ny+nz,1);  
  ];
  df = df + ccost';  % As in MINOS, the linear cost row is additive wrt
                     % any nonlinear cost.
end

return

% Currently fmincon won't use the Hessian for medium scale problems;
% when it does, or when the LargeScale methods support g, geq, A, Aeq,
% the following must be fixed for mixed poly/pwl costs.

%% ---- evaluate cost Hessian -----
if nargout > 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 = [pgbas:qgend]';
  d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2]);
end

return;

