error regarding wgdx and rgdx

Hi, please someone tell me regarding integration of MATLAB code with GAMS, i want to implement DCOPF. I receive this error in MATLAB:Error using wgdx:Amount of full format .val data exceeds number of UELs. Please try to resolve it, is my wgdx and rgdx functions working?
This is my MATLAB code:
clear all
close all
clc
%inputs
mpc=load (‘grid_data_24bus.mat’);
profile_L=mpc.WD(:,3);
profile_gen_rand=[1.0057;1.027;0.98;0.967;0.9837;0.971;1.001;1.04063;1.0129;0.961;0.989;0.9554;1.000129;0.9937;1.0497;1.0316;0.9986;1.03944;0.963;0.989;1.0427;1.041749;1.0214;1.0118];
n_time = 24; % time intervals
DeltaT = 24/n_time; % granularity
n_bus = size(mpc.BusData,1); % Number of buses
n_brnc = size(mpc.branches,1); % Number of branches
n_gen = size(mpc.Gen,1); % num of generators
slack = 13; % slack choice
VS = 1; % default value for slack voltage
V_up=1.1; % Voltage limits
V_lo=0.9; % Voltage limits
% Time profiles
mpc.BusData = repmat (mpc.BusData, [1,1,n_time]);
mpc.Gen = repmat (mpc.Gen, [1,1,n_time]);
for t=1:n_time
mpc.BusData(:,2,t)=profile_L(t).*mpc.BusData(:,2,t);% column 2 have Pd
for g=1:size(mpc.Gen,1)
mpc.Gen(g,3:4,t)=profile_gen_rand(t).*( mpc.Gen(g,3:4,t));%pmax and pmin
end
end
mpc_pre=mpc; %save mpcstruct
%%% ESS allocation
SOC_max = 100; % battery size [MWh]
SOC_0 = 0.2*SOC_max/mpc.Sbase; % initial SoC [per unit]
eta_c = 0.95; % charge efficiency
eta_d = 0.9; % discharge efficiency
VOLL = 10000; % Value of Loss of Load
VWC = 50; % Value of loss of wind ($/MWh)
% Define Wind Capacity and Generation Profiles for Specific Buses
% Define Wind Capacity and Generation Profiles for Specific Buses
Wcap = zeros(n_bus, 1); % Wind capacity array for all buses
Wcap(8) = 200; % Bus 8: 200 MW
Wcap(19) = 150; % Bus 19: 150 MW
Wcap(21) = 100; % Bus 21: 100 MW
% Wind data over 24 hours (from ‘w’ column)
wind_profile = mpc.WD(:, 2); % 24-hour wind profile
% Wind Generation Matrix (Buses x Time)
wind_gen = Wcap * wind_profile’; % Corrected dimensions

% sending data to GAMS
% Time index for GAMS file
TT=cell(1,n_time);
for tt=1:n_time
TT{tt}=[‘t’,num2str(tt)];
end
y.name=’t’;
y.type=’set’;
y.uels=TT;

% Nodes
bus = linspace(1,n_bus,n_bus)’;
Bus.name=’bus’;
Bus.type=’set’;
Bus.uels=num2cell(bus)’;

% Slack bus
Slack.name=’slack’;
Slack.type=’set’;
Slack.val=slack;
Slack.uels=Bus.uels;

% Generators
GG=cell(1,n_gen);
for gg=1:n_gen
GG{gg}=[‘g’,num2str(gg)];
end
gen.name=’Gen’;
gen.type=’set’;
gen.uels=GG;

% S_base
S_base.name=’sbase’;
S_base.type=’parameter’;
S_base.val=mpc.Sbase;
S_base.form=’full’;
S_base.dim=0;

% VOLL
VoLL.name=’VOLL’;
VoLL.type=’parameter’;
VoLL.val=VOLL;
VoLL.form=’full’;
VoLL.dim=0;

% Value wind curtailment
Vwc.name=’VWC’;
Vwc.type=’parameter’;
Vwc.val=VWC;
Vwc.form=’full’;
Vwc.dim=0;

% Charge efficiency []
Eta_c.name=’eta_c’;
Eta_c.type=’parameter’;
Eta_c.val=eta_c;
Eta_c.form=’full’;
Eta_c.dim=0;

% Discharge efficiency []
Eta_d.name=’eta_d’;
Eta_d.type=’parameter’;
Eta_d.val=eta_d;
Eta_d.form=’full’;
Eta_d.dim=0;

