Dear Gamze
Some years ago I needed to minimize the deviation of the power exchange
(active and reactive) between areas in a power system, calculation the
deviation from the start point. I think it will useful for your problem.
The main idea was to maintain the power exchange between areas as closer is
possible to the start point. When other variables had changed, load and
generation.
This what Professor Zimmerman propose to me at that time:
"So, I don't see a way to add the costs you want directly to the current
model. However, I do have an idea.
How about this? To simplify the explanation, let's suppose each area
interchange can be defined as a summation of Pf (power injected into the
"from" end) of a subset of branches (the tie lines). Then for each tie line
do the following:
1. If bus f is the current "from" bus for the branch, add a new bus f' and
connect the tie line to bus f' instead of bus f.
2. Add a new dummy generator at bus f and another at bus f', to represent
the flow in the branch.
3. Add an equality constraint between the voltage angle at bus f and at bus
f'.
4. Add an equality constraint between the voltage magnitude at bus f and at
bus f'.
5. Add an equality constraint between the power injected by the new
generator at bus f and the negative of the power injected by the new
generator at bus f' (one constraint for real power and one for reactive
power).
The standard OPF solution for this new network model should match that of
the original, but now you have optimization variables that are equal to the
flows in the tie lines and you can put costs on them. You can use N and
rhat in the user-defined cost model to compute the deviation from the
scheduled interchange and put a cost on that.
Would that give you what you are looking for?"
I implemented the function "min_power_exchange.m
<http://electrica.uc3m.es/alvaro/mat/min_power_exchange.m>" and solved my
problem. But this function has a bug and I didn't fixed yet. The function
generate island systems by area, I assume, if there is no reference bus in
in each island system it crash. One solution, put a reference bus en each
subsystem and the function will work.
I hope this approach will be useful for your problem.
Best regards,
Álvaro Jaramillo Duque
function [results, success] = min_power_exchange(mpc, mpopt)
%MIN_POWER_EXCHANGE Minimization of the power exchanged across the tie lines
% [RESULTS, SUCCESS] = RUNOPF(CASEDATA, MPOPT)
%
% Minimize the deviation of the power exchange (active and reactive)
% between areas. Deviation computed from the start point of the
% system.
%
% Inputs:
% CASEDATA : MATPOWER case struct, the struct must to come from a
% soleved case with flows in the branchs.
% See also CASEFORMAT and LOADCASE.
% MPOPT : MATPOWER options vector to override default options
% can be used to specify the solution algorithm, output options
% termination tolerances, and more (see also MPOPTION).
%
% Outputs (all are optional):
% RESULTS : results struct, with the following fields:
% (all fields from the input MATPOWER case, i.e. bus, branch,
% gen, etc., but with solved voltages, power flows, etc.)
% order - info used in external <-> internal data conversion
% et - elapsed time in seconds
% success - success flag, 1 = succeeded, 0 = failed
% (additional OPF fields, see OPF for details)
% SUCCESS : the success flag can additionally be returned as
% a second output argument
%
% Calling syntax options:
% results = min_power_exchange(casedata);
% results = min_power_exchange(casedata, mpopt);
%
% Example:
% mpc = loadcase('case30.m');
% mpopt = mpoption('VERBOSE',0,'OUT_ALL',0);
% mpc = runopf(mpc,mpopt);
% define_constants;
% %% create map of external bus numbers to bus indices
% i2e = mpc.bus(:, BUS_I);
% e2i = sparse(max(i2e), 1);
% e2i(i2e) = (1:size(mpc.bus, 1))';
% %% ties lines
% ties = find(mpc.bus(e2i(mpc.branch(:, F_BUS)), BUS_AREA) ~=
mpc.bus(e2i(mpc.branch(:, T_BUS)), BUS_AREA));
% % ties_idx = [mpc.branch(ties, F_BUS) mpc.branch(ties, T_BUS)];
% flow_t0 = [mpc.branch(ties, PF) mpc.branch(ties, QF)];
% % Reduction of the load at bus 26
% mpc.bus(26,[PD QD]) = 0.8*mpc.bus(26,[PD QD]);
% results = min_power_exchange(mpc,mpopt);
% flow_tf = [results.branch(ties, PF) results.branch(ties, QF)];
% % difference of the tie flows
% flow_t0 - flow_tf
% MATPOWER
% $Id: min_power_exchange.m
% by Ray Zimmerman, PSERC Cornell
% and Álvaro Jaramillo Duque, UC3M
% Copyright (c) 2008-2011 by Power System Engineering Research Center
(PSERC)
% Copyright (c) 2011 by UC3M
%
% This file is part of MATPOWER.
% See http://www.pserc.cornell.edu/matpower/ for more info.
%
% MATPOWER is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published
% by the Free Software Foundation, either version 3 of the License,
% or (at your option) any later version.
%
% MATPOWER is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
%
% Additional permission under GNU GPL version 3 section 7
%
% If you modify MATPOWER, or any covered work, to interface with
% other modules (such as MATLAB code and MEX-files) available in a
% MATLAB(R) or comparable environment containing parts covered
% under other licensing terms, the licensors of MATPOWER grant
% you additional permission to convey the resulting work.
% define named indices into bus, gen, mpc.branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
if ~(size(mpc.branch,2) >= 17)
error('The system must to come from a solved case with flows in the
branch’s')
end
if nargin == 1
mpopt = mpoption; %% default options
end
%% data dimensions
nb = size(mpc.bus, 1); %% number of buses
nl = size(mpc.branch, 1); %% number of branches
ng = size(mpc.gen, 1); %% number of dispatchable injections
%% create map of external bus numbers to bus indices
i2e = mpc.bus(:, BUS_I);
e2i = sparse(max(i2e), 1);
e2i(i2e) = (1:size(mpc.bus, 1))';
%% parameters
ties = find(mpc.bus(e2i(mpc.branch(:, F_BUS)), BUS_AREA) ~= ...
mpc.bus(e2i(mpc.branch(:, T_BUS)), BUS_AREA));
% ties lines Inter-ties
nt = length(ties);
% save the original generator cost function
gencost_original = mpc.gencost;
% for each tie line add a virtual bus f'
f.index = e2i(mpc.branch(ties, F_BUS));
f.bus_i = mpc.bus(f.index,BUS_I);
mpc.bus = ([mpc.bus; mpc.bus(f.index,:)]);
% fixing the f' bus index
f_v.index = (nb+1):(nb+nt);
f_v.bus_i = [(max(mpc.bus(:, BUS_I))+1):(max(mpc.bus(:, BUS_I))+nt) ]';
mpc.bus( f_v.index , BUS_I ) = f_v.bus_i;
% new bus number
nb_n = nb + nt;
% setting some values of the virtual bus f'
mpc.bus(f_v.index, [PD, QD, GS, BS] ) =0;
mpc.bus(f_v.index, BUS_TYPE ) = 2;
% setting some values of the bus f
mpc.bus(f.index, [BUS_TYPE] ) = 2;
% for each tie line, change the connection of the f bus to the virtual f' bus
mpc.branch(ties, F_BUS) = f_v.bus_i;
% Add a new dummy generator at bus f and another at virtual f' bus, to
represent
% the flow in the mpc.branch
dummy_g = zeros(nt, size(mpc.gen,2)); % initialize variable
dummy_g(:,GEN_BUS) = f.bus_i; % index of the f bus for each tie line
dummy_g(:,PG ) = -mpc.branch(ties,PF); % negative flow, generated flow is
entering in the node
dummy_g(:,QG ) = -mpc.branch(ties,QF); %
dummy_g(:,QMAX ) = mpc.branch(ties,RATE_A);
dummy_g(:,QMIN ) = -mpc.branch(ties,RATE_A);
dummy_g(:,VG ) = mpc.bus(f.index,VM); % voltage magnitude p.u.
dummy_g(:,MBASE ) = mpc.baseMVA;
dummy_g(:,GEN_STATUS ) = ones(nt,1);
dummy_g(:,PMAX ) = mpc.branch(ties,RATE_A);
dummy_g(:,PMIN ) = -mpc.branch(ties,RATE_A);
% generators at f buses
mpc.gen = [mpc.gen; dummy_g ];
% generators at f' buses
dummy_g(:,GEN_BUS) = f_v.bus_i; % index of the f' bus for each tie line
dummy_g(:,PG ) = -dummy_g(:,PG ); % negative flow respect to the generator
connected in f
dummy_g(:,QG ) = -dummy_g(:,QG );
mpc.gen = [mpc.gen; dummy_g ];
% new generator number
ng_n = ng+nt*2;
% zero function cost
mpc.gencost = ones(ng_n,1)*[2, 0, 0, 2, 0, 0];
gen.f = [ zeros(ng,1) ; (1:nt)' ; zeros(nt,1) ];
gen.f_v = [ zeros(ng,1) ; zeros(nt,1) ; (1:nt)' ];
index_gen = NaN(nt,2);
cooples = NaN(nt,2);
for i=1:nt
index_gen(i,:) = [ find(gen.f==i) , find(gen.f_v==i) ];
cooples(i,:) = [ mpc.gen(index_gen(i,1), 1) , mpc.gen(index_gen(i,2),
1) ];
end
% positions x = [Va Vm Pg Qg]'
thbas = 1; thend = thbas+nb_n-1; % Va
vbas = thend+1; vend = vbas +nb_n-1; % Vm
pgbas = vend+1; pgend = pgbas+ng_n-1; % Pg
qgbas = pgend+1; qgend = qgbas+ng_n-1; % Qg
nxyz = 2*nb_n + 2*ng_n;
% V == V'
A_V = sparse( [1:nt 1:nt], [ f.index' f_v.index] , [ ones(1,nt) ,-ones(1,nt)]
, nt , nb_n);
%PQ == -PQ'
A_PQg = sparse( [1:nt 1:nt], index_gen(:) , [ ones(1,nt) ones(1,nt)] , nt ,
ng_n);
% Va == Va'
A_Va = [A_V sparse(nt, nb_n + ng_n*2)];
% Vm == Vm'
A_Vm = [ sparse(nt,nb_n) A_V sparse(nt,ng_n*2)];
%Pg == -Pg'
A_Pg = [ sparse(nt,nb_n*2) A_PQg sparse(nt,ng_n)];
%Qg == -Qg'
A_Qg = [ sparse(nt, nb_n*2 + ng_n) A_PQg ];
% equality restrictions
mpc.A = [ A_Va ; A_Vm ; A_Pg ; A_Qg ];
%spy(mpc.A); % to see the matrix structure, is OK!
%limits
mpc.l = zeros(nt*4 ,1);
mpc.u = zeros(nt*4 ,1);
% restrictions for the virtual generators Figure 5-2 in the MATPOWER Manual
mpc.N = sparse((1:nt*2)', [((pgend-nt+1):pgend)';((qgend-nt+1):qgend)'] ,
mpc.baseMVA * ones(nt*2,1), nt*2, nxyz);
mpc.fparm = ones(nt*2,1) * [ 2 0 0 1 ]; % [d rh k m]
% using line flow as start point
mpc.fparm(:,2) = [ mpc.gen(index_gen(:,2),PG) ;mpc.gen(index_gen(:,2),QG)]; %
rh
mpc.H = sparse(0,0);
mpc.Cw = ones(nt*2,1);
% do not print
mpopt_o = mpopt;
mpopt = mpoption(mpopt,'VERBOSE',0,'OUT_ALL',0);
%% run OPF with user-defined cost, quadratic term
[results, success] = runopf(mpc, mpopt);
%restart the original options
mpopt = mpopt_o;
results.bus = results.bus(1:nb,:); % remove the virtual f' buses
results.gen = results.gen(1:ng,:); % remove the virtual generators
results.branch = results.branch;
results.branch(ties, F_BUS) = f.bus_i; % fixing id's
results.gencost = gencost_original;
% remove the extra terms
clear results.A results.l results.u results.N results.fparm results.H
results.Cw
if success
printpf(results, 1, mpopt);
end
if nargout == 0
clear results
end
