Re: Ramp AGC

[Stephen Suffian]

I believe I've seen a similar question asked before. See this link
 which
indicates that RAMP_10 is used for reserve calculations, RAMP_30 for
dispatch scheduling, and RAMP_AGC Is not currently used by matpower
according to this link
.

On Thu, Mar 8, 2018 at 6:51 PM, Bainan Xia  wrote:

> Hi,
>
>
>
> I’m running a huge power system simulation using MOST recently. I’ve got
> some information on ‘ramp rate (% of the rated capacity per min)’ for each
> type of generator from a paper. I would like to add these constraints into
> my case. I’ve found a field in mpc.gen, called ‘RAMP_AGC’, which is ramp
> rate for load following/AGC (MW/min). I’m wondering whether that is the
> right place to fill. Or the following columns ‘RAMP_10’ and ‘RAMP_30’ are
> the right choices. I’m also a little confused on the difference between
> ‘RAMP_10’ and ‘RAMP_30’. What does the word ‘reserve’ mean in the
> documentation.
>
>
>
> Thanks for the help J
>
>
>
> Best,
>
> Bainan
>

Storage on Redispatch

[Stephen Suffian]

Hey everyone,

I have a question regarding how storage is calculated for each time step
during an economic dispatch run. I am attempting to simulate a 24-hour
economic dispatch followed by hourly re-dispatch of units with an updated
load forecast. With an identical forecast, I would assume that the
redispatch would be identical to the initial dispatch. This is in fact the
case with the normal generators (generators in mpc) and also any wind
generators that are added.

However, when I add hydro (as energy-limiting storage units), it seems to
shift around which time periods are using different levels of storage power
during the re-dispatch, despite all other conditions being identical. This
is problematic, because I calculate an added operational cost to the grid
between the initial dispatch and the redispatch, despite the same load
being met.

The following is a simple example with a single generator, a single storage
unit, a flat load, running on a no-bus network, with no terminal targets.
To reiterate, I ideally would want the levels of power dispatched from
storage to be identical between the initial 4-hour dispatch and each
subsequent re-dispatch, but that doesn't seem to be the case below.

Any help understanding better how the storage outputs are calculated, or
whether I am doing something wrong in this code or otherwise, would be
greatly appreciated.

Thanks!
Stephen Suffian

define_constants
this_load_profile.values = [125; 125; 125; 125];
nt = size(this_load_profile.values,1);

warning('off','all')
mpopt = mpoption('verbose', 0);
mpopt = mpoption(mpopt, 'most.dc_model', 0);% use model with no network
mpopt = mpoption(mpopt, 'most.storage.cyclic', 0);
mpopt = mpoption(mpopt, 'most.storage.terminal_target',0);

mpc_orig.version = '2';
mpc_orig.baseMVA = 100;
mpc_orig.bus = [ 1300001101351
1.050.95;];
mpc_orig.gen = [
1125025-2511001200000
0000025060000;
1-200000110010-60000
0000050050000;
];
mpc_orig.branch = [];
mpc_orig.gencost = [
20030   1  0;
200301000;
];
xgd_table.colnames = {};
xgd_table.data = [];

%Storage
storage.gen=[ 1000011001150 00
00000151500];
storage.xgd_table.colnames = {};
storage.xgd_table.data = [];
storage.sd_table.colnames = {
'InitialStorage', ...
'InitialStorageLowerBound', ...
'InitialStorageUpperBound', ...
'InitialStorageCost', ...
'TerminalStoragePrice', ...
'MinStorageLevel', ...
'MaxStorageLevel', ...
'OutEff', ...
'InEff', ...
'LossFactor', ...
'rho', ...

'ExpectedTerminalStorageMin',...

'ExpectedTerminalStorageMax',...
};
storage.sd_table.data=[3535350003511
00  0 0];

