Home > matpower4.1 > lpopf_solver.m

lpopf_solver

PURPOSE ^

------------------------------ deprecated ------------------------------

SYNOPSIS ^

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

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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.22 2011/06/16 17:46:37 cvs 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 %% split l <= A*x <= u into less than, equal to, greater than, and
0104 %% doubly-bounded sets
0105 ieq = find( abs(u-l) <= eps );          %% equality
0106 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0107 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0108 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0109 Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0110 bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
0111 Afeq = A(ieq, :);
0112 bfeq = u(ieq);
0113 
0114 %% build admittance matrices
0115 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0116 
0117 %% find branches with flow limits
0118 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0119 nl2 = length(il);           %% number of constrained lines
0120 
0121 %% set mpopt defaults
0122 mpopt(15)  = 2 * nb + length(bfeq); %% set number of equality constraints
0123 if mpopt(19) == 0           %% CONSTR_MAX_IT
0124   mpopt(19) = 150 + 2*nb;
0125 end
0126 
0127 %% run load flow to get starting point
0128 [x, success_lf] = LPeqslvr(x0, om2, Ybus, Yf, Yt, Afeq, bfeq, Af, bf, mpopt);
0129 if success_lf ~= 1
0130   error('Sorry, cannot find a starting point using power flow, please check data!');
0131 end
0132 
0133 %% set step size
0134 cstep = 0;
0135 if ny > 0
0136   PgQg = [gen(:, PG); gen(:, QG)];
0137   ipwl = find(gencost(:, MODEL) == PW_LINEAR);  %% piece-wise linear costs
0138   Cp = totcost(gencost(ipwl, :), PgQg(ipwl));
0139   cstep = max(abs(Cp));
0140   if cstep < 1.0e6, cstep = 1.0e6; end
0141 end
0142 step0 = ones(size(x0));
0143 step0(vv.i1.Va:vv.iN.Va) = 2;
0144 step0(vv.i1.Vm:vv.iN.Vm) = 1;
0145 step0(vv.i1.Pg:vv.iN.Pg) = 0.6;
0146 step0(vv.i1.Qg:vv.iN.Qg) = 0.3;
0147 if ny > 0
0148   step0(vv.i1.y:vv.iN.y)   = cstep;
0149 end
0150 
0151 %% run optimization
0152 [x, lambda, success] = LPconstr('fun_copf', x0, mpopt, step0, [], [], 'grad_copf', ...
0153                        'LPeqslvr', om2, Ybus, Yf(il,:), Yt(il,:), Afeq, bfeq, Af, bf, mpopt, il);
0154 info = success;
0155 
0156 %% get final objective function value & constraint values
0157 [f, g] = feval('fun_copf', x, om2, Ybus, Yf(il,:), Yt(il,:), Afeq, bfeq, Af, bf, mpopt, il);
0158 
0159 %% update solution data
0160 Va = x(vv.i1.Va:vv.iN.Va);
0161 Vm = x(vv.i1.Vm:vv.iN.Vm);
0162 Pg = x(vv.i1.Pg:vv.iN.Pg);
0163 Qg = x(vv.i1.Qg:vv.iN.Qg);
0164 V = Vm .* exp(1j*Va);
0165 
0166 %%-----  calculate return values  -----
0167 %% update voltages & generator outputs
0168 bus(:, VA) = Va * 180/pi;
0169 bus(:, VM) = Vm;
0170 gen(:, PG) = Pg * baseMVA;
0171 gen(:, QG) = Qg * baseMVA;
0172 gen(:, VG) = Vm(gen(:, GEN_BUS));
0173 
0174 %% compute branch flows
0175 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0176 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0177 branch(:, PF) = real(Sf) * baseMVA;
0178 branch(:, QF) = imag(Sf) * baseMVA;
0179 branch(:, PT) = real(St) * baseMVA;
0180 branch(:, QT) = imag(St) * baseMVA;
0181 
0182 %% package up results
0183 nA = length(u);
0184 neq = length(ieq);
0185 nlt = length(ilt);
0186 ngt = length(igt);
0187 nbx = length(ibx);
0188 
0189 %% extract multipliers (lambda is ordered as Pmis, Qmis, Afeq, Sf, St, Af)
0190 %% nonlinear constraints
0191 inln = [(1:2*nb), (1:2*nl2) + 2*nb+neq];
0192 kl = find(lambda(inln) < 0);
0193 ku = find(lambda(inln) > 0);
0194 nl_mu_l = zeros(2*(nb+nl2), 1);
0195 nl_mu_u = zeros(2*(nb+nl2), 1);
0196 nl_mu_l(kl) = -lambda(inln(kl));
0197 nl_mu_u(ku) =  lambda(inln(ku));
0198 
0199 %% linear constraints
0200 ilin = [(1:neq)+2*nb, (1:(nlt+ngt+2*nbx)) + 2*nb+neq+2*nl2];
0201 kl = find(lambda(ilin(1:neq)) < 0);
0202 ku = find(lambda(ilin(1:neq)) > 0);
0203 
0204 mu_l = zeros(nA, 1);
0205 mu_l(ieq) = -lambda(ilin(1:neq));
0206 mu_l(ieq(ku)) = 0;
0207 mu_l(igt) = lambda(ilin(neq+nlt+(1:ngt)));
0208 mu_l(ibx) = lambda(ilin(neq+nlt+ngt+nbx+(1:nbx)));
0209 
0210 mu_u = zeros(nA, 1);
0211 mu_u(ieq) = lambda(ilin(1:neq));
0212 mu_u(ieq(kl)) = 0;
0213 mu_u(ilt) = lambda(ilin(neq+(1:nlt)));
0214 mu_u(ibx) = lambda(ilin(neq+nlt+ngt+(1:nbx)));
0215 
0216 %% variable bounds
0217 muLB = mu_l(ll.i1.varlims:ll.iN.varlims);
0218 muUB = mu_u(ll.i1.varlims:ll.iN.varlims);
0219 mu_l(ll.i1.varlims:ll.iN.varlims) = [];
0220 mu_u(ll.i1.varlims:ll.iN.varlims) = [];
0221 
0222 %% line constraint is actually on square of limit
0223 %% so we must fix multipliers
0224 muSf = zeros(nl, 1);
0225 muSt = zeros(nl, 1);
0226 muSf(il) = 2 * nl_mu_u((1:nl2)+2*nb    ) .* branch(il, RATE_A) / baseMVA;
0227 muSt(il) = 2 * nl_mu_u((1:nl2)+2*nb+nl2) .* branch(il, RATE_A) / baseMVA;
0228 
0229 %% resize mu for nonlinear constraints
0230 nl_mu_l = [nl_mu_l(1:2*nb); zeros(2*nl, 1)];
0231 nl_mu_u = [nl_mu_u(1:2*nb); muSf; muSt];
0232 
0233 %% update Lagrange multipliers
0234 bus(:, MU_VMAX)  = muUB(vv.i1.Vm:vv.iN.Vm);
0235 bus(:, MU_VMIN)  = muLB(vv.i1.Vm:vv.iN.Vm);
0236 gen(:, MU_PMAX)  = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0237 gen(:, MU_PMIN)  = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0238 gen(:, MU_QMAX)  = muUB(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0239 gen(:, MU_QMIN)  = muLB(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0240 bus(:, LAM_P)    = (nl_mu_u(nn.i1.Pmis:nn.iN.Pmis) - nl_mu_l(nn.i1.Pmis:nn.iN.Pmis)) / baseMVA;
0241 bus(:, LAM_Q)    = (nl_mu_u(nn.i1.Qmis:nn.iN.Qmis) - nl_mu_l(nn.i1.Qmis:nn.iN.Qmis)) / baseMVA;
0242 branch(:, MU_SF) = muSf / baseMVA;
0243 branch(:, MU_ST) = muSt / baseMVA;
0244 
0245 mu = struct( ...
0246   'var', struct('l', muLB, 'u', muUB), ...
0247   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0248   'lin', struct('l', mu_l, 'u', mu_u) );
0249 
0250 results = mpc;
0251 [results.bus, results.branch, results.gen, ...
0252     results.om, results.x, results.mu, results.f] = ...
0253         deal(bus, branch, gen, om, x, mu, f);
0254 
0255 pimul = [ ...
0256   results.mu.nln.l - results.mu.nln.u;
0257   results.mu.lin.l - results.mu.lin.u;
0258   -ones(ny>0, 1);
0259   results.mu.var.l - results.mu.var.u;
0260 ];
0261 raw = struct('xr', x, 'pimul', pimul, 'info', info);

Generated on Mon 26-Jan-2015 15:00:13 by m2html © 2005