function [f, g] = fun_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)%FUN_CCV  Evaluates objective function & constraints for OPF.%   [f, g] = fun_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)%   MATPOWER%   $Id: fun_ccv.m,v 1.6 2004/09/08 12:37:35 ray Exp $%   by 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;%% constantj = sqrt(-1);%% generator infoon = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?%% 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 on%% check for costs for Qg[pcost, qcost] = pqcost(gencost, size(gen, 1), on);if isempty(qcost)       %% set number of cost variables    ncv = ng;           %% only Cpelse    ncv = 2 * ng;       %% Cp and Cqend%% set up indexing for xj1 = 1;         j2  = npv;              %% j1:j2    - V angle of pv busesj3 = j2 + 1;    j4  = j2 + npq;         %% j3:j4    - V angle of pq 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 generatorsj11 = j10 + 1;  j12 = j10 + ng;         %% j11:j12  - Cp, cost of Pg%% grab Pg & Qg and their costsPg = x(j7:j8);                              %% active generation in p.u.Qg = x(j9:j10);                             %% reactive generation in p.u.Cp = x(j11:j12);                            %% active costs in $/hrif ncv > ng             %% no free VArs    j13 = j12 + 1;  j14 = j12 + ng;         %% j13:j14  - Cq, cost of Qg    Cq = x(j13:j14);                        %% reactive costs in $/hrend%%----- evaluate objective function -----%% put Pg & Qg back in gengen(on, PG) = Pg * baseMVA;                 %% active generation in MWgen(on, QG) = Qg * baseMVA;                 %% reactive generation in MVAr%% compute objective valueif ncv > ng             %% no free VArs    f = sum([Cp; Cq]);else                    %% free VArs    f = sum(Cp);end%%----- 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    br = find(branch(:, BR_STATUS));    Sf = 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 generator cost constraints ( costfcn @ Pg - Cp , etc.)    Qcc = [];    nsegs = pcost(:, N) - 1;            %% number of cost constraints for each gen    ncc = sum(nsegs);                   %% total number of cost constraints    Pcc = zeros(ncc, 1);    for i = 1:ng        xx = pcost(i,       COST:2:( COST + 2*(nsegs(i))    ))';        yy = pcost(i,   (COST+1):2:( COST + 2*(nsegs(i)) + 1))';        i1 = 1:nsegs(i);        i2 = 2:(nsegs(i) + 1);        m = (yy(i2) - yy(i1)) ./ (xx(i2) - xx(i1));        b = yy(i1) - m .* xx(i1);        Pcc(sum(nsegs(1:(i-1))) + [1:nsegs(i)]) = ...                            m .* gen(on(i), PG) + b - Cp(i);    end    if ncv > ng             %% no free VArs        nsegs = qcost(:, N) - 1;            %% number of cost constraints for each gen        ncc = sum(nsegs);                   %% total number of cost constraints        Qcc = zeros(ncc, 1);        for i = 1:ng            xx = qcost(i,       COST:2:( COST + 2*(nsegs(i))    ))';            yy = qcost(i,   (COST+1):2:( COST + 2*(nsegs(i)) + 1))';            i1 = 1:nsegs(i);            i2 = 2:(nsegs(i) + 1);            m = (yy(i2) - yy(i1)) ./ (xx(i2) - xx(i1));            b = yy(i1) - m .* xx(i1);            Qcc(sum(nsegs(1:(i-1))) + [1:nsegs(i)]) = ...                                m .* gen(on(i), QG) + b - Cq(i);        end    end        %% compute line flow constraints    if mpopt(24) == 1   %% branch active power limits        flow_limit = [            real(Sf) - branch(br, RATE_A) / baseMVA;    %% from bus            real(St) - branch(br, RATE_A) / baseMVA;    %% to bus        ];    else                %% branch apparent power limits        flow_limit = [            abs(Sf) - branch(br, RATE_A) / baseMVA;     %% from bus            abs(St) - branch(br, RATE_A) / baseMVA;     %% to bus        ];    end        %% 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) - Vm;                  %% lower voltage limit for var V        Vm - bus(:, VMAX);                  %% upper voltage limit for var V        gen(on, PMIN) / baseMVA - Pg;       %% lower generator P limit        Pg - gen(on, PMAX) / baseMVA;       %% upper generator P limit        gen(on, QMIN) / baseMVA - Qg;       %% lower generator Q limit        Qg - gen(on, QMAX) / baseMVA;       %% upper generator Q limit        %% inequality constraints (line flow limits)        flow_limit;                %% inequality cost constraints        Pcc;        Qcc;    ];endreturn;