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 parametersnb = 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 & QgPg = x(pgbas:pgend);                %% active generation in p.u.Qg = x(qgbas:qgend);                %% reactive generation in p.u.%% put Pg & Qg back in gengen(:, PG) = Pg * baseMVA;          %% active generation in MWgen(:, 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 costsipwl = find(gencost(:, MODEL) == PW_LINEAR);     % pw_lin MW and MVAr costsxx = [ gen(:, PG); gen(:, QG)];                  % if ~isempty(ipol)  f = sum( totcost(gencost(ipol, :), xx(ipol)) );  %% cost of poly P or Qelse  f = 0;endif 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.endreturn% 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]);endreturn;