Home > matpower5.1 > dcopf_solver.m

dcopf_solver

PURPOSE ^

DCOPF_SOLVER Solves a DC optimal power flow.

SYNOPSIS ^

function [results, success, raw] = dcopf_solver(om, mpopt)

DESCRIPTION ^

DCOPF_SOLVER  Solves a DC optimal power flow.

   [RESULTS, SUCCESS, RAW] = DCOPF_SOLVER(OM, MPOPT)

   Inputs are an OPF model object and a MATPOWER options struct.

   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.

   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
   branch, gen, gencost fields, along with the following additional
   fields:
       .order      see 'help ext2int' for details of this field
       .x          final value of optimization variables (internal order)
       .f          final objective function value
       .mu         shadow prices on ...
           .var
               .l  lower bounds on variables
               .u  upper bounds on variables
           .lin
               .l  lower bounds on linear constraints
               .u  upper bounds on linear constraints

   SUCCESS     1 if solver converged successfully, 0 otherwise

   RAW         raw output in form returned by MINOS
       .xr     final value of optimization variables
       .pimul  constraint multipliers
       .info   solver specific termination code
       .output solver specific output information

   See also OPF, QPS_MATPOWER.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005