for index = 1: 2
mpc = mpc_orig;
xgd = loadxgendata(xgd_table, mpc_orig);
profiles = getprofiles(this_load_profile);
profiles(1).values=profiles(1).values(index:end);
[iess, mpc, xgd, sd] = addstorage(storage, mpc, xgd);
remaining_nt = nt-index+1;
mdi = loadmd(mpc, remaining_nt, xgd, sd, [], profiles);
mdo_re = most(mdi, mpopt);
mdo_re.results.ExpectedDispatch
end


%%OUTPUT (The third row should be identical (so if the first ans is [0 15
15 5], the second should be [15 15 5], as it is the re-dispatch at time t=2.

ans =

  125.  110.  110.  120.
 -125. -125. -125. -125.
 0   15.   15.5.


ans =

  120.  110.  110.
 -125. -125. -125.
5.   15.   15.

Advice Regarding UCED

[Stephen Suffian]

I am attempting to look at the value of reduced forecasting error by
running a deterministic 24-hour UC using MOST, followed by a re-dispatch
and each hour (where it runs the same UCED but with the correct demand
value at the first time point in the redispatch, and then the forecasted
value for the rest fo the timepoints). I am currently running into a
problem where the re-dispatch doesn't have the correct initial conditions.

For each redispatch, I was able to assign the correct ramping constraints
at t=0 using profiles that adjust the PMAX and PMIN values of each
generator. However, I then ran into a problem where the re-dispatch is
choosing to de-commit generators at t=0 that were previously committed. I
used the profiles again to enforce the units that are already committed by
changing the CommitKey in the xgd table. I believe there may have been a
bug in apply_profile.m that I found that prevents people from updating the
xgd table with a profile, but I think I found the solution and submitted it
as an issue on github.

My end goal is to look at the added system cost of forecasting error, much
like was done in this paper <http://ieeexplore.ieee.org/document/5428351/>,
where they conduct hourly dispatch using MATPOWER opf. My current issue is
in regards to adding reserves. I don't want to add fixed zonal reserves, as
that won't help with the forecasting error. Instead, I'd like to add some
form of ramping reserves. I've played with the
PositiveLoadFollowingQuantity in the xgd table, but that seems to provide
an additional ramping constraint on each generator rather than guarantee a
level of ramping reserve.

I seem to have functioning code for the small tutorial 3 bus test grid, but
when I move up to the iEEE 30 bus network, I am still running into several
issues where I either see no difference in cost or there is significant
load shedding.

Does anyone have any advice in regards to how to compare an imperfect UCED
with hourly dispatch vs. a perfect UCED in Matpower MOST? Perhaps I am
going about this all wrong. Any advice in terms of generally how to conduct
that comparison would be greatly appreciated.

Thanks!
Stephen Suffian

Re: Non-optimized UCED Depending on VOLL

[Stephen Suffian]

Apologies, please ignore this email. I needed to add the cplex options into
mpopts. When I add the following, it works as intended:

mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);   %% dual simplex
mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.mipgap', 0);
mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.absmipgap', 0);
mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);

On Mon, Sep 25, 2017 at 11:37 AM, Stephen Suffian <ssuff...@villanova.edu>
wrote:

