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 vector.
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,v 1.39 2011/11/11 15:42:46 cvs Exp $
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 verbose = mpopt(31);    %% VERBOSE
0078 alg     = mpopt(26);    %% OPF_ALG_DC
0079 
0080 %% default solver
0081 if alg == 0
0082     if have_fcn('cplex')        %% use CPLEX by default, if available
0083         alg = 500;
0084     elseif have_fcn('mosek')    %% if not, then MOSEK, if available
0085         alg = 600;
0086     elseif have_fcn('gurobi')   %% if not, then Gurobi, if available
0087         alg = 700;
0088     elseif have_fcn('bpmpd')    %% if not, then BPMPD_MEX, if available
0089         alg = 100;
0090     elseif have_fcn('quadprog') %% if not, then Optimization Tbx, if available
0091         alg = 300;
0092     else                        %% otherwise MIPS
0093         alg = 200;
0094     end
0095 end
0096 
0097 %% unpack data
0098 mpc = get_mpc(om);
0099 [baseMVA, bus, gen, branch, gencost] = ...
0100     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0101 cp = get_cost_params(om);
0102 [N, H, Cw] = deal(cp.N, cp.H, cp.Cw);
0103 fparm = [cp.dd cp.rh cp.kk cp.mm];
0104 Bf = userdata(om, 'Bf');
0105 Pfinj = userdata(om, 'Pfinj');
0106 [vv, ll] = get_idx(om);
0107 
0108 %% problem dimensions
0109 ipol = find(gencost(:, MODEL) == POLYNOMIAL); %% polynomial costs
0110 ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0111 nb = size(bus, 1);          %% number of buses
0112 nl = size(branch, 1);       %% number of branches
0113 nw = size(N, 1);            %% number of general cost vars, w
0114 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0115 nxyz = getN(om, 'var');     %% total number of control vars of all types
0116 
0117 %% linear constraints & variable bounds
0118 [A, l, u] = linear_constraints(om);
0119 [x0, xmin, xmax] = getv(om);
0120 
0121 %% set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X
0122 %% where X = [x;y;z]. First set up as quadratic function of w,
0123 %% f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We
0124 %% will be building on the (optionally present) user supplied parameters.
0125 
0126 %% piece-wise linear costs
0127 any_pwl = (ny > 0);
0128 if any_pwl
0129     Npwl = sparse(ones(ny,1), vv.i1.y:vv.iN.y, 1, 1, nxyz);     %% sum of y vars
0130     Hpwl = 0;
0131     Cpwl = 1;
0132     fparm_pwl = [1 0 0 1];
0133 else
0134     Npwl = sparse(0, nxyz);
0135     Hpwl = [];
0136     Cpwl = [];
0137     fparm_pwl = [];
0138 end
0139 
0140 %% quadratic costs
0141 npol = length(ipol);
0142 if any(find(gencost(ipol, NCOST) > 3))
0143     error('DC opf cannot handle polynomial costs with higher than quadratic order.');
0144 end
0145 iqdr = find(gencost(ipol, NCOST) == 3);
0146 ilin = find(gencost(ipol, NCOST) == 2);
0147 polycf = zeros(npol, 3);                            %% quadratic coeffs for Pg
0148 if ~isempty(iqdr)
0149   polycf(iqdr, :)   = gencost(ipol(iqdr), COST:COST+2);
0150 end
0151 polycf(ilin, 2:3) = gencost(ipol(ilin), COST:COST+1);
0152 polycf = polycf * diag([ baseMVA^2 baseMVA 1]);     %% convert to p.u.
0153 Npol = sparse(1:npol, vv.i1.Pg-1+ipol, 1, npol, nxyz);         %% Pg vars
0154 Hpol = sparse(1:npol, 1:npol, 2*polycf(:, 1), npol, npol);
0155 Cpol = polycf(:, 2);
0156 fparm_pol = ones(npol,1) * [ 1 0 0 1 ];
0157 
0158 %% combine with user costs
0159 NN = [ Npwl; Npol; N ];
0160 HHw = [ Hpwl, sparse(any_pwl, npol+nw);
0161         sparse(npol, any_pwl), Hpol, sparse(npol, nw);
0162         sparse(nw, any_pwl+npol), H   ];
0163 CCw = [Cpwl; Cpol; Cw];
0164 ffparm = [ fparm_pwl; fparm_pol; fparm ];
0165 
0166 %% transform quadratic coefficients for w into coefficients for X
0167 nnw = any_pwl+npol+nw;
0168 M   = sparse(1:nnw, 1:nnw, ffparm(:, 4), nnw, nnw);
0169 MR  = M * ffparm(:, 2);
0170 HMR = HHw * MR;
0171 MN  = M * NN;
0172 HH = MN' * HHw * MN;
0173 CC = full(MN' * (CCw - HMR));
0174 C0 = 1/2 * MR' * HMR + sum(polycf(:, 3));   %% constant term of cost
0175 
0176 %% set up input for QP solver
0177 opt = struct('alg', alg, 'verbose', verbose);
0178 switch alg
0179     case {200, 250}
0180         %% try to select an interior initial point
0181         Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0182 
0183         lb = xmin; ub = xmax;
0184         lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0185         ub(xmax ==  Inf) =  1e10;
0186         x0 = (lb + ub) / 2;
0187         x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0188         if ny > 0
0189             ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0190             c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0191             x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0192         end
0193 
0194         %% set up options
0195         feastol = mpopt(81);    %% PDIPM_FEASTOL
0196         gradtol = mpopt(82);    %% PDIPM_GRADTOL
0197         comptol = mpopt(83);    %% PDIPM_COMPTOL
0198         costtol = mpopt(84);    %% PDIPM_COSTTOL
0199         max_it  = mpopt(85);    %% PDIPM_MAX_IT
0200         max_red = mpopt(86);    %% SCPDIPM_RED_IT
0201         if feastol == 0
0202             feastol = mpopt(16);    %% = OPF_VIOLATION by default
0203         end
0204         opt.mips_opt = struct(  'feastol', feastol, ...
0205                                 'gradtol', gradtol, ...
0206                                 'comptol', comptol, ...
0207                                 'costtol', costtol, ...
0208                                 'max_it', max_it, ...
0209                                 'max_red', max_red, ...
0210                                 'cost_mult', 1  );
0211     case 400
0212         opt.ipopt_opt = ipopt_options([], mpopt);
0213     case 500
0214         opt.cplex_opt = cplex_options([], mpopt);
0215     case 600
0216         opt.mosek_opt = mosek_options([], mpopt);
0217     case 700
0218         opt.grb_opt = gurobi_options([], mpopt);
0219 end
0220 
0221 %%-----  run opf  -----
0222 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0223 success = (info == 1);
0224 
0225 %%-----  calculate return values  -----
0226 if ~any(isnan(x))
0227     %% update solution data
0228     Va = x(vv.i1.Va:vv.iN.Va);
0229     Pg = x(vv.i1.Pg:vv.iN.Pg);
0230     f = f + C0;
0231 
0232     %% update voltages & generator outputs
0233     bus(:, VA) = Va * 180/pi;
0234     gen(:, PG) = Pg * baseMVA;
0235 
0236     %% compute branch flows
0237     branch(:, [QF, QT]) = zeros(nl, 2);
0238     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0239     branch(:, PT) = -branch(:, PF);
0240 end
0241 
0242 %% package up results
0243 mu_l = lambda.mu_l;
0244 mu_u = lambda.mu_u;
0245 muLB = lambda.lower;
0246 muUB = lambda.upper;
0247 
0248 %% update Lagrange multipliers
0249 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0250 bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
0251 gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
0252 branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
0253 bus(:, LAM_P)       = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
0254 branch(il, MU_SF)   = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0255 branch(il, MU_ST)   = mu_u(ll.i1.Pt:ll.iN.Pt) / baseMVA;
0256 gen(:, MU_PMIN)     = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0257 gen(:, MU_PMAX)     = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0258 pimul = [
0259   mu_l - mu_u;
0260  -ones(ny>0, 1);    %% dummy entry corresponding to linear cost row in A (in MINOS)
0261   muLB - muUB
0262 ];
0263 
0264 mu = struct( ...
0265   'var', struct('l', muLB, 'u', muUB), ...
0266   'lin', struct('l', mu_l, 'u', mu_u) );
0267 
0268 results = mpc;
0269 [results.bus, results.branch, results.gen, ...
0270     results.om, results.x, results.mu, results.f] = ...
0271         deal(bus, branch, gen, om, x, mu, f);
0272 
0273 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', output);
dcopf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
