Multiperiod deterministic problem, UC does not change

[yadong-zhang]

I solved the problem in https://github.com/MATPOWER/most/issues/33, but encountered new problems, can you please give a hint?

1. UC doesn’t change when “RAMP_30" for thermal gens is set less than half of the “PMAX”, which actually manifests no ramp constraints (Line 84 -- 90);
2. UC doesn’t change when thermal gen “StartUp" and “ShutDown" cost is larger than zero, even very small (Line 108 -- 120);
3. UC doesn't change when adding reserve in the computation, whatever reserve requirement is set (Line 142 -- 151);
4. The “gencost” of dispatchable loads seems to impact SCUC computation, small value (e.g., 10) will lead to unreasonable results (e.g., load is not fully satisfied even if there is still generation capacity). But everything is correct when I use large value (e.g., 1e6).

Thanks a lot for your kind help!!

function [mpc, xgd_table, wind_UC] = GetFiles()

% Generate necessary files for running MOST

%% Define constants
define_constants;

%% Read bus and zone info
clc;

% Read bus info
load_buses = readmatrix('zones/load_bus.csv');
gen_buses = readmatrix('zones/gen_bus.csv');
wind_buses = readmatrix('zones/wind_bus.csv');
thermal_buses = readmatrix('zones/thermal_buses.csv');

zone1_buses = readmatrix('zones/zone1_bus.csv');
zone2_buses = readmatrix('zones/zone2_bus.csv');
zone3_buses = readmatrix('zones/zone3_bus.csv');

zone1_loads = readmatrix('zones/zone1_load_bus.csv');
zone2_loads = readmatrix('zones/zone2_load_bus.csv');
zone3_loads = readmatrix('zones/zone3_load_bus.csv');

zone1_gens = readmatrix('zones/zone1_gen_bus.csv');
zone2_gens = readmatrix('zones/zone2_gen_bus.csv');
zone3_gens = readmatrix('zones/zone3_gen_bus.csv');

zone1_winds = readmatrix('zones/zone1_wind_bus.csv');
zone2_winds = readmatrix('zones/zone2_wind_bus.csv');
zone3_winds = readmatrix('zones/zone3_wind_bus.csv');

zone1_thermals = readmatrix('zones/zone1_thermal_buses.csv');
zone2_thermals = readmatrix('zones/zone2_thermal_buses.csv');
zone3_thermals = readmatrix('zones/zone3_thermal_buses.csv');

% Set numbers
num_zones = 3;
num_buses = 118;
num_gens = size(gen_buses, 1);
num_loads = size(load_buses, 1);
num_winds = size(wind_buses, 1);
num_thermals = size(thermal_buses, 1);

num_zone1_buses = size(zone1_buses, 1);
num_zone2_buses = size(zone2_buses, 1);
num_zone3_buses = size(zone3_buses, 1);

num_zone1_laods = size(zone1_loads, 1);
num_zone2_laods = size(zone2_loads, 1);
num_zone3_laods = size(zone3_loads, 1);

num_zone1_gens = size(zone1_gens, 1);
num_zone2_gens = size(zone2_gens, 1);
num_zone3_gens = size(zone3_gens, 1);

num_zone1_winds = size(zone1_winds, 1);
num_zone2_winds = size(zone2_winds, 1);
num_zone3_winds = size(zone3_winds, 1);

num_zone1_thermals = size(zone1_thermals, 1);
num_zone2_thermals = size(zone2_thermals, 1);
num_zone3_thermals = size(zone3_thermals, 1);

%% Set power grid case file
clc;

% Read casefile
casefile = PowerGrid(num_zones);     
mpc = loadcase(casefile);

% Set mpc.bus 
mpc.bus(load_buses, [PD QD]) = 0;  

% Set mpc.gen
mpc.gen(1:num_gens, [PG QG QMAX QMIN ...
                     PC1 PC2 ...
                     QC1MIN QC1MAX ...
                     QC2MIN QC2MAX ...
                     RAMP_Q APF]) = 0;     % Set unnecessary values as 0

mpc.gen(1:num_gens, PMAX) = 200;           %%%%%%%%% Max gen capacity
mpc.gen(1:num_gens, PMIN) = 0.1;           %%%%%%%%% Min gen capacity (cannot be zero)
mpc.gen(1:num_gens, RAMP_10) = 0; 
mpc.gen(1:num_gens, RAMP_AGC) = 0; 
mpc.gen(1:num_gens, RAMP_30) = 1e6;         %%%%%%%%% Ramp up/dowm limit within 30 mins
% Dispatchable loads
mpc.gen(num_gens+1:end, [PG QG QMAX QMIN ...
                         PC1 PC2 QC1MIN ...
                         QC1MAX QC2MIN ...
                         QC2MAX RAMP_Q APF]) = 0;
