error regarding wgdx and rgdx

18 次查看（过去 30 天）
Ayesha
Ayesha 2024-10-1
回答： Pooja Kumari 2024-10-8
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;

请先登录，再回答此问题。

回答（1 个）

Pooja Kumari
Pooja Kumari 2024-10-8
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 :-

类别

在 Help Center 和 File Exchange 中查找有关 Modeling 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by