QCQP_OPF Builds a quadratically-constrained quadratic program.
[nVAR, nEQ, nINEQ, C, c, A, a, B, b] = QCQP_OPF(CASEDATA,MODEL)
Inputs (all are optional):
CASEDATA : either a MATPOWER case struct or a string containing
the name of the file with the case data (default is 'case9')
(see also CASEFORMAT and LOADCASE)
MODEL : number equal to 0, 1, or 2 to indicate whether the output
matrices are desired to be complex (default value 0), Hermitian
(value 1), or real symmetric (value 2).
Outputs (all are optional):
nVAR : number of variables
nEQ : number of equality constraints
nINEQ : number of inequality constraints
C : square matrix of size nVAR defining the coefficients in the
cost function
c : real number defining the constant term in the cost function
A : cell array of square matrices, each of size nVAR, defining the
coefficients in the equality constraints
a : column vector of size nEQ defining the constant terms in the
equality constraints
B : cell array of square matrices, each of size nVAR, defining the
coefficients in the inequality constraints
b : column vector of size nINEQ defining the constant terms in the
inequality constraints
S : two row matrix which contains the sparsity pattern of the
quadratically-constrained quadratic program, that is to say the set
of indexes i and j between 1 and nVAR such that either x_i x_j,
x_i x_j', or x_i' x_j has a nonzero coefficient in the objective or
constraint functions. (The apostrophe stands for complex
conjugate).
Examples:
[nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9,0);
[nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9,1);
[nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9,2);
[nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case89pegase);
[nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9241pegase);
The optimal power flow problem can be viewed as an instance of
quadratically-constrained quadratic programming. In order for this to
be true, we consider the objective function of the optimal power flow
problem to be a linear function of active power. Higher degree terms
are discarded from the objective function. Moreover, current line flow
constraints are enforced instead of apparent line flow constraints in
order to have quadratic constraints only. The optimal power flow
problem remains non-convex despite the slightly simplified framework
we consider.
The output of this code defines the problem that consists in solving
for a column vector variable x of size nVAR with the aim to
minimize
x' * C * x + c
subject to nEQ equality constraints:
x' * A{k} * x = a(k) , k = 1,...,nEQ,
and subject to nINEQ inequality constraints
x' * B{k} * x <= b(k) , k = 1,...,nINEQ,
where the apostrophe stands for conjugate transpose.
If MODEL == 0 (default value), then
1) x, a, and b are complex vectors;
2) C is a Hermitian matrix;
3) A{1}, ..., A{nEQ}, B{1}, ..., and B{nINEQ} are complex matrices;
4) for k = 1,...,nINEQ, the inequality x' * B{k} * x <= b(k) is defined
by real(x' * B{k} * x) <= real(b(k)) and imag(x' * B{k} * x) <=
imag(b(k));
5) x corresponds to the complex voltages at each bus.
If MODEL == 1, then
1) x is complex vector;
2) a and b are real vectors;
3) C, A{1}, ..., A{nEQ}, B{1}, ..., and B{nINEQ} are Hermitian
matrices;
4) x corresponds to the complex voltages at each bus.
If MODEL == 2, then
1) x, a, and b are real vectors;
2) C, A{1}, ..., A{nEQ}, B{1}, ..., and B{nINEQ} are real symmetric
matrices;
3) x corresponds to the real parts of the complex voltages at each bus
followed by the imaginary parts of the complex voltages at each bus.
When publishing results based on this code, please cite:
C. Josz, S. Fliscounakis, J. Maeght, and P. Panciatici, "AC Power Flow
Data in MATPOWER and QCQP Format: iTesla, RTE Snapshots, and PEGASE"
https://arxiv.org/abs/1603.01533
Contacts:
Cédric Josz, Stéphane Fliscounakis, Jean Maeght, Patrick Panciatici
firstname.lastname@rte-france.com
Réseau de Transport d'Electricité (French Transmission System Operator)
Département Expertise Système, Immeuble "Le Colbert"
9 rue de la Porte de Buc, 78000 Versailles Cedex, France
March 4th, 20160001 function [nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(casedata,model) 0002 %QCQP_OPF Builds a quadratically-constrained quadratic program. 0003 % [nVAR, nEQ, nINEQ, C, c, A, a, B, b] = QCQP_OPF(CASEDATA,MODEL) 0004 % 0005 % Inputs (all are optional): 0006 % CASEDATA : either a MATPOWER case struct or a string containing 0007 % the name of the file with the case data (default is 'case9') 0008 % (see also CASEFORMAT and LOADCASE) 0009 % MODEL : number equal to 0, 1, or 2 to indicate whether the output 0010 % matrices are desired to be complex (default value 0), Hermitian 0011 % (value 1), or real symmetric (value 2). 0012 % 0013 % Outputs (all are optional): 0014 % nVAR : number of variables 0015 % nEQ : number of equality constraints 0016 % nINEQ : number of inequality constraints 0017 % C : square matrix of size nVAR defining the coefficients in the 0018 % cost function 0019 % c : real number defining the constant term in the cost function 0020 % A : cell array of square matrices, each of size nVAR, defining the 0021 % coefficients in the equality constraints 0022 % a : column vector of size nEQ defining the constant terms in the 0023 % equality constraints 0024 % B : cell array of square matrices, each of size nVAR, defining the 0025 % coefficients in the inequality constraints 0026 % b : column vector of size nINEQ defining the constant terms in the 0027 % inequality constraints 0028 % S : two row matrix which contains the sparsity pattern of the 0029 % quadratically-constrained quadratic program, that is to say the set 0030 % of indexes i and j between 1 and nVAR such that either x_i x_j, 0031 % x_i x_j', or x_i' x_j has a nonzero coefficient in the objective or 0032 % constraint functions. (The apostrophe stands for complex 0033 % conjugate). 0034 % 0035 % Examples: 0036 % [nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9,0); 0037 % [nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9,1); 0038 % [nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9,2); 0039 % [nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case89pegase); 0040 % [nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(case9241pegase); 0041 % 0042 % The optimal power flow problem can be viewed as an instance of 0043 % quadratically-constrained quadratic programming. In order for this to 0044 % be true, we consider the objective function of the optimal power flow 0045 % problem to be a linear function of active power. Higher degree terms 0046 % are discarded from the objective function. Moreover, current line flow 0047 % constraints are enforced instead of apparent line flow constraints in 0048 % order to have quadratic constraints only. The optimal power flow 0049 % problem remains non-convex despite the slightly simplified framework 0050 % we consider. 0051 % 0052 % The output of this code defines the problem that consists in solving 0053 % for a column vector variable x of size nVAR with the aim to 0054 % 0055 % minimize 0056 % 0057 % x' * C * x + c 0058 % 0059 % subject to nEQ equality constraints: 0060 % 0061 % x' * A{k} * x = a(k) , k = 1,...,nEQ, 0062 % 0063 % and subject to nINEQ inequality constraints 0064 % 0065 % x' * B{k} * x <= b(k) , k = 1,...,nINEQ, 0066 % 0067 % where the apostrophe stands for conjugate transpose. 0068 % 0069 % If MODEL == 0 (default value), then 0070 % 1) x, a, and b are complex vectors; 0071 % 2) C is a Hermitian matrix; 0072 % 3) A{1}, ..., A{nEQ}, B{1}, ..., and B{nINEQ} are complex matrices; 0073 % 4) for k = 1,...,nINEQ, the inequality x' * B{k} * x <= b(k) is defined 0074 % by real(x' * B{k} * x) <= real(b(k)) and imag(x' * B{k} * x) <= 0075 % imag(b(k)); 0076 % 5) x corresponds to the complex voltages at each bus. 0077 % 0078 % If MODEL == 1, then 0079 % 1) x is complex vector; 0080 % 2) a and b are real vectors; 0081 % 3) C, A{1}, ..., A{nEQ}, B{1}, ..., and B{nINEQ} are Hermitian 0082 % matrices; 0083 % 4) x corresponds to the complex voltages at each bus. 0084 % 0085 % If MODEL == 2, then 0086 % 1) x, a, and b are real vectors; 0087 % 2) C, A{1}, ..., A{nEQ}, B{1}, ..., and B{nINEQ} are real symmetric 0088 % matrices; 0089 % 3) x corresponds to the real parts of the complex voltages at each bus 0090 % followed by the imaginary parts of the complex voltages at each bus. 0091 % 0092 % When publishing results based on this code, please cite: 0093 % 0094 % C. Josz, S. Fliscounakis, J. Maeght, and P. Panciatici, "AC Power Flow 0095 % Data in MATPOWER and QCQP Format: iTesla, RTE Snapshots, and PEGASE" 0096 % https://arxiv.org/abs/1603.01533 0097 % 0098 % Contacts: 0099 % Cédric Josz, Stéphane Fliscounakis, Jean Maeght, Patrick Panciatici 0100 % firstname.lastname@rte-france.com 0101 % Réseau de Transport d'Electricité (French Transmission System Operator) 0102 % Département Expertise Système, Immeuble "Le Colbert" 0103 % 9 rue de la Porte de Buc, 78000 Versailles Cedex, France 0104 % 0105 % March 4th, 2016 0106 0107 % MATPOWER 0108 % Copyright (c) 2016, Power Systems Engineering Research Center (PSERC) 0109 % by Cedric Josz, Jean Maeght, Stephane Fliscounakis, and Patrick Panciatici 0110 % 0111 % This file is part of MATPOWER Extras. 0112 % Covered by the 3-clause BSD License (see LICENSE file for details). 0113 % See https://github.com/MATPOWER/matpower-extras for more info. 0114 0115 %% default arguments 0116 if nargin < 2 0117 model = 0; % default output are complex matrices 0118 if nargin < 1 0119 casedata = 'case9'; % default data file is 'case9.m' 0120 end 0121 end 0122 0123 %% compute admittance matrix 0124 casedata = loadcase(casedata); 0125 mpc = ext2int(casedata); 0126 [Ybus, Yf, Yt] = makeYbus(mpc); 0127 if isfield(mpc,'gencost') == 0 0128 error('Reference to non-existent field ''gencost''.'); 0129 end 0130 if size(mpc.gen,1) ~= size(mpc.gencost) 0131 error('mpc.gen/mpc.gencost: gen and gencost must have the same number of rows'); 0132 end 0133 Ybus = mpc.baseMVA*Ybus; % conversion from per unit to physical unit scaling 0134 % All constraints are in physical units apart from current flow constraints 0135 % which are in per units for better conditioning. 0136 0137 %% define named indices 0138 define_constants; 0139 PQbus = mpc.bus(mpc.bus(:,BUS_TYPE) == 1,BUS_I); % PQ bus numbers 0140 nPQbus = length(PQbus); % number of PQ buses 0141 PVbus = mpc.bus(mpc.bus(:,BUS_TYPE) ~= 1,BUS_I); % PV bus numbers 0142 nPVbus = length(PVbus); % number of PV buses 0143 LFbound = find(mpc.branch(:,RATE_A)>0); % numbers of branches with flow bounds 0144 nLFbound = length(LFbound); % number of lines with flow bounds 0145 0146 %% build cost matrix 0147 % linear function of active power 0148 if sum( mpc.gencost(:,MODEL) ~= 2*ones(size(mpc.gen,1),1) ) 0149 error('mpc.gencost: the objective must be a polynomial and cannot contain piecewise linear terms'); 0150 end 0151 nbus = size(mpc.bus,1); % number of buses 0152 nVAR = nbus; % number of variables 0153 costs = zeros(nbus,1); 0154 costs(mpc.gen(:,GEN_BUS)) = mpc.gencost(:,end-1); % extracting linear cost coefficients that multiply active power in optimal power flow objective function 0155 P = sparse(diag(costs)*Ybus); 0156 C = (P'+P)/2; 0157 c = sum(costs.*mpc.bus(:,PD)) + sum(mpc.gencost(:,end)); % extracting constant costs in optimal power flow objective function 0158 0159 %% build equality matrix and vector 0160 % power balance equations 0161 nEQ = nPQbus; % number of constraints 0162 A = cell(nEQ,1); % initializing cell array 0163 a = zeros(nEQ,1); % initializing vector 0164 for k = 1:nPQbus 0165 num = PQbus(k); % bus number 0166 % The next two lines encode 0167 % V(num)*I(num)' = I' e_num e_num' V = ... 0168 % V' ( Ybus' e_num e_num' ) V = - Pdem(num) - 1i*Qdem(num) 0169 % where e_num is the column vector of size nVAR that contains only one 0170 % nonzero element in position num equal to 1. Matrix Ybus' e_num e_num' 0171 % is equal to Ybus' where all columns except column k are set to zero. 0172 A{k} = sparse(1:nVAR, num, Ybus(num,:)', nVAR, nVAR); 0173 a(k) = -mpc.bus(num,PD) - 1i*mpc.bus(num,QD); 0174 end 0175 0176 %% build inequality matrix and vector 0177 0178 nINEQ = 2*nbus+2*nLFbound+2*nPVbus; 0179 B = cell(nINEQ,1); 0180 b = zeros(nINEQ,1); 0181 0182 % voltage magnitude bounds 0183 for k = 1:nbus 0184 % The next three lines encode 0185 % Vmin(k)^2 <= |V(k)|^2 <= Vmax(k)^2 0186 B{2*k-1} = sparse(k,k,1,nbus,nbus); 0187 B{2*k} = sparse(k,k,-1,nbus,nbus); 0188 b(2*k-1:2*k) = [ mpc.bus(k,VMAX).^2 ; -mpc.bus(k,VMIN).^2 ] ; 0189 end 0190 0191 % current flow bounds 0192 count = 2*nbus; 0193 for k = 1:nLFbound 0194 num = LFbound(k); % branch number 0195 yf = Yf(num,:); % If = Yf * V so If(num) = yf * V where "num" is a branch 0196 yt = Yt(num,:); % It = Yt * V so It(num) = yt * V where "num" is a branch 0197 % The next four lines encode: 0198 % |If(num)|^2 = V' * ( yf'*yf ) * V <= Imax(num)^2 0199 % |It(num)|^2 = V' * ( yt'*yt ) * V <= Imax(num)^2 0200 % Per unit scaling is used for better conditioning. 0201 B{count+2*k-1} = sparse(yf'*yf); 0202 B{count+2*k} = sparse(yt'*yt); 0203 b(count+(2*k-1:2*k)) = [ (mpc.branch(num,RATE_A)/mpc.baseMVA).^2 ; ... 0204 (mpc.branch(num,RATE_A)/mpc.baseMVA).^2 ]; 0205 end 0206 0207 % power generation bounds 0208 count = count+2*nLFbound; 0209 for k = 1:nPVbus 0210 num = PVbus(k); % bus number 0211 % The next six lines encode 0212 % Smin(num) - Sdem(num) <= V(num)*I(num)' <= Smax(num) - Sdem(num) 0213 % where Smin and Smax are lower and upper bounds on complex power 0214 % generation, and where Sdem is the complex power demand. 0215 B{count+2*k-1} = sparse(1:nVAR, num, Ybus(num,:)', nVAR, nVAR); 0216 B{count+2*k} = sparse(1:nVAR, num, -Ybus(num,:)', nVAR, nVAR); 0217 mult_gen = find(mpc.gen(:,GEN_BUS) == num); 0218 b(count+(2*k-1:2*k)) = ... 0219 [ sum(mpc.gen(mult_gen,PMAX))-mpc.bus(num,PD) + 1i*( sum(mpc.gen(mult_gen,QMAX))-mpc.bus(num,QD) ) ; ... % upper bound on complex power generation 0220 -( sum(mpc.gen(mult_gen,PMIN))-mpc.bus(num,PD) + 1i*( sum(mpc.gen(mult_gen,QMIN))-mpc.bus(num,QD) ) ) ]; % lower bound on complex power generation 0221 end 0222 0223 % Compute sparsity pattern 0224 S = unique(sort(mpc.branch(:,1:2),2),'rows'); 0225 0226 %% construct Hermitian output (if MODEL == 1) 0227 if model == 1 0228 0229 % equality constraints 0230 nEQ = 2*nEQ; 0231 HA = cell(nEQ,1); 0232 ha = zeros(nEQ,1); 0233 for k = 1:nEQ/2 0234 HA{2*k-1} = (A{k}+A{k}')/2; 0235 HA{2*k} = (A{k}-A{k}')/(2*1i); 0236 ha(2*k-1) = real(a(k)); 0237 ha(2*k) = imag(a(k)); 0238 end 0239 A = HA; 0240 a = ha; 0241 0242 % inequality constraints 0243 nINEQ = 2*nbus+2*nLFbound+2*(2*nPVbus); 0244 HB = cell(nINEQ,1); 0245 hb = zeros(nINEQ,1); 0246 for k = 1:(2*nbus+2*nLFbound) 0247 HB{k} = B{k}; 0248 hb(k) = b(k); 0249 end 0250 for k = (2*nbus+2*nLFbound+1):(2*nbus+2*nLFbound+2*nPVbus) 0251 HB{-2*nbus-2*nLFbound+2*k-1} = (B{k}+B{k}')/2; 0252 HB{-2*nbus-2*nLFbound+2*k} = (B{k}-B{k}')/(2*1i); 0253 hb(-2*nbus-2*nLFbound+2*k-1) = real(b(k)); 0254 hb(-2*nbus-2*nLFbound+2*k) = imag(b(k)); 0255 end 0256 B = HB; 0257 b = hb; 0258 0259 end 0260 0261 %% construct real symmetric output (if MODEL == 2) 0262 if model == 2 0263 0264 % objective function 0265 C = [real(C) -imag(C); imag(C) real(C)]; 0266 0267 % equality constraints 0268 nEQ = 2*nEQ; 0269 RA = cell(nEQ,1); 0270 ra = zeros(nEQ,1); 0271 for k = 1:nEQ/2 0272 RA{2*k-1} = (A{k}+A{k}')/2; 0273 RA{2*k-1} = [real(RA{2*k-1}) -imag(RA{2*k-1}); ... 0274 imag(RA{2*k-1}) real(RA{2*k-1})]; 0275 RA{2*k} = (A{k}-A{k}')/(2*1i); 0276 RA{2*k} = [real(RA{2*k}) -imag(RA{2*k}); ... 0277 imag(RA{2*k}) real(RA{2*k})]; 0278 ra(2*k-1) = real(a(k)); 0279 ra(2*k) = imag(a(k)); 0280 end 0281 A = RA; 0282 a = ra; 0283 0284 % inequality constraints 0285 nINEQ = 2*nbus+2*nLFbound+2*(2*nPVbus); 0286 RB = cell(nINEQ,1); 0287 rb = zeros(nINEQ,1); 0288 for k = 1:(2*nbus+2*nLFbound) 0289 RB{k} = [real(B{k}) -imag(B{k}); imag(B{k}) real(B{k})]; 0290 rb(k) = b(k); 0291 end 0292 count = -2*nbus-2*nLFbound; 0293 for k = (2*nbus+2*nLFbound+1):(2*nbus+2*nLFbound+2*nPVbus) 0294 RB{count+2*k-1} = (B{k}+B{k}')/2; 0295 RB{count+2*k-1} = [real(RB{count+2*k-1}) -imag(RB{count+2*k-1}); ... 0296 imag(RB{count+2*k-1}) real(RB{count+2*k-1})]; 0297 RB{count+2*k} = (B{k}-B{k}')/(2*1i); 0298 RB{count+2*k} = [real(RB{count+2*k}) -imag(RB{count+2*k}); ... 0299 imag(RB{count+2*k}) real(RB{count+2*k})]; 0300 rb(count+2*k-1) = real(b(k)); 0301 rb(count+2*k) = imag(b(k)); 0302 end 0303 B = RB; 0304 b = rb; 0305 0306 % sparsity pattern 0307 T = zeros(size(S)); 0308 T(:,2) = nVAR; 0309 U = zeros(size(S)); 0310 U(:,1) = nVAR; 0311 S = [S; S+T; S+U; S+T+U]; 0312 0313 % number of variables 0314 nVAR = 2*nVAR; 0315 0316 end 0317 0318 %% run matlab interior point solver 0319 % to do so, uncomment the following section: 0320 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0321 % mpc.gencost(:,COST) = 0; 0322 % results = runopf(mpc,mpoption('opf.ac.solver','MIPS','opf.flow_lim',... 0323 % 'I','out.suppress_detail','1')); 0324 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%