> I am running a very simple 3 plant grid simulation, and am not sure
> whether I have come across a bug due to a rounding error somewhere in the
> code, or a misunderstanding on my part in terms of calculations.
>
> I am running it in DC mode with no branch constraints, no reserve
> constraints, and the same load for all hours.
>
> The problem relates to the relative size of the VOLL and the cost of the
> generator units, which gives a very different dispatch from a slight change
> in VOLL despite no loss of load taking place.
>
> If I run the following simulation for 1 time period, I get what I believe
> to be the optimal solution, a dispatch of [12,12,6] however if I run it for
> 2 time periods, I get a sub-optimal solution for both periods [15,15,0]. I
> have found the a threshold of where the difference in 1 dollar changes
> whether or not this problem shows up, even though in the simulation no load
> is lost. You can see the difference using the code below and changing the
> VOLL value from 30014 (no problem) to 30015 (yields suboptimal solution).
>
> If someone could run this code and let me know whether this is a bug or a
> misunderstanding, that would be greatly appreciated!
>
> clear
> define_constants;
> VOLL=30014 %30015
> mpc.version = '2';
> %%-  Power Flow Data  -%%
> %% system MVA base
> mpc.baseMVA = 100;
>
> %% bus data
> %bus_itypePdQdGsBsareaVmVa
> baseKVzoneVmaxVmin
> mpc.bus = [
> 13000011013511.050.95;
> ];
>
> %% generator data
> %busPgQgQmaxQminVgmBasestatusPmax
> PminPc1Pc2Qc1minQc1maxQc2minQc2maxramp_agc
> ramp_10ramp_30ramp_qapf
> mpc.gen = [
> 10000  1100 1   20  0   0
> 00   0   0   00   000   0;
> 1000 0  1100 1   35  0   0
> 000000000   0;
> 10000  1100 1   25  0   0
> 000000000   0;
> 1000 0  1100 1   0   -300
> 000000000   0;
> ];
>
> mpc.branch = [];
>
> mpc.gencost = [
> 200310   0;
> 200310   0;
> 200320   0;
> 20030VOLL0;
> ];
>
> warning('off','all')
> mpopt = mpoption('verbose', 0);
> mpopt = mpoption(mpopt, 'most.dc_model', 0);
> loadprofile = struct( ...
> 'type', 'mpcData', ...
> 'table', CT_TLOAD, ...
> 'rows', 0, ...
> 'col', CT_LOAD_ALL_PQ, ...
> 'chgtype', CT_REP, ...
> 'values', [] );
>
> xgd_table.colnames = {
> 'CommitKey', ...
> 'CommitSched', ...
> };
> xgd_table.data = [
> 1   1;
> 1   1;
> 1   1;
> 2   1 ;
> ];
>
> loadprofile.values(:, 1, 1) = [30,30];
> profiles = getprofiles(loadprofile);
> nt = size(profiles(1).values, 1);   % number of periods
> xgd = loadxgendata(xgd_table, mpc);
>
> %This runs the dispatch twice: once for two hours, then for only one hour.
> for index = 1 : nt
> mpc = loadcase(mpc);
> remaining_nt = nt-index+1;
> profiles(1).values=profiles(1).values(index:end);
> mdi = loadmd(mpc, remaining_nt, xgd, [], [], profiles);
> mdo = most(mdi, mpopt);
> mdo.results.ExpectedDispatch
> end
>
>

Hydropower in MATPOWER MOST

[Stephen Suffian]

Ray and all,

I was wondering if anyone has experience in modeling hydro commitment in
MOST. I understand that a comprehensive hydropower commitment incorporates
medium and long-term solutions.

However, I am hoping to build a simplified model that possibly assumes
energy constraints on the hydro at the daily or weekly level, and then
looks to solve hourly commitment based on these constraints. Have other
people incorporated hydro in some way? Perhaps there is a way to consider
it a form of storage that can only be externally recharged (through
rainfall)? Do I develop a profile in order to have a time-dependent cost
function? Any help on the matter would be greatly appreciated!

Stephen Suffian

Re: Adding Wind Generation to Profiles

[Stephen Suffian]

Check out the 'help get_profile'. Much like you pass the previous profiles
object on the load profile line, you can pass the previous profiles object
on the other lines:

profiles = getprofiles('wind_profile1',iwind);
profiles = getprofiles('wind_profile2',profiles,iwind);
profiles = getprofiles('load_profile_30',profiles)

See if that works!

On Thu, Sep 21, 2017 at 7:24 PM, Joshua Sebben <joshkeep...@gmail.com>
wrote:

