Home > matpower6.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 %   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 http://www.pserc.cornell.edu/matpower/ 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 = get_mpc(om);
0057 [baseMVA, bus, gen, branch, gencost] = ...
0058     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0059 cp = get_cost_params(om);
0060 [N, H, Cw] = deal(cp.N, cp.H, cp.Cw);
0061 fparm = [cp.dd cp.rh cp.kk cp.mm];
0062 Bf = userdata(om, 'Bf');
0063 Pfinj = userdata(om, 'Pfinj');
0064 [vv, ll] = get_idx(om);
0065 
0066 %% problem dimensions
0067 ipol = find(gencost(:, MODEL) == POLYNOMIAL); %% polynomial costs
0068 ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0069 nb = size(bus, 1);          %% number of buses
0070 nl = size(branch, 1);       %% number of branches
0071 nw = size(N, 1);            %% number of general cost vars, w
0072 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0073 nxyz = getN(om, 'var');     %% total number of control vars of all types
0074 
0075 %% linear constraints & variable bounds
0076 [A, l, u] = linear_constraints(om);
0077 [x0, xmin, xmax] = getv(om);
0078 
0079 %% set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X
0080 %% where X = [x;y;z]. First set up as quadratic function of w,
0081 %% f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We
0082 %% will be building on the (optionally present) user supplied parameters.
0083 
0084 %% piece-wise linear costs
0085 any_pwl = (ny > 0);
0086 if any_pwl
0087     Npwl = sparse(ones(ny,1), vv.i1.y:vv.iN.y, 1, 1, nxyz);     %% sum of y vars
0088     Hpwl = 0;
0089     Cpwl = 1;
0090     fparm_pwl = [1 0 0 1];
0091 else
0092     Npwl = sparse(0, nxyz);
0093     Hpwl = [];
0094     Cpwl = [];
0095     fparm_pwl = [];
0096 end
0097 
0098 %% quadratic costs
0099 npol = length(ipol);
0100 if any(find(gencost(ipol, NCOST) > 3))
0101     error('DC opf cannot handle polynomial costs with higher than quadratic order.');
0102 end
0103 iqdr = find(gencost(ipol, NCOST) == 3);
0104 ilin = find(gencost(ipol, NCOST) == 2);
0105 polycf = zeros(npol, 3);                            %% quadratic coeffs for Pg
0106 if ~isempty(iqdr)
0107   polycf(iqdr, :)   = gencost(ipol(iqdr), COST:COST+2);
0108 end
0109 polycf(ilin, 2:3) = gencost(ipol(ilin), COST:COST+1);
0110 polycf = polycf * diag([ baseMVA^2 baseMVA 1]);     %% convert to p.u.
0111 Npol = sparse(1:npol, vv.i1.Pg-1+ipol, 1, npol, nxyz);         %% Pg vars
0112 Hpol = sparse(1:npol, 1:npol, 2*polycf(:, 1), npol, npol);
0113 Cpol = polycf(:, 2);
0114 fparm_pol = ones(npol,1) * [ 1 0 0 1 ];
0115 
0116 %% combine with user costs
0117 NN = [ Npwl; Npol; N ];
0118 HHw = [ Hpwl, sparse(any_pwl, npol+nw);
0119         sparse(npol, any_pwl), Hpol, sparse(npol, nw);
0120         sparse(nw, any_pwl+npol), H   ];
0121 CCw = [Cpwl; Cpol; Cw];
0122 ffparm = [ fparm_pwl; fparm_pol; fparm ];
0123 
0124 %% transform quadratic coefficients for w into coefficients for X
0125 nnw = any_pwl+npol+nw;
0126 M   = sparse(1:nnw, 1:nnw, ffparm(:, 4), nnw, nnw);
0127 MR  = M * ffparm(:, 2);
0128 HMR = HHw * MR;
0129 MN  = M * NN;
0130 HH = MN' * HHw * MN;
0131 CC = full(MN' * (CCw - HMR));
0132 C0 = 1/2 * MR' * HMR + sum(polycf(:, 3));   %% constant term of cost
0133 
0134 %% options
0135 if isempty(HH) || ~any(any(HH))
0136     model = 'LP';
0137 else
0138     model = 'QP';
0139 end
0140 opt = mpopt2qpopt(mpopt, model);
0141 
0142 %% try to select an interior initial point, unless requested not to
0143 if mpopt.opf.init_from_mpc ~= 1 && ...
0144         (strcmp(opt.alg, 'MIPS') || strcmp(opt.alg, 'IPOPT'))
0145     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0146 
0147     lb = xmin; ub = xmax;
0148     lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0149     ub(xmax ==  Inf) =  1e10;
0150     x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
0151     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0152     x0(k) = xmax(k) - 1;                    %% set just below upper bound
0153     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0154     x0(k) = xmin(k) + 1;                    %% set just above lower bound
0155     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0156     if ny > 0
0157         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0158         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0159         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0160     end
0161 end
0162 
0163 %%-----  run opf  -----
0164 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0165 success = (info == 1);
0166 
0167 %%-----  calculate return values  -----
0168 if ~any(isnan(x))
0169     %% update solution data
0170     Va = x(vv.i1.Va:vv.iN.Va);
0171     Pg = x(vv.i1.Pg:vv.iN.Pg);
0172     f = f + C0;
0173 
0174     %% update voltages & generator outputs
0175     bus(:, VM) = ones(nb, 1);
0176     bus(:, VA) = Va * 180/pi;
0177     gen(:, PG) = Pg * baseMVA;
0178 
0179     %% compute branch flows
0180     branch(:, [QF, QT]) = zeros(nl, 2);
0181     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0182     branch(:, PT) = -branch(:, PF);
0183 end
0184 
0185 %% package up results
0186 mu_l = lambda.mu_l;
0187 mu_u = lambda.mu_u;
0188 muLB = lambda.lower;
0189 muUB = lambda.upper;
0190 
0191 %% update Lagrange multipliers
0192 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0193 bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
0194 gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
0195 branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
0196 bus(:, LAM_P)       = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
0197 branch(il, MU_SF)   = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0198 branch(il, MU_ST)   = mu_l(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0199 gen(:, MU_PMIN)     = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0200 gen(:, MU_PMAX)     = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0201 pimul = [
0202   mu_l - mu_u;
0203  -ones(ny>0, 1);    %% dummy entry corresponding to linear cost row in A (in MINOS)
0204   muLB - muUB
0205 ];
0206 
0207 mu = struct( ...
0208   'var', struct('l', muLB, 'u', muUB), ...
0209   'lin', struct('l', mu_l, 'u', mu_u) );
0210 
0211 results = mpc;
0212 [results.bus, results.branch, results.gen, ...
0213     results.om, results.x, results.mu, results.f] = ...
0214         deal(bus, branch, gen, om, x, mu, f);
0215 
0216 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', output);
Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005