Home > matpower6.0 > extras > misc > qcqp_opf.m

qcqp_opf

PURPOSE ^

QCQP_OPF Builds a quadratically-constrained quadratic program.

SYNOPSIS ^

function [nVAR, nEQ, nINEQ, C, c, A, a, B, b, S] = qcqp_opf(casedata,model)

DESCRIPTION ^

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"
     http://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, 2016

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 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 %     http://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.
0112 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0113 %   See http://www.pserc.cornell.edu/matpower/ 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005