From what I can tell both kse and ksdensity are used to get a smoothly defined probability distribution function from vector of finite data. Is there a way to get MATLAB to produce the actual function as the output instead of finite set of input/output vectors? I need to be able to call on the kernal density estimator at an arbitrary point in my domain and I’d prefer not to just interpolate between the values provided. Is there a functionality of either of these two that gives the true function/estimator or do I need to write my own code to accomplish this?From what I can tell both kse and ksdensity are used to get a smoothly defined probability distribution function from vector of finite data. Is there a way to get MATLAB to produce the actual function as the output instead of finite set of input/output vectors? I need to be able to call on the kernal density estimator at an arbitrary point in my domain and I’d prefer not to just interpolate between the values provided. Is there a functionality of either of these two that gives the true function/estimator or do I need to write my own code to accomplish this? From what I can tell both kse and ksdensity are used to get a smoothly defined probability distribution function from vector of finite data. Is there a way to get MATLAB to produce the actual function as the output instead of finite set of input/output vectors? I need to be able to call on the kernal density estimator at an arbitrary point in my domain and I’d prefer not to just interpolate between the values provided. Is there a functionality of either of these two that gives the true function/estimator or do I need to write my own code to accomplish this? statistics MATLAB Answers — New Questions

​

basiaclly whenever im typing code in the command window and the editor there is never auto completion options just today they stopped showing up. is there a way to turn them back? i checked preferences and under the auto completion tab everything is checked i dont know the issue.basiaclly whenever im typing code in the command window and the editor there is never auto completion options just today they stopped showing up. is there a way to turn them back? i checked preferences and under the auto completion tab everything is checked i dont know the issue. basiaclly whenever im typing code in the command window and the editor there is never auto completion options just today they stopped showing up. is there a way to turn them back? i checked preferences and under the auto completion tab everything is checked i dont know the issue. auto completion MATLAB Answers — New Questions

​

Can someone give me a steer on the best way to remove random speckle from a monochrome image?
I have a series of 16-bit monochrome images which show the pattern of light uniformity on a flat surface at different wavelengths. The images are stored in a 3-d array such that the first two dimensions (row, column) are the image pixel values and the third dimension (‘page’) is the wavelength value.
The images are ‘noisy’. That is, if you chose a page of the array and plotted its 2-d surface, it would not be smooth but show small peaks and troughs either side of a mean value at every pixel location. I want to smooth out these peaks and troughs.
The imgaussfilt function applies a Gaussian smoothing filter to images. There is also a method using the polyfitn function on file exchange. What is the difference in outcome between these approaches?
I tried comparing the two using the code below (my histogram comparison method is crude!). Visually, method 2 looks best but I don’t know whether it’s a fair test. Also, I can select the polyfitn polynomial order that works best only when I know the original noise-less signal. How should I select the order when the signal is unknown?
%Make up 3 dummy test images 64 x 80
rows = 64;
columns = 80;
y=1:rows;
x=1:columns;
x1=(2-exp(x./1000)).*50000;
y1= 1 – (y.*0.008-0.4).^2; % profile at wavelength 1
y2= 1 – (y.*0.009-0.4).^2; % profile at wavelength 2
y3= 1 – (y.*0.01-0.4).^2; % profile at wavelength 3
y1=y1′;y2=y2′;y3=y3′;
zSignal= cat (3, y1*x1, y2*x1, y3*x1);

% Add some noise to images
z= zSignal + 350.*rand(rows,columns,3); % randn more realistic?
z= round(z); %

% SMOOTHING METHOD 1 (Try blurring the image on page 1)
image1 = uint16(squeeze(z(:,:,1))); % create the dummy image
zBlur = double(imgaussfilt(image1,2));
figure;
histogram (zBlur-zSignal(:,:,1)); % check for absolute error #1

