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 <[email protected]> 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_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; > ]; > > %% 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 0 0 0 0 1 100 1 20 0 0 > 0 0 0 0 0 0 0 0 0 0; > 1 0 0 0 0 1 100 1 35 0 0 > 0 0 0 0 0 0 0 0 0 0; > 1 0 0 0 0 1 100 1 25 0 0 > 0 0 0 0 0 0 0 0 0 0; > 1 0 0 0 0 1 100 1 0 -30 0 > 0 0 0 0 0 0 0 0 0 0; > ]; > > mpc.branch = []; > > mpc.gencost = [ > 2 0 0 3 1 0 0; > 2 0 0 3 1 0 0; > 2 0 0 3 2 0 0; > 2 0 0 3 0 VOLL 0; > ]; > > 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 > >
