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:\Users\user\Desktop\matpower7.1\data\dcopf_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;

Answer 1

[Pooja Kumari]

Hi Ayesha,

I understand that you are facing error using wgdx while trying to implement DCOPF using integration of MATLAB code with GAMS.

The error message indicates a mismatch between the data structure you are passing from MATLAB to GAMS. This usually means the data dimensions or UELs do not align with what is expected in your GAMS model.

A possible workaround is matching the dimensions of .val data matrix with the expected format in GAMS.

Check for DLL and Mex file dependencies :-

• GAMS distributes its function in a set of compiled DLL and MEX files. It is possible to use compiled function (format mexw64 for example) to your compiled application, as long as you add them to the application when you package it.
• If the function depends from some standard DLL, in theory you should not have to add these DLL to the application. The compiled function should be able to find these DLL in the common folder in a process as described in the following page: https://docs.microsoft.com/en-us/windows/desktop/dlls/dynamic-link-library-search-order#standard-search-order-for-desktop-applications