Description of dcopf

0001 function [results, success, raw] = dcopf_solver(om, mpopt)
0002 %DCOPF_SOLVER  Solves a DC optimal power flow.
0003 %
0004 %   [RESULTS, SUCCESS, RAW] = DCOPF_SOLVER(OM, MPOPT)
0005 %
0006 %   Inputs are an OPF model object and a MATPOWER options struct.
0007 %
0008 %   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.
0009 %
0010 %   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
0011 %   branch, gen, gencost fields, along with the following additional
0012 %   fields:
0013 %       .order      see 'help ext2int' for details of this field
0014 %       .x          final value of optimization variables (internal order)
0015 %       .f          final objective function value
0016 %       .mu         shadow prices on ...
0017 %           .var
0018 %               .l  lower bounds on variables
0019 %               .u  upper bounds on variables
0020 %           .lin
0021 %               .l  lower bounds on linear constraints
0022 %               .u  upper bounds on linear constraints
0023 %
0024 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0025 %
0026 %   RAW         raw output in form returned by MINOS
0027 %       .xr     final value of optimization variables
0028 %       .pimul  constraint multipliers
0029 %       .info   solver specific termination code
0030 %       .output solver specific output information
0031 %
0032 %   See also OPF, QPS_MATPOWER.
0033 
0034 %   MATPOWER
0035 %   Copyright (c) 2000-2015 by Power System Engineering Research Center (PSERC)
0036 %   by Ray Zimmerman, PSERC Cornell
0037 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0038 %
0039 %   $Id: dcopf_solver.m 2644 2015-03-11 19:34:22Z ray $
0040 %
0041 %   This file is part of MATPOWER.
0042 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0043 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0044 
0045 %%----- initialization -----
0046 %% define named indices into data matrices
0047 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0048     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0049 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0050     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0051     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0052 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0053     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0054     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0055 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0056 
0057 %% options
0058 alg = upper(mpopt.opf.dc.solver);
0059 
0060 %% unpack data
0061 mpc = get_mpc(om);
0062 [baseMVA, bus, gen, branch, gencost] = ...
0063     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0064 cp = get_cost_params(om);
0065 [N, H, Cw] = deal(cp.N, cp.H, cp.Cw);
0066 fparm = [cp.dd cp.rh cp.kk cp.mm];
0067 Bf = userdata(om, 'Bf');
0068 Pfinj = userdata(om, 'Pfinj');
0069 [vv, ll] = get_idx(om);
0070 
0071 %% problem dimensions
0072 ipol = find(gencost(:, MODEL) == POLYNOMIAL); %% polynomial costs
0073 ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0074 nb = size(bus, 1);          %% number of buses
0075 nl = size(branch, 1);       %% number of branches
0076 nw = size(N, 1);            %% number of general cost vars, w
0077 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0078 nxyz = getN(om, 'var');     %% total number of control vars of all types
0079 
0080 %% linear constraints & variable bounds
0081 [A, l, u] = linear_constraints(om);
0082 [x0, xmin, xmax] = getv(om);
0083 
0084 %% set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X
0085 %% where X = [x;y;z]. First set up as quadratic function of w,
0086 %% f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We
0087 %% will be building on the (optionally present) user supplied parameters.
0088 
0089 %% piece-wise linear costs
0090 any_pwl = (ny > 0);
0091 if any_pwl
0092     Npwl = sparse(ones(ny,1), vv.i1.y:vv.iN.y, 1, 1, nxyz);     %% sum of y vars
0093     Hpwl = 0;
0094     Cpwl = 1;
0095     fparm_pwl = [1 0 0 1];
0096 else
0097     Npwl = sparse(0, nxyz);
0098     Hpwl = [];
0099     Cpwl = [];
0100     fparm_pwl = [];
0101 end
0102 
0103 %% quadratic costs
0104 npol = length(ipol);
0105 if any(find(gencost(ipol, NCOST) > 3))
0106     error('DC opf cannot handle polynomial costs with higher than quadratic order.');
0107 end
0108 iqdr = find(gencost(ipol, NCOST) == 3);
0109 ilin = find(gencost(ipol, NCOST) == 2);
0110 polycf = zeros(npol, 3);                            %% quadratic coeffs for Pg
0111 if ~isempty(iqdr)
0112   polycf(iqdr, :)   = gencost(ipol(iqdr), COST:COST+2);
0113 end
0114 polycf(ilin, 2:3) = gencost(ipol(ilin), COST:COST+1);
0115 polycf = polycf * diag([ baseMVA^2 baseMVA 1]);     %% convert to p.u.
0116 Npol = sparse(1:npol, vv.i1.Pg-1+ipol, 1, npol, nxyz);         %% Pg vars
0117 Hpol = sparse(1:npol, 1:npol, 2*polycf(:, 1), npol, npol);
0118 Cpol = polycf(:, 2);
0119 fparm_pol = ones(npol,1) * [ 1 0 0 1 ];
0120 
0121 %% combine with user costs
0122 NN = [ Npwl; Npol; N ];
0123 HHw = [ Hpwl, sparse(any_pwl, npol+nw);
0124         sparse(npol, any_pwl), Hpol, sparse(npol, nw);
0125         sparse(nw, any_pwl+npol), H   ];
0126 CCw = [Cpwl; Cpol; Cw];
0127 ffparm = [ fparm_pwl; fparm_pol; fparm ];
0128 
0129 %% transform quadratic coefficients for w into coefficients for X
0130 nnw = any_pwl+npol+nw;
0131 M   = sparse(1:nnw, 1:nnw, ffparm(:, 4), nnw, nnw);
0132 MR  = M * ffparm(:, 2);
0133 HMR = HHw * MR;
0134 MN  = M * NN;
0135 HH = MN' * HHw * MN;
0136 CC = full(MN' * (CCw - HMR));
0137 C0 = 1/2 * MR' * HMR + sum(polycf(:, 3));   %% constant term of cost
0138 
0139 %% default solver
0140 if strcmp(alg, 'DEFAULT')
0141     if have_fcn('cplex')        %% use CPLEX by default, if available
0142         alg = 'CPLEX';
0143     elseif have_fcn('gurobi')   %% if not, then Gurobi, if available
0144         alg = 'GUROBI';
0145     elseif have_fcn('mosek')    %% if not, then MOSEK, if available
0146         alg = 'MOSEK';
0147     elseif have_fcn('bpmpd')    %% if not, then BPMPD_MEX, if available
0148         alg = 'BPMPD';
0149     elseif have_fcn('quadprog') %% if not, then Optimization Tbx, if available
0150         alg = 'OT';
0151     elseif (isempty(HH) || ~any(any(HH))) && have_fcn('glpk') %% if not, and
0152         alg = 'GLPK';           %% prob is LP (not QP), then GLPK, if available
0153     else                        %% otherwise MIPS
0154         alg = 'MIPS';
0155     end
0156 end
0157 
0158 %% set up input for QP solver
0159 opt = struct('alg', alg, 'verbose', mpopt.verbose);
0160 switch alg
0161     case {'MIPS', 'IPOPT'}
0162         %% try to select an interior initial point
0163         if mpopt.opf.init_from_mpc ~= 1
0164             Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0165 
0166             lb = xmin; ub = xmax;
0167             lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0168             ub(xmax ==  Inf) =  1e10;
0169             x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
0170             k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0171             x0(k) = xmax(k) - 1;                    %% set just below upper bound
0172             k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0173             x0(k) = xmin(k) + 1;                    %% set just above lower bound
0174             x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0175             if ny > 0
0176                 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0177                 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0178                 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0179             end
0180         end
0181         switch alg
0182             case 'MIPS'
0183                 %% set up options
0184                 opt.mips_opt = mpopt.mips;
0185                 opt.mips_opt.verbose = mpopt.verbose;
0186                 if opt.mips_opt.feastol == 0
0187                     opt.mips_opt.feastol = mpopt.opf.violation;  %% = MPOPT.opf.violation by default
0188                 end
0189                 if ~isfield(opt.mips_opt, 'cost_mult') || isempty(opt.mips_opt.cost_mult)
0190                     opt.mips_opt.cost_mult = 1;
0191                 end
0192             case 'IPOPT'
0193                 opt.ipopt_opt = ipopt_options([], mpopt);
0194         end
0195     case 'CLP'
0196         opt.clp_opt = clp_options([], mpopt);
0197     case 'CPLEX'
0198         opt.cplex_opt = cplex_options([], mpopt);
0199     case 'GLPK'
0200         opt.glpk_opt = glpk_options([], mpopt);
0201     case 'GUROBI'
0202         opt.grb_opt = gurobi_options([], mpopt);
0203     case 'MOSEK'
0204         opt.mosek_opt = mosek_options([], mpopt);
0205     case 'OT'
0206         if isfield(mpopt, 'linprog') && ~isempty(mpopt.linprog)
0207             opt.linprog_opt = mpopt.linprog;
0208         end
0209         if isfield(mpopt, 'quadprog') && ~isempty(mpopt.quadprog)
0210             opt.quadprog_opt = mpopt.quadprog;
0211         end
0212 end
0213 
0214 %%-----  run opf  -----
0215 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0216 success = (info == 1);
0217 
0218 %%-----  calculate return values  -----
0219 if ~any(isnan(x))
0220     %% update solution data
0221     Va = x(vv.i1.Va:vv.iN.Va);
0222     Pg = x(vv.i1.Pg:vv.iN.Pg);
0223     f = f + C0;
0224 
0225     %% update voltages & generator outputs
0226     bus(:, VM) = ones(nb, 1);
0227     bus(:, VA) = Va * 180/pi;
0228     gen(:, PG) = Pg * baseMVA;
0229 
0230     %% compute branch flows
0231     branch(:, [QF, QT]) = zeros(nl, 2);
0232     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0233     branch(:, PT) = -branch(:, PF);
0234 end
0235 
0236 %% package up results
0237 mu_l = lambda.mu_l;
0238 mu_u = lambda.mu_u;
0239 muLB = lambda.lower;
0240 muUB = lambda.upper;
0241 
0242 %% update Lagrange multipliers
0243 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0244 bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
0245 gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
0246 branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
0247 bus(:, LAM_P)       = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
0248 branch(il, MU_SF)   = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0249 branch(il, MU_ST)   = mu_l(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0250 gen(:, MU_PMIN)     = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0251 gen(:, MU_PMAX)     = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0252 pimul = [
0253   mu_l - mu_u;
0254  -ones(ny>0, 1);    %% dummy entry corresponding to linear cost row in A (in MINOS)
0255   muLB - muUB
0256 ];
0257 
0258 mu = struct( ...
0259   'var', struct('l', muLB, 'u', muUB), ...
0260   'lin', struct('l', mu_l, 'u', mu_u) );
0261 
0262 results = mpc;
0263 [results.bus, results.branch, results.gen, ...
0264     results.om, results.x, results.mu, results.f] = ...
0265         deal(bus, branch, gen, om, x, mu, f);
0266 
0267 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', output);
dcopf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