% Maximum battery SoC [MWh]
socmax.name=’SOCmax’;
socmax.type=’parameter’;
socmax.val=SOC_max;
socmax.form=’full’;
socmax.dim=0;

% Initial battery SoC [0-1]
socini.name=’SOC0′;
socini.type=’parameter’;
socini.val=SOC_0;
socini.form=’full’;
socini.dim=0;

% Wind and load power profiles []
%WD(t,*) wind-demand daily variation patterns: ‘w’ wind ‘d’ demand
wd.name=’WD’;
wd.type=’parameter’;
wd.val=mpc.WD(:,2:3);
wd.form=’full’; %missing indices are zero
wd.dim=2;
wd.uels={TT,{‘w’,’d’}};

% Branches data
%branches(bus,node,*) ‘r(pu)’ ‘x(pu)’ ‘b(pu)’ ‘limit(MVA)’
branches = [];
for ii=1:size(branches,1)
branches=vertcat(branches, zeros(4,4));
branches((ii-1)*4+1:ii*4,:)=….
[ branches(ii,1) branches(ii,2) 1 branches(ii,3); … from to ‘r’ value
branches(ii,1) branches(ii,2) 2 branches(ii,4); … from to ‘x’ value
branches(ii,1) branches(ii,2) 3 branches(ii,5); … from to ‘b’ value
branches(ii,1) branches(ii,2) 4 branches(ii,6) ]; … from to ‘limit’ value
end
Branches.name=’branches’;
Branches.type=’parameter’;
Branches.val=mpc.Branches;
Branches.form=’sparse’; %missing indices are zero
Branches.uels={Bus.uels, Bus.uels, {‘r’,’x’,’b’,’limit’}};
Branches.dim=3;

% Gen bus location
gb.name=’GB’;
gb.type=’set’;
gb.val=mpc.Gen(:,1:2); %GB(Gen,bus) generating buses
gb.val(:,1) = 1:n_gen;
gb.form=’sparse’;
gb.dim=2;
gb.uels={gen.uels,{1:24}}; % Covering all 24 buses in the system

% Gen data
GenData=[mpc.Gen(:,5)’; mpc.Gen(:,4)’; mpc.Gen(:,3)’; mpc.Gen(:,8)’; mpc.Gen(:,9)’]’; %GenData GD(Gen,*) ‘b’ ‘pmin’ ‘pmax’ ‘RU’ ‘RD’ -> ‘b’ t.c. Gen_cost = a*Pg^2 + b*Pg + c
Gen_Data.name=’GD’;
Gen_Data.type=’parameter’;
Gen_Data.val=GenData;
Gen_Data.form=’full’;
Gen_Data.dim=2;
Gen_Data.uels={GG,{‘b’,’pmin’,’pmax’,’RU’,’RD’}};
% Define the WCap parameter
WCap.name = ‘Wcap’; % Parameter name
WCap.type = ‘parameter’; % Define it as a parameter
WCap.val = Wcap; % Assign wind capacity values from Wcap
WCap.form = ‘full’; % Missing indices are assumed zero
WCap.dim = 1; % One-dimensional parameter (indexed by buses)
WCap.uels = {Bus.uels}; % Link the parameter to the bus indices
% Bus data
Bus_Data.name=’BusData’;
Bus_Data.type=’parameter’;
Bus_Data.val=[mpc.BusData(:,1) ones(size(mpc.BusData,1),1) mpc.BusData(:,2) mpc.BusData(:,1) 2*ones(size(mpc.BusData,1),1) mpc.BusData(:,3)];
Bus_Data.val=reshape(Bus_Data.val’, [size(Bus_Data.val’,1)/2,size(Bus_Data.val’,2)*2]);
Bus_Data.val=Bus_Data.val’;
Bus_Data.form=’sparse’;
Bus_Data.dim=2;
Bus_Data.uels={Bus.uels,{‘Pd’,’Qd’}};

%% GAMS calculation
GAMSfolder=[];
wgdxName=’dcopf_bess’;

% write gdx file as input for GAMS model
wgdx(wgdxName,y,gen,Bus,Slack,S_base,VoLL,Vwc,Eta_c,Eta_d,socmax,socini,wd,gb,Gen_Data,WCap,Bus_Data,Branches);
% gdxInfo DC_OPF %to print and check gdx input file

