Home > matpower7.1 > lib > nlpopf_solver.m

nlpopf_solver

PURPOSE ^

NLPOPF_SOLVER Solves AC optimal power flow using MP-Opt-Model.

SYNOPSIS ^

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

DESCRIPTION ^

NLPOPF_SOLVER  Solves AC optimal power flow using MP-Opt-Model.

   [RESULTS, SUCCESS, RAW] = NLPOPF_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    (deprecated) 2*nb+2*nl - Pmis, Qmis, Sf, St
               .l  lower bounds on nonlinear constraints
               .u  upper bounds on nonlinear constraints
           .nle    nonlinear equality constraints
           .nli    nonlinear inequality 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] = nlpopf_solver(om, mpopt)
0002 %NLPOPF_SOLVER  Solves AC optimal power flow using MP-Opt-Model.
0003 %
0004 %   [RESULTS, SUCCESS, RAW] = NLPOPF_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    (deprecated) 2*nb+2*nl - Pmis, Qmis, Sf, St
0021 %               .l  lower bounds on nonlinear constraints
0022 %               .u  upper bounds on nonlinear constraints
0023 %           .nle    nonlinear equality constraints
0024 %           .nli    nonlinear inequality constraints
0025 %           .lin
0026 %               .l  lower bounds on linear constraints
0027 %               .u  upper bounds on linear constraints
0028 %
0029 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0030 %
0031 %   RAW         raw output in form returned by MINOS
0032 %       .xr     final value of optimization variables
0033 %       .pimul  constraint multipliers
0034 %       .info   solver specific termination code
0035 %       .output solver specific output information
0036 %
0037 %   See also OPF, MIPS.
0038 
0039 %   MATPOWER
0040 %   Copyright (c) 2000-2020, Power Systems Engineering Research Center (PSERC)
0041 %   by Ray Zimmerman, PSERC Cornell
0042 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0043 %
0044 %   This file is part of MATPOWER.
0045 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0046 %   See https://matpower.org 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 % mips_opt = mpopt.mips;
0062 % mips_opt.verbose = mpopt.verbose;
0063 % if mips_opt.feastol == 0
0064 %     mips_opt.feastol = mpopt.opf.violation;  %% = MPOPT.opf.violation by default
0065 % end
0066 % if ~isfield(mips_opt, 'cost_mult') || isempty(mips_opt.cost_mult)
0067 %     mips_opt.cost_mult = 1e-4;
0068 % end
0069 % opt = struct('mips_opt', mips_opt);
0070 
0071 %% unpack data
0072 mpc = om.get_mpc();
0073 [baseMVA, bus, gen, branch, gencost] = ...
0074     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0075 [vv, ll, nne, nni] = om.get_idx();
0076 
0077 %% problem dimensions
0078 nb = size(bus, 1);          %% number of buses
0079 nl = size(branch, 1);       %% number of branches
0080 ny = om.getN('var', 'y');   %% number of piece-wise linear costs
0081 
0082 %% options
0083 model = om.problem_type();
0084 opt = mpopt2nlpopt(mpopt, model);
0085 
0086 %% try to select an interior initial point, unless requested not to
0087 if mpopt.opf.start < 2
0088     [x0, xmin, xmax] = om.params_var();     %% init var & bounds
0089     s = 1;                      %% set init point inside bounds by s
0090     lb = xmin; ub = xmax;
0091     lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
0092     ub(xmax ==  Inf) =  1e10;
0093     x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
0094     k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
0095     x0(k) = xmax(k) - s;                    %% set just below upper bound
0096     k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
0097     x0(k) = xmin(k) + s;                    %% set just above lower bound
0098     Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0099     Vmax = min(bus(:, VMAX), 1.5);
0100     Vmin = max(bus(:, VMIN), 0.5);
0101     Vm = (Vmax + Vmin) / 2;
0102     if mpopt.opf.v_cartesian
0103         V = Vm * exp(1j*Varefs(1));
0104         x0(vv.i1.Vr:vv.iN.Vr) = real(V);
0105         x0(vv.i1.Vi:vv.iN.Vi) = imag(V);
0106     else
0107         x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
0108         x0(vv.i1.Vm:vv.iN.Vm) = Vm;         %% voltage magnitudes
0109         if ny > 0
0110             ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0111             c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
0112             x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0113         end
0114     end
0115     opt.x0 = x0;
0116 end
0117 
0118 %% find branches with flow limits
0119 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0120 nl2 = length(il);           %% number of constrained lines
0121 
0122 %%-----  run opf  -----
0123 [x, f, eflag, output, lambda] = om.solve(opt);
0124 success = (eflag > 0);
0125 
0126 %% update solution data
0127 if mpopt.opf.v_cartesian
0128     Vi = x(vv.i1.Vi:vv.iN.Vi);
0129     Vr = x(vv.i1.Vr:vv.iN.Vr);
0130     V = Vr + 1j*Vi;
0131     Va = angle(V);
0132     Vm = abs(V);
0133 else
0134     Va = x(vv.i1.Va:vv.iN.Va);
0135     Vm = x(vv.i1.Vm:vv.iN.Vm);
0136     V = Vm .* exp(1j*Va);
0137 end
0138 Pg = x(vv.i1.Pg:vv.iN.Pg);
0139 Qg = x(vv.i1.Qg:vv.iN.Qg);
0140 
0141 %%-----  calculate return values  -----
0142 %% update voltages & generator outputs
0143 bus(:, VA) = Va * 180/pi;
0144 bus(:, VM) = Vm;
0145 gen(:, PG) = Pg * baseMVA;
0146 gen(:, QG) = Qg * baseMVA;
0147 gen(:, VG) = Vm(gen(:, GEN_BUS));
0148 
0149 %% compute branch flows
0150 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);    %% build admittance matrices
0151 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);   %% cplx pwr at "from" bus, p.u.
0152 St = V(branch(:, T_BUS)) .* conj(Yt * V);   %% cplx pwr at "to" bus, p.u.
0153 branch(:, PF) = real(Sf) * baseMVA;
0154 branch(:, QF) = imag(Sf) * baseMVA;
0155 branch(:, PT) = real(St) * baseMVA;
0156 branch(:, QT) = imag(St) * baseMVA;
0157 
0158 %% line constraint is typically on square of limit
0159 %% so we must fix multipliers
0160 muSf = zeros(nl, 1);
0161 muSt = zeros(nl, 1);
0162 if ~isempty(il)
0163     if upper(mpopt.opf.flow_lim(1)) == 'P'
0164         muSf(il) = lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf);
0165         muSt(il) = lambda.ineqnonlin(nni.i1.St:nni.iN.St);
0166     else
0167         muSf(il) = 2 * lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf) .* branch(il, RATE_A) / baseMVA;
0168         muSt(il) = 2 * lambda.ineqnonlin(nni.i1.St:nni.iN.St) .* branch(il, RATE_A) / baseMVA;
0169     end
0170 end
0171 
0172 %% update Lagrange multipliers
0173 if mpopt.opf.v_cartesian
0174     if om.userdata.veq
0175         lam = lambda.eqnonlin(nne.i1.Veq:nne.iN.Veq);
0176         mu_Vmax = zeros(size(lam));
0177         mu_Vmin = zeros(size(lam));
0178         mu_Vmax(lam > 0) =  lam(lam > 0);
0179         mu_Vmin(lam < 0) = -lam(lam < 0);
0180         bus(om.userdata.veq, MU_VMAX) = mu_Vmax;
0181         bus(om.userdata.veq, MU_VMIN) = mu_Vmin;
0182     end
0183     bus(om.userdata.viq, MU_VMAX) = lambda.ineqnonlin(nni.i1.Vmax:nni.iN.Vmax);
0184     bus(om.userdata.viq, MU_VMIN) = lambda.ineqnonlin(nni.i1.Vmin:nni.iN.Vmin);
0185 else
0186     bus(:, MU_VMAX)  = lambda.upper(vv.i1.Vm:vv.iN.Vm);
0187     bus(:, MU_VMIN)  = lambda.lower(vv.i1.Vm:vv.iN.Vm);
0188 end
0189 gen(:, MU_PMAX)  = lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0190 gen(:, MU_PMIN)  = lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0191 gen(:, MU_QMAX)  = lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0192 gen(:, MU_QMIN)  = lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0193 if mpopt.opf.current_balance
0194     %% convert current balance shadow prices to equivalent lamP and lamQ
0195     %% P + jQ = (Vr + jVi) * (M - jN)
0196     %% M = (Vr P + Vi Q) / (Vr^2 + Vi^2)
0197     %% N = (Vi P - Vr Q) / (Vr^2 + Vi^2)
0198     %% lamP = df/dP = df/dM * dM/dP + df/dN + dN/dP
0199     %% lamQ = df/dQ = df/dM * dM/dQ + df/dN + dN/dQ
0200     VV = V ./ (V .* conj(V));   %% V / Vm^2
0201     VVr = real(VV);
0202     VVi = imag(VV);
0203     lamM = lambda.eqnonlin(nne.i1.rImis:nne.iN.rImis);
0204     lamN = lambda.eqnonlin(nne.i1.iImis:nne.iN.iImis);
0205     bus(:, LAM_P) = (VVr.*lamM + VVi.*lamN) / baseMVA;
0206     bus(:, LAM_Q) = (VVi.*lamM - VVr.*lamN) / baseMVA;
0207 else
0208     bus(:, LAM_P) = lambda.eqnonlin(nne.i1.Pmis:nne.iN.Pmis) / baseMVA;
0209     bus(:, LAM_Q) = lambda.eqnonlin(nne.i1.Qmis:nne.iN.Qmis) / baseMVA;
0210 end
0211 branch(:, MU_SF) = muSf / baseMVA;
0212 branch(:, MU_ST) = muSt / baseMVA;
0213 
0214 %% package up results
0215 nlnN = 2*nb + 2*nl;     %% because muSf and muSt are nl x 1, not nl2 x 1
0216 
0217 %% extract multipliers for nonlinear constraints
0218 kl = find(lambda.eqnonlin(1:2*nb) < 0);
0219 ku = find(lambda.eqnonlin(1:2*nb) > 0);
0220 nl_mu_l = zeros(nlnN, 1);
0221 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0222 nl_mu_l(kl) = -lambda.eqnonlin(kl);
0223 nl_mu_u(ku) =  lambda.eqnonlin(ku);
0224 
0225 if isfield(lambda, 'ineqnonlin')
0226     lam_nli = lambda.ineqnonlin;
0227 else
0228     lam_nli = [];
0229 end
0230 
0231 mu = struct( ...
0232   'var', struct('l', lambda.lower, 'u', lambda.upper), ...
0233   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0234   'nle', lambda.eqnonlin, ...
0235   'nli', lam_nli, ...
0236   'lin', struct('l', lambda.mu_l, 'u', lambda.mu_u) );
0237 
0238 results = mpc;
0239 [results.bus, results.branch, results.gen, ...
0240     results.om, results.x, results.mu, results.f] = ...
0241         deal(bus, branch, gen, om, x, mu, f);
0242 
0243 pimul = [ ...
0244   results.mu.nln.l - results.mu.nln.u;
0245   results.mu.lin.l - results.mu.lin.u;
0246   -ones(ny>0, 1);
0247   results.mu.var.l - results.mu.var.u;
0248 ];
0249 raw = struct('xr', x, 'pimul', pimul, 'info', eflag, 'output', output);

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005