Home > matpower6.0 > mipsopf_solver.m

mipsopf_solver

PURPOSE ^

MIPSOPF_SOLVER Solves AC optimal power flow using MIPS.

SYNOPSIS ^

function [results, success, raw] = mipsopf_solver(om, mpopt)

DESCRIPTION ^

MIPSOPF_SOLVER  Solves AC optimal power flow using MIPS.

   [RESULTS, SUCCESS, RAW] = MIPSOPF_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
           .nln
               .l  lower bounds on nonlinear constraints
               .u  upper bounds on nonlinear constraints
           .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, MIPS.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [results, success, raw] = mipsopf_solver(om, mpopt)
0002 %MIPSOPF_SOLVER  Solves AC optimal power flow using MIPS.
0003 %
0004 %   [RESULTS, SUCCESS, RAW] = MIPSOPF_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 %           .nln
0021 %               .l  lower bounds on nonlinear constraints
0022 %               .u  upper bounds on nonlinear constraints
0023 %           .lin
0024 %               .l  lower bounds on linear constraints
0025 %               .u  upper bounds on linear constraints
0026 %
0027 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0028 %
0029 %   RAW         raw output in form returned by MINOS
0030 %       .xr     final value of optimization variables
0031 %       .pimul  constraint multipliers
0032 %       .info   solver specific termination code
0033 %       .output solver specific output information
0034 %
0035 %   See also OPF, MIPS.
0036 
0037 %   MATPOWER
0038 %   Copyright (c) 2000-2016, Power Systems Engineering Research Center (PSERC)
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0041 %
0042 %   This file is part of MATPOWER.
0043 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0044 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0045 
0046 %%----- initialization -----
0047 %% define named indices into data matrices
0048 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0049     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0050 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0051     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0052     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0053 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0054     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0055     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0056 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0057 
0058 %% options
0059 opt = mpopt.mips;
0060 opt.verbose = mpopt.verbose;
0061 if opt.feastol == 0
0062     opt.feastol = mpopt.opf.violation;  %% = MPOPT.opf.violation by default
0063 end
0064 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult)
0065     opt.cost_mult = 1e-4;
0066 end
0067 
0068 %% unpack data
0069 mpc = get_mpc(om);
0070 [baseMVA, bus, gen, branch, gencost] = ...
0071     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0072 [vv, ll, nn] = get_idx(om);
0073 
0074 %% problem dimensions
0075 nb = size(bus, 1);          %% number of buses
0076 nl = size(branch, 1);       %% number of branches
0077 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0078 
0079 %% linear constraints
0080 [A, l, u] = linear_constraints(om);
0081 
0082 %% bounds on optimization vars
0083 [x0, xmin, xmax] = getv(om);
0084 
0085 %% build admittance matrices
0086 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0087 
0088 %% try to select an interior initial point
0089 if mpopt.opf.init_from_mpc ~= 1
0090     ll = xmin; uu = xmax;
0091     ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0092     uu(xmax ==  Inf) =  1e10;
0093     x0 = (ll + uu) / 2;         %% set x0 mid-way between bounds
0094     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0095     x0(k) = xmax(k) - 1;                    %% set just below upper bound
0096     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0097     x0(k) = xmin(k) + 1;                    %% set just above lower bound
0098     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0099     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0100     if ny > 0
0101         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0102     %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0103     %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0104         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0105         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0106     %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0107     end
0108 end
0109 
0110 %% find branches with flow limits
0111 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0112 nl2 = length(il);           %% number of constrained lines
0113 
0114 %%-----  run opf  -----
0115 f_fcn = @(x)opf_costfcn(x, om);
0116 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0117 hess_fcn = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0118 [x, f, info, Output, Lambda] = ...
0119   mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0120 success = (info > 0);
0121 
0122 %% update solution data
0123 Va = x(vv.i1.Va:vv.iN.Va);
0124 Vm = x(vv.i1.Vm:vv.iN.Vm);
0125 Pg = x(vv.i1.Pg:vv.iN.Pg);
0126 Qg = x(vv.i1.Qg:vv.iN.Qg);
0127 V = Vm .* exp(1j*Va);
0128 
0129 %%-----  calculate return values  -----
0130 %% update voltages & generator outputs
0131 bus(:, VA) = Va * 180/pi;
0132 bus(:, VM) = Vm;
0133 gen(:, PG) = Pg * baseMVA;
0134 gen(:, QG) = Qg * baseMVA;
0135 gen(:, VG) = Vm(gen(:, GEN_BUS));
0136 
0137 %% compute branch flows
0138 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0139 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0140 branch(:, PF) = real(Sf) * baseMVA;
0141 branch(:, QF) = imag(Sf) * baseMVA;
0142 branch(:, PT) = real(St) * baseMVA;
0143 branch(:, QT) = imag(St) * baseMVA;
0144 
0145 %% line constraint is actually on square of limit
0146 %% so we must fix multipliers
0147 muSf = zeros(nl, 1);
0148 muSt = zeros(nl, 1);
0149 if ~isempty(il)
0150     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0151     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0152 end
0153 
0154 %% update Lagrange multipliers
0155 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0156 bus(:, MU_VMIN)  = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0157 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0158 gen(:, MU_PMIN)  = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0159 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0160 gen(:, MU_QMIN)  = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0161 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0162 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0163 branch(:, MU_SF) = muSf / baseMVA;
0164 branch(:, MU_ST) = muSt / baseMVA;
0165 
0166 %% package up results
0167 nlnN = getN(om, 'nln');
0168 
0169 %% extract multipliers for nonlinear constraints
0170 kl = find(Lambda.eqnonlin < 0);
0171 ku = find(Lambda.eqnonlin > 0);
0172 nl_mu_l = zeros(nlnN, 1);
0173 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0174 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0175 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0176 
0177 mu = struct( ...
0178   'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0179   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0180   'lin', struct('l', Lambda.mu_l, 'u', Lambda.mu_u) );
0181 
0182 results = mpc;
0183 [results.bus, results.branch, results.gen, ...
0184     results.om, results.x, results.mu, results.f] = ...
0185         deal(bus, branch, gen, om, x, mu, f);
0186 
0187 pimul = [ ...
0188   results.mu.nln.l - results.mu.nln.u;
0189   results.mu.lin.l - results.mu.lin.u;
0190   -ones(ny>0, 1);
0191   results.mu.var.l - results.mu.var.u;
0192 ];
0193 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);
Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005