> Stephen,
>
> Thank you for your reply.
>
> casefile = 'case30_modified_MedPen';
> mpc = loadcase(casefile);
> xgd = loadxgendata('xgd_uc_30_MedPen', mpc);
> [iwind, mpc, xgd] = addwind('wind_uc_30_MedPen', mpc, xgd);
> profiles = getprofiles('wind_profile1', iwind);
> profiles = getprofiles('wind_profile2', iwind);
> profiles = getprofiles('wind_profile3', iwind);
> profiles = getprofiles('wind_profile4', iwind);
> profiles = getprofiles('wind_profile5', iwind);
> profiles = getprofiles('wind_profile6', iwind);
> profiles = getprofiles('wind_profile7', iwind);
> profiles = getprofiles('load_profile_30', profiles);
> nt = size(profiles(1).values, 1);   % number of periods
>
> I created the 7 wind profiles with each of the row counts changed to 1 to
> 7 for their respective generators.
>
> This would then change my main code assuming to the snippet above?
> This again also has the problem of only changing the last generator in the
> list as each of the previous wind profile changes are overwritten by the
> new getprofile?
>
> Regards,
> Josh
>
> On 22 September 2017 at 04:02, Stephen Suffian <ssuff...@villanova.edu>
> wrote:
>
>> I also couldn't figure it out, but I got around it by adding a separate
>> profile for each generator and changing the rows value. So if a wind
>> generator is in row 5, you would have the struct look like this below.
>>
>>
>> windprofile = struct( ...
>> 'type', 'mpcData', ...
>> 'table', CT_TGEN, ...
>> 'rows', 5, ...
>> 'col', PMAX, ...
>> 'chgtype', CT_REL, ...
>> 'values', [] );
>>
>> And then for row 6 you would add another profile with the struct looking
>> tlike this:
>>
>>
>> windprofile = struct( ...
>> 'type', 'mpcData', ...
>> 'table', CT_TGEN, ...
>> 'rows', 6, ...
>> 'col', PMAX, ...
>> 'chgtype', CT_REL, ...
>>     'values', [] );
>>
>> It is a bit of a work around, but I believe it should work (I did
>> something similar to set the Pmax for each conventional generator in a
>> profile, so I imagine it will work the same).
>>
>> On Thu, Sep 21, 2017 at 1:52 PM, Stephen Suffian <
>> stephen.suff...@gmail.com> wrote:
>>
>>> I also couldn't figure it out, but I got around it by adding a separate
>>> profile for each generator and changing the rows value. So if a wind
>>> generator is in row 5, you would have the struct look like this below.
>>>
>>>
>>> windprofile = struct( ...
>>> 'type', 'mpcData', ...
>>> 'table', CT_TGEN, ...
>>> 'rows', 5, ...
>>> 'col', PMAX, ...
>>> 'chgtype', CT_REL, ...
>>> 'values', [] );
>>>
>>> And then for row 6 you would add another profile with the struct looking
>>> tlike this:
>>>
>>>
>>> windprofile = struct( ...
>>> 'type', 'mpcData', ...
>>> 'table', CT_TGEN, ...
>>> 'rows', 6, ...
>>> 'col', PMAX, ...
>>> 'chgtype', CT_REL, ...
>>> 'values', [] );
>>>
>>> It is a bit of a work around, but I believe it should work (I did
>>> something similar to set the Pmax for each conventional generator in a
>>> profile, so I imagine it will work the same).
>>>
>>> On Thu, Sep 21, 2017 at 5:39 AM, Joshua Sebben <joshkeep...@gmail.com>
>>> wrote:
>>>
>>>> By the way,
>>>>
>>>>
>>>>
>>>> I have also tried setting the row count to [1 2 3 4 5 6 7] for my 7
>>>> extra generators that I want to add, however when I run my code I get an
>>>> error:
>>>>
>>>>
>>>>
>>>> Error using apply_profile (line 148)
>>>>
>>>> apply_profile: third dimension of profile.values should match length of
>>>> pro=
>>>>
>>>> file.rows
>>>>
>>>>
>>>>
>>>> Error in loadmd (line 508)
>>>>
>>>> optab =3D apply_profile(profiles(p), optab);
>>>>
>>>>
>>>>
>>>&g

