lpopf_solver

------------------------------  deprecated  ------------------------------
   OPF solvers based on LPCONSTR to be removed in a future version.
--------------------------------------------------------------------------
LPOPF_SOLVER  Solves AC optimal power flow using succesive LPs.

   [RESULTS, SUCCESS, RAW] = LPOPF_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

   See also OPF, LPCONSTR, LPEQSLVR, FUN_COPF, GRAD_COPF.

0001 function [results, success, raw] = lpopf_solver(om, mpopt)
0002 %------------------------------  deprecated  ------------------------------
0003 %   OPF solvers based on LPCONSTR to be removed in a future version.
0004 %--------------------------------------------------------------------------
0005 %LPOPF_SOLVER  Solves AC optimal power flow using succesive LPs.
0006 %
0007 %   [RESULTS, SUCCESS, RAW] = LPOPF_SOLVER(OM, MPOPT)
0008 %
0009 %   Inputs are an OPF model object and a MATPOWER options vector.
0010 %
0011 %   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.
0012 %
0013 %   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
0014 %   branch, gen, gencost fields, along with the following additional
0015 %   fields:
0016 %       .order      see 'help ext2int' for details of this field
0017 %       .x          final value of optimization variables (internal order)
0018 %       .f          final objective function value
0019 %       .mu         shadow prices on ...
0020 %           .var
0021 %               .l  lower bounds on variables
0022 %               .u  upper bounds on variables
0023 %           .nln
0024 %               .l  lower bounds on nonlinear constraints
0025 %               .u  upper bounds on nonlinear constraints
0026 %           .lin
0027 %               .l  lower bounds on linear constraints
0028 %               .u  upper bounds on linear constraints
0029 %
0030 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0031 %
0032 %   RAW         raw output in form returned by MINOS
0033 %       .xr     final value of optimization variables
0034 %       .pimul  constraint multipliers
0035 %       .info   solver specific termination code
0036 %
0037 %   See also OPF, LPCONSTR, LPEQSLVR, FUN_COPF, GRAD_COPF.
0038 
0039 %   MATPOWER
0040 %   $Id: lpopf_solver.m,v 1.21 2010/06/09 14:56:58 ray Exp $
0041 %   by Ray Zimmerman, PSERC Cornell
0042 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0043 %   Copyright (c) 2000-2010 by Power System Engineering Research Center (PSERC)
0044 %
0045 %   This file is part of MATPOWER.
0046 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0047 %
0048 %   MATPOWER is free software: you can redistribute it and/or modify
0049 %   it under the terms of the GNU General Public License as published
0050 %   by the Free Software Foundation, either version 3 of the License,
0051 %   or (at your option) any later version.
0052 %
0053 %   MATPOWER is distributed in the hope that it will be useful,
0054 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0055 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0056 %   GNU General Public License for more details.
0057 %
0058 %   You should have received a copy of the GNU General Public License
0059 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0060 %
0061 %   Additional permission under GNU GPL version 3 section 7
0062 %
0063 %   If you modify MATPOWER, or any covered work, to interface with
0064 %   other modules (such as MATLAB code and MEX-files) available in a
0065 %   MATLAB(R) or comparable environment containing parts covered
0066 %   under other licensing terms, the licensors of MATPOWER grant
0067 %   you additional permission to convey the resulting work.
0068 
0069 %%----- initialization -----
0070 %% define named indices into data matrices
0071 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0072     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0073 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0074     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0075     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0076 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0077     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0078     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0079 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0080 
0081 %% unpack data
0082 mpc = get_mpc(om);
0083 [baseMVA, bus, gen, branch, gencost] = ...
0084     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0085 vv = get_idx(om);
0086 
0087 %% problem dimensions
0088 nb = size(bus, 1);          %% number of buses
0089 nl = size(branch, 1);       %% number of branches
0090 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0091 
0092 %% bounds on optimization vars
0093 [x0, LB, UB] = getv(om);
0094 
0095 %% linear constraints
0096 %% WORKAROUND: Add bounds on all vars to A, l, u
0097 nxyz = length(x0);
0098 om2 = om;
0099 om2 = add_constraints(om2, 'varlims', speye(nxyz, nxyz), LB, UB);
0100 [vv, ll, nn] = get_idx(om2);
0101 [A, l, u] = linear_constraints(om2);
0102 
0103 %% so, can we do anything good about lambda initialization?
0104 if all(bus(:, LAM_P) == 0)
0105   bus(:, LAM_P) = (10)*ones(nb, 1);
0106 end
0107 
0108 %% split l <= A*x <= u into less than, equal to, greater than, and
0109 %% doubly-bounded sets
0110 ieq = find( abs(u-l) <= eps );          %% equality
0111 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0112 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0113 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0114 Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0115 bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
0116 Afeq = A(ieq, :);
0117 bfeq = u(ieq);
0118 
0119 %% build admittance matrices
0120 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0121 
0122 %% find branches with flow limits
0123 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0124 nl2 = length(il);           %% number of constrained lines
0125 
0126 %% set mpopt defaults
0127 mpopt(15)  = 2 * nb + length(bfeq); %% set number of equality constraints
0128 if mpopt(19) == 0           %% CONSTR_MAX_IT
0129   mpopt(19) = 150 + 2*nb;
0130 end
0131 
0132 %% run load flow to get starting point
0133 [x, success_lf] = LPeqslvr(x0, om2, Ybus, Yf, Yt, Afeq, bfeq, Af, bf, mpopt);
0134 if success_lf ~= 1
0135   error('Sorry, cannot find a starting point using power flow, please check data!');
0136 end
0137 
0138 %% set step size
0139 cstep = 0;
0140 if ny > 0
0141   PgQg = [gen(:, PG); gen(:, QG)];
0142   ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0143   Cp = totcost(gencost(ipwl, :), PgQg(ipwl));
0144   cstep = max(abs(Cp));
0145   if cstep < 1.0e6, cstep = 1.0e6; end
0146 end
0147 step0 = ones(size(x0));
0148 step0(vv.i1.Va:vv.iN.Va) = 2;
0149 step0(vv.i1.Vm:vv.iN.Vm) = 1;
0150 step0(vv.i1.Pg:vv.iN.Pg) = 0.6;
0151 step0(vv.i1.Qg:vv.iN.Qg) = 0.3;
0152 if ny > 0
0153   step0(vv.i1.y:vv.iN.y)   = cstep;
0154 end
0155 
0156 %% run optimization
0157 [x, lambda, success] = LPconstr('fun_copf', x0, mpopt, step0, [], [], 'grad_copf', ...
0158                        'LPeqslvr', om2, Ybus, Yf(il,:), Yt(il,:), Afeq, bfeq, Af, bf, mpopt, il);
0159 info = success;
0160 
0161 %% get final objective function value & constraint values
0162 [f, g] = feval('fun_copf', x, om2, Ybus, Yf(il,:), Yt(il,:), Afeq, bfeq, Af, bf, mpopt, il);
0163 
0164 %% update solution data
0165 Va = x(vv.i1.Va:vv.iN.Va);
0166 Vm = x(vv.i1.Vm:vv.iN.Vm);
0167 Pg = x(vv.i1.Pg:vv.iN.Pg);
0168 Qg = x(vv.i1.Qg:vv.iN.Qg);
0169 V = Vm .* exp(1j*Va);
0170 
0171 %%-----  calculate return values  -----
0172 %% update voltages & generator outputs
0173 bus(:, VA) = Va * 180/pi;
0174 bus(:, VM) = Vm;
0175 gen(:, PG) = Pg * baseMVA;
0176 gen(:, QG) = Qg * baseMVA;
0177 gen(:, VG) = Vm(gen(:, GEN_BUS));
0178 
0179 %% compute branch flows
0180 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0181 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0182 branch(:, PF) = real(Sf) * baseMVA;
0183 branch(:, QF) = imag(Sf) * baseMVA;
0184 branch(:, PT) = real(St) * baseMVA;
0185 branch(:, QT) = imag(St) * baseMVA;
0186 
0187 %% package up results
0188 nA = length(u);
0189 neq = length(ieq);
0190 nlt = length(ilt);
0191 ngt = length(igt);
0192 nbx = length(ibx);
0193 
0194 %% extract multipliers (lambda is ordered as Pmis, Qmis, Afeq, Sf, St, Af)
0195 %% nonlinear constraints
0196 inln = [(1:2*nb), (1:2*nl2) + 2*nb+neq];
0197 kl = find(lambda(inln) < 0);
0198 ku = find(lambda(inln) > 0);
0199 nl_mu_l = zeros(2*(nb+nl2), 1);
0200 nl_mu_u = zeros(2*(nb+nl2), 1);
0201 nl_mu_l(kl) = -lambda(inln(kl));
0202 nl_mu_u(ku) =  lambda(inln(ku));
0203 
0204 %% linear constraints
0205 ilin = [(1:neq)+2*nb, (1:(nlt+ngt+2*nbx)) + 2*nb+neq+2*nl2];
0206 kl = find(lambda(ilin(1:neq)) < 0);
0207 ku = find(lambda(ilin(1:neq)) > 0);
0208 
0209 mu_l = zeros(nA, 1);
0210 mu_l(ieq) = -lambda(ilin(1:neq));
0211 mu_l(ieq(ku)) = 0;
0212 mu_l(igt) = lambda(ilin(neq+nlt+(1:ngt)));
0213 mu_l(ibx) = lambda(ilin(neq+nlt+ngt+nbx+(1:nbx)));
0214 
0215 mu_u = zeros(nA, 1);
0216 mu_u(ieq) = lambda(ilin(1:neq));
0217 mu_u(ieq(kl)) = 0;
0218 mu_u(ilt) = lambda(ilin(neq+(1:nlt)));
0219 mu_u(ibx) = lambda(ilin(neq+nlt+ngt+(1:nbx)));
0220 
0221 %% variable bounds
0222 muLB = mu_l(ll.i1.varlims:ll.iN.varlims);
0223 muUB = mu_u(ll.i1.varlims:ll.iN.varlims);
0224 mu_l(ll.i1.varlims:ll.iN.varlims) = [];
0225 mu_u(ll.i1.varlims:ll.iN.varlims) = [];
0226 
0227 %% line constraint is actually on square of limit
0228 %% so we must fix multipliers
0229 muSf = zeros(nl, 1);
0230 muSt = zeros(nl, 1);
0231 muSf(il) = 2 * nl_mu_u((1:nl2)+2*nb    ) .* branch(il, RATE_A) / baseMVA;
0232 muSt(il) = 2 * nl_mu_u((1:nl2)+2*nb+nl2) .* branch(il, RATE_A) / baseMVA;
0233 
0234 %% resize mu for nonlinear constraints
0235 nl_mu_l = [nl_mu_l(1:2*nb); zeros(2*nl, 1)];
0236 nl_mu_u = [nl_mu_u(1:2*nb); muSf; muSt];
0237 
0238 %% update Lagrange multipliers
0239 bus(:, MU_VMAX)  = muUB(vv.i1.Vm:vv.iN.Vm);
0240 bus(:, MU_VMIN)  = muLB(vv.i1.Vm:vv.iN.Vm);
0241 gen(:, MU_PMAX)  = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0242 gen(:, MU_PMIN)  = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0243 gen(:, MU_QMAX)  = muUB(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0244 gen(:, MU_QMIN)  = muLB(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0245 bus(:, LAM_P)    = (nl_mu_u(nn.i1.Pmis:nn.iN.Pmis) - nl_mu_l(nn.i1.Pmis:nn.iN.Pmis)) / baseMVA;
0246 bus(:, LAM_Q)    = (nl_mu_u(nn.i1.Qmis:nn.iN.Qmis) - nl_mu_l(nn.i1.Qmis:nn.iN.Qmis)) / baseMVA;
0247 branch(:, MU_SF) = muSf / baseMVA;
0248 branch(:, MU_ST) = muSt / baseMVA;
0249 
0250 mu = struct( ...
0251   'var', struct('l', muLB, 'u', muUB), ...
0252   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0253   'lin', struct('l', mu_l, 'u', mu_u) );
0254 
0255 results = mpc;
0256 [results.bus, results.branch, results.gen, ...
0257     results.om, results.x, results.mu, results.f] = ...
0258         deal(bus, branch, gen, om, x, mu, f);
0259 
0260 pimul = [ ...
0261   results.mu.nln.l - results.mu.nln.u;
0262   results.mu.lin.l - results.mu.lin.u;
0263   -ones(ny>0, 1);
0264   results.mu.var.l - results.mu.var.u;
0265 ];
0266 raw = struct('xr', x, 'pimul', pimul, 'info', info);