%SMOOTHING METHOD 2 (Create a model zPred for each page)
% (with thanks to ImageAnalyst)
zPred = zeros (64,80,3);
[X,Y] = meshgrid(1:columns, 1:rows);
x1d = reshape(X, numel(X), 1);
y1d = reshape(Y, numel(Y), 1);
zPred = zeros (64,80,3);
for page = 1:3
image = uint16(squeeze(z(:,:,page))); % create the dummy image
z1d = double(reshape(image, numel(image), 1));
polynomialOrder = 4;
p = polyfitn([x1d y1d], z1d, polynomialOrder);
zg = polyvaln(p, [x1d y1d]);
zPred(:,:,page) = reshape(zg, [rows columns]);
end
figure;
histogram (zPred(:,:,1)-zSignal(:,:,1)); % check for error #2
return;Can someone give me a steer on the best way to remove random speckle from a monochrome image?
I have a series of 16-bit monochrome images which show the pattern of light uniformity on a flat surface at different wavelengths. The images are stored in a 3-d array such that the first two dimensions (row, column) are the image pixel values and the third dimension (‘page’) is the wavelength value.
The images are ‘noisy’. That is, if you chose a page of the array and plotted its 2-d surface, it would not be smooth but show small peaks and troughs either side of a mean value at every pixel location. I want to smooth out these peaks and troughs.
The imgaussfilt function applies a Gaussian smoothing filter to images. There is also a method using the polyfitn function on file exchange. What is the difference in outcome between these approaches?
I tried comparing the two using the code below (my histogram comparison method is crude!). Visually, method 2 looks best but I don’t know whether it’s a fair test. Also, I can select the polyfitn polynomial order that works best only when I know the original noise-less signal. How should I select the order when the signal is unknown?
%Make up 3 dummy test images 64 x 80
rows = 64;
columns = 80;
y=1:rows;
x=1:columns;
x1=(2-exp(x./1000)).*50000;
y1= 1 – (y.*0.008-0.4).^2; % profile at wavelength 1
y2= 1 – (y.*0.009-0.4).^2; % profile at wavelength 2
y3= 1 – (y.*0.01-0.4).^2; % profile at wavelength 3
y1=y1′;y2=y2′;y3=y3′;
zSignal= cat (3, y1*x1, y2*x1, y3*x1);

% Add some noise to images
z= zSignal + 350.*rand(rows,columns,3); % randn more realistic?
z= round(z); %

% SMOOTHING METHOD 1 (Try blurring the image on page 1)
image1 = uint16(squeeze(z(:,:,1))); % create the dummy image
zBlur = double(imgaussfilt(image1,2));
figure;
histogram (zBlur-zSignal(:,:,1)); % check for absolute error #1

%SMOOTHING METHOD 2 (Create a model zPred for each page)
% (with thanks to ImageAnalyst)
zPred = zeros (64,80,3);
[X,Y] = meshgrid(1:columns, 1:rows);
x1d = reshape(X, numel(X), 1);
y1d = reshape(Y, numel(Y), 1);
zPred = zeros (64,80,3);
for page = 1:3
image = uint16(squeeze(z(:,:,page))); % create the dummy image
z1d = double(reshape(image, numel(image), 1));
polynomialOrder = 4;
p = polyfitn([x1d y1d], z1d, polynomialOrder);
zg = polyvaln(p, [x1d y1d]);
zPred(:,:,page) = reshape(zg, [rows columns]);
end
figure;
histogram (zPred(:,:,1)-zSignal(:,:,1)); % check for error #2
return; Can someone give me a steer on the best way to remove random speckle from a monochrome image?
I have a series of 16-bit monochrome images which show the pattern of light uniformity on a flat surface at different wavelengths. The images are stored in a 3-d array such that the first two dimensions (row, column) are the image pixel values and the third dimension (‘page’) is the wavelength value.
The images are ‘noisy’. That is, if you chose a page of the array and plotted its 2-d surface, it would not be smooth but show small peaks and troughs either side of a mean value at every pixel location. I want to smooth out these peaks and troughs.
The imgaussfilt function applies a Gaussian smoothing filter to images. There is also a method using the polyfitn function on file exchange. What is the difference in outcome between these approaches?
I tried comparing the two using the code below (my histogram comparison method is crude!). Visually, method 2 looks best but I don’t know whether it’s a fair test. Also, I can select the polyfitn polynomial order that works best only when I know the original noise-less signal. How should I select the order when the signal is unknown?
%Make up 3 dummy test images 64 x 80
rows = 64;
columns = 80;
y=1:rows;
x=1:columns;
x1=(2-exp(x./1000)).*50000;
y1= 1 – (y.*0.008-0.4).^2; % profile at wavelength 1
y2= 1 – (y.*0.009-0.4).^2; % profile at wavelength 2
y3= 1 – (y.*0.01-0.4).^2; % profile at wavelength 3
y1=y1′;y2=y2′;y3=y3′;
zSignal= cat (3, y1*x1, y2*x1, y3*x1);

