FixedReserve constraint as a percentage of some unit's power (Pg)

[lordleoo]

The fixedReserve feature in MOST and MATPOWER allows setting a fixed reserve requirement in power units (i.e. in MW).
The system can be divided into zones, and a fixed requirement can be set for each zone.
Of course, the system can be declared as one zone, and one fixed reserve value (in MW) is required.

However, going through the literature, i came across many rules-of-thumb for setting the reserve, such as:

1. reserve must be larger than the output of the largest online unit
2. you want to operate some renewable energy source at 5% below its maximum (available) output
3. reserve must be larger than a percentage of a dispatchable/elastic load.

it took me less than an hour to implement this in MOST. but i haven't gone through MATPOWER files to do it.

in the 'mpc.reserve' structure, i'm just going to define a new field called pweights
pweights has as many rows as zones
pweights has as many columns as the number of generators
that is:
size(mpc.reserves.pweights) = [nrz , ng]
where
nrz = size(mpc.reserves.zones,1);
ng = size(mpc.gen,1);

Each element in 'pweights' has a coefficient to be multiplied by Pg of a unit.
for example,
if i have one zone, and i want the reserve to be at least 50% of the power output of the third generator, then:
pweights = [0 0 0.5 0]; % assuming i have 4 generators

The required changes in most.m are just 2 lines of code. I'll submit this as a pull-request later. my thesis is due in 3 weeks.
the mathematical formula of the old constraint is:
sum_i(igr) { R(i(igr)} >= fixed_reserve_requirement
the fixed reserve requirement is called Rreq in most.m
the new mathematical formula allows both, fixed reserve (in MW) and a percentage of a unit's output
sum_i(igr) { R(i(igr)} >= fixed_reserve_requirement + coefficient * Pg
rearrange to have all variables on one side
sum_i(igr) { R(i(igr)} - coefficient * Pg>= fixed_reserve_requirement

usually i'd set the fixed_reserve_requirement in this case to 0.
but this way, the constraint is generic. you can apply either type or both.

code changes:
at line 1042 of most.m
{
% original block of code:

          A = r.zones(:, r.igr);
          l = r.req / mpc.baseMVA;
          vs = struct('name', {'R'}, 'idx', {{t,j,k}});
          om.add_lin_constraint('Rreq', {t,j,k}, A, l, [], vs);

}

% new block of code

if ~isfield(r,'pweights')
r.pweights = zeros(nrz,ng); % for backward-compatibility
end
assert(size(r.pweights,2)==ng,'Number of columns of reserve.pweights mst be = n_gen');
A = [r.zones(:, r.igr)]; % P is already in pu, no need to divide by baseMVA
l = r.req / mpc.baseMVA;
vs = struct('name', {'R','Pg'}, 'idx', {{t,j,k}, {t,j,k}});
om.add_lin_constraint('Rreq', {t,j,k}, [A, -r.pweights], l, [], vs);

% ================================================
a few notes:

• let ng be the number of generators

• let igr be the indices of generators which are willing to provide reserve. if all units provide reserve, then igr=1:ng;

• to implement the reserve criterion that: total reserve >= max(Pg(:))

  • you will need to define ng zones. that is, create r.zones which has ng rows. (such that size(reserves.zones,1) = ng)
  • define pweights with ng rows; that is: size(reserves.pweights,1) = ng'
    and remember to:
  • set the coefficient in r.pweights corresponding to dispatchable loads to zero, (column of dispatchable load should be all zeros)
  • and don't create a new zone for dispatchable loads. that is, delete the row corresponding to the dispatchable load

  now define these:

mpc.reserves.zones = ones(ng,igr);
mpc.reserves.pweights = eye(ng,ng);
%{
% if you had dispatchable load, do this:
which_dispatchable = isload(mpc.gen);
mpc.reserves.zones(:,which_dispatchable)=0;
mpc.reserves.zones(which_dispatchable,:)=[ ];
mpc.reserves.pweights(:,which_dispatchable)=0;
%}
mpc.reserves.req = zeros(size(mpc.reserves.zones,1),1); % or some margin of your liking

• this code here assumes that the reserve from generator (i) is included in the total reserve to replace it. to exclude a generator from contributing reserve to replace its own self, remove its entry from r.zones. for example, remove generator 1 from the total reserve required to replace generator 1

r.zones(1,1)=0;
%in general
r.zones = r.zones - eye(size(r.zones));

be aware that in small systems (used for tutorials and simple examples), the number of generators is small. requiring all generators to be able to replace on full generator might make the problem infeasible.

• if you want your fixed_reserve to be a percentage of some fixed load, you don't need to use the pweights method. just grab the sum of all loads from mpc.bus(:,PD) and throw it in reserves.qty
• if you want your fixed_reserve to be a percentage of some dispatchable load, then remember to use negative weight for the load in pweights, because dispatchable load is negative

% ========================================================================
% MOST example 3; see how fixed_resreves work

% ADD THESE LINES AT THE END OF EX_CASE3B
%{
ng=size(mpc.gen,1);
which_load = isload(mpc.gen);
igr=sum(~which_load);
mpc.reserves.zones = ones(ng,igr);
mpc.reserves.zones(which_load,:)=[];
mpc.reserves.zones(:,which_load)=0;
mpc.reserves.pweights = eye(size(mpc.reserves.zones,1),ng);
mpc.reserves.req = zeros(size(mpc.reserves.zones,1),1);
% 
%}
clc; clear all; define_constants;
% mpc = loadcase('ex_case3_2zones');
mpc = loadcase('ex_case3b');
% mpopt = mpoption(mpopt, 'most.dc_model', 0); % mpopt must have already been defined before. this line will throw an error
mpopt = mpoption('most.dc_model', 0); % use model with no network
mdi = loadmd(mpc);
mdi.FixedReserves = mpc.reserves; % include fixed zonal reserves
% It is NOT enough to have only mpc.reserves defined and most.fixed_res=1
% you must also defined mdi.FixedReserves

mdo = most(mdi, mpopt);
r2 = mdo.flow.mpc;
Pg2 = r2.gen(:, PG);        % active generation
lam2 = r2.bus(:, LAM_P);    % nodal energy price
R2 = r2.reserves.R;         % reserve quantity
prc2 = r2.reserves.prc;     % reserve price
ms = most_summary(mdo);
display(R2);
disp(['Total reserve is ',num2str(sum(R2(:)))]);

% ========================================================================
% this is the case file now. notice the definition of reserve.pweights

function mpc = ex_case3_2zones
%   MOST
%   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
1     3     0     0     0     0     1     1     0   135     1  1.05  0.95
2     2     0     0     0     0     1     1     0   135     1  1.05  0.95
3     2     0     0     0     0     1     1     0   135     1  1.05  0.95
];

%% generator data
%bus Pg    Qg   Qmax  Qmin   Vg	mBase status Pmax	Pmin    Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
1   125     0    25   -25     1   100     1   200    60     0     0     0     0     0     0     0   250   250     0     0
1   125     0    25   -25     1   100     1   200    65     0     0     0     0     0     0     0   250   250     0     0
2   200     0    50   -50     1   100     1   500    60     0     0     0     0     0     0     0   600   600     0     0
3  -450     0     0     0     1   100     1     0  -450     0     0     0     0     0     0     0   500   500     0     0
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
1      2  0.005   0.01      0    300    300    300      0      0      1   -360    360
1      3  0.005   0.01      0    240    240    240      0      0      1   -360    360
2      3  0.005   0.01      0    300    300    300      0      0      1   -360    360
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
2      0      0      2     25      0
2    200    200      2     30      0
2   3000    600      2     40      0
2      0      0      2   1000      0
];

%%-----  Reserve Data  -----%%
%% (same order as gens, but skipping any gen that does not belong to any zone)
mpc.reserves.cost  = [	1;	3;	5;	];

%% OPTIONAL max reserve quantities for each gen that belongs to at least 1 zone
%% (same order as gens, but skipping any gen that does not belong to any zone)
mpc.reserves.qty   = [	100;	100;	200;	];

switch 1
case 0
disp('Two zones of reserve, requirement is a fixed value in MW');
mpc.reserves.zones = [
	1	0	0	0;
    0   1   1   0;
];
% number of columns  of pweights must be the same as the number of generators
mpc.reserves.pweights = [
0	0	0       0; % no relative constraint on reserve in this zone
0   0   0.75    0; % reserve in zone 1 must be > 0.75 Pg of generator 3 + fixed margin of 50 (see below)
];
mpc.reserves.req   = [50; 50];

case 1 % reserve is larger than max(Pg(:))
disp('Reserve requirement is: larger than max(Pg)');
mpc.reserves.zones = ones(3);
mpc.reserves.pweights = eye(3,4);
mpc.reserves.req = zeros(3,1);

case 2 % reserve is > 25% of dispatchable load
disp('Two zones of reserve, requirement is 25% of total generation ( = load)');
mpc.reserves.zones = [1	1	1	0];
mpc.reserves.pweights = [0 0 0 -0.25]; %notice the negative
mpc.reserves.req   = [0];
otherwise
error('No can do');
end

[rdzman]

I like this! Nice idea.

Pull requests are welcome, but I think I'd prefer to start by adding it to MATPOWER where it'll be easier to test and document.

[lordleoo]

I'm glad you liked it. I hope this makes up for my other issue earlier today.

For MOST manual, you can refer to the following papers and make the following argument:

This feature is particularly useful with renewable energy sources. According to [1-4], operating a wind turbine at 5% below the available (maximum) power from wind contributes significant cost savings on primary reserve.
i know that this margin we're talking about is a fraction of PMAX, and it can be defined in mpc.reserve.req as a fixed number; but for a generic one-time solution that applies to all time steps (t) and all scenarios (j), you can just enter it as
pweight(:,iwind) = (1/0.95)-1=0.0526;

[1] Ye Wang ; Herman Bayem ; Maria Giralt-Devant ; Vera Silva ; Xavier Guillaud ; Bruno Francois, " Methods for Assessing Available Wind Primary Power Reserve," IEEE {Trans.} on Sustainable Energy, vol. 6, no. 1, pp. 272-280, January 2015.
[2] E. Saiz-Marin, J. Garcia-Gonzalez, J. Barquin, and E. Lobato, “Economic assessment of the participation of wind generation in the secondary regulation market,” IEEE Trans. Power Syst., vol. 27, no. 2, pp. 866–874, May 2012.
[3] M. Hedayati-Mehdiabadi, K. W. Hedman, and J. Zhang, “Reserve policy optimization for scheduling wind energy and reserve,” IEEE Trans. Power Syst., vol. 33, no. 1, pp. 19–31, Jan. 2018.

and if i may cite my ownself:
[4] B. Mohandes, M. S. E. Moursi, N. Hatziargyriou and S. E. Khatib, "A Review of Power System Flexibility With High Penetration of Renewables," in IEEE Transactions on Power Systems, vol. 34, no. 4, pp. 3140-3155, July 2019. doi: 10.1109/TPWRS.2019.2897727

[lordleoo]

Hi. I'm back. I've implemented this in MATPOWER.
i should've done it earlier, it's just 4 lines of code

in MATPOWER library, toggle_reserves.m
replace the line: om.add_lin_constraint('Rreq', r.zones(:, igr), lreq, [], {'R'}); %the original
replace with:

if ~isfield(r,'pweights')
    nrz = size(r.zones,1);
    r.pweights = zeros(nrz,ng); % remember to check how & where 'r' was created, and whether its fields are transferred from another structure
  end
  assert(size(r.pweights,2)==ng,'Number of columns of reserve.pweights must be = n_gen');
  om.add_lin_constraint('Rreq', [r.zones(:, igr), -r.pweights], lreq, [], {'R','Pg'});

I'm not good with github, not sure how to submit pull-requests. i'll try.
@rdzman Dr. Rayman, can we discuss preparting the documentation for this. what is your preferred style, etc.

[rdzman]

Thanks, again. I really do like this idea.

I've opened MATPOWER issue #104 to track the inclusion of this feature in MATPOWER, so we can move the discussion there until it's in MATPOWER, and use this space for subsequent inclusion into MOST.