function [g, geq, dg, dgeq] = consfmin(x, baseMVA, bus, gen, gencost, branch, areas, Ybus, Yf, Yt, mpopt, parms, ccost)%CONSFMIN  Evaluates nonlinear constraints and their Jacobian for OPF.%   [g, geq, dg, dgeq] = consfmin(x, baseMVA, bus, gen, gencost, ...%                                  branch, areas, Ybus, Yf, Yt, mpopt)%   MATPOWER%   $Id: consfmin.m,v 1.4 2004/09/07 18:27:37 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 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;%% constantj = sqrt(-1);%% generator infoon = find(gen(:, GEN_STATUS) > 0);    %% which generators are on?gbus = gen(on, GEN_BUS);        %% what buses are they at?%% sizes of thingsnb = size(bus, 1);nl = size(branch, 1);%npv  = length(pv);%npq  = length(pq);ng = length(on);            %% number of generators that are turned onbr = find(branch(:, BR_STATUS));nbr = length(br);%% set up indexing for xj1 = 1;      j4 = nb;             %% j1:j4 - V angle of all busesj5 = j4 + 1;  j6  = j4 + nb;      %% j5:j6  - V mag of all busesj7 = j6 + 1;  j8  = j6 + ng;      %% j7:j8  - P of generatorsj9 = j8 + 1;  j10  = j8 + ng;      %% j9:j10  - Q of generatorsnb = 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(j7:j8);                %% active generation in p.u.Qg = x(j9:j10);                %% reactive generation in p.u.%% put Pg & Qg back in gengen(on, PG) = Pg * baseMVA;          %% active generation in MWgen(on, QG) = Qg * baseMVA;          %% reactive generation in MVAr %% rebuild SbusSbus = makeSbus(baseMVA, bus, gen);    %% net injected power in p.u.%% ----- evaluate constraints -----%% reconstruct VVa = zeros(nb, 1);Va = x(j1:j4);Vm = x(j5:j6);V = Vm .* exp(j * Va);%% evaluate power flow equationsmis = V .* conj(Ybus * V) - Sbus;%% compute branch power flowsSf = V(branch(br, F_BUS)) .* conj(Yf(br, :) * V);  %% complex power injected at "from" bus (p.u.)St = V(branch(br, T_BUS)) .* conj(Yt(br, :) * V);  %% complex power injected at "to" bus (p.u.)%% compute constraint function values%% First the equality constraints (power flow)geq = [ real(mis);              %% active power mismatch for all buses        imag(mis) ];            %% reactive power mismatch for all buses%% Then the inequality constraints (branch limits)if mpopt(24) == 1   %% P limit (Pan Wei)  g = [ real(Sf) - branch(br, RATE_A)/baseMVA;  %% branch apparent power limits (from bus)        real(St) - branch(br, RATE_A)/baseMVA ]; %% branch apparent power limits (to bus)else                %% |S| limit  g = [ abs(Sf) - branch(br, RATE_A)/baseMVA;  %% branch apparent power limits (from bus)        abs(St) - branch(br, RATE_A)/baseMVA ]; %% branch apparent power limits (to bus)end%%----- evaluate partials of constraints -----if nargout > 2  %% 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);  %% construct Jacobian of equality constraints (load flows) and transpose it  dgeq = [    %% equality constraints    real(dSbus_dVa), real(dSbus_dVm), ...          real(dSbus_dPg), real(dSbus_dQg), sparse(nb, ny+nz);   %% P mismatch    imag(dSbus_dVa), imag(dSbus_dVm), ...          imag(dSbus_dPg), imag(dSbus_dQg), sparse(nb, ny+nz);   %% Q mismatch   ]';   %% construct Jacobian of inequality constraints (branch limits)  %% and transpose it so fmincon likes it  if mpopt(24) == 1     %% P limit (Pan Wei)    dg = [      real(dSf_dVa(br,:)), real(dSf_dVm(br,:)), sparse(nbr,2*ng+ny+nz);		%% Pf limit      real(dSt_dVa(br,:)), real(dSt_dVm(br,:)), sparse(nbr,2*ng+ny+nz);		%% Pt limit    ]';  else                  %% |S| limit    dg = [      dAf_dVa(br,:), dAf_dVm(br,:), sparse(nbr,2*ng+ny+nz);  %% |Sf| limit      dAt_dVa(br,:), dAt_dVm(br,:), sparse(nbr,2*ng+ny+nz);  %% |St| limit    ]';  end    %% the following lines can be removed if/when fmincon  %% supports sparse matrices for non-linearly constrained problems  dgeq = full(dgeq);  dg = full(dg);endreturn;