% Add some noise to images
z= zSignal + 350.*rand(rows,columns,3); % randn more realistic?
z= round(z); %

% SMOOTHING METHOD 1 (Try blurring the image on page 1)
image1 = uint16(squeeze(z(:,:,1))); % create the dummy image
zBlur = double(imgaussfilt(image1,2));
figure;
histogram (zBlur-zSignal(:,:,1)); % check for absolute error #1

%SMOOTHING METHOD 2 (Create a model zPred for each page)
% (with thanks to ImageAnalyst)
zPred = zeros (64,80,3);
[X,Y] = meshgrid(1:columns, 1:rows);
x1d = reshape(X, numel(X), 1);
y1d = reshape(Y, numel(Y), 1);
zPred = zeros (64,80,3);
for page = 1:3
image = uint16(squeeze(z(:,:,page))); % create the dummy image
z1d = double(reshape(image, numel(image), 1));
polynomialOrder = 4;
p = polyfitn([x1d y1d], z1d, polynomialOrder);
zg = polyvaln(p, [x1d y1d]);
zPred(:,:,page) = reshape(zg, [rows columns]);
end
figure;
histogram (zPred(:,:,1)-zSignal(:,:,1)); % check for error #2
return; polyfitn, polyvaln, imgaussfilt, image MATLAB Answers — New Questions

​

Hello all, I currently have the code below that runs calculations and exports a matrix for my .csv files. Right now this is assuming my data is a clean 3Xwhatever matrix in my .csv files. Within my .csv files I only want to look at data in columns starting at B10, C10, and D10 and then running ot the end of those columns. Can matlab point to certain data within a .csv? I am assumning I need to tweak the M = readmatrix(fullfile(fn(ii).folder,fn(ii).name)) line or M = readmatrix(fullfile(fn(ii).folder,fn(ii).name)) line.

Any help is much appreciated! Thank you.

% appropriate dir() call that returns info
% about the files you want to process:
fn = dir(‘*.csv’); % this call returns info about .csv files in the current directory;
% you may need to modify it to work for your file locations
% (see dir documentation)
% number of files:
N_files = numel(fn);
% pre-allocate results matrix (one row per file, 3 columns):
results = zeros(N_files,3);
% read and process each file:
for ii = 1:N
% read the file:
M = readmatrix(fullfile(fn(ii).folder,fn(ii).name));
% process the file’s data:
S = sum(abs(diff(M,1,1)),1);
% store the result:
results(ii,:) = S;
end
% write the results file (can be located anywhere):
writematrix(results,’results.csv’)Hello all, I currently have the code below that runs calculations and exports a matrix for my .csv files. Right now this is assuming my data is a clean 3Xwhatever matrix in my .csv files. Within my .csv files I only want to look at data in columns starting at B10, C10, and D10 and then running ot the end of those columns. Can matlab point to certain data within a .csv? I am assumning I need to tweak the M = readmatrix(fullfile(fn(ii).folder,fn(ii).name)) line or M = readmatrix(fullfile(fn(ii).folder,fn(ii).name)) line.

Any help is much appreciated! Thank you.

% appropriate dir() call that returns info
% about the files you want to process:
fn = dir(‘*.csv’); % this call returns info about .csv files in the current directory;
% you may need to modify it to work for your file locations
% (see dir documentation)
% number of files:
N_files = numel(fn);
% pre-allocate results matrix (one row per file, 3 columns):
results = zeros(N_files,3);
% read and process each file:
for ii = 1:N
% read the file:
M = readmatrix(fullfile(fn(ii).folder,fn(ii).name));
% process the file’s data:
S = sum(abs(diff(M,1,1)),1);
% store the result:
results(ii,:) = S;
end
% write the results file (can be located anywhere):
writematrix(results,’results.csv’) Hello all, I currently have the code below that runs calculations and exports a matrix for my .csv files. Right now this is assuming my data is a clean 3Xwhatever matrix in my .csv files. Within my .csv files I only want to look at data in columns starting at B10, C10, and D10 and then running ot the end of those columns. Can matlab point to certain data within a .csv? I am assumning I need to tweak the M = readmatrix(fullfile(fn(ii).folder,fn(ii).name)) line or M = readmatrix(fullfile(fn(ii).folder,fn(ii).name)) line.

