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.37 2010/12/16 20:59:57 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('mosek')        %% use MOSEK by default if available
0083         alg = 600;
0084     elseif have_fcn('bpmpd')    %% if not, then BPMPD_MEX if available
0085         alg = 100;
0086     elseif have_fcn('cplex')    %% if not, then CPLEX if available
0087         alg = 500;
0088     else                        %% otherwise MIPS
0089         alg = 200;
0090     end
0091 end
0092 
0093 %% unpack data
0094 mpc = get_mpc(om);
0095 [baseMVA, bus, gen, branch, gencost] = ...
0096     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0097 cp = get_cost_params(om);
0098 [N, H, Cw] = deal(cp.N, cp.H, cp.Cw);
0099 fparm = [cp.dd cp.rh cp.kk cp.mm];
0100 Bf = userdata(om, 'Bf');
0101 Pfinj = userdata(om, 'Pfinj');
0102 [vv, ll] = get_idx(om);
0103 
0104 %% problem dimensions
0105 ipol = find(gencost(:, MODEL) == POLYNOMIAL); %% polynomial costs
0106 ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0107 nb = size(bus, 1);          %% number of buses
0108 nl = size(branch, 1);       %% number of branches
0109 nw = size(N, 1);            %% number of general cost vars, w
0110 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0111 nxyz = getN(om, 'var');     %% total number of control vars of all types
0112 
0113 %% linear constraints & variable bounds
0114 [A, l, u] = linear_constraints(om);
0115 [x0, xmin, xmax] = getv(om);
0116 
0117 %% set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X
0118 %% where X = [x;y;z]. First set up as quadratic function of w,
0119 %% f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We
0120 %% will be building on the (optionally present) user supplied parameters.
0121 
0122 %% piece-wise linear costs
0123 any_pwl = (ny > 0);
0124 if any_pwl
0125     Npwl = sparse(ones(ny,1), vv.i1.y:vv.iN.y, 1, 1, nxyz);     %% sum of y vars
0126     Hpwl = 0;
0127     Cpwl = 1;
0128     fparm_pwl = [1 0 0 1];
0129 else
0130     Npwl = sparse(0, nxyz);
0131     Hpwl = [];
0132     Cpwl = [];
0133     fparm_pwl = [];
0134 end
0135 
0136 %% quadratic costs
0137 npol = length(ipol);
0138 if any(find(gencost(ipol, NCOST) > 3))
0139     error('DC opf cannot handle polynomial costs with higher than quadratic order.');
0140 end
0141 iqdr = find(gencost(ipol, NCOST) == 3);
0142 ilin = find(gencost(ipol, NCOST) == 2);
0143 polycf = zeros(npol, 3);                            %% quadratic coeffs for Pg
0144 if ~isempty(iqdr)
0145   polycf(iqdr, :)   = gencost(ipol(iqdr), COST:COST+2);
0146 end
0147 polycf(ilin, 2:3) = gencost(ipol(ilin), COST:COST+1);
0148 polycf = polycf * diag([ baseMVA^2 baseMVA 1]);     %% convert to p.u.
0149 Npol = sparse(1:npol, vv.i1.Pg-1+ipol, 1, npol, nxyz);         %% Pg vars
0150 Hpol = sparse(1:npol, 1:npol, 2*polycf(:, 1), npol, npol);
0151 Cpol = polycf(:, 2);
0152 fparm_pol = ones(npol,1) * [ 1 0 0 1 ];
0153 
0154 %% combine with user costs
0155 NN = [ Npwl; Npol; N ];
0156 HHw = [ Hpwl, sparse(any_pwl, npol+nw);
0157         sparse(npol, any_pwl), Hpol, sparse(npol, nw);
0158         sparse(nw, any_pwl+npol), H   ];
0159 CCw = [Cpwl; Cpol; Cw];
0160 ffparm = [ fparm_pwl; fparm_pol; fparm ];
0161 
0162 %% transform quadratic coefficients for w into coefficients for X
0163 nnw = any_pwl+npol+nw;
0164 M   = sparse(1:nnw, 1:nnw, ffparm(:, 4), nnw, nnw);
0165 MR  = M * ffparm(:, 2);
0166 HMR = HHw * MR;
0167 MN  = M * NN;
0168 HH = MN' * HHw * MN;
0169 CC = full(MN' * (CCw - HMR));
0170 C0 = 1/2 * MR' * HMR + sum(polycf(:, 3));   %% constant term of cost
0171 
0172 %% set up input for QP solver
0173 opt = struct('alg', alg, 'verbose', verbose);
0174 switch alg
0175     case {200, 250}
0176         %% try to select an interior initial point
0177         Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0178 
0179         lb = xmin; ub = xmax;
0180         lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0181         ub(xmax ==  Inf) =  1e10;
0182         x0 = (lb + ub) / 2;
0183         x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0184         if ny > 0
0185             ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0186             c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0187             x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0188         end
0189 
0190         %% set up options
0191         feastol = mpopt(81);    %% PDIPM_FEASTOL
0192         gradtol = mpopt(82);    %% PDIPM_GRADTOL
0193         comptol = mpopt(83);    %% PDIPM_COMPTOL
0194         costtol = mpopt(84);    %% PDIPM_COSTTOL
0195         max_it  = mpopt(85);    %% PDIPM_MAX_IT
0196         max_red = mpopt(86);    %% SCPDIPM_RED_IT
0197         if feastol == 0
0198             feastol = mpopt(16);    %% = OPF_VIOLATION by default
0199         end
0200         opt.mips_opt = struct(  'feastol', feastol, ...
0201                                 'gradtol', gradtol, ...
0202                                 'comptol', comptol, ...
0203                                 'costtol', costtol, ...
0204                                 'max_it', max_it, ...
0205                                 'max_red', max_red, ...
0206                                 'cost_mult', 1  );
0207     case 400
0208         opt.ipopt_opt = ipopt_options([], mpopt);
0209     case 500
0210         opt.cplex_opt = cplex_options([], mpopt);
0211     case 600
0212         opt.mosek_opt = mosek_options([], mpopt);
0213 end
0214 
0215 %%-----  run opf  -----
0216 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0217 success = (info == 1);
0218 
0219 %%-----  calculate return values  -----
0220 if ~any(isnan(x))
0221     %% update solution data
0222     Va = x(vv.i1.Va:vv.iN.Va);
0223     Pg = x(vv.i1.Pg:vv.iN.Pg);
0224     f = f + C0;
0225 
0226     %% update voltages & generator outputs
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_u(ll.i1.Pt:ll.iN.Pt) / 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
