function [buso, geno, brancho, f, success, info, et] = dcopf(baseMVA, ...                   bus, gen, branch, areas, gencost, mpopt)%DCOPF  Solves a DC optimal power flow.%%   [bus, gen, branch, f, success] = dcopf(casefile, mpopt)%%   [bus, gen, branch, f, success] = dcopf(baseMVA, bus, gen, branch, areas, ...%                                    gencost, mpopt)%%   [bus, gen, branch, f, success, info, et] = dcopf(casefile)%%   The data for the problem can be specified in one of 3 ways: (1) the name of%   a case file which defines the data matrices baseMVA, bus, gen, branch,%   areas and gencost, (2) a struct containing the data matrices as fields, or%   (3) the data matrices themselves.%%   The optional mpopt vector specifies MATPOWER options. Type 'help mpoption'%   for details and default values.%%   The solved case is returned in the data matrices, bus, gen and branch. Also,%   returned are the final objective function value (f) and a flag which is%   true if the algorithm was successful in finding a solution (success).%   Additional optional return values are an algorithm specific return status%   (info), elapsed time in seconds (et).%   MATPOWER%   $Id: dcopf.m,v 1.11 2005/01/14 17:20:13 ray Exp $%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales%   and Ray Zimmerman, PSERC Cornell%   Copyright (c) 1996-2005 by Power System Engineering Research Center (PSERC)%   See http://www.pserc.cornell.edu/matpower/ for more info.%%----- initialization -----%% Sort argumentsif isstr(baseMVA) | isstruct(baseMVA)  casefile = baseMVA;  if nargin == 1    mpopt = mpoption;  elseif nargin == 2    mpopt = bus;  else    error('dcopf.m: Incorrect input parameter order, number or type');  end  [baseMVA, bus, gen, branch, areas, gencost] = loadcase(casefile);else  if nargin == 6    mpopt = mpoption;  elseif nargin ~= 7    error('dcopf.m: Incorrect input parameter order, number or type');  endend%% optionsmpopt(10) = 1;  % force DC treatmentverbose = mpopt(31);npts = mpopt(14);       %% number of points to evaluate when converting                        %% polynomials to piece-wise linear%% define constantsj = sqrt(-1);%% define named indices into data 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;[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;[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...    GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, N, COST] = idx_cost;% If tables do not have multiplier/extra columns, append zero colsif size(bus,2) < MU_VMIN  bus = [bus zeros(size(bus,1),MU_VMIN-size(bus,2)) ];endif size(gen,2) < MU_QMIN  gen = [ gen zeros(size(gen,1),MU_QMIN-size(gen,2)) ];endif size(branch,2) < MU_ST  branch = [ branch zeros(size(branch,1),MU_ST-size(branch,2)) ];end% Filter out inactive generators and branches; save original bus & branchcomgen = find(gen(:,GEN_STATUS) > 0);offgen = find(gen(:,GEN_STATUS) <= 0);onbranch  = find(branch(:,BR_STATUS) ~= 0);offbranch = find(branch(:,BR_STATUS) == 0);genorg = gen;branchorg = branch;ng = size(gen,1);         % original size(gen), at least temporallygen   = gen(comgen, :);branch = branch(onbranch, :);if size(gencost,1) == ng  gencost = gencost(comgen, :);else  gencost = gencost( [comgen; comgen+ng], :);end% Renumber buses consecutively[i2e, bus, gen, branch, areas] = ext2int(bus, gen, branch, areas);%% get bus index lists of each type of bus[ref, pv, pq] = bustypes(bus, gen);%% build B matrices and phase shift injections[B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);%%-----  check/convert costs, set default algorithm  -----%% get cost model, check consistencymodel = gencost(:, MODEL);i_pwln = find(model == PW_LINEAR);i_poly = find(model == POLYNOMIAL);if any(i_pwln) & any(i_poly) & verbose    fprintf('not all generators use same cost model, all will be converted to piece-wise linear\n');endif any(i_pwln) | any(find(gencost(:, N) > 3))    formulation = 2;    %% use piecewise linear formulationelse    formulation = 1;    %% use polynomial cost formulation    if any(find(gencost(:, N) ~= 3))        error('DC opf with polynomial costs can only handle quadratic costs.');    endend%% convert polynomials to piece-wise linear (if necessary)if any(i_poly) & formulation == 2    if verbose        fprintf('converting from polynomial to piece-wise linear cost model\n');    end    [pcost, qcost] = pqcost(gencost, size(gen, 1));    i_poly = find(pcost(:, MODEL) == POLYNOMIAL);    tmp = poly2pwl(pcost(i_poly, :), gen(i_poly, PMIN), gen(i_poly, PMAX), npts);    pcost(i_poly, 1:size(tmp,2)) = tmp;    if ~isempty(qcost)        i_poly = find(qcost(:, MODEL) == POLYNOMIAL);        tmp = poly2pwl(qcost(i_poly, :), gen(i_poly, QMIN), gen(i_poly, QMAX), npts);        qcost(i_poly, 1:size(tmp,2)) = tmp;    end    gencost = [pcost; qcost];end%%-----  run opf  -----%% start timert0 = clock;%% gen 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);ng = length(on);                        %% number of generators that are turned on%% initial stateVa  = bus(:, VA) * (pi/180);Pg  = gen(on, PG) / baseMVA;%% check for costs for Qg[pcost, qcost] = pqcost(gencost, size(gen, 1), on);%% set up x along with indexingj1 = 1;         j2  = nb;               %% j1:j2    - bus V anglesj3 = j2 + 1;    j4  = j2 + ng;          %% j3:j4    - P of generatorsj5 = j4 + 1;    j6  = j4 + ng;          %% j5:j6    - Cp, cost of Pgif formulation == 2             %% piece-wise linear costs    Cp = totcost(pcost, Pg * baseMVA);else                            %% polynomial costs    Cp = [];endnc = length(Cp);                        %% number of cost variables (ng or 0)x = [Va; Pg; Cp];%% set up constraints and indexing where,   AA * x <= bb% i0 = 1;                           %% 1 - voltage angle referencei1 = 2;         i2 = nb + 1;        %% i1:i2 - P mismatch, all busesi3 = i2 + 1;    i4 = i2 + ng;       %% i3:i4 - Pmin, gen busesi5 = i4 + 1;    i6 = i4 + ng;       %% i5:i6 - Pmax, gen busesi7 = i6 + 1;    i8 = i6 + nl;       %% i7:i8 - |Pf| line limiti9 = i8 + 1;    i10 = i8 + nl;      %% i9:i10 - |Pt| line limitif formulation == 2             %% piece-wise linear costs    %% compute cost constraints [ Cp >= m * Pg + b ] => [ m * Pg - Cp <= -b ]    nsegs = pcost(:, N) - 1;            %% number of cost constraints for each gen    ncc = sum(nsegs);                   %% total number of cost constraints    Acc = sparse(ncc, nb+ng+nc);    bcc = 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))';        k1 = 1:nsegs(i);        k2 = 2:(nsegs(i) + 1);        m = (yy(k2) - yy(k1)) ./ (xx(k2) - xx(k1));        b = yy(k1) - m .* xx(k1);        Acc(sum(nsegs(1:(i-1))) + [1:nsegs(i)], nb+i)       = m * baseMVA;        Acc(sum(nsegs(1:(i-1))) + [1:nsegs(i)], nb+ng+i)    = -ones(nsegs(i), 1);        bcc(sum(nsegs(1:(i-1))) + [1:nsegs(i)])             = -b;    endelse                            %% polynomial costs    Acc = sparse(0, nb+ng+nc);    bcc = zeros(0, 1);endAA = [        sparse(1, ref, 1, 1, nb+ng+nc);                     %% reference angle        B,  -sparse(gen(on, GEN_BUS), 1:ng, ones(ng, 1), nb, ng), sparse(nb, nc);                                                            %% real power flow eqns        sparse(ng, nb), -speye(ng, ng), sparse(ng, nc);     %% lower limit on Pg        sparse(ng, nb), speye(ng, ng), sparse(ng, nc);      %% upper limit on Pg        Bf, sparse(nl, ng+nc);                              %% flow limit on Pf        -Bf, sparse(nl, ng+nc);                             %% flow limit on Pt        Acc;                                                %% cost constraints];bb = [        Va(ref);                                            %% reference angle        -(bus(:, PD) + bus(:, GS)) / baseMVA - Pbusinj;     %% real power flow eqns        -gen(on, PMIN) / baseMVA;                           %% lower limit on Pg        gen(on, PMAX) / baseMVA;                            %% upper limit on Pg        branch(:, RATE_A) / baseMVA - Pfinj;                %% flow limit on Pf        branch(:, RATE_A) / baseMVA + Pfinj;                %% flow limit on Pt        bcc;                                                %% cost constraints];%% set up objective function of the form:  0.5 * x'*H*x + c'*xif formulation == 2             %% piece-wise linear costs    H = sparse(nb+ng+nc, nb+ng+nc);    c = [   zeros(nb+ng, 1);              ones(nc, 1)     ];else                            %% polynomial costs    polycf = pcost(:, COST:COST+2) * diag([ baseMVA^2 baseMVA 1]);  %% coeffs for Pg in p.u.    H = sparse(j3:j4, j3:j4, 2*polycf(:, 1), nb+ng, nb+ng );    c = [   zeros(nb, 1);              polycf(:, 2)    ];end%% run QP solvermpopt(15)   = nb + 1;               %% set number of equality constraintsif verbose > 1                      %% print QP progress for verbose levels 2 & 3    qpverbose = 1;else    qpverbose = -1;endif ~have_fcn('sparse_qp') | mpopt(51) == 0 %% don't use sparse matrices    AA = full(AA);    H = full(H);end[x, lambda, how, success] = mp_qp(H, c, AA, bb, [], [], x, mpopt(15), qpverbose);info = success;%% update solution dataVa = x(j1:j2);Pg = x(j3:j4);f = sum(totcost(pcost, Pg * baseMVA));%%-----  calculate return values  -----%% update voltages & generator outputsbus(:, VM) = ones(nb, 1);bus(:, VA) = Va * 180 / pi;gen(on, PG) = Pg * baseMVA;%% compute flows etc.branch(:, [QF, QT]) = zeros(nl, 2);branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;branch(:, PT) = -branch(:, PF);%% update lambda's and mu'sbus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);bus(:, LAM_P)       = lambda(i1:i2) / baseMVA;gen(on, MU_PMIN)    = lambda(i3:i4) / baseMVA;gen(on, MU_PMAX)    = lambda(i5:i6) / baseMVA;branch(:, MU_SF)    = lambda(i7:i8) / baseMVA;branch(:, MU_ST)    = lambda(i9:i10) / baseMVA;% convert to original external bus ordering[bus, gen, branch, areas] = int2ext(i2e, bus, gen, branch, areas);% Now create output matrices with all lines, all generators, committed and% non-committedgeno = genorg;brancho = branchorg;geno(comgen, : ) = gen;brancho(onbranch, :)  = branch;% And zero out appropriate fields of non-comitted generators and linestmp = zeros(length(offgen), 1);geno(offgen, PG) = tmp;geno(offgen, QG) = tmp;geno(offgen, MU_PMAX) = tmp;geno(offgen, MU_PMIN) = tmp;tmp = zeros(length(offbranch), 1);brancho(offbranch, PF) = tmp;brancho(offbranch, QF) = tmp;brancho(offbranch, PT) = tmp;brancho(offbranch, QT) = tmp;brancho(offbranch, MU_SF) = tmp;brancho(offbranch, MU_ST) = tmp;%% compute elapsed timeet = etime(clock, t0);if (nargout == 0) & (info > 0)  printpf(baseMVA, bus, geno, brancho, f, success, et, 1, mpopt);endif nargout, buso = bus; endreturn;