------------------------------ deprecated ------------------------------
OPF solvers based on CONSTR & LPCONSTR to be removed in future version.
--------------------------------------------------------------------------
FUN_COPF Evaluates objective function and constraints for OPF.
[F, G] = FUN_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT)
[F, G] = FUN_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT, IL)
Objective and constraint evaluation function for AC optimal power
flow, suitable for use with CONSTR.
Inputs:
X : optimization vector
OM : OPF model object
YBUS : bus admittance matrix
YF : admittance matrix for "from" end of constrained branches
YT : admittance matrix for "to" end of constrained branches
AFEQ : sparse A matrix for linear equality constraints
BFEQ : right hand side vector for linear equality constraints
AF : sparse A matrix for linear inequality constraints
BF : right hand side vector for linear inequality constraints
MPOPF : MATPOWER options vector
IL : (optional) vector of branch indices corresponding to
branches with flow limits (all others are assumed to be
unconstrained). The default is [1:nl] (all branches).
YF and YT contain only the rows corresponding to IL.
Outputs:
F : value of objective function
G : vector of constraint values, in the following order
nonlinear equality (power balance)
linear equality
nonlinear inequality (flow limits)
linear inequality
The flow limits (limit^2 - flow^2) can be apparent power
real power or current, depending on value of OPF_FLOW_LIM
in MPOPT (only for constrained lines)0001 function [f, g] = fun_copf(x, om, Ybus, Yf, Yt, Afeq, bfeq, Af, bf, mpopt, il) 0002 %------------------------------ deprecated ------------------------------ 0003 % OPF solvers based on CONSTR & LPCONSTR to be removed in future version. 0004 %-------------------------------------------------------------------------- 0005 %FUN_COPF Evaluates objective function and constraints for OPF. 0006 % [F, G] = FUN_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT) 0007 % [F, G] = FUN_COPF(X, OM, YBUS, YF, YT, AFEQ, BFEQ, AF, BF, MPOPT, IL) 0008 % 0009 % Objective and constraint evaluation function for AC optimal power 0010 % flow, suitable for use with CONSTR. 0011 % 0012 % Inputs: 0013 % X : optimization vector 0014 % OM : OPF model object 0015 % YBUS : bus admittance matrix 0016 % YF : admittance matrix for "from" end of constrained branches 0017 % YT : admittance matrix for "to" end of constrained branches 0018 % AFEQ : sparse A matrix for linear equality constraints 0019 % BFEQ : right hand side vector for linear equality constraints 0020 % AF : sparse A matrix for linear inequality constraints 0021 % BF : right hand side vector for linear inequality constraints 0022 % MPOPF : MATPOWER options vector 0023 % IL : (optional) vector of branch indices corresponding to 0024 % branches with flow limits (all others are assumed to be 0025 % unconstrained). The default is [1:nl] (all branches). 0026 % YF and YT contain only the rows corresponding to IL. 0027 % 0028 % Outputs: 0029 % F : value of objective function 0030 % G : vector of constraint values, in the following order 0031 % nonlinear equality (power balance) 0032 % linear equality 0033 % nonlinear inequality (flow limits) 0034 % linear inequality 0035 % The flow limits (limit^2 - flow^2) can be apparent power 0036 % real power or current, depending on value of OPF_FLOW_LIM 0037 % in MPOPT (only for constrained lines) 0038 0039 % MATPOWER 0040 % $Id: fun_copf.m,v 1.17 2010/06/09 14:56:58 ray Exp $ 0041 % by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales 0042 % and Ray Zimmerman, PSERC Cornell 0043 % Copyright (c) 1996-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 %%----- initialize ----- 0070 %% define named indices into data matrices 0071 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0072 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0073 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0074 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0075 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0076 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0077 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost; 0078 0079 %% default args 0080 if nargin < 11 0081 il = []; 0082 end 0083 0084 %% unpack data 0085 mpc = get_mpc(om); 0086 [baseMVA, bus, gen, branch, gencost] = ... 0087 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost); 0088 cp = get_cost_params(om); 0089 [N, Cw, H, dd, rh, kk, mm] = deal(cp.N, cp.Cw, cp.H, cp.dd, ... 0090 cp.rh, cp.kk, cp.mm); 0091 vv = get_idx(om); 0092 0093 %% problem dimensions 0094 ny = getN(om, 'var', 'y'); %% number of piece-wise linear costs 0095 nxyz = length(x); %% total number of control vars of all types 0096 0097 %% set default constrained lines 0098 if isempty(il) 0099 nl = size(branch, 1); %% number of branches 0100 il = (1:nl); %% all lines have limits by default 0101 end 0102 0103 %% grab Pg & Qg 0104 Pg = x(vv.i1.Pg:vv.iN.Pg); %% active generation in p.u. 0105 Qg = x(vv.i1.Qg:vv.iN.Qg); %% reactive generation in p.u. 0106 0107 %% put Pg & Qg back in gen 0108 gen(:, PG) = Pg * baseMVA; %% active generation in MW 0109 gen(:, QG) = Qg * baseMVA; %% reactive generation in MVAr 0110 0111 %%----- evaluate objective function ----- 0112 %% polynomial cost of P and Q 0113 % use totcost only on polynomial cost; in the minimization problem 0114 % formulation, pwl cost is the sum of the y variables. 0115 ipol = find(gencost(:, MODEL) == POLYNOMIAL); %% poly MW and MVAr costs 0116 xx = [ gen(:, PG); gen(:, QG)]; 0117 if ~isempty(ipol) 0118 f = sum( totcost(gencost(ipol, :), xx(ipol)) ); %% cost of poly P or Q 0119 else 0120 f = 0; 0121 end 0122 0123 %% piecewise linear cost of P and Q 0124 if ny > 0 0125 ccost = full(sparse(ones(1,ny), vv.i1.y:vv.iN.y, ones(1,ny), 1, nxyz)); 0126 f = f + ccost * x; 0127 end 0128 0129 %% generalized cost term 0130 if ~isempty(N) 0131 nw = size(N, 1); 0132 r = N * x - rh; %% Nx - rhat 0133 iLT = find(r < -kk); %% below dead zone 0134 iEQ = find(r == 0 & kk == 0); %% dead zone doesn't exist 0135 iGT = find(r > kk); %% above dead zone 0136 iND = [iLT; iEQ; iGT]; %% rows that are Not in the Dead region 0137 iL = find(dd == 1); %% rows using linear function 0138 iQ = find(dd == 2); %% rows using quadratic function 0139 LL = sparse(iL, iL, 1, nw, nw); 0140 QQ = sparse(iQ, iQ, 1, nw, nw); 0141 kbar = sparse(iND, iND, [ ones(length(iLT), 1); 0142 zeros(length(iEQ), 1); 0143 -ones(length(iGT), 1)], nw, nw) * kk; 0144 rr = r + kbar; %% apply non-dead zone shift 0145 M = sparse(iND, iND, mm(iND), nw, nw); %% dead zone or scale 0146 diagrr = sparse(1:nw, 1:nw, rr, nw, nw); 0147 0148 %% linear rows multiplied by rr(i), quadratic rows by rr(i)^2 0149 w = M * (LL + QQ * diagrr) * rr; 0150 0151 f = f + (w' * H * w) / 2 + Cw' * w; 0152 end 0153 0154 %% ----- evaluate constraints ----- 0155 %% reconstruct V 0156 Va = x(vv.i1.Va:vv.iN.Va); 0157 Vm = x(vv.i1.Vm:vv.iN.Vm); 0158 V = Vm .* exp(1j * Va); 0159 0160 %% rebuild Sbus 0161 Sbus = makeSbus(baseMVA, bus, gen); %% net injected power in p.u. 0162 0163 %% evaluate power flow equations 0164 mis = V .* conj(Ybus * V) - Sbus; 0165 0166 %%----- evaluate constraint function values ----- 0167 %% first, the equality constraints (power flow) 0168 geq = [ real(mis); %% active power mismatch for all buses 0169 imag(mis) ]; %% reactive power mismatch for all buses 0170 0171 %% then, the inequality constraints (branch flow limits) 0172 flow_max = (branch(il, RATE_A)/baseMVA).^2; 0173 flow_max(flow_max == 0) = Inf; 0174 if mpopt(24) == 2 %% current magnitude limit, |I| 0175 If = Yf * V; 0176 It = Yt * V; 0177 gineq = [ If .* conj(If) - flow_max; %% branch current limits (from bus) 0178 It .* conj(It) - flow_max ]; %% branch current limits (to bus) 0179 else 0180 %% compute branch power flows 0181 Sf = V(branch(il, F_BUS)) .* conj(Yf * V); %% complex power injected at "from" bus (p.u.) 0182 St = V(branch(il, T_BUS)) .* conj(Yt * V); %% complex power injected at "to" bus (p.u.) 0183 if mpopt(24) == 1 %% active power limit, P (Pan Wei) 0184 gineq = [ real(Sf).^2 - flow_max; %% branch real power limits (from bus) 0185 real(St).^2 - flow_max ]; %% branch real power limits (to bus) 0186 else %% apparent power limit, |S| 0187 gineq = [ Sf .* conj(Sf) - flow_max; %% branch apparent power limits (from bus) 0188 St .* conj(St) - flow_max ]; %% branch apparent power limits (to bus) 0189 end 0190 end 0191 0192 g = [geq; Afeq * x - bfeq; gineq; Af * x - bf];