Re: Adding Wind Generation to Profiles

[Stephen Suffian]

I also couldn't figure it out, but I got around it by adding a separate
profile for each generator and changing the rows value. So if a wind
generator is in row 5, you would have the struct look like this below.


windprofile = struct( ...
'type', 'mpcData', ...
'table', CT_TGEN, ...
'rows', 5, ...
'col', PMAX, ...
'chgtype', CT_REL, ...
'values', [] );

And then for row 6 you would add another profile with the struct looking
tlike this:


windprofile = struct( ...
'type', 'mpcData', ...
'table', CT_TGEN, ...
'rows', 6, ...
'col', PMAX, ...
'chgtype', CT_REL, ...
'values', [] );

It is a bit of a work around, but I believe it should work (I did something
similar to set the Pmax for each conventional generator in a profile, so I
imagine it will work the same).

On Thu, Sep 21, 2017 at 1:52 PM, Stephen Suffian <stephen.suff...@gmail.com>
wrote:

> I also couldn't figure it out, but I got around it by adding a separate
> profile for each generator and changing the rows value. So if a wind
> generator is in row 5, you would have the struct look like this below.
>
>
> windprofile = struct( ...
> 'type', 'mpcData', ...
> 'table', CT_TGEN, ...
> 'rows', 5, ...
> 'col', PMAX, ...
> 'chgtype', CT_REL, ...
> 'values', [] );
>
> And then for row 6 you would add another profile with the struct looking
> tlike this:
>
>
> windprofile = struct( ...
> 'type', 'mpcData', ...
> 'table', CT_TGEN, ...
> 'rows', 6, ...
> 'col', PMAX, ...
> 'chgtype', CT_REL, ...
> 'values', [] );
>
> It is a bit of a work around, but I believe it should work (I did
> something similar to set the Pmax for each conventional generator in a
> profile, so I imagine it will work the same).
>
> On Thu, Sep 21, 2017 at 5:39 AM, Joshua Sebben <joshkeep...@gmail.com>
> wrote:
>
>> By the way,
>>
>>
>>
>> I have also tried setting the row count to [1 2 3 4 5 6 7] for my 7 extra
>> generators that I want to add, however when I run my code I get an error:
>>
>>
>>
>> Error using apply_profile (line 148)
>>
>> apply_profile: third dimension of profile.values should match length of
>> pro=
>>
>> file.rows
>>
>>
>>
>> Error in loadmd (line 508)
>>
>> optab =3D apply_profile(profiles(p), optab);
>>
>>
>>
>> Error in Test (line 33)
>>
>> mdi =3D loadmd(mpc, transmat, xgd, [], [], profiles);
>>
>>
>> Regards,
>>
>> Josh
>>
>> On 20 September 2017 at 22:17, Joshua Sebben <joshkeep...@gmail.com>
>> wrote:
>>
>>> Currently working on adding wind generator units to my model. I am
>>> working off the 30 bus example case in MOST.  I am trying to run the
>>> following code
>>>
>>> casefile = 'case30';
>>> mpc = loadcase(casefile);
>>> xgd = loadxgendata('xgd_uc', mpc);
>>> [iwind, mpc, xgd] = addwind('wind_uc_30', mpc, xgd);
>>> profiles = getprofiles('wind_profile', iwind);
>>> profiles = getprofiles('load_profile', profiles);
>>> nt = size(profiles(1).values, 1);   % number of periods
>>>
>>> %%-  Full Transition Probabilities  -
>>> transmat = ex_transmat(nt);
>>> mdi = loadmd(mpc, transmat, xgd, [], [], profiles);
>>> mdo = most(mdi, mpopt);
>>> if verbose
>>> ms = most_summary(mdo);
>>> end
>>>
>>> However the added generator units of which there are 7 of them in
>>> wind_uc_30 don't seem to get the profile set to them in the output. Instead
>>> only the first Generator unit gets the profile while the rest of them are
>>> set at PMAX.
>>>
>>> windprofile = struct( ...
>>> 'type', 'mpcData', ...
>>> 'table', CT_TGEN, ...
>>> 'rows', 1, ...
>>> 'col', PMAX, ...
>>> 'chgtype', CT_REL, ...
>>> 'values', [] );
>>>
>>>  windprofile.values(:, :, 1) = [
>>>   0.80;
>>>   0.65;
>>>   0.60;
>>>   0.82;
>>>   1.00;
>>>   0.70;
>>>   0.50;
>>>   0.85;
>>>   1.00;
>>>   1.10;
>>>   1.06;
>>>   0.95;
>>>
>>> Above is a snippet from the wind_profile.  I am assuming it has
>>> something to do with the row count, however, I haven't been able to find a
>>> solution.
>>>
>>> Could I please get some help. to ensure all the generator units receive
>>> the load profile.
>>>
>>> Thankyou,
>>>
>>> --
>>> J.Sebben
>>>
>>
>>
>>
>> --
>> J.Sebben
>>
>
>

Spinning Reserves for UC

[Stephen Suffian]

Matpower List,

I am trying to run the determininstic unit commitment example (example 6)
but include a spinning reserve requirement. I am using the spinning reserve
example data that is used in example 1, but I get the following error:

Index exceeds matrix dimensions

Error in most (line 594)
  r = mdi.FixedReserves(t,j,k)


I imagine that means I need to incorporate a time component into the
reserve data. I tried a few different ways to update the matrices below but
was unable to get past this error. Is there a different way to incorporate
reserves in multiperiod problems?

Thanks!

%%-  Reserve Data  -%%
%% reserve zones, element i, j is 1 if gen j is in zone i, 0 otherwise
mpc.reserves.zones = [
1110;
];

%% reserve requirements for each zone in MW
mpc.reserves.req   = 150;

%% reserve costs in $/MW 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.cost  = [5;5;21;];
% mpc.reserves.cost  = [5;5;16.25;];
% mpc.reserves.cost  = [0;0;11.25;];
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;];