Any help is much appreciated! Thank you.

% appropriate dir() call that returns info
% about the files you want to process:
fn = dir(‘*.csv’); % this call returns info about .csv files in the current directory;
% you may need to modify it to work for your file locations
% (see dir documentation)
% number of files:
N_files = numel(fn);
% pre-allocate results matrix (one row per file, 3 columns):
results = zeros(N_files,3);
% read and process each file:
for ii = 1:N
% read the file:
M = readmatrix(fullfile(fn(ii).folder,fn(ii).name));
% process the file’s data:
S = sum(abs(diff(M,1,1)),1);
% store the result:
results(ii,:) = S;
end
% write the results file (can be located anywhere):
writematrix(results,’results.csv’) .csv, matrices, matrix MATLAB Answers — New Questions

​

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

% 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

% 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; power system MATLAB Answers — New Questions

​

In Matlab OOP, I wanna use a function like mustBeEqualSize(a,b) from the Examples on function property validation to validate a property:
classdef foo
properties
a
b {mustBeEqualSize(a,b)}
end
end
However, this gives me an error.In Matlab OOP, I wanna use a function like mustBeEqualSize(a,b) from the Examples on function property validation to validate a property:
classdef foo
properties
a
b {mustBeEqualSize(a,b)}
end
end
However, this gives me an error. In Matlab OOP, I wanna use a function like mustBeEqualSize(a,b) from the Examples on function property validation to validate a property:
classdef foo
properties
a
b {mustBeEqualSize(a,b)}
end
end
However, this gives me an error. matlab, oop MATLAB Answers — New Questions

​

I got this error "Nonlinear constraint function is undefined at initial point. Fmincon cannot continue." when i ran a simulation to remove offset using extended kalman filter "EKF" in quadruple system. The target is to control level of water in the lower tank h1 and h2 by manipulating v1 and v2. check below for the code;
NMPC CODE:
clc, clear all
%% Simulink Model
% To run this example, Simulink(R) is required.
if ~mpcchecktoolboxinstalled(‘simulink’)
disp(‘Simulink is required to run this example.’)
return
end
%%initialize NMPC object and set properties
nx = 8;
ny = 2;
%nu = 2;
%’UD’, [5,6,7,8,9,10]
nlobj = nlmpc(nx, ny, ‘MV’, [1,2], ‘UD’, [3,4,5,6]);
%The prediction model sample time is the same as the controller sample time
Ts = 3; %Ts = 2
nlobj.Ts = Ts;
%Prediction and control horizon
nlobj.PredictionHorizon = 20;
nlobj.ControlHorizon = [1];
%weights on output and manipulated variables
nlobj.Weights.OutputVariables = [1000,1000]; %0.1
nlobj.Weights.ManipulatedVariables = [10,10]; %0.05;
%nlobj.Weights.ManipulatedVariablesRate = [10,10];
%constraints on manipulated variables
nlobj.MV(1).RateMin = -1;
nlobj.MV(2).RateMin = -1;
nlobj.MV(1).RateMax = 1;
nlobj.Mv(2).RateMax = 1;
nlobj.MV(1).Min = 2.7;
nlobj.MV(2).Min = 2.7;
nlobj.MV(1).Max = 3.3;
nlobj.MV(1).Max = 3.3;
%constraint on output variables to keep them within reasonable bounds
nlobj.OV(1).Min = 11.16;
nlobj.OV(2).Min = 11.43;
%nlobj.OV(3).Min = 2.5;
%nlobj.OV(4).Min = 2.5;
nlobj.OV(1).Max = 13.64;
nlobj.OV(2).Max = 13.97;
%nlobj.OV(3).Max = 8;
%nlobj.OV(4).Max = 8;
%% specify the nonlinear state and output functions
nlobj.Model.StateFcn = ‘four_tankStateFcnCT’;
nlobj.Model.OutputFcn = ‘four_tankOutputFcn’;
SolverOptions.Algorithm = ‘sqp’;
%% validate the nonlinear MPC object with initial conditions
h0 = [12.44;13.17;4.73;4.986;0;0;0;0]; %[0.0033;920.86; 0.0104; 0.0085; 0.0118; 798];
u0 = [3.15;3.15];
validateFcns(nlobj, h0, u0);

