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 %   $Id: dcopf_solver.m 2439 2014-12-05 14:33:30Z ray $
0036 %   by Ray Zimmerman, PSERC Cornell
0037 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0038 %   Copyright (c) 2000-2010 by Power System Engineering Research Center (PSERC)
0039 %
0040 %   This file is part of MATPOWER.
0041 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0042 %
0043 %   MATPOWER is free software: you can redistribute it and/or modify
0044 %   it under the terms of the GNU General Public License as published
0045 %   by the Free Software Foundation, either version 3 of the License,
0046 %   or (at your option) any later version.
0047 %
0048 %   MATPOWER is distributed in the hope that it will be useful,
0049 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0050 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0051 %   GNU General Public License for more details.
0052 %
0053 %   You should have received a copy of the GNU General Public License
0054 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0055 %
0056 %   Additional permission under GNU GPL version 3 section 7
0057 %
0058 %   If you modify MATPOWER, or any covered work, to interface with
0059 %   other modules (such as MATLAB code and MEX-files) available in a
0060 %   MATLAB(R) or comparable environment containing parts covered
0061 %   under other licensing terms, the licensors of MATPOWER grant
0062 %   you additional permission to convey the resulting work.
0063 
0064 %%----- initialization -----
0065 %% define named indices into data matrices
0066 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0067     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0068 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0069     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0070     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0071 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0072     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0073     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0074 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0075 
0076 %% options
0077 alg = upper(mpopt.opf.dc.solver);
0078 
0079 %% unpack data
0080 mpc = get_mpc(om);
0081 [baseMVA, bus, gen, branch, gencost] = ...
0082     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0083 cp = get_cost_params(om);
0084 [N, H, Cw] = deal(cp.N, cp.H, cp.Cw);
0085 fparm = [cp.dd cp.rh cp.kk cp.mm];
0086 Bf = userdata(om, 'Bf');
0087 Pfinj = userdata(om, 'Pfinj');
0088 [vv, ll] = get_idx(om);
0089 
0090 %% problem dimensions
0091 ipol = find(gencost(:, MODEL) == POLYNOMIAL); %% polynomial costs
0092 ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0093 nb = size(bus, 1);          %% number of buses
0094 nl = size(branch, 1);       %% number of branches
0095 nw = size(N, 1);            %% number of general cost vars, w
0096 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0097 nxyz = getN(om, 'var');     %% total number of control vars of all types
0098 
0099 %% linear constraints & variable bounds
0100 [A, l, u] = linear_constraints(om);
0101 [x0, xmin, xmax] = getv(om);
0102 
0103 %% set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X
0104 %% where X = [x;y;z]. First set up as quadratic function of w,
0105 %% f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We
0106 %% will be building on the (optionally present) user supplied parameters.
0107 
0108 %% piece-wise linear costs
0109 any_pwl = (ny > 0);
0110 if any_pwl
0111     Npwl = sparse(ones(ny,1), vv.i1.y:vv.iN.y, 1, 1, nxyz);     %% sum of y vars
0112     Hpwl = 0;
0113     Cpwl = 1;
0114     fparm_pwl = [1 0 0 1];
0115 else
0116     Npwl = sparse(0, nxyz);
0117     Hpwl = [];
0118     Cpwl = [];
0119     fparm_pwl = [];
0120 end
0121 
0122 %% quadratic costs
0123 npol = length(ipol);
0124 if any(find(gencost(ipol, NCOST) > 3))
0125     error('DC opf cannot handle polynomial costs with higher than quadratic order.');
0126 end
0127 iqdr = find(gencost(ipol, NCOST) == 3);
0128 ilin = find(gencost(ipol, NCOST) == 2);
0129 polycf = zeros(npol, 3);                            %% quadratic coeffs for Pg
0130 if ~isempty(iqdr)
0131   polycf(iqdr, :)   = gencost(ipol(iqdr), COST:COST+2);
0132 end
0133 polycf(ilin, 2:3) = gencost(ipol(ilin), COST:COST+1);
0134 polycf = polycf * diag([ baseMVA^2 baseMVA 1]);     %% convert to p.u.
0135 Npol = sparse(1:npol, vv.i1.Pg-1+ipol, 1, npol, nxyz);         %% Pg vars
0136 Hpol = sparse(1:npol, 1:npol, 2*polycf(:, 1), npol, npol);
0137 Cpol = polycf(:, 2);
0138 fparm_pol = ones(npol,1) * [ 1 0 0 1 ];
0139 
0140 %% combine with user costs
0141 NN = [ Npwl; Npol; N ];
0142 HHw = [ Hpwl, sparse(any_pwl, npol+nw);
0143         sparse(npol, any_pwl), Hpol, sparse(npol, nw);
0144         sparse(nw, any_pwl+npol), H   ];
0145 CCw = [Cpwl; Cpol; Cw];
0146 ffparm = [ fparm_pwl; fparm_pol; fparm ];
0147 
0148 %% transform quadratic coefficients for w into coefficients for X
0149 nnw = any_pwl+npol+nw;
0150 M   = sparse(1:nnw, 1:nnw, ffparm(:, 4), nnw, nnw);
0151 MR  = M * ffparm(:, 2);
0152 HMR = HHw * MR;
0153 MN  = M * NN;
0154 HH = MN' * HHw * MN;
0155 CC = full(MN' * (CCw - HMR));
0156 C0 = 1/2 * MR' * HMR + sum(polycf(:, 3));   %% constant term of cost
0157 
0158 %% default solver
0159 if strcmp(alg, 'DEFAULT')
0160     if have_fcn('cplex')        %% use CPLEX by default, if available
0161         alg = 'CPLEX';
0162     elseif have_fcn('gurobi')   %% if not, then Gurobi, if available
0163         alg = 'GUROBI';
0164     elseif have_fcn('mosek')    %% if not, then MOSEK, if available
0165         alg = 'MOSEK';
0166     elseif have_fcn('bpmpd')    %% if not, then BPMPD_MEX, if available
0167         alg = 'BPMPD';
0168     elseif have_fcn('quadprog') %% if not, then Optimization Tbx, if available
0169         alg = 'OT';
0170     elseif (isempty(HH) || ~any(any(HH))) && have_fcn('glpk') %% if not, and
0171         alg = 'GLPK';           %% prob is LP (not QP), then GLPK, if available
0172     else                        %% otherwise MIPS
0173         alg = 'MIPS';
0174     end
0175 end
0176 
0177 %% set up input for QP solver
0178 opt = struct('alg', alg, 'verbose', mpopt.verbose);
0179 switch alg
0180     case {'MIPS', 'IPOPT'}
0181         %% try to select an interior initial point
0182         if mpopt.opf.init_from_mpc ~= 1
0183             Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0184 
0185             lb = xmin; ub = xmax;
0186             lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0187             ub(xmax ==  Inf) =  1e10;
0188             x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
0189             k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0190             x0(k) = xmax(k) - 1;                    %% set just below upper bound
0191             k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0192             x0(k) = xmin(k) + 1;                    %% set just above lower bound
0193             x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0194             if ny > 0
0195                 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0196                 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0197                 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0198             end
0199         end
0200         switch alg
0201             case 'MIPS'
0202                 %% set up options
0203                 opt.mips_opt = mpopt.mips;
0204                 opt.mips_opt.verbose = mpopt.verbose;
0205                 if opt.mips_opt.feastol == 0
0206                     opt.mips_opt.feastol = mpopt.opf.violation;  %% = MPOPT.opf.violation by default
0207                 end
0208                 if ~isfield(opt.mips_opt, 'cost_mult') || isempty(opt.mips_opt.cost_mult)
0209                     opt.mips_opt.cost_mult = 1;
0210                 end
0211             case 'IPOPT'
0212                 opt.ipopt_opt = ipopt_options([], mpopt);
0213         end
0214     case 'CPLEX'
0215         opt.cplex_opt = cplex_options([], mpopt);
0216     case 'GLPK'
0217         opt.glpk_opt = glpk_options([], mpopt);
0218     case 'GUROBI'
0219         opt.grb_opt = gurobi_options([], mpopt);
0220     case 'MOSEK'
0221         opt.mosek_opt = mosek_options([], mpopt);
0222     case 'OT'
0223         if isfield(mpopt, 'linprog') && ~isempty(mpopt.linprog)
0224             opt.linprog_opt = mpopt.linprog;
0225         end
0226         if isfield(mpopt, 'quadprog') && ~isempty(mpopt.quadprog)
0227             opt.quadprog_opt = mpopt.quadprog;
0228         end
0229 end
0230 
0231 %%-----  run opf  -----
0232 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0233 success = (info == 1);
0234 
0235 %%-----  calculate return values  -----
0236 if ~any(isnan(x))
0237     %% update solution data
0238     Va = x(vv.i1.Va:vv.iN.Va);
0239     Pg = x(vv.i1.Pg:vv.iN.Pg);
0240     f = f + C0;
0241 
0242     %% update voltages & generator outputs
0243     bus(:, VM) = ones(nb, 1);
0244     bus(:, VA) = Va * 180/pi;
0245     gen(:, PG) = Pg * baseMVA;
0246 
0247     %% compute branch flows
0248     branch(:, [QF, QT]) = zeros(nl, 2);
0249     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0250     branch(:, PT) = -branch(:, PF);
0251 end
0252 
0253 %% package up results
0254 mu_l = lambda.mu_l;
0255 mu_u = lambda.mu_u;
0256 muLB = lambda.lower;
0257 muUB = lambda.upper;
0258 
0259 %% update Lagrange multipliers
0260 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0261 bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
0262 gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
0263 branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
0264 bus(:, LAM_P)       = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
0265 branch(il, MU_SF)   = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0266 branch(il, MU_ST)   = mu_l(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0267 gen(:, MU_PMIN)     = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0268 gen(:, MU_PMAX)     = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0269 pimul = [
0270   mu_l - mu_u;
0271  -ones(ny>0, 1);    %% dummy entry corresponding to linear cost row in A (in MINOS)
0272   muLB - muUB
0273 ];
0274 
0275 mu = struct( ...
0276   'var', struct('l', muLB, 'u', muUB), ...
0277   'lin', struct('l', mu_l, 'u', mu_u) );
0278 
0279 results = mpc;
0280 [results.bus, results.branch, results.gen, ...
0281     results.om, results.x, results.mu, results.f] = ...
0282         deal(bus, branch, gen, om, x, mu, f);
0283 
0284 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', output);
dcopf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