CPLEX Issues

[Stephen Suffian]

Dear All,

I've recently installed the academic version of the CPLEX studio (12.7.1),
and am able to run cplexlp and cplexqp commands from within Octave, but the
have_fcn('cplex') returns 0 and when I try to use it as an mpoption, I
receive the following error:

mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);   %% dual simplex

error: mpoption: 'cplex' is not a valid option name

Does anyone know if this is due to my improper installation of CPLEX or
whether there is something matpower related that is preventing me from
using CPLEX?

I looked in have_fcn() and it is searching for a mex file which doesn't
exist in my cplex folder. Instead, there are .m and .p files as well as a
.mexa64. But the which command used in have_fcn doesn't return mexa64
files.

Any help would be greatly appreciated. Thanks!
Stephen Suffian

MOST Question regarding Quadratic Gen Costs

[Stephen Suffian]

Hello all,

I have been able to succesfully run MOST for several days, but I recently
came across a problem when attempting to load cost data into mpc.gencost. I
My original gencost was below:

mpc.gencost = [
20022000;
20023000;
20025000;
200210000;
];

When I changed it to quadratic, it works as long as the first term is 0.

mpc.gencost = [
2003   02000;
2003   03000;
2003   05000;
2003   010000;
];

However, when I change the first term to a non-zero value:

mpc.gencost = [
2003   0.0042000;
2003   03000;
2003   05000;
2003   010000;
];

I get the following error.

error: mpopt2qpopt: Sorry, no solver available for MIQP models
error: called from
mpopt2qpopt at line 92 column 13
most at line 2065 column 14
test at line 41 column 5

I am using the glpk solver. Does anyone know if this is a solver issue, a
human error (on my part) issue, or otherwise?

Thanks!
Stephen Suffian