Hi all,
Background:
I am working on a project that involves simulating a system using MATPOWER's
MOST function. The system consists of a single bus with a PV generator, a load,
and a battery storage. The objective is to optimize the operation of the
battery based on the PV generation profile and the load.
System Configuration:
In my case setup, the distributed load is represented as a generator with
negative power values. To enforce a strict demand, I've set the Pg, PMAX, and
PMIN parameters of this "load generator" to be equal. This ensures that the
load demand remains constant and is not subjected to optimization. I am running
dcopf model:
mpopt = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'opf.dc.solver', 'MIPS');
mpopt = mpoption(mpopt, 'most.solver', 'DEFAULT');
Problem:
I've encountered a peculiar issue with the MOST optimization. Whenever I set
the PV generation (`PVg`) to be less than half of the load (`0.5*load`), the
optimization fails, and many of the variables from 'most.results' are missing.
This happens despite the battery storage and discharge rate being set high
enough to provide the excess power being demanded. This behavior persists
regardless of the specific values I use for the load and PVg. The error is
consistent: as long as PVg is less than half of the load, the error occurs.
Specifically, the error I receive is:
------------------------------------------------
Unrecognized field name "Rrp".
Error in most_summary (line 62)
Rup = mdo.results.Rrp;
------------------------------------------------
Furthermore, I've observed that the error does not occur when the PV generation
is equal to or greater than half of the load. I suspect it might be related to
some constraint that prevents the battery generation from exceeding the PV
generation, but I haven't been able to pinpoint the exact cause.
Question:
Has anyone encountered a similar issue or can provide insights into why this
might be happening? Are there specific constraints or parameters in MOST that
could be causing this behavior? Any guidance or suggestions would be greatly
appreciated.
Thank you in advance for your assistance!
I have pasted my code below:
1. Main program
2. Case_file
3. PV generation profile funtion ("wind" just refers to "PV". Couldnt be
bothered changing all of the native functions from "wind" to "solar")
4. storage data function
1)% Main program to run DCOPF for no line connection.
clear all;
close all;
%% set up options
verbose = 3;
mpopt = mpoption('verbose', verbose);
mpopt = mpoption(mpopt, 'out.all', 1); % Print all results
mpopt = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'opf.dc.solver', 'MIPS');
mpopt = mpoption(mpopt, 'most.solver', 'DEFAULT');
mpopt = mpoption(mpopt, 'out.lim.all', 2);
mpopt = mpoption(mpopt, 'out.force', 1);
mpopt = mpoption(mpopt, 'out.suppress_detail', 0);
mpopt = mpoption(mpopt, 'most.security_constraints', 0);
%% Load case study
casefile = 'casefarm_noline';
mpc = loadcase(casefile);
%% XGEN data
xgd_table.colnames = {
'CommitKey', ...
'CommitSched', ...
};
xgd_table.data = [
2 1;
];
xgd = loadxgendata(xgd_table, mpc);
%% ----- + Solar -----
% 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
wind.gen = [
1 0 0 50 -50 1 1 1 46 -2 0 0 0 0 0 0 20 20 20 0 0; % Keep ramp_10 & ramp_30
large toensure initial demand deifned in profile can be fully supplied
];
% generator cost data
wind.gencost = [
2 0 0 2 1 0; % PV Load
];
% xGenData
wind.xgd_table.colnames = {
'CommitKey', 'CommitSched'
};
wind.xgd_table.data = [
2 1;
];
[iwind, mpc, xgd] = addwind(wind, mpc, xgd);
profiles = getprofiles('wind_profile_farm_noline', iwind);
nt = size(profiles(1).values, 1); % number of periods
%% ----- + storage -----
% mpopt = mpoption(mpopt, 'most.storage.cyclic', 0);
[iess, mpc, xgd, sd] = addstorage('storage_farm_noline', mpc, xgd);
%% simulation
mdi = loadmd(mpc, nt, xgd, sd, [], profiles);
mdo = most(mdi, mpopt); % Run simulation
ms = most_summary(mdo); % print Summary of results
2)
function mpc = casefarm_noline
%% MATPOWER Case Format : Version 2
mpc.version = '2';
%%----- Power Flow Data -----%%
%% system MVA base
mpc.baseMVA = 1;
baseKV = 22;
baseMVA = mpc.baseMVA;
%% 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 baseKV 1 1.1 0.9;
];
%% generator data
Pmax = -40; % convert kw value to MW
Pmax = Pmax / 1000; % convert kw value to MW
Pmin = Pmax;
%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 Pmax 0 50 -50 1 1 1 Pmax Pmin 0 0 0 0 0 0 0 60 60 0 0; % load. Keep
PMAX=Pmin, so it acts as a strictly defined Pd.
];
%% branch data (none for single bus)
mpc.branch = [];
%%----- OPF Data -----%%
%% generator cost data
mpc.gencost = [
2 0 0 2 1 0; % Slack Load
];
function windprofile = wind_profile_farm_noline
%WIND_PROFILE_FARM_NOLINE Example wind profile data
%% define constants
[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;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
CT_MODCOST_X] = idx_ct;
%------- Values are in MW --------
windprofile = struct( ...
'type', 'mpcData', ...
'table', CT_TGEN, ...
'rows', 0, ...
'col', PMAX, ...
'chgtype', CT_REP, ...
'values', [] );
windprofile.values(:, :, 1) = [
550;
450;
300;
19; % If load = 20kw and i set this to 19kw, then the dcopf wont converge
on 4th time period
]/1000; % /1000 to convert kW values to MW
function storage = storage_farm_noline(mpc)
%%----- storage -----
ecap = 15; % energy capacity (MWh)
pcap = ecap*0.25; % power capacity (MW). Max charge/discharge rate
scost = 1; % cost/value of initial/residual stored energy
SOC_IC = 0.5*ecap; % initial storage
%% 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
storage.gen = [
1 0 0 0 0 1 1 1 pcap -pcap 0 0 0 0 0 0 0 20 20 0 0;
]; % 4th and 3rd last columns specify the charge/discharge rate
%%----- OPF Data -----%%
%% generator cost data
% 1 startup shutdown n x1 y1 ... xn yn
% 2 startup shutdown n c(n-1) ... c0
storage.gencost = [
2 0 0 2 1 0;
];
%% xGenData
storage.xgd_table.colnames = {
'CommitKey', ...
'CommitSched', ...
};
storage.xgd_table.data = [
2 1;
];
%% StorageData
storage.sd_table.OutEff = 1;
storage.sd_table.InEff = 1;
storage.sd_table.LossFactor = 0;
storage.sd_table.rho = 1;
storage.sd_table.colnames = {
'InitialStorage', ... % MWh
'InitialStorageLowerBound', ...
'InitialStorageUpperBound', ...
'InitialStorageCost', ...
'TerminalStoragePrice', ...
'MinStorageLevel', ... % MWh
'MaxStorageLevel', ... % MWh
'OutEff', ...
'InEff', ...
'LossFactor', ...
'rho', ...
};
Eff = 1; % battery efficiency (typically between 0.85 to 0.95)
storage.sd_table.data = [
SOC_IC 7 ecap scost scost 0 ecap Eff Eff 0 1;
];