mpc.gen(num_gens+1:end, VG) = 1;           % Set Vg
mpc.gen(num_gens+1:end, 7) = 100;          % Set mBase = 100
mpc.gen(num_gens+1:end, 8) = 1;            % Set status = 1
mpc.gen(num_gens+1:end, PMAX) = 0;         %%%%%%%%% Must be zero
mpc.gen(num_gens+1:end, PMIN) = -1000;        %%%%%%%%% Must be < 0
mpc.gen(num_gens+1:end, RAMP_AGC) = 0;
mpc.gen(num_gens+1:end, RAMP_10) = 0;
mpc.gen(num_gens+1:end, RAMP_30) = 1e6;    % Set unlimited ramp up/down limit

% Set mpc.branch

% Set mpc.gencost
% Gens
mpc.gencost(1:num_gens, 1) = 2;             % Set cost model
mpc.gencost(1:num_gens, 2) = 0;          %%%%%%%%%%%%%%%%%%%%%% StartUp cost, cannot be too big
mpc.gencost(1:num_gens, 3) = 0;          %%%%%%%%%%%%%%%%%%%%%% ShutDown cost, cannot be too big
mpc.gencost(1:num_gens, 4) = 2;             % Set num of params
mpc.gencost(1:num_gens, 5) = 1.5;           % Set cost coeff
mpc.gencost(1:num_gens, 6) = 0;
% Dispatchable loads
mpc.gencost(num_gens+1:end, 1) = 2;
mpc.gencost(num_gens+1:end, [2 3]) = 0;     % Set StartUp and StutDown cost = 0
mpc.gencost(num_gens+1:end, 4) = 2;         % Set number of parameters
mpc.gencost(num_gens+1:end, 5) = 1e6;       %%%%%%%%%%%%%%%%%%%%%% Set as +Inf
mpc.gencost(num_gens+1:end, 6) = 0;

% Set mpc.reserves (wind gens have no reserve)
temp1 = zeros(1, num_gens);
bidx = ismember(gen_buses, zone1_thermals);
temp1(bidx) = 1;

temp2 = zeros(1, num_gens);
bidx = ismember(gen_buses, zone2_thermals);
temp2(bidx) = 1;

temp3 = zeros(1, num_gens);
bidx = ismember(gen_buses, zone3_thermals);
temp3(bidx) = 1;

temp4 = zeros(1, num_loads);    % Dispatchable loads

mpc.reserves.zones(1, :) = cat(2, temp1, temp4);
mpc.reserves.zones(2, :) = cat(2, temp2, temp4);
mpc.reserves.zones(3, :) = cat(2, temp3, temp4);

% reserve requirements for each zone in MW
mpc.reserves.req = [10; 50; 10];             % Should be consistent with "num_zones"

% 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 = 0*ones(num_thermals, 1);    % Skip winds and loads

% 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 = 200*ones(num_thermals, 1);

%% Set XGD file
clc;

% Read xgd table
xgd_table = XGDTable;

% Gens
xgd_table.data(1:num_gens, 1) = 1;          % Set thermal gens UC = commitable
xgd_table.data(1:num_gens, 2) = 1;      
xgd_table.data(1:num_gens, [3 4]) = 2;      %%%%%%%%%%%%%%%%%%%%%% Set thermal gens MinUp and MinDown
xgd_table.data(1:num_gens, 5) = 0;          % PositiveActiveReservePrice
xgd_table.data(1:num_gens, 6) = 0;          % PositiveActiveReserveQuantity
xgd_table.data(1:num_gens, 7) = 0;          % NegativeActiveReservePrice
xgd_table.data(1:num_gens, 8) = 0;          % NegativeActiveReserveQuantity
xgd_table.data(1:num_gens, 9) = 0;          % PositiveActiveDeltaPrice
xgd_table.data(1:num_gens, 10) = 0;         % NegativeActiveDeltaPrice
xgd_table.data(1:num_gens, 11) = 0;         % PositiveLoadFollowReservePrice
xgd_table.data(1:num_gens, 12) = 1e6;       % PositiveLoadFollowReserveQuantity (already specified by RAMP_30)
xgd_table.data(1:num_gens, 13) = 0;         % NegativeLoadFollowReservePrice
xgd_table.data(1:num_gens, 14) = 1e6;       % NegativeLoadFollowReserveQuantity (already specified by RAMP_30)

%%%%%%%%%%%% Set UC = off for thermals at winds %%%%%%%%%%%%%%%%%%
temp1 = ismember(gen_buses, wind_buses);    
temp2 = zeros(num_loads, 1);
bidx = logical(cat(1, temp1, temp2));
xgd_table.data(bidx, 1) = -1;

