check_feasibility

CHECK_FEASIBILITY  Returns the maximum constraint violation in p.u.
   V = CHECK_FEASIBILITY(MPC)
   [V, F, HN, GN, AL, AU, XL, XU] = CHECK_FEASIBILITY(MPC, MPOPT)

   Not thoroughly tested.

0001 function [v, f, hn, gn, Al, Au, xl, xu] = check_feasibility(mpc, mpopt)
0002 %CHECK_FEASIBILITY  Returns the maximum constraint violation in p.u.
0003 %   V = CHECK_FEASIBILITY(MPC)
0004 %   [V, F, HN, GN, AL, AU, XL, XU] = CHECK_FEASIBILITY(MPC, MPOPT)
0005 %
0006 %   Not thoroughly tested.
0007 
0008 %   by Ray Zimmerman, PSERC Cornell
0009 %   Copyright (c) 2010 by Power System Engineering Research Center (PSERC)
0010 %
0011 %   This file is part of MATPOWER.
0012 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0013 %
0014 %   MATPOWER is free software: you can redistribute it and/or modify
0015 %   it under the terms of the GNU General Public License as published
0016 %   by the Free Software Foundation, either version 3 of the License,
0017 %   or (at your option) any later version.
0018 %
0019 %   MATPOWER is distributed in the hope that it will be useful,
0020 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0021 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0022 %   GNU General Public License for more details.
0023 %
0024 %   You should have received a copy of the GNU General Public License
0025 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0026 %
0027 %   Additional permission under GNU GPL version 3 section 7
0028 %
0029 %   If you modify MATPOWER, or any covered work, to interface with
0030 %   other modules (such as MATLAB code and MEX-files) available in a
0031 %   MATLAB(R) or comparable environment containing parts covered
0032 %   under other licensing terms, the licensors of MATPOWER grant
0033 %   you additional permission to convey the resulting work.
0034 
0035 %%----- initialization -----
0036 %% define named indices into data matrices
0037 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0038     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0039 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0040     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0041     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0042 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0043     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0044     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0045 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0046 
0047 if nargin < 2
0048     mpopt = mpoption;
0049 end
0050 
0051 [mpc, mpopt] = opf_args(mpc, mpopt);
0052 mpc = ext2int(mpc);
0053 om = opf_setup(mpc, mpopt);
0054 om = build_cost_params(om);
0055 
0056 %% unpack data
0057 mpc = get_mpc(om);
0058 [vv, ll, nn] = get_idx(om);
0059 
0060 %% problem dimensions
0061 nb = size(mpc.bus, 1);      %% number of buses
0062 nl = size(mpc.branch, 1);   %% number of branches
0063 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0064 
0065 %% linear constraints
0066 [A, l, u] = linear_constraints(om);
0067 
0068 %% bounds on optimization vars
0069 [x, xmin, xmax] = getv(om);
0070 
0071 %% build admittance matrices
0072 [Ybus, Yf, Yt] = makeYbus(mpc.baseMVA, mpc.bus, mpc.branch);
0073 
0074 %% set y variables (ONLY IMPLEMENTED FOR PG)
0075 if ny > 0
0076     ipwl = find(mpc.gencost(:, MODEL) == PW_LINEAR);
0077     ig = vv.i1.Pg:vv.iN.Pg;
0078     x(vv.i1.y:vv.iN.y) = totcost(mpc.gencost(ipwl, :), x(ig(ipwl)) * mpc.baseMVA);
0079 end
0080 
0081 %% find branches with flow limits
0082 il = find(mpc.branch(:, RATE_A) ~= 0 & mpc.branch(:, RATE_A) < 1e10);
0083 nl2 = length(il);           %% number of constrained lines
0084 
0085 %%-----  run opf  -----
0086 f_fcn = @(x)opf_costfcn(x, om);
0087 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0088 hess_fcn = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0089 
0090 [f, df, d2f] = f_fcn(x);
0091 [hn, gn, dhn, dgn] = gh_fcn(x);
0092 
0093 Ax = A * x;
0094 Au = Ax - u;
0095 Al = l - Ax;
0096 
0097 xu = x - xmax;
0098 xl = xmin - x;
0099 
0100 v = max([hn; abs(gn); Au; Al; xu; xl]);