-
Notifications
You must be signed in to change notification settings - Fork 21
/
Copy pathrun_nBus.m
74 lines (54 loc) · 1.82 KB
/
run_nBus.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
% This script runs compares MPEC-AMPL to NFXP for various sample sizes
clear all
%add path to the functions
currdir=pwd;
addpath([currdir filesep 'matlab_include']);
addpath([currdir filesep 'nfxp_Rust87']);
addpath([currdir filesep 'nBus']);
addpath('..');
cd nBus
% paths for knitro and ampl
param.ampl_command='ampl';
%Need to have correct env variable for knitro to be able to run under ampl!!!
%setenv('DYLD_LIBRARY_PATH','/usr/local/knitro900/lib') %Fedor
setenv('DYLD_LIBRARY_PATH','/usr/local/knitro/lib') %Bertel
% sys switches
param.figure=0; %should figures be drawn
% parameters
param.MC=1; %number of MC simulated samples
param.nT = 120;
param.N = 175;
param.beta=0.9999;
param.RC=11.7257;
param.thetaCost=2.45569; %cost corresponding to John's grid size
param.thetaProbs=[0.0937 0.4475 0.4459 0.0127 0.0002]';
%MAIN LOOP over dimensions of the grid
param.out_suff='_nBus';
param.dontAMPL=0; % 1 to skip AMPL, only compute X0
param.dontNFXP=0; % 1 to skip AMPL, only compute X0
param.beta=0.9999;
nBus=(100:100:10000);
mvec=zeros(size(nBus));
for iBus=1:numel(nBus)
param.nBus = nBus(iBus);
%number of non-zeros in rows of transition probability matrix
param.M=numel(param.thetaProbs);
mvec(iBus)=5;
%solve the model
myEV=trueEV(param);
%simulate data (all MC samples)
[data.MC_dt,data.MC_xt,data.MC_dx]=simdata(param,myEV);
%run AMPL solver (and create X0)
ampl{iBus}=MPEC_AMPL(param,data);
%NFXP
if ~param.dontNFXP
nfxp{iBus}=run_nfxp_nBus(param,data,ampl{iBus}.X0);
end
% Result table (appended on each iteration)
% RTable(param,[dimvec;mvec],ampl,nfxp,1);
RTable(param,nBus,ampl,nfxp,0); %no plot
save('Results_nBus','ampl','nBus', 'mvec', 'param', 'nfxp');
end % end of loop over sample size
cd ..
for iBus=1:numel(nBus); disp([ampl{iBus}.runtime ampl{iBus}.runtime_matlab ]); end
plot_cpu_vs_samplesize