Description of mipsopf

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