four_tankStateFcnCT CODE:
function dhdt = four_tankStateFcnCT(h, u)
% Define system parameters
g = 981; % gravitational constant (m/s^2)

% Tank areas (in cm^2)
A1 = 28;
A2 = 32;
A3 = 28;
A4 = 32;

% Outlet areas (in cm^2)
a1 = 0.071;
a2 = 0.057;
a3 = 0.071;
a4 = 0.057;

% Flow distribution coefficients
gamma1 = 0.43;
gamma2 = 0.34;

% Pump constants (flow per voltage)
k1 = 3.14;
k2 = 3.29;

% State variables (water levels in each tank)
h1 = h(1);
h2 = h(2);
h3 = h(3);
h4 = h(4);

% Inputs (pump voltages)
v1 = u(1);
v2 = u(2);
dh1 = h1;
dh2 = h2;
dh3 = h3;
dh4 = h4;
dhdt = zeros(8,1);
% Equations from the journal
dhdt(1) = -a1/A1 * sqrt(2 * g * h1) + a3/A1 * sqrt(2 * g * h3) + (gamma1 * k1 * v1) / A1 + dh1;
dhdt(2) = -a2/A2 * sqrt(2 * g * h2) + a4/A2 * sqrt(2 * g * h4) + (gamma2 * k2 * v2) / A2 + dh2;
dhdt(3) = -a3/A3 * sqrt(2 * g * h3) + ((1 – gamma2) * k2 * v2) / A3 + dh3;
dhdt(4) = -a4/A4 * sqrt(2 * g * h4) + ((1 – gamma1) * k1 * v1) / A4 + dh4;
dhdt(5) = u(3);
dhdt(6) = u(4);
dhdt(7) = u(5);
dhdt(8) = u(6);

four_tankStateFcnDT CODE:
function hk1 = four_tankStateFcnDT(hk,uk)
hstep = 1;
hk1 = hk(:);
Nsteps = 2; % Number of integration time steps for Euler method
uk1 = [uk(:);0;0;0;0]; % The manipulated input + white noise of size 2
for i = 1:Nsteps
hk1 = hk1 + hstep*four_tankStateFcnCT(hk1,uk1);
end

four_tankMeasFcn CODE:
function y = four_tankMeasFcn(h)

y = [h(1) ;
h(2)] ;

