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