dcopf_solver

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.

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-2016, Power Systems Engineering Research Center (PSERC)
0036 %   by Ray Zimmerman, PSERC Cornell
0037 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0038 %
0039 %   This file is part of MATPOWER.
0040 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0041 %   See https://matpower.org for more info.
0042 
0043 %%----- initialization -----
0044 %% define named indices into data matrices
0045 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0046     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0047 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0048     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0049     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0050 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0051     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0052     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0053 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0054 
0055 %% unpack data
0056 mpc = om.get_mpc();
0057 [baseMVA, bus, gen, branch, gencost] = ...
0058     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0059 cp = om.get_cost_params();
0060 Bf = om.get_userdata('Bf');
0061 Pfinj = om.get_userdata('Pfinj');
0062 [vv, ll] = om.get_idx();
0063 
0064 %% problem dimensions
0065 nb = size(bus, 1);          %% number of buses
0066 nl = size(branch, 1);       %% number of branches
0067 ny = om.getN('var', 'y');   %% number of piece-wise linear costs
0068 
0069 %% linear constraints & variable bounds
0070 [A, l, u] = om.params_lin_constraint();
0071 [x0, xmin, xmax] = om.params_var();
0072 
0073 %% get cost parameters
0074 [HH, CC, C0] = om.params_quad_cost();
0075 
0076 %% options
0077 if isempty(HH) || ~any(any(HH))
0078     model = 'LP';
0079 else
0080     model = 'QP';
0081 end
0082 opt = mpopt2qpopt(mpopt, model);
0083 
0084 %% try to select an interior initial point, unless requested not to
0085 if mpopt.opf.start < 2 && ...
0086         (strcmp(opt.alg, 'MIPS') || strcmp(opt.alg, 'IPOPT'))
0087     s = 1;                      %% set init point inside bounds by s
0088     lb = xmin; ub = xmax;
0089     lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0090     ub(xmax ==  Inf) =  1e10;
0091     x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
0092     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0093     x0(k) = xmax(k) - s;                    %% set just below upper bound
0094     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0095     x0(k) = xmin(k) + s;                    %% set just above lower bound
0096     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0097     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0098     if ny > 0
0099         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0100         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0101         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0102     end
0103 end
0104 
0105 %%-----  run opf  -----
0106 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0107 success = (info == 1);
0108 
0109 %%-----  calculate return values  -----
0110 if ~any(isnan(x))
0111     %% update solution data
0112     Va = x(vv.i1.Va:vv.iN.Va);
0113     Pg = x(vv.i1.Pg:vv.iN.Pg);
0114     f = f + C0;
0115 
0116     %% update voltages & generator outputs
0117     bus(:, VM) = ones(nb, 1);
0118     bus(:, VA) = Va * 180/pi;
0119     gen(:, PG) = Pg * baseMVA;
0120 
0121     %% compute branch flows
0122     branch(:, [QF, QT]) = zeros(nl, 2);
0123     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0124     branch(:, PT) = -branch(:, PF);
0125 end
0126 
0127 %% package up results
0128 mu_l = lambda.mu_l;
0129 mu_u = lambda.mu_u;
0130 muLB = lambda.lower;
0131 muUB = lambda.upper;
0132 
0133 %% update Lagrange multipliers
0134 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0135 bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
0136 gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
0137 branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
0138 bus(:, LAM_P)       = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
0139 branch(il, MU_SF)   = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0140 branch(il, MU_ST)   = mu_l(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0141 gen(:, MU_PMIN)     = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0142 gen(:, MU_PMAX)     = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0143 pimul = [
0144   mu_l - mu_u;
0145  -ones(ny>0, 1);    %% dummy entry corresponding to linear cost row in A (in MINOS)
0146   muLB - muUB
0147 ];
0148 
0149 mu = struct( ...
0150   'var', struct('l', muLB, 'u', muUB), ...
0151   'lin', struct('l', mu_l, 'u', mu_u) );
0152 
0153 results = mpc;
0154 [results.bus, results.branch, results.gen, ...
0155     results.om, results.x, results.mu, results.f] = ...
0156         deal(bus, branch, gen, om, x, mu, f);
0157 
0158 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', output);