mipsopf_solver

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 vector.

   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.

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 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 %           .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 %   $Id: mipsopf_solver.m,v 1.13 2010/06/09 14:56:58 ray Exp $
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %   Copyright (c) 2000-2010 by Power System Engineering Research Center (PSERC)
0042 %
0043 %   This file is part of MATPOWER.
0044 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0045 %
0046 %   MATPOWER is free software: you can redistribute it and/or modify
0047 %   it under the terms of the GNU General Public License as published
0048 %   by the Free Software Foundation, either version 3 of the License,
0049 %   or (at your option) any later version.
0050 %
0051 %   MATPOWER is distributed in the hope that it will be useful,
0052 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0053 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0054 %   GNU General Public License for more details.
0055 %
0056 %   You should have received a copy of the GNU General Public License
0057 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0058 %
0059 %   Additional permission under GNU GPL version 3 section 7
0060 %
0061 %   If you modify MATPOWER, or any covered work, to interface with
0062 %   other modules (such as MATLAB code and MEX-files) available in a
0063 %   MATLAB(R) or comparable environment containing parts covered
0064 %   under other licensing terms, the licensors of MATPOWER grant
0065 %   you additional permission to convey the resulting work.
0066 
0067 %%----- initialization -----
0068 %% define named indices into data matrices
0069 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0070     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0071 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0072     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0073     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0074 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0075     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0076     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0077 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0078 
0079 %% options
0080 verbose = mpopt(31);    %% VERBOSE
0081 feastol = mpopt(81);    %% PDIPM_FEASTOL
0082 gradtol = mpopt(82);    %% PDIPM_GRADTOL
0083 comptol = mpopt(83);    %% PDIPM_COMPTOL
0084 costtol = mpopt(84);    %% PDIPM_COSTTOL
0085 max_it  = mpopt(85);    %% PDIPM_MAX_IT
0086 max_red = mpopt(86);    %% SCPDIPM_RED_IT
0087 step_control = (mpopt(11) == 565);  %% OPF_ALG == 565, MIPS-sc
0088 if feastol == 0
0089     feastol = mpopt(16);    %% = OPF_VIOLATION by default
0090 end
0091 opt = struct(   'feastol', feastol, ...
0092                 'gradtol', gradtol, ...
0093                 'comptol', comptol, ...
0094                 'costtol', costtol, ...
0095                 'max_it', max_it, ...
0096                 'max_red', max_red, ...
0097                 'step_control', step_control, ...
0098                 'cost_mult', 1e-4, ...
0099                 'verbose', verbose  );
0100 
0101 %% unpack data
0102 mpc = get_mpc(om);
0103 [baseMVA, bus, gen, branch, gencost] = ...
0104     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0105 [vv, ll, nn] = get_idx(om);
0106 
0107 %% problem dimensions
0108 nb = size(bus, 1);          %% number of buses
0109 nl = size(branch, 1);       %% number of branches
0110 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0111 
0112 %% linear constraints
0113 [A, l, u] = linear_constraints(om);
0114 
0115 %% bounds on optimization vars
0116 [x0, xmin, xmax] = getv(om);
0117 
0118 %% build admittance matrices
0119 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0120 
0121 %% try to select an interior initial point
0122 ll = xmin; uu = xmax;
0123 ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0124 uu(xmax ==  Inf) =  1e10;
0125 x0 = (ll + uu) / 2;
0126 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0127 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0128 if ny > 0
0129     ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0130 %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0131 %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0132     c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0133     x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0134 %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0135 end
0136 
0137 %% find branches with flow limits
0138 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0139 nl2 = length(il);           %% number of constrained lines
0140 
0141 %%-----  run opf  -----
0142 f_fcn = @(x)opf_costfcn(x, om);
0143 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0144 hess_fcn = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0145 [x, f, info, Output, Lambda] = ...
0146   mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0147 success = (info > 0);
0148 
0149 %% update solution data
0150 Va = x(vv.i1.Va:vv.iN.Va);
0151 Vm = x(vv.i1.Vm:vv.iN.Vm);
0152 Pg = x(vv.i1.Pg:vv.iN.Pg);
0153 Qg = x(vv.i1.Qg:vv.iN.Qg);
0154 V = Vm .* exp(1j*Va);
0155 
0156 %%-----  calculate return values  -----
0157 %% update voltages & generator outputs
0158 bus(:, VA) = Va * 180/pi;
0159 bus(:, VM) = Vm;
0160 gen(:, PG) = Pg * baseMVA;
0161 gen(:, QG) = Qg * baseMVA;
0162 gen(:, VG) = Vm(gen(:, GEN_BUS));
0163 
0164 %% compute branch flows
0165 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0166 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0167 branch(:, PF) = real(Sf) * baseMVA;
0168 branch(:, QF) = imag(Sf) * baseMVA;
0169 branch(:, PT) = real(St) * baseMVA;
0170 branch(:, QT) = imag(St) * baseMVA;
0171 
0172 %% line constraint is actually on square of limit
0173 %% so we must fix multipliers
0174 muSf = zeros(nl, 1);
0175 muSt = zeros(nl, 1);
0176 if ~isempty(il)
0177     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0178     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0179 end
0180 
0181 %% update Lagrange multipliers
0182 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0183 bus(:, MU_VMIN)  = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0184 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0185 gen(:, MU_PMIN)  = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0186 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0187 gen(:, MU_QMIN)  = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0188 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0189 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0190 branch(:, MU_SF) = muSf / baseMVA;
0191 branch(:, MU_ST) = muSt / baseMVA;
0192 
0193 %% package up results
0194 nlnN = getN(om, 'nln');
0195 
0196 %% extract multipliers for nonlinear constraints
0197 kl = find(Lambda.eqnonlin < 0);
0198 ku = find(Lambda.eqnonlin > 0);
0199 nl_mu_l = zeros(nlnN, 1);
0200 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0201 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0202 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0203 
0204 mu = struct( ...
0205   'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0206   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0207   'lin', struct('l', Lambda.mu_l, 'u', Lambda.mu_u) );
0208 
0209 results = mpc;
0210 [results.bus, results.branch, results.gen, ...
0211     results.om, results.x, results.mu, results.f] = ...
0212         deal(bus, branch, gen, om, x, mu, f);
0213 
0214 pimul = [ ...
0215   results.mu.nln.l - results.mu.nln.u;
0216   results.mu.lin.l - results.mu.lin.u;
0217   -ones(ny>0, 1);
0218   results.mu.var.l - results.mu.var.u;
0219 ];
0220 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);