Home > matpower5.0 > 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 %   $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);

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005