clear all
start_v = ones(5,1);
% read in the raw data for the network planning model
small5data

Xs.name = 'X';
Xs.ts   = 'all branches';
Xs.val  = [ linspace(1,6,6) ]';
Xs.type = 'set';
Xs.form = 'sparse';
Xs.uels = {'L1', 'L2','L3','L4','L5','L6'};

Ks.name = 'K';
Ks.ts   = 'new branches';
Ks.val  = [ linspace(3,3,6) ]';
Ks.type = 'set';
Ks.form = 'sparse';
Ks.uels = {'L4','L5','L6'};

Ys.name = 'Y';
Ys.ts   = 'all nodes';
Ys.val  = [ linspace(1,nN,nN) ]';
Ys.type = 'set';
Ys.form = 'sparse';
Ys.uels = {'N1', 'N2','N3','N4','N5'};

ALs.name = 'AL';
ALs.ts   = 'alternatives in network';
ALs.val  = [ linspace(1,kN,kN) ]';
ALs.type = 'set';
ALs.form = 'sparse';
ALs.uels = {'AL1','AL2','AL3','AL4'};

loads.name = 'load';
loads.ts   = 'mean active power on each node';
loads.type = 'parameter';
loads.val = current;
loads.form = 'full';
loads.dim = 1;
loads.uels = Ys.uels;

ress.name = 'res';
ress.ts   = 'resistance of all lines';
ress.type = 'parameter';
ress.val = Line_res;
ress.form = 'full';
ress.dim = 2;
ress.uels{2} = Xs.uels;
ress.uels{1} = ALs.uels;   % the row

curts.name = 'current_max';
curts.ts   = 'capacities of all lines';
curts.type = 'parameter';
curts.val = Line_cap;
curts.form = 'full';
curts.dim = 2;
curts.uels{2} = Xs.uels;
curts.uels{1} = ALs.uels;

capexs.name = 'c_cost';
capexs.ts   = 'capex of the lines';
capexs.type = 'parameter';
capexs.val = Line_capex;
capexs.form = 'full';
capexs.dim = 2;
capexs.uels{2} = Xs.uels;
capexs.uels{1} = ALs.uels;

opexs.name = 'o_cost';
opexs.ts   = 'opex of the lines';
opexs.type = 'parameter';
opexs.val = Line_opex;
opexs.form = 'full';
opexs.dim = 2;
opexs.uels{2} = Xs.uels;
opexs.uels{1} = ALs.uels;

incids.name = 'sf';
incids.ts   = 'node-branch incidence';
incids.type = 'parameter';
incids.val = -S;
incids.form = 'full';
incids.dim = 2;
incids.uels{2} = Xs.uels;
incids.uels{1} = Ys.uels;    

volls.name = 'voltage_min';
volls.ts   = 'the lower voltage boundary';
volls.type = 'parameter';
volls.val = Vmin';
volls.form = 'full';
volls.dim = 1;
volls.uels = Ys.uels;

volus.name = 'voltage_max';
volus.ts   = 'the upper voltage boundary';
volus.type = 'parameter';
volus.val = Vmax';
volus.form = 'full';
volus.dim = 1;
volus.uels = Ys.uels;

ms.name = 'M';
ms.ts   = 'a very big number';
ms.type = 'parameter';
ms.val = M;
ms.form = 'full';
ms.dim = 0;

Ibs.name = 'Ibase';
Ibs.ts   = 'base value for current';
Ibs.type = 'parameter';
Ibs.val = Ibase;
Ibs.form = 'full';
Ibs.dim = 0;

wgdx ('smalldat5', Xs, Ys, Ks, ALs,loads,ress,capexs, opexs,curts,incids,volls,volus,ms,Ibs);

gams_output = 'std';
gams('smallmat5');
solGDX = 'small5sol.gdx';

% now get the outputs from the GDX file produced by the GAMS run
os = struct('name','modelStat','form','sparse');
o = rgdx (solGDX, os);
modelstat = o.val;

os.name = 'solveStat';
o = rgdx (solGDX, os);
solvestat = o.val;

os.name = 'c_total';
os.field = 'l';
o = rgdx (solGDX, os);
c_total = o.val;

os.name = 'Pim';
os.field = 'l';
o = rgdx (solGDX, os);
Pim = o.val;

os.name = 'I';
o = rgdx (solGDX, os);
I = o.val;

os.name = 'V';
o = rgdx (solGDX, os);
V = o.val;

oos = struct('name','modelStat','form','full');

oos.name = 'dec_l';
oos.uels{1} = Xs.uels;
oos.uels{2} = ALs.uels;
o = rgdx (solGDX, oos);
dec_l = o.val

fprintf('modelstat = %d (1 = optimal solution)\n', modelstat);
fprintf('solvestat = %d (1 = normal completion)\n\n', solvestat);
fprintf('optimal value = %g\n\n', c_total);
%fprintf('optimal decision = %g\n\n', dec_l);