% run Gams file
gams([ ‘dcopf_bess’ ‘ lo=2’ ]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% read the gdx output file from GAMS
solGDX=[GAMSfolder ‘dcopf_bess_Solution.gdx’];
rs = struct(‘name’,’returnStat’,’compress’,’true’);
stats = rgdx(solGDX,rs);

gdxInfo DC_OPF_Solution %print gdx results

%matout file contains:
%OF, Pij(bus,node,t), delta(bus,t),
%Pg(Gen,t), lsh (bus,t), Pw(bus,t), pwc(bus,t),
%SOC(bus,t), Pd(bus,t), Pc(bus,t),NESS(bus);

% print of the stats of the GAMS run
disp([‘Best possible: ‘,num2str(stats.val(1,2))]);
disp([‘Objective: ‘,num2str(stats.val(2,2))]);
disp([‘Optimality gap: ‘,num2str((stats.val(2,2)-stats.val(1,2))/(stats.val(1,2))*100),’ %’]);
disp([‘Time: ‘,num2str(stats.val(3,2)), ‘ s’]);
disp([‘Model status: ‘,num2str(stats.val(4,2)), ‘(1: optimal, 8: integer solution model found)’]);
disp(‘Done!’)
disp(‘ ‘)

%% Results

% reading of the results related to charge power
P_c = struct(‘name’,’Pc’,’form’,’sparse’,’compress’,’true’);
sP_c=rgdx(solGDX,P_c);
P_c = zeros(n_bus,n_time);
for ijk=1:length(sP_c.val(:,3))
P_c(sP_c.val(ijk,1),sP_c.val(ijk,2)) = sP_c.val(ijk,3);
end

% reading of the results related to discharge power
P_d = struct(‘name’,’Pd’,’form’,’sparse’,’compress’,’true’);
sP_d=rgdx(solGDX,P_d);
P_d = zeros(n_bus,n_time);
for ijk=1:length(sP_d.val(:,3))
P_d(sP_d.val(ijk,1),sP_d.val(ijk,2)) = sP_d.val(ijk,3);
end

% reading of the results related to the State of Charge
SoC_batt=struct(‘name’,’SOC’,’form’,’full’,’compress’,’true’);
sSoC_B=rgdx(solGDX,SoC_batt);
SoC = zeros(n_bus,n_time);
SoC = sSoC_B.val;

% reading of the results related to number of batteries
NESS=struct(‘name’,’NESS’,’form’,’full’,’compress’,’true’);
sNESS=rgdx(solGDX,NESS);
idx = str2double(sNESS.uels{1});
N_ESS = zeros(n_bus,1);
N_ESS(idx) = sNESS.val;
figure;
stairs(N_ESS)
xlabel(‘bus number’)
ylabel(‘Number of ESS’)
ylim([0 max(N_ESS)+0.2])

% reading of the results related to generators
P_gen=struct(‘name’,’Pg’,’form’,’full’,’compress’,’true’);
sP_gen=rgdx(solGDX,P_gen);
P_gen = zeros(n_gen,n_time);
P_gen = sP_gen.val;

% wind curtailment result
P_cut=struct(‘name’,’pwc’,’form’,’full’,’compress’,’true’);
sP_Cut=rgdx(solGDX,P_cut);
P_cur = zeros(n_bus,n_time);
P_cur = sP_Cut.val;

% load shedding result
P_sh=struct(‘name’,’lsh’,’form’,’full’,’compress’,’true’);
sP_sh=rgdx(solGDX,P_sh);
Psh = zeros(n_bus,n_time);
Psh = sP_sh.val;
figure;
plot(P_cur)
hold on
plot(Psh)
hold off
xlabel(‘time [h]’)
ylabel(‘Wind curtailment and load shedding’)
legend(‘Curtailment’,’Load shedding’)

% Pij flows result
P_ij = struct(‘name’,’Pij’,’form’,’sparse’,’compress’,’true’);
sP_ij=rgdx(solGDX,P_ij);
P_ij = zeros(n_bus, n_bus, n_time); % If Pij is bus-to-bus over time
for ijk=1:length(sP_ij.val(:,3))
P_ij(sP_ij.val(ijk,1), sP_ij.val(ijk,2), sP_ij.val(ijk,3)) = sP_ij.val(ijk,4);
end
% Plotting power flow between bus i and bus j over time
i = 1; % Sending bus
j = 2; % Receiving bus
figure;
stairs(1:n_time, squeeze(P_ij(i,j,:)));
xlabel(‘Time’);
ylabel([‘Power flow from bus ‘, num2str(i), ‘ to bus ‘, num2str(j)]);
title(‘Power Flow between Buses over Time’);
This is integration with GAMS to see its output:
$set matout1 "’dcopf_bess_Solution.gdx’,returnStat,t,OF,Pd,Pc,SOC,Pij,Pg,lsh,pwc"
Sets
t time
DeltaT granularity
bus buses
Gen generators
GB(Gen,bus) generators’ bus location
slack(bus) slack bus;
alias(bus,node);
scalar
sbase base power
V_up upper voltage limit
V_lo lower voltage limit
eta_c charge efficiency
eta_d discharge efficiency
VWC Value of Wind curtailment
VOLL Value of Loss of Load;
parameters
SOCmax battery size [MWh]
SOC0 initial SoC [per unit]
GD(Gen,*) Generator data
branches(bus,node,*) r(pu) x(pu) b(pu) limit(MVA)
WD(t,*) Wind-demand daily variation patterns: ‘w’ wind ‘d’ demand
Wcap(bus) Wind capacities
BusData(bus,*) Demands of each bus in MW;
$gdxin "C:UsersuserDesktopmatpower7.1datadcopf_bess.gdx"
$load t Gen sbase bus slack VOLL VWC eta_c eta_d SOCmax SOC0 WD GD GB WCap BusData branches
$gdxin
Scalars Nmax;
branches(bus,node,’x’)$(branches(bus,node,’x’)=0)=branches(node,bus,’x’);
branches(bus,node,’limit’)$(branches(bus,node,’limit’)=0)=branches(node,bus,’limit’);
branches(bus,node,’bij’)$(conex(bus,node))=1/branches(node,bus,’x’);
Parameter conex(bus,node) Bus connectivity matrix;
conex(bus,node)$(conex(node,bus))=1;
conex(bus,node)$(branches(bus,node,’limit’) and branches(node,bus,’limit’)) =1;
Variables OF, Pij(bus,node,t),Pg(Gen,t),delta(bus,t), pw(bus,t),lsh(bus,t),SOC(bus,t),pwc(bus,t),Pd(bus,t),Pc(bus,t);
Integer variable NESS(bus);
Equations objfnc, balance,const1,const2,const3,const4,constESS,constESS2,constESS3(bus,t),constESS4(bus,t),constESS5(bus,t),constESS6(bus,t);
balance(bus,t)..sum(node$(conex(node,bus)),Pij(bus,node,t))=e=sum(Gen$(GB(Gen,bus)),Pg(Gen,t))+lsh(bus,t)$busdata(bus,’Pd’)+pw(bus,t)$Wcap(bus)-Pc(bus,t)$SOCmax(bus)+Pd(bus,t)$SOCmax(bus)-WD(t,’d’)*busdata(bus,’Pd’)/sbase;
const1(bus,node,t)$conex(bus,node)..Pij(bus,node,t)=e=branches(bus,node,’bij’)*(delta(bus,t)-delta(node,t));
const2(Gen,t)..Pg(Gen,t+1)-Pg(Gen,t)=l=Gendata(Gen,’RU’)/sbase;
const3(Gen,t)..Pg(Gen,t-1)-Pg(Gen,t)=l=Gendata(Gen,’RD’)/sbase;
const4(bus,t)$Wcap(bus)..pwc(bus,t)=e=WD(t,’w’)*Wcap(bus)/sbase-pw(bus,t);
constESS(bus,t)$SOCmax(bus)..SOC(bus,t)=e=(0.2*NESS(bus)*SOCmax/sbase)$(ord(t)=1)+SOC(bus,t-1)$(ord (t)>1)+(Pc(bus,t)$SOCmax(bus)*eta_c)-(Pd(bus,t)$SOCmax(bus)/eta_d);
constESS2(bus,t) .. SOC(bus,t)=l=NESS(bus)*SOCmax/sbase ;
*P_max charge
constESS3(bus,t) .. Pc(bus,t)=l=0.2*NESS(bus)*SOCmax/sbase ;
*P_max discharge
constESS4(bus,t) .. Pd(bus,t)=l=0.2*NESS(bus)*SOCmax/sbase ;
*limit on total ESSs
constESS5 .. sum(bus,NESS(bus))=l=Nmax;
*final SoC imposition (SoC0=0.2*SOCmax)
constESS6(bus) .. SOC(bus,’t24′)=e=0.2*NESS(bus)*SOCmax/sbase ;
objfnc..OF=e=sum((bus,Gen,t)$GB(Gen,bus),Pg(Gen,t)*Gendata(Gen,’b’)*sbase)+sum((bus,t),VOLL*lsh(bus,t)*sbase$(busdata(bus,’Pd’))+VOLW*pwc(bus,t)*Sbase$(Wcap(bus)));
Model load_flow/all/;
Pij.lo(bus,node,t)$(conex(bus,node))=-1*branches(bus,node,’limit’)/sbase;
Pij.up(bus,node,t)$(conex(bus,node))=1*branches(bus,node,’limit’)/sbase;
Pg.lo(Gen,t)=GD(Gen,’Pmin’)/sbase;
Pg.up(Gen,t)=GD(Gen,’Pmax’)/sbase;
delta.lo(bus,t)=-pi/2;
delta.up(bus,t)=pi/2;
delta.fx(slack,t)=0;
Pc.lo(bus,t)=0;
Pd.lo(bus,t)=0;
SOC.lo(bus,t)=0;
NESS.up(bus)=5;
Nmax=110;
pw.lo(bus,t)=0;
pw.up(bus,t)=WD(t,’w’)*Wcap(bus)/Sbase;
pwc.lo(bus,t)=0;
pwc.up(bus,t)=WD(t,’w’)*Wcap(bus)/Sbase;
lsh.lo(bus,t)=0;
lsh.up(bus,t)=WD(t,’d’)*busdata(bus,’Pd’)/Sbase;
Solve load_flow using MIP minimizing OF;

set stat /bestpossible,objective,tgams,tsolver,modelstat,solvestat/;
parameter returnStat(stat);
returnStat(‘bestpossible’) = loadflow.objest;
returnStat(‘objective’) = loadflow.objval;
returnStat(‘tsolver’) = loadflow.resusd;
returnStat(‘modelstat’) = loadflow.modelstat;
returnStat(‘solvestat’) = loadflow.solvestat;
execute_unload %matout1%
display Pg.l,Pij.l,delta.l, pwc.l,lsh.l,Pc.l,Pd.l;Hi, please someone tell me regarding integration of MATLAB code with GAMS, i want to implement DCOPF. I receive this error in MATLAB:Error using wgdx:Amount of full format .val data exceeds number of UELs. Please try to resolve it, is my wgdx and rgdx functions working?
This is my MATLAB code:
clear all
close all
clc
%inputs
mpc=load (‘grid_data_24bus.mat’);
profile_L=mpc.WD(:,3);
profile_gen_rand=[1.0057;1.027;0.98;0.967;0.9837;0.971;1.001;1.04063;1.0129;0.961;0.989;0.9554;1.000129;0.9937;1.0497;1.0316;0.9986;1.03944;0.963;0.989;1.0427;1.041749;1.0214;1.0118];
n_time = 24; % time intervals
DeltaT = 24/n_time; % granularity
n_bus = size(mpc.BusData,1); % Number of buses
n_brnc = size(mpc.branches,1); % Number of branches
n_gen = size(mpc.Gen,1); % num of generators
slack = 13; % slack choice
VS = 1; % default value for slack voltage
V_up=1.1; % Voltage limits
V_lo=0.9; % Voltage limits
% Time profiles
mpc.BusData = repmat (mpc.BusData, [1,1,n_time]);
mpc.Gen = repmat (mpc.Gen, [1,1,n_time]);
for t=1:n_time
mpc.BusData(:,2,t)=profile_L(t).*mpc.BusData(:,2,t);% column 2 have Pd
for g=1:size(mpc.Gen,1)
mpc.Gen(g,3:4,t)=profile_gen_rand(t).*( mpc.Gen(g,3:4,t));%pmax and pmin
end
end
mpc_pre=mpc; %save mpcstruct
%%% ESS allocation
SOC_max = 100; % battery size [MWh]
SOC_0 = 0.2*SOC_max/mpc.Sbase; % initial SoC [per unit]
eta_c = 0.95; % charge efficiency
eta_d = 0.9; % discharge efficiency
VOLL = 10000; % Value of Loss of Load
VWC = 50; % Value of loss of wind ($/MWh)
% Define Wind Capacity and Generation Profiles for Specific Buses
% Define Wind Capacity and Generation Profiles for Specific Buses
Wcap = zeros(n_bus, 1); % Wind capacity array for all buses
Wcap(8) = 200; % Bus 8: 200 MW
Wcap(19) = 150; % Bus 19: 150 MW
Wcap(21) = 100; % Bus 21: 100 MW
% Wind data over 24 hours (from ‘w’ column)
wind_profile = mpc.WD(:, 2); % 24-hour wind profile
% Wind Generation Matrix (Buses x Time)
wind_gen = Wcap * wind_profile’; % Corrected dimensions