four_tankOutputFcn CODE:
function y = four_tankOutputFcn(h,u)
y = [h(1) ;
h(2)] ;I got this error "Nonlinear constraint function is undefined at initial point. Fmincon cannot continue." when i ran a simulation to remove offset using extended kalman filter "EKF" in quadruple system. The target is to control level of water in the lower tank h1 and h2 by manipulating v1 and v2. check below for the code;
NMPC CODE:
clc, clear all
%% Simulink Model
% To run this example, Simulink(R) is required.
if ~mpcchecktoolboxinstalled(‘simulink’)
disp(‘Simulink is required to run this example.’)
return
end
%%initialize NMPC object and set properties
nx = 8;
ny = 2;
%nu = 2;
%’UD’, [5,6,7,8,9,10]
nlobj = nlmpc(nx, ny, ‘MV’, [1,2], ‘UD’, [3,4,5,6]);
%The prediction model sample time is the same as the controller sample time
Ts = 3; %Ts = 2
nlobj.Ts = Ts;
%Prediction and control horizon
nlobj.PredictionHorizon = 20;
nlobj.ControlHorizon = [1];
%weights on output and manipulated variables
nlobj.Weights.OutputVariables = [1000,1000]; %0.1
nlobj.Weights.ManipulatedVariables = [10,10]; %0.05;
%nlobj.Weights.ManipulatedVariablesRate = [10,10];
%constraints on manipulated variables
nlobj.MV(1).RateMin = -1;
nlobj.MV(2).RateMin = -1;
nlobj.MV(1).RateMax = 1;
nlobj.Mv(2).RateMax = 1;
nlobj.MV(1).Min = 2.7;
nlobj.MV(2).Min = 2.7;
nlobj.MV(1).Max = 3.3;
nlobj.MV(1).Max = 3.3;
%constraint on output variables to keep them within reasonable bounds
nlobj.OV(1).Min = 11.16;
nlobj.OV(2).Min = 11.43;
%nlobj.OV(3).Min = 2.5;
%nlobj.OV(4).Min = 2.5;
nlobj.OV(1).Max = 13.64;
nlobj.OV(2).Max = 13.97;
%nlobj.OV(3).Max = 8;
%nlobj.OV(4).Max = 8;
%% specify the nonlinear state and output functions
nlobj.Model.StateFcn = ‘four_tankStateFcnCT’;
nlobj.Model.OutputFcn = ‘four_tankOutputFcn’;
SolverOptions.Algorithm = ‘sqp’;
%% validate the nonlinear MPC object with initial conditions
h0 = [12.44;13.17;4.73;4.986;0;0;0;0]; %[0.0033;920.86; 0.0104; 0.0085; 0.0118; 798];
u0 = [3.15;3.15];
validateFcns(nlobj, h0, u0);

four_tankStateFcnCT CODE:
function dhdt = four_tankStateFcnCT(h, u)
% Define system parameters
g = 981; % gravitational constant (m/s^2)

% Tank areas (in cm^2)
A1 = 28;
A2 = 32;
A3 = 28;
A4 = 32;

% Outlet areas (in cm^2)
a1 = 0.071;
a2 = 0.057;
a3 = 0.071;
a4 = 0.057;

% Flow distribution coefficients
gamma1 = 0.43;
gamma2 = 0.34;

% Pump constants (flow per voltage)
k1 = 3.14;
k2 = 3.29;

% State variables (water levels in each tank)
h1 = h(1);
h2 = h(2);
h3 = h(3);
h4 = h(4);

% Inputs (pump voltages)
v1 = u(1);
v2 = u(2);
dh1 = h1;
dh2 = h2;
dh3 = h3;
dh4 = h4;
dhdt = zeros(8,1);
% Equations from the journal
dhdt(1) = -a1/A1 * sqrt(2 * g * h1) + a3/A1 * sqrt(2 * g * h3) + (gamma1 * k1 * v1) / A1 + dh1;
dhdt(2) = -a2/A2 * sqrt(2 * g * h2) + a4/A2 * sqrt(2 * g * h4) + (gamma2 * k2 * v2) / A2 + dh2;
dhdt(3) = -a3/A3 * sqrt(2 * g * h3) + ((1 – gamma2) * k2 * v2) / A3 + dh3;
dhdt(4) = -a4/A4 * sqrt(2 * g * h4) + ((1 – gamma1) * k1 * v1) / A4 + dh4;
dhdt(5) = u(3);
dhdt(6) = u(4);
dhdt(7) = u(5);
dhdt(8) = u(6);

four_tankStateFcnDT CODE:
function hk1 = four_tankStateFcnDT(hk,uk)
hstep = 1;
hk1 = hk(:);
Nsteps = 2; % Number of integration time steps for Euler method
uk1 = [uk(:);0;0;0;0]; % The manipulated input + white noise of size 2
for i = 1:Nsteps
hk1 = hk1 + hstep*four_tankStateFcnCT(hk1,uk1);
end

four_tankMeasFcn CODE:
function y = four_tankMeasFcn(h)

y = [h(1) ;
h(2)] ;