% Dispatchable loads
xgd_table.data(num_gens+1:end, 1) = 2;      % Set dispatchable loads UC = 2
xgd_table.data(num_gens+1:end, 2) = 1;
xgd_table.data(num_gens+1:end, [3 4]) = 1;  % Set dispatchable loads MinUp/MinDown
xgd_table.data(num_gens+1:end, 5) = 0;      % PositiveActiveReservePrice
xgd_table.data(num_gens+1:end, 6) = 0;      % PositiveActiveReserveQuantity
xgd_table.data(num_gens+1:end, 7) = 0;      % NegativeActiveReservePrice
xgd_table.data(num_gens+1:end, 8) = 0;      % NegativeActiveReserveQuantity
xgd_table.data(num_gens+1:end, 9) = 0;      % PositiveActiveDeltaPrice
xgd_table.data(num_gens+1:end, 10) = 0;     % NegativeActiveDeltaPrice
xgd_table.data(num_gens+1:end, 11) = 0;     % PositiveLoadFollowReservePrice
xgd_table.data(num_gens+1:end, 12) = 1e6;   % PositiveLoadFollowReserveQuantity
xgd_table.data(num_gens+1:end, 13) = 0;     % NegativeLoadFollowReservePrice (set as unlimited)
xgd_table.data(num_gens+1:end, 14) = 1e6;   % NegativeLoadFollowReserveQuantity (set as unlimited)

%% Set WindUC file
clc;

% Read wind UC file
wind_UC = WindUC;

% Gens
wind_UC.gen(:, PMAX) = 100;         %%%%%%%%%%%%%%%%%%% This value will be replaced by stochastic winds
wind_UC.gen(:, PMIN) = 0; 
wind_UC.gen(:, RAMP_AGC) = 0;       
wind_UC.gen(:, RAMP_10) = 0;
wind_UC.gen(:, RAMP_30) = 1e6;      % Set unlimited ramp rate
% xgd_table
wind_UC.xgd_table.data(:, 1) = 2;   %%%%%%%%%%%%%%%%%%% Set UC = on all the time
wind_UC.xgd_table.data(:, 2) = 1;
wind_UC.xgd_table.data(:, 3) = 0;   
wind_UC.xgd_table.data(:, 4) = 0;   % Set cost = 0 to make sure wind is always preferably used
wind_UC.xgd_table.data(:, 5:end) = 0;  % Set all = 0
wind_UC.xgd_table.data(:, 12) = 1e6;   % Set unlimited reserve and load follow capacity
wind_UC.xgd_table.data(:, 14) = 1e6;   % Set unlimited reserve and load follow capacity

end

[rdzman]

1. Are any of the ramps from one period to the next in the original solution exceeding twice the RAMP_30 values you specify? If not, it simply means that none of the ramp constraints are binding. Also, a binding ramp constraint will affect the dispatch, but may not necessarily modify the UC.
2. Does your original UC solution have any units with startup and shutdown events? If so, if you add a large enough startup/shutdown cost on one of these units, it should force the unit to be off/on for the whole planning horizon, assuming that is feasible.
3. I assume you are talking about fixed zonal reserves. If you have any units that are offline in the original commitment, it should be possible to increase a reserve requirement enough to get them committed. Just make sure that the units of interest are included in the zone whose reserve requirement you are increasing.
4. A small "cost" on a dispatchable load means that there is very little value to serving it, or that it can be curtailed at a low price. If the LMP at the bus is higher, it will be curtailed.

[yadong-zhang]

1. When RAMP_30 is set as unlimited (e.g., 1e6), UC will change over time. However, when setting RAMP_30 as <= half of PMAX, UC doesn't change. By UC does not change, I mean UC for all generators are always 1 for all time steps. For the former, there are cases with significant ramp up or down, while for the later, all gens are turned on and many of them are just working at PMIN, which is a very small value, e.g., 0.1.
2. Yes, there are units with startup and shutdown events. But I tried tiny cost (e.g., 1e-4) for startup and shutdown, it turned out that UC was always on for all generators for all time steps. However, everything was correct when I sued a smaller value, e.g., 1e-6. Perhaps I should have mentioned (which is also in the code), that I use the same price for all generators.
3. Yes, that is guaranteed.
4. Does that mean I should always use high cost value for dispatchable loads, e.g., 1e3, 1e6? I also notice that it is related to the cost of generators, for example, if cost for gens is 0.01 then 10 is fine for dispatchable loads, but if it is 10 for the former then the later requires >1e3.

[rdzman]

Ok, I misunderstood what you meant by "UC does not change". I thought you meant the entire commitment schedule over the horizon was identical for the run with ramp limits and without. But you mean that all generators are on for all time periods. If you expect to see changing commitment statuses (i.e. startups and shutdowns) as you move through the planning horizon, you will need a changing load (if dispatchable, with a "cost" that is higher than the generation units being committed, so it's not more economical to curtail load than commit new units), variation in costs between generators, so the expensive ones are only used in peak hours.

1. So yes, you would expect ramp constraints to make shutdown/startup events less likely, since they generally increase the ramping requirements of other units.
2. Adding startup and shutdown costs will also make those events less likely ... i.e. more likely to see units staying on for all periods.
4. Unless you are modeling true dispatchable loads (in which case you use their actual cost of curtailment), yes, you would use a high value, something that represents the value of lost load. Think of the curtailment of a dispatchable load as if it were a generator competing with other generators. The lower cost ones will be selected first by the optimization.

[yadong-zhang]

I see! I should definitely try PMAX and gencost with some variation instead of the same. Thanks a lot for the help!