Home > matpower5.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 %   $Id: mipsopf_solver.m 2438 2014-12-05 14:32:42Z ray $
0039 %   by Ray Zimmerman, PSERC Cornell
0040 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %   Copyright (c) 2000-2013 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 opt = mpopt.mips;
0081 opt.verbose = mpopt.verbose;
0082 if opt.feastol == 0
0083     opt.feastol = mpopt.opf.violation;  %% = MPOPT.opf.violation by default
0084 end
0085 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult)
0086     opt.cost_mult = 1e-4;
0087 end
0088 
0089 %% unpack data
0090 mpc = get_mpc(om);
0091 [baseMVA, bus, gen, branch, gencost] = ...
0092     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0093 [vv, ll, nn] = get_idx(om);
0094 
0095 %% problem dimensions
0096 nb = size(bus, 1);          %% number of buses
0097 nl = size(branch, 1);       %% number of branches
0098 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0099 
0100 %% linear constraints
0101 [A, l, u] = linear_constraints(om);
0102 
0103 %% bounds on optimization vars
0104 [x0, xmin, xmax] = getv(om);
0105 
0106 %% build admittance matrices
0107 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0108 
0109 %% try to select an interior initial point
0110 if mpopt.opf.init_from_mpc ~= 1
0111     ll = xmin; uu = xmax;
0112     ll(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0113     uu(xmax ==  Inf) =  1e10;
0114     x0 = (ll + uu) / 2;         %% set x0 mid-way between bounds
0115     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0116     x0(k) = xmax(k) - 1;                    %% set just below upper bound
0117     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0118     x0(k) = xmin(k) + 1;                    %% set just above lower bound
0119     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0120     x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0121     if ny > 0
0122         ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0123     %     PQ = [gen(:, PMAX); gen(:, QMAX)];
0124     %     c = totcost(gencost(ipwl, :), PQ(ipwl));
0125         c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0126         x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0127     %     x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c);
0128     end
0129 end
0130 
0131 %% find branches with flow limits
0132 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0133 nl2 = length(il);           %% number of constrained lines
0134 
0135 %%-----  run opf  -----
0136 f_fcn = @(x)opf_costfcn(x, om);
0137 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0138 hess_fcn = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0139 [x, f, info, Output, Lambda] = ...
0140   mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0141 success = (info > 0);
0142 
0143 %% update solution data
0144 Va = x(vv.i1.Va:vv.iN.Va);
0145 Vm = x(vv.i1.Vm:vv.iN.Vm);
0146 Pg = x(vv.i1.Pg:vv.iN.Pg);
0147 Qg = x(vv.i1.Qg:vv.iN.Qg);
0148 V = Vm .* exp(1j*Va);
0149 
0150 %%-----  calculate return values  -----
0151 %% update voltages & generator outputs
0152 bus(:, VA) = Va * 180/pi;
0153 bus(:, VM) = Vm;
0154 gen(:, PG) = Pg * baseMVA;
0155 gen(:, QG) = Qg * baseMVA;
0156 gen(:, VG) = Vm(gen(:, GEN_BUS));
0157 
0158 %% compute branch flows
0159 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0160 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0161 branch(:, PF) = real(Sf) * baseMVA;
0162 branch(:, QF) = imag(Sf) * baseMVA;
0163 branch(:, PT) = real(St) * baseMVA;
0164 branch(:, QT) = imag(St) * baseMVA;
0165 
0166 %% line constraint is actually on square of limit
0167 %% so we must fix multipliers
0168 muSf = zeros(nl, 1);
0169 muSt = zeros(nl, 1);
0170 if ~isempty(il)
0171     muSf(il) = 2 * Lambda.ineqnonlin(1:nl2)       .* branch(il, RATE_A) / baseMVA;
0172     muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0173 end
0174 
0175 %% update Lagrange multipliers
0176 bus(:, MU_VMAX)  = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0177 bus(:, MU_VMIN)  = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0178 gen(:, MU_PMAX)  = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0179 gen(:, MU_PMIN)  = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0180 gen(:, MU_QMAX)  = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0181 gen(:, MU_QMIN)  = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0182 bus(:, LAM_P)    = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0183 bus(:, LAM_Q)    = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0184 branch(:, MU_SF) = muSf / baseMVA;
0185 branch(:, MU_ST) = muSt / baseMVA;
0186 
0187 %% package up results
0188 nlnN = getN(om, 'nln');
0189 
0190 %% extract multipliers for nonlinear constraints
0191 kl = find(Lambda.eqnonlin < 0);
0192 ku = find(Lambda.eqnonlin > 0);
0193 nl_mu_l = zeros(nlnN, 1);
0194 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0195 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0196 nl_mu_u(ku) =  Lambda.eqnonlin(ku);
0197 
0198 mu = struct( ...
0199   'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0200   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0201   'lin', struct('l', Lambda.mu_l, 'u', Lambda.mu_u) );
0202 
0203 results = mpc;
0204 [results.bus, results.branch, results.gen, ...
0205     results.om, results.x, results.mu, results.f] = ...
0206         deal(bus, branch, gen, om, x, mu, f);
0207 
0208 pimul = [ ...
0209   results.mu.nln.l - results.mu.nln.u;
0210   results.mu.lin.l - results.mu.lin.u;
0211   -ones(ny>0, 1);
0212   results.mu.var.l - results.mu.var.u;
0213 ];
0214 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005