four_tankOutputFcn CODE:
function y = four_tankOutputFcn(h,u)
y = [h(1) ;
h(2)] ; I got this error "Nonlinear constraint function is undefined at initial point. Fmincon cannot continue." when i ran a simulation to remove offset using extended kalman filter "EKF" in quadruple system. The target is to control level of water in the lower tank h1 and h2 by manipulating v1 and v2. check below for the code;
NMPC CODE:
clc, clear all
%% Simulink Model
% To run this example, Simulink(R) is required.
if ~mpcchecktoolboxinstalled(‘simulink’)
disp(‘Simulink is required to run this example.’)
return
end
%%initialize NMPC object and set properties
nx = 8;
ny = 2;
%nu = 2;
%’UD’, [5,6,7,8,9,10]
nlobj = nlmpc(nx, ny, ‘MV’, [1,2], ‘UD’, [3,4,5,6]);
%The prediction model sample time is the same as the controller sample time
Ts = 3; %Ts = 2
nlobj.Ts = Ts;
%Prediction and control horizon
nlobj.PredictionHorizon = 20;
nlobj.ControlHorizon = [1];
%weights on output and manipulated variables
nlobj.Weights.OutputVariables = [1000,1000]; %0.1
nlobj.Weights.ManipulatedVariables = [10,10]; %0.05;
%nlobj.Weights.ManipulatedVariablesRate = [10,10];
%constraints on manipulated variables
nlobj.MV(1).RateMin = -1;
nlobj.MV(2).RateMin = -1;
nlobj.MV(1).RateMax = 1;
nlobj.Mv(2).RateMax = 1;
nlobj.MV(1).Min = 2.7;
nlobj.MV(2).Min = 2.7;
nlobj.MV(1).Max = 3.3;
nlobj.MV(1).Max = 3.3;
%constraint on output variables to keep them within reasonable bounds
nlobj.OV(1).Min = 11.16;
nlobj.OV(2).Min = 11.43;
%nlobj.OV(3).Min = 2.5;
%nlobj.OV(4).Min = 2.5;
nlobj.OV(1).Max = 13.64;
nlobj.OV(2).Max = 13.97;
%nlobj.OV(3).Max = 8;
%nlobj.OV(4).Max = 8;
%% specify the nonlinear state and output functions
nlobj.Model.StateFcn = ‘four_tankStateFcnCT’;
nlobj.Model.OutputFcn = ‘four_tankOutputFcn’;
SolverOptions.Algorithm = ‘sqp’;
%% validate the nonlinear MPC object with initial conditions
h0 = [12.44;13.17;4.73;4.986;0;0;0;0]; %[0.0033;920.86; 0.0104; 0.0085; 0.0118; 798];
u0 = [3.15;3.15];
validateFcns(nlobj, h0, u0);

four_tankStateFcnCT CODE:
function dhdt = four_tankStateFcnCT(h, u)
% Define system parameters
g = 981; % gravitational constant (m/s^2)

% Tank areas (in cm^2)
A1 = 28;
A2 = 32;
A3 = 28;
A4 = 32;

% Outlet areas (in cm^2)
a1 = 0.071;
a2 = 0.057;
a3 = 0.071;
a4 = 0.057;

% Flow distribution coefficients
gamma1 = 0.43;
gamma2 = 0.34;

% Pump constants (flow per voltage)
k1 = 3.14;
k2 = 3.29;

% State variables (water levels in each tank)
h1 = h(1);
h2 = h(2);
h3 = h(3);
h4 = h(4);

% Inputs (pump voltages)
v1 = u(1);
v2 = u(2);
dh1 = h1;
dh2 = h2;
dh3 = h3;
dh4 = h4;
dhdt = zeros(8,1);
% Equations from the journal
dhdt(1) = -a1/A1 * sqrt(2 * g * h1) + a3/A1 * sqrt(2 * g * h3) + (gamma1 * k1 * v1) / A1 + dh1;
dhdt(2) = -a2/A2 * sqrt(2 * g * h2) + a4/A2 * sqrt(2 * g * h4) + (gamma2 * k2 * v2) / A2 + dh2;
dhdt(3) = -a3/A3 * sqrt(2 * g * h3) + ((1 – gamma2) * k2 * v2) / A3 + dh3;
dhdt(4) = -a4/A4 * sqrt(2 * g * h4) + ((1 – gamma1) * k1 * v1) / A4 + dh4;
dhdt(5) = u(3);
dhdt(6) = u(4);
dhdt(7) = u(5);
dhdt(8) = u(6);

four_tankStateFcnDT CODE:
function hk1 = four_tankStateFcnDT(hk,uk)
hstep = 1;
hk1 = hk(:);
Nsteps = 2; % Number of integration time steps for Euler method
uk1 = [uk(:);0;0;0;0]; % The manipulated input + white noise of size 2
for i = 1:Nsteps
hk1 = hk1 + hstep*four_tankStateFcnCT(hk1,uk1);
end

four_tankMeasFcn CODE:
function y = four_tankMeasFcn(h)

y = [h(1) ;
h(2)] ;

four_tankOutputFcn CODE:
function y = four_tankOutputFcn(h,u)
y = [h(1) ;
h(2)] ; optimization, fmincon MATLAB Answers — New Questions

​

I’m studing the example "CreateConfigureAndSimulate5GNRNodesExample.mlx".
In this example you can add mobility to a UE using the method:
addMobility(ue,BoundaryShape="rectangle")
But, if I have coordinates from my cenario and UE by along the time (a list of diferent lon,lat coods from user monitoring, for example) it is possible to simulate the mobility using this data?I’m studing the example "CreateConfigureAndSimulate5GNRNodesExample.mlx".
In this example you can add mobility to a UE using the method:
addMobility(ue,BoundaryShape="rectangle")
But, if I have coordinates from my cenario and UE by along the time (a list of diferent lon,lat coods from user monitoring, for example) it is possible to simulate the mobility using this data? I’m studing the example "CreateConfigureAndSimulate5GNRNodesExample.mlx".
In this example you can add mobility to a UE using the method:
addMobility(ue,BoundaryShape="rectangle")
But, if I have coordinates from my cenario and UE by along the time (a list of diferent lon,lat coods from user monitoring, for example) it is possible to simulate the mobility using this data? ue, mobility, 5g MATLAB Answers — New Questions

​

I am trying to write a script to run a series of test cases and would like to open Test Manager as little as possible. I have found that you can programmatically create a test case and specify the type as either ‘baseline’,’equivalence’, or ‘simulation’. However, I have not been able to find a type to create a test case for real-time simulation without having to open Test Manager.I am trying to write a script to run a series of test cases and would like to open Test Manager as little as possible. I have found that you can programmatically create a test case and specify the type as either ‘baseline’,’equivalence’, or ‘simulation’. However, I have not been able to find a type to create a test case for real-time simulation without having to open Test Manager. I am trying to write a script to run a series of test cases and would like to open Test Manager as little as possible. I have found that you can programmatically create a test case and specify the type as either ‘baseline’,’equivalence’, or ‘simulation’. However, I have not been able to find a type to create a test case for real-time simulation without having to open Test Manager. simulink test, real time test, test manager, test case MATLAB Answers — New Questions

​

I have pictures of cracks in soil taken in an outdoor experiment which I am looking to quanitfy. My ultimate goal is to segment the cracks from the intact soil peds to calculate width and then skeletonise to get length, intersection angles etc. As the photos are taken outdoors when i try to segment the cracks from the intact soil peds there is a lot of noise due to the uneven lighting conditions and highly pitted surface caused by rain action. Can anyone help with problem and suggest how to best formulate a code to address a high degree of hetereogenity between images?I have pictures of cracks in soil taken in an outdoor experiment which I am looking to quanitfy. My ultimate goal is to segment the cracks from the intact soil peds to calculate width and then skeletonise to get length, intersection angles etc. As the photos are taken outdoors when i try to segment the cracks from the intact soil peds there is a lot of noise due to the uneven lighting conditions and highly pitted surface caused by rain action. Can anyone help with problem and suggest how to best formulate a code to address a high degree of hetereogenity between images? I have pictures of cracks in soil taken in an outdoor experiment which I am looking to quanitfy. My ultimate goal is to segment the cracks from the intact soil peds to calculate width and then skeletonise to get length, intersection angles etc. As the photos are taken outdoors when i try to segment the cracks from the intact soil peds there is a lot of noise due to the uneven lighting conditions and highly pitted surface caused by rain action. Can anyone help with problem and suggest how to best formulate a code to address a high degree of hetereogenity between images? cracks, image segmentation, image processing, uneven lighting MATLAB Answers — New Questions

​