Email: helpdesk@telkomuniversity.ac.id

This Portal for internal use only!

  • My Download
  • Checkout
Application Package Repository Telkom University
All Categories

All Categories

  • Visual Paradigm
  • IBM
  • Adobe
  • Google
  • Matlab
  • Microsoft
    • Microsoft Apps
    • Analytics
    • AI + Machine Learning
    • Compute
    • Database
    • Developer Tools
    • Internet Of Things
    • Learning Services
    • Middleware System
    • Networking
    • Operating System
    • Productivity Tools
    • Security
    • VLS
      • Windows
      • Office
  • Opensource
  • Wordpress
    • Plugin WP
    • Themes WP
  • Others

Search

0 Wishlist

Cart

Categories
  • Microsoft
    • Microsoft Apps
    • Office
    • Operating System
    • VLS
    • Developer Tools
    • Productivity Tools
    • Database
    • AI + Machine Learning
    • Middleware System
    • Learning Services
    • Analytics
    • Networking
    • Compute
    • Security
    • Internet Of Things
  • Adobe
  • Matlab
  • Google
  • Visual Paradigm
  • WordPress
    • Plugin WP
    • Themes WP
  • Opensource
  • Others
More Categories Less Categories
  • Get Pack
    • Product Category
    • Simple Product
    • Grouped Product
    • Variable Product
    • External Product
  • My Account
    • Download
    • Cart
    • Checkout
    • Login
  • About Us
    • Contact
    • Forum
    • Frequently Questions
    • Privacy Policy
  • Forum
    • News
      • Category
      • News Tag

iconTicket Service Desk

  • My Download
  • Checkout
Application Package Repository Telkom University
All Categories

All Categories

  • Visual Paradigm
  • IBM
  • Adobe
  • Google
  • Matlab
  • Microsoft
    • Microsoft Apps
    • Analytics
    • AI + Machine Learning
    • Compute
    • Database
    • Developer Tools
    • Internet Of Things
    • Learning Services
    • Middleware System
    • Networking
    • Operating System
    • Productivity Tools
    • Security
    • VLS
      • Windows
      • Office
  • Opensource
  • Wordpress
    • Plugin WP
    • Themes WP
  • Others

Search

0 Wishlist

Cart

Menu
  • Home
    • Download Application Package Repository Telkom University
    • Application Package Repository Telkom University
    • Download Official License Telkom University
    • Download Installer Application Pack
    • Product Category
    • Simple Product
    • Grouped Product
    • Variable Product
    • External Product
  • All Pack
    • Microsoft
      • Operating System
      • Productivity Tools
      • Developer Tools
      • Database
      • AI + Machine Learning
      • Middleware System
      • Networking
      • Compute
      • Security
      • Analytics
      • Internet Of Things
      • Learning Services
    • Microsoft Apps
      • VLS
    • Adobe
    • Matlab
    • WordPress
      • Themes WP
      • Plugin WP
    • Google
    • Opensource
    • Others
  • My account
    • Download
    • Get Pack
    • Cart
    • Checkout
  • News
    • Category
    • News Tag
  • Forum
  • About Us
    • Privacy Policy
    • Frequently Questions
    • Contact
Home/News

Category: News

Error setting up  UAV Toolbox support for PX4 Autopilots
Matlab News

Error setting up UAV Toolbox support for PX4 Autopilots

PuTI / 2025-06-22

Hi everyone. When setting up UAV Toolbox support for PX4 Autopilots, I created succesfully the wsl environment and cloned the git repository with 1.14.0 version of PX4 source code. The problem arises when trying to run the next part, which is the ‘bash ./ubuntu.sh’. I reach to an error that says the following. Does anyone know what may be happening and how should I fix this.
Installing PX4 Python3 dependencies
error: externally-managed-environment

× This environment is externally managed
╰─> To install Python packages system-wide, try apt install
python3-xyz, where xyz is the package you are trying to
install.

If you wish to install a non-Debian-packaged Python package,
create a virtual environment using python3 -m venv path/to/venv.
Then use path/to/venv/bin/python and path/to/venv/bin/pip. Make
sure you have python3-full installed.

If you wish to install a non-Debian packaged Python application,
it may be easiest to use pipx install xyz, which will manage a
virtual environment for you. Make sure you have pipx installed.

See /usr/share/doc/python3.12/README.venv for more information.Hi everyone. When setting up UAV Toolbox support for PX4 Autopilots, I created succesfully the wsl environment and cloned the git repository with 1.14.0 version of PX4 source code. The problem arises when trying to run the next part, which is the ‘bash ./ubuntu.sh’. I reach to an error that says the following. Does anyone know what may be happening and how should I fix this.
Installing PX4 Python3 dependencies
error: externally-managed-environment

× This environment is externally managed
╰─> To install Python packages system-wide, try apt install
python3-xyz, where xyz is the package you are trying to
install.

If you wish to install a non-Debian-packaged Python package,
create a virtual environment using python3 -m venv path/to/venv.
Then use path/to/venv/bin/python and path/to/venv/bin/pip. Make
sure you have python3-full installed.

If you wish to install a non-Debian packaged Python application,
it may be easiest to use pipx install xyz, which will manage a
virtual environment for you. Make sure you have pipx installed.

See /usr/share/doc/python3.12/README.venv for more information. Hi everyone. When setting up UAV Toolbox support for PX4 Autopilots, I created succesfully the wsl environment and cloned the git repository with 1.14.0 version of PX4 source code. The problem arises when trying to run the next part, which is the ‘bash ./ubuntu.sh’. I reach to an error that says the following. Does anyone know what may be happening and how should I fix this.
Installing PX4 Python3 dependencies
error: externally-managed-environment

× This environment is externally managed
╰─> To install Python packages system-wide, try apt install
python3-xyz, where xyz is the package you are trying to
install.

If you wish to install a non-Debian-packaged Python package,
create a virtual environment using python3 -m venv path/to/venv.
Then use path/to/venv/bin/python and path/to/venv/bin/pip. Make
sure you have python3-full installed.

If you wish to install a non-Debian packaged Python application,
it may be easiest to use pipx install xyz, which will manage a
virtual environment for you. Make sure you have pipx installed.

See /usr/share/doc/python3.12/README.venv for more information. px4 toolbox MATLAB Answers — New Questions

​

Can not see the variable when double clicks them on workspace
Matlab News

Can not see the variable when double clicks them on workspace

PuTI / 2025-06-22

As the attached picture shows , when I double click a variable on workspace I cannot see the exact number. Thanks for any one who can helps.As the attached picture shows , when I double click a variable on workspace I cannot see the exact number. Thanks for any one who can helps. As the attached picture shows , when I double click a variable on workspace I cannot see the exact number. Thanks for any one who can helps. workspace, variable, variable display MATLAB Answers — New Questions

​

The sine wave comes out as flat as flat line .
Matlab News

The sine wave comes out as flat as flat line .

PuTI / 2025-06-22

When I do the sine wave to scope normally its ok but when i put it into a system it come out as a flat line and the system output is also not the right one .When I do the sine wave to scope normally its ok but when i put it into a system it come out as a flat line and the system output is also not the right one . When I do the sine wave to scope normally its ok but when i put it into a system it come out as a flat line and the system output is also not the right one . sine, wave, scope MATLAB Answers — New Questions

​

modeling an advanced energy management system for electric vehicles using MATLAB/Simulink
Matlab News

modeling an advanced energy management system for electric vehicles using MATLAB/Simulink

PuTI / 2025-06-22

I’m working on modeling an advanced energy management system for electric vehicles using MATLAB/Simulink. My project involves integrating batteries and supercapacitors to manage energy efficiently, where supercapacitors handle high energy demands during acceleration and store energy recovered from braking. Additionally, I’m incorporating a backup power source using fuel and exploring the integration of renewable energy sources like solar power. I’m facing challenges in connecting and configuring these components correctly in Simulink. Can someone guide me on the best practices for setting up this hybrid system and suggest specific blocks or configurations that would work efficiently for this application?I’m working on modeling an advanced energy management system for electric vehicles using MATLAB/Simulink. My project involves integrating batteries and supercapacitors to manage energy efficiently, where supercapacitors handle high energy demands during acceleration and store energy recovered from braking. Additionally, I’m incorporating a backup power source using fuel and exploring the integration of renewable energy sources like solar power. I’m facing challenges in connecting and configuring these components correctly in Simulink. Can someone guide me on the best practices for setting up this hybrid system and suggest specific blocks or configurations that would work efficiently for this application? I’m working on modeling an advanced energy management system for electric vehicles using MATLAB/Simulink. My project involves integrating batteries and supercapacitors to manage energy efficiently, where supercapacitors handle high energy demands during acceleration and store energy recovered from braking. Additionally, I’m incorporating a backup power source using fuel and exploring the integration of renewable energy sources like solar power. I’m facing challenges in connecting and configuring these components correctly in Simulink. Can someone guide me on the best practices for setting up this hybrid system and suggest specific blocks or configurations that would work efficiently for this application? all MATLAB Answers — New Questions

​

Precision lost when combining Int32 integers with single precision numerical numbers
Matlab News

Precision lost when combining Int32 integers with single precision numerical numbers

PuTI / 2025-06-22

I have a column data A composed of Int32 numbers, and another column data B composed of single precision numbers. When I try to put them into one array C, my single precision numbers were botchered into integers.
C = [A, B];
Why is Matlab set up this way? Due to the loss of precision, my final calcualted values are way off. It took me quite some time to find out this is the reason.I have a column data A composed of Int32 numbers, and another column data B composed of single precision numbers. When I try to put them into one array C, my single precision numbers were botchered into integers.
C = [A, B];
Why is Matlab set up this way? Due to the loss of precision, my final calcualted values are way off. It took me quite some time to find out this is the reason. I have a column data A composed of Int32 numbers, and another column data B composed of single precision numbers. When I try to put them into one array C, my single precision numbers were botchered into integers.
C = [A, B];
Why is Matlab set up this way? Due to the loss of precision, my final calcualted values are way off. It took me quite some time to find out this is the reason. int32, array MATLAB Answers — New Questions

​

Waterfall diagram and fft for a vibration of an electric motor
Matlab News

Waterfall diagram and fft for a vibration of an electric motor

PuTI / 2025-06-21

Hello Everyone,

I have a data set: exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.mat
That consist of a test going slowyly from 0 to 3000 RPM of a motor, from which i want to do a waterfall diagramm to check for resonances, but i cant manage to do so. can somebody help?

I have an own code i tried to do but i dont hink its correct. here it is:

clear all
close all
% % % For WorkStation dSpace recorders!
file_name = ‘exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated’;
s=load(file_name);
%%% Select window type for FFT
% a=no window, b=Rectangular, c=Hann, d=Hamming, e=Flattop, f=Blackman-Harris, g=Nuttall, h=Chebyshev
winType = ‘b’;
length=length(s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data);
x=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data(1:length);
%%% Important base recording params
Fs=20000; %%% sampling frequency [Hz]
T=1/Fs;
%%% important base FFT params
fft_cycles = 2 ;
speed_aver_window = 1; %%% in sec
speed_aver_wind_points = speed_aver_window / T;
multipleORHz = 1; %%% 1 – multiply, 2 is Hz
%%% main indexis
index_speed_kistler = 1;
index_torque_an = 2;
index_torque_an_filt = 3;
index_vibr = 4;
index_vibr_filt = 5;
index_iq_act = 7;
index_id_act = 6;
index_speed_sw = 8;
index_speed_sew = 9;
%%% Define staret point of FFT:
time_start_fft =40;
point_start_fft = round(time_start_fft / T);
%%%% Indexes for analysis
index_main_speed = index_speed_kistler;
index_main_torque = index_torque_an_filt;
index_main_vibr = index_vibr;
speed_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_speed).Data(1:length);
torque_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_torque).Data(1:length);
vibr_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_vibr).Data(1:length);
%%% Extraction of currents
Id_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_id_act).Data(1:length);
Iq_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_iq_act).Data(1:length);
%%% Extraction of speed and torque
SEW_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sew).Data(1:length);
Kistler_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_kistler).Data(1:length);
Kistler_torque_an=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an).Data(1:length);
Kistler_torque_an_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an_filt).Data(1:length);
N_pres=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sw).Data(1:length);
vibr=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr).Data(1:length);
vibr_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr_filt).Data(1:length);
%%% plot results of the test
set(gcf,’color’,’white’)
ax1=subplot(5,1,1);
plot(x,Kistler_torque_an,’r’, x,Kistler_torque_an_filt,’b’);
title(‘Torque Kistler’);
xlabel(‘Time [s]’)
ylabel(‘Torque [Nm]’);
grid on
ax2=subplot(5,1,2);
plot(x,N_pres,’r’,x,SEW_speed,’b’,x,Kistler_speed,’g’);
title(‘Speed’);
xlabel(‘Time [s]’)
ylabel(‘Speed [rpm]’);
legend(‘Sensor’,’SEW’,’Kistler’);
grid on
ax3=subplot(5,1,3);
plot(x,vibr,’.- r’, x,vibr_filt,’b’);
title(‘Vibrations’);
xlabel(‘Time [s]’)
ylabel(‘Acceleration [g]’);
grid on
ax4=subplot(5,1,4);
plot(x,Id_act_1,’r’);
title(‘Id act vs ref’);
xlabel(‘Time [s]’)
ylabel(‘Id [A]’);
grid on
ax5=subplot(5,1,5);
plot(x,Iq_act_1,’r’);
title(‘Iq act vs ref’);
xlabel(‘Time [s]’)
ylabel(‘Iq [A]’);
grid on
linkaxes([ax1,ax2,ax3,ax4,ax5],’x’);
%%% Making window to analyze FFT
speed_main_rpm_aver = abs(mean(speed_main(point_start_fft:(point_start_fft + speed_aver_wind_points))));
speed_rps = abs(speed_main_rpm_aver / 60);
period = 1 / speed_rps;
window_sec = fft_cycles * period;
window_poins = round(window_sec * Fs);
%%% Select and build desired window
switch lower(winType)
case ‘a’ % No window (raw data)
win = ones(window_poins,1);
case ‘b’ % Rectangular
win = rectwin(window_poins);
case ‘c’ % Hann
win = hann(window_poins);
case ‘d’ % Hamming
win = hamming(window_poins);
case ‘e’ % Flattop
win = flattopwin(window_poins);
case ‘f’ % Blackman-Harris
win = blackmanharris(window_poins);
case ‘g’ % Nuttall
win = nuttallwin(window_poins);
case ‘h’ % Chebyshev
win = chebwin(window_poins, 60); % 60 dB sidelobe suppression
otherwise
error(‘Unknown winType "%s"’, winType);
end
% Compute mean values for DC removal
torque_main_aver = mean(torque_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
vibr_main_aver = mean(vibr_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
% Preallocate arrays
time_wind = zeros(window_poins, 1);
torque_wind = zeros(window_poins, 1);
vibr_wind = zeros(window_poins, 1);
speed_wind = zeros(window_poins, 1);
% Create windowed signals
j = 1;
for i = point_start_fft:(point_start_fft + window_poins – 1)
time_wind(j) = x(i);
torque_wind(j) = (torque_main(i) – torque_main_aver) * win(j);
vibr_wind(j) = (vibr_main(i) – vibr_main_aver) * win(j);
speed_wind(j) = speed_main(i);
j = j + 1;
end
figure
set(gcf,’color’,’white’)
bx1=subplot(3,1,1);
plot(time_wind,torque_wind,’r’);
title(‘Torque window’);
xlabel(‘Time [s]’)
ylabel(‘Torque [Nm]’);
grid on
bx2=subplot(3,1,2);
plot(time_wind,vibr_wind,’r’);
title(‘Acceleretion window’);
xlabel(‘Time [s]’)
ylabel(‘Acceleration [g]’);
grid on
bx3=subplot(3,1,3);
plot(time_wind,speed_wind,’r’);
title(‘Speed window’);
xlabel(‘Time [s]’)
ylabel(‘Speed [rpm]’);
grid on
linkaxes([bx1,bx2,bx3],’x’);
%%% FFT of Vibrations
L=window_poins;
t = (0:L-1)*T; % Time vector
Y1=fft(vibr_wind,L);
P2 = abs(Y1/L).^2;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
if multipleORHz == 1
freq_ref = speed_main_rpm_aver/60;
freq_plot_lim = 100;
elseif multipleORHz == 2
freq_ref = 1;
freq_plot_lim = 4000;
end
f1 = (Fs*(0:(L/2))/L)/(freq_ref);
figure
set(gcf,’color’,’white’)
bar(f1,P1)
title(‘Single-Sided Amplitude Spectrum of X(t)’)
xlabel(‘f (multiple of mechanical frequency)’)
xlim([0 freq_plot_lim])
ylabel(‘Absolute value of Harmonic VIBRATIONS [g]’)
% — Define cutoff RPM for data selection —
rpm_min = 0;
rpm_max = 4200;
% Find indices corresponding to rpm_min and rpm_max in time vector x
ind_min = find(speed_main >= rpm_min, 1, ‘first’);
ind_max = find(speed_main <= rpm_max, 1, ‘last’);
% Cut signals based on these indices
time_cut = x(ind_min:ind_max);
rpm_cut = speed_main(ind_min:ind_max);
vibr_cut = vibr_main(ind_min:ind_max);
if any(rpm_cut <= 0)
warning(‘RPM data contains non-positive values. Fixing…’);
rpm_cut(rpm_cut <= 0) = NaN; % Or interpolate, or remove
% Option 1: interpolate missing values (if feasible)
rpm_cut = fillmissing(rpm_cut, ‘linear’);
% Option 2: truncate all rows with NaNs (if that’s acceptable)
validIdx = ~isnan(rpm_cut);
vibr_cut = vibr_cut(validIdx);
rpm_cut = rpm_cut(validIdx);
end
% Sampling frequency from original code
Fs = 4000; % Hz
% % —- Generate order map and waterfall plot —-
% Assuming you have the rpmordermap function in your path.
% If not, I can help to replace it with an equivalent.
% Frequency-based RPM map
[map, freq, rpm_axis, time_map, res] = rpmordermap(vibr_cut, Fs, rpm_cut, 2, …
‘Scale’, ‘dB’, ‘Window’, ‘hann’, ‘Amplitude’, ‘rms’);
[fr, rp] = meshgrid(freq, rpm_axis);
% Waterfall plot in frequency
figure;
waterfall(fr, rp, map’);
view(6, -60);
grid on;
xlabel(‘Frequency [Hz]’);
ylabel(‘RPM’);
zlabel(‘Amplitude [dB]’);
title(‘Waterfall Plot (Frequency Map)’);
% — FFT over speed steps (similar to your reference code) —
% Define speed steps (adjust based on your rpm range)
speed_steps = rpm_min:100:rpm_max;
% Define FFT window length in seconds (e.g., 0.5 s)
fft_window_sec = 0.5;
fft_window_points = round(fft_window_sec * Fs);
% Initialize figure for FFT results
for i = 1:length(speed_steps)
% Find start time index closest to current speed step
idx_start = find(rpm_cut >= speed_steps(i), 1, ‘first’);
if isempty(idx_start) || (idx_start + fft_window_points – 1) > length(vibr_cut)
continue; % Skip if index invalid or window exceeds data length
end
% Extract segment
segment = vibr_cut(idx_start : idx_start + fft_window_points – 1);
time_segment = time_cut(idx_start : idx_start + fft_window_points – 1);
% Remove DC
segment = segment – mean(segment);
% Apply Hann window
winvec = hann(length(segment));
% FFT
L = length(segment);
Y = fft(segment .* winvec);
P2 = abs(Y / L);
P1 = P2(1 : floor(L/2) + 1);
P1(2:end-1) = 2 * P1(2:end-1);
f = Fs * (0:(L/2)) / L;
% Plot time-domain and FFT for this step
subplot(2,1,1);
plot(time_segment, segment);
title(sprintf(‘Vibration at Speed Step: %d RPM’, speed_steps(i)));
xlabel(‘Time [s]’);
ylabel(‘Acceleration [g]’);
grid on;
hold on;
subplot(2,1,2);
bar(f, P1);
title(sprintf(‘FFT Spectrum at Speed Step: %d RPM’, speed_steps(i)));
xlabel(‘Frequency [Hz]’);
ylabel(‘Amplitude’);
grid on;
hold on;
end

code should be correct until the part of the water fall.

also the index of the data is no longer 1 to 9 but 1 to 3 being speed-raw vibration- filtered vibration

I would appreciete if somebody could help me make a code i can use with different data sets, maybe some have 1 to 6 and with the option of ploting the sigbals at the beginning and the fft how in the code.

Thanks in advancedHello Everyone,

I have a data set: exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.mat
That consist of a test going slowyly from 0 to 3000 RPM of a motor, from which i want to do a waterfall diagramm to check for resonances, but i cant manage to do so. can somebody help?

I have an own code i tried to do but i dont hink its correct. here it is:

clear all
close all
% % % For WorkStation dSpace recorders!
file_name = ‘exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated’;
s=load(file_name);
%%% Select window type for FFT
% a=no window, b=Rectangular, c=Hann, d=Hamming, e=Flattop, f=Blackman-Harris, g=Nuttall, h=Chebyshev
winType = ‘b’;
length=length(s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data);
x=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data(1:length);
%%% Important base recording params
Fs=20000; %%% sampling frequency [Hz]
T=1/Fs;
%%% important base FFT params
fft_cycles = 2 ;
speed_aver_window = 1; %%% in sec
speed_aver_wind_points = speed_aver_window / T;
multipleORHz = 1; %%% 1 – multiply, 2 is Hz
%%% main indexis
index_speed_kistler = 1;
index_torque_an = 2;
index_torque_an_filt = 3;
index_vibr = 4;
index_vibr_filt = 5;
index_iq_act = 7;
index_id_act = 6;
index_speed_sw = 8;
index_speed_sew = 9;
%%% Define staret point of FFT:
time_start_fft =40;
point_start_fft = round(time_start_fft / T);
%%%% Indexes for analysis
index_main_speed = index_speed_kistler;
index_main_torque = index_torque_an_filt;
index_main_vibr = index_vibr;
speed_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_speed).Data(1:length);
torque_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_torque).Data(1:length);
vibr_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_vibr).Data(1:length);
%%% Extraction of currents
Id_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_id_act).Data(1:length);
Iq_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_iq_act).Data(1:length);
%%% Extraction of speed and torque
SEW_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sew).Data(1:length);
Kistler_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_kistler).Data(1:length);
Kistler_torque_an=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an).Data(1:length);
Kistler_torque_an_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an_filt).Data(1:length);
N_pres=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sw).Data(1:length);
vibr=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr).Data(1:length);
vibr_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr_filt).Data(1:length);
%%% plot results of the test
set(gcf,’color’,’white’)
ax1=subplot(5,1,1);
plot(x,Kistler_torque_an,’r’, x,Kistler_torque_an_filt,’b’);
title(‘Torque Kistler’);
xlabel(‘Time [s]’)
ylabel(‘Torque [Nm]’);
grid on
ax2=subplot(5,1,2);
plot(x,N_pres,’r’,x,SEW_speed,’b’,x,Kistler_speed,’g’);
title(‘Speed’);
xlabel(‘Time [s]’)
ylabel(‘Speed [rpm]’);
legend(‘Sensor’,’SEW’,’Kistler’);
grid on
ax3=subplot(5,1,3);
plot(x,vibr,’.- r’, x,vibr_filt,’b’);
title(‘Vibrations’);
xlabel(‘Time [s]’)
ylabel(‘Acceleration [g]’);
grid on
ax4=subplot(5,1,4);
plot(x,Id_act_1,’r’);
title(‘Id act vs ref’);
xlabel(‘Time [s]’)
ylabel(‘Id [A]’);
grid on
ax5=subplot(5,1,5);
plot(x,Iq_act_1,’r’);
title(‘Iq act vs ref’);
xlabel(‘Time [s]’)
ylabel(‘Iq [A]’);
grid on
linkaxes([ax1,ax2,ax3,ax4,ax5],’x’);
%%% Making window to analyze FFT
speed_main_rpm_aver = abs(mean(speed_main(point_start_fft:(point_start_fft + speed_aver_wind_points))));
speed_rps = abs(speed_main_rpm_aver / 60);
period = 1 / speed_rps;
window_sec = fft_cycles * period;
window_poins = round(window_sec * Fs);
%%% Select and build desired window
switch lower(winType)
case ‘a’ % No window (raw data)
win = ones(window_poins,1);
case ‘b’ % Rectangular
win = rectwin(window_poins);
case ‘c’ % Hann
win = hann(window_poins);
case ‘d’ % Hamming
win = hamming(window_poins);
case ‘e’ % Flattop
win = flattopwin(window_poins);
case ‘f’ % Blackman-Harris
win = blackmanharris(window_poins);
case ‘g’ % Nuttall
win = nuttallwin(window_poins);
case ‘h’ % Chebyshev
win = chebwin(window_poins, 60); % 60 dB sidelobe suppression
otherwise
error(‘Unknown winType "%s"’, winType);
end
% Compute mean values for DC removal
torque_main_aver = mean(torque_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
vibr_main_aver = mean(vibr_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
% Preallocate arrays
time_wind = zeros(window_poins, 1);
torque_wind = zeros(window_poins, 1);
vibr_wind = zeros(window_poins, 1);
speed_wind = zeros(window_poins, 1);
% Create windowed signals
j = 1;
for i = point_start_fft:(point_start_fft + window_poins – 1)
time_wind(j) = x(i);
torque_wind(j) = (torque_main(i) – torque_main_aver) * win(j);
vibr_wind(j) = (vibr_main(i) – vibr_main_aver) * win(j);
speed_wind(j) = speed_main(i);
j = j + 1;
end
figure
set(gcf,’color’,’white’)
bx1=subplot(3,1,1);
plot(time_wind,torque_wind,’r’);
title(‘Torque window’);
xlabel(‘Time [s]’)
ylabel(‘Torque [Nm]’);
grid on
bx2=subplot(3,1,2);
plot(time_wind,vibr_wind,’r’);
title(‘Acceleretion window’);
xlabel(‘Time [s]’)
ylabel(‘Acceleration [g]’);
grid on
bx3=subplot(3,1,3);
plot(time_wind,speed_wind,’r’);
title(‘Speed window’);
xlabel(‘Time [s]’)
ylabel(‘Speed [rpm]’);
grid on
linkaxes([bx1,bx2,bx3],’x’);
%%% FFT of Vibrations
L=window_poins;
t = (0:L-1)*T; % Time vector
Y1=fft(vibr_wind,L);
P2 = abs(Y1/L).^2;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
if multipleORHz == 1
freq_ref = speed_main_rpm_aver/60;
freq_plot_lim = 100;
elseif multipleORHz == 2
freq_ref = 1;
freq_plot_lim = 4000;
end
f1 = (Fs*(0:(L/2))/L)/(freq_ref);
figure
set(gcf,’color’,’white’)
bar(f1,P1)
title(‘Single-Sided Amplitude Spectrum of X(t)’)
xlabel(‘f (multiple of mechanical frequency)’)
xlim([0 freq_plot_lim])
ylabel(‘Absolute value of Harmonic VIBRATIONS [g]’)
% — Define cutoff RPM for data selection —
rpm_min = 0;
rpm_max = 4200;
% Find indices corresponding to rpm_min and rpm_max in time vector x
ind_min = find(speed_main >= rpm_min, 1, ‘first’);
ind_max = find(speed_main <= rpm_max, 1, ‘last’);
% Cut signals based on these indices
time_cut = x(ind_min:ind_max);
rpm_cut = speed_main(ind_min:ind_max);
vibr_cut = vibr_main(ind_min:ind_max);
if any(rpm_cut <= 0)
warning(‘RPM data contains non-positive values. Fixing…’);
rpm_cut(rpm_cut <= 0) = NaN; % Or interpolate, or remove
% Option 1: interpolate missing values (if feasible)
rpm_cut = fillmissing(rpm_cut, ‘linear’);
% Option 2: truncate all rows with NaNs (if that’s acceptable)
validIdx = ~isnan(rpm_cut);
vibr_cut = vibr_cut(validIdx);
rpm_cut = rpm_cut(validIdx);
end
% Sampling frequency from original code
Fs = 4000; % Hz
% % —- Generate order map and waterfall plot —-
% Assuming you have the rpmordermap function in your path.
% If not, I can help to replace it with an equivalent.
% Frequency-based RPM map
[map, freq, rpm_axis, time_map, res] = rpmordermap(vibr_cut, Fs, rpm_cut, 2, …
‘Scale’, ‘dB’, ‘Window’, ‘hann’, ‘Amplitude’, ‘rms’);
[fr, rp] = meshgrid(freq, rpm_axis);
% Waterfall plot in frequency
figure;
waterfall(fr, rp, map’);
view(6, -60);
grid on;
xlabel(‘Frequency [Hz]’);
ylabel(‘RPM’);
zlabel(‘Amplitude [dB]’);
title(‘Waterfall Plot (Frequency Map)’);
% — FFT over speed steps (similar to your reference code) —
% Define speed steps (adjust based on your rpm range)
speed_steps = rpm_min:100:rpm_max;
% Define FFT window length in seconds (e.g., 0.5 s)
fft_window_sec = 0.5;
fft_window_points = round(fft_window_sec * Fs);
% Initialize figure for FFT results
for i = 1:length(speed_steps)
% Find start time index closest to current speed step
idx_start = find(rpm_cut >= speed_steps(i), 1, ‘first’);
if isempty(idx_start) || (idx_start + fft_window_points – 1) > length(vibr_cut)
continue; % Skip if index invalid or window exceeds data length
end
% Extract segment
segment = vibr_cut(idx_start : idx_start + fft_window_points – 1);
time_segment = time_cut(idx_start : idx_start + fft_window_points – 1);
% Remove DC
segment = segment – mean(segment);
% Apply Hann window
winvec = hann(length(segment));
% FFT
L = length(segment);
Y = fft(segment .* winvec);
P2 = abs(Y / L);
P1 = P2(1 : floor(L/2) + 1);
P1(2:end-1) = 2 * P1(2:end-1);
f = Fs * (0:(L/2)) / L;
% Plot time-domain and FFT for this step
subplot(2,1,1);
plot(time_segment, segment);
title(sprintf(‘Vibration at Speed Step: %d RPM’, speed_steps(i)));
xlabel(‘Time [s]’);
ylabel(‘Acceleration [g]’);
grid on;
hold on;
subplot(2,1,2);
bar(f, P1);
title(sprintf(‘FFT Spectrum at Speed Step: %d RPM’, speed_steps(i)));
xlabel(‘Frequency [Hz]’);
ylabel(‘Amplitude’);
grid on;
hold on;
end

code should be correct until the part of the water fall.

also the index of the data is no longer 1 to 9 but 1 to 3 being speed-raw vibration- filtered vibration

I would appreciete if somebody could help me make a code i can use with different data sets, maybe some have 1 to 6 and with the option of ploting the sigbals at the beginning and the fft how in the code.

Thanks in advanced Hello Everyone,

I have a data set: exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.mat
That consist of a test going slowyly from 0 to 3000 RPM of a motor, from which i want to do a waterfall diagramm to check for resonances, but i cant manage to do so. can somebody help?

I have an own code i tried to do but i dont hink its correct. here it is:

clear all
close all
% % % For WorkStation dSpace recorders!
file_name = ‘exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated’;
s=load(file_name);
%%% Select window type for FFT
% a=no window, b=Rectangular, c=Hann, d=Hamming, e=Flattop, f=Blackman-Harris, g=Nuttall, h=Chebyshev
winType = ‘b’;
length=length(s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data);
x=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.X(1).Data(1:length);
%%% Important base recording params
Fs=20000; %%% sampling frequency [Hz]
T=1/Fs;
%%% important base FFT params
fft_cycles = 2 ;
speed_aver_window = 1; %%% in sec
speed_aver_wind_points = speed_aver_window / T;
multipleORHz = 1; %%% 1 – multiply, 2 is Hz
%%% main indexis
index_speed_kistler = 1;
index_torque_an = 2;
index_torque_an_filt = 3;
index_vibr = 4;
index_vibr_filt = 5;
index_iq_act = 7;
index_id_act = 6;
index_speed_sw = 8;
index_speed_sew = 9;
%%% Define staret point of FFT:
time_start_fft =40;
point_start_fft = round(time_start_fft / T);
%%%% Indexes for analysis
index_main_speed = index_speed_kistler;
index_main_torque = index_torque_an_filt;
index_main_vibr = index_vibr;
speed_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_speed).Data(1:length);
torque_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_torque).Data(1:length);
vibr_main = s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_main_vibr).Data(1:length);
%%% Extraction of currents
Id_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_id_act).Data(1:length);
Iq_act_1=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_iq_act).Data(1:length);
%%% Extraction of speed and torque
SEW_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sew).Data(1:length);
Kistler_speed=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_kistler).Data(1:length);
Kistler_torque_an=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an).Data(1:length);
Kistler_torque_an_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_torque_an_filt).Data(1:length);
N_pres=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_speed_sw).Data(1:length);
vibr=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr).Data(1:length);
vibr_filt=s.exp1_acoustic_cpm150A0_ASQ_N7_30x100rpmx2sec_20kHz_2kHz_LPF_automated.Y(index_vibr_filt).Data(1:length);
%%% plot results of the test
set(gcf,’color’,’white’)
ax1=subplot(5,1,1);
plot(x,Kistler_torque_an,’r’, x,Kistler_torque_an_filt,’b’);
title(‘Torque Kistler’);
xlabel(‘Time [s]’)
ylabel(‘Torque [Nm]’);
grid on
ax2=subplot(5,1,2);
plot(x,N_pres,’r’,x,SEW_speed,’b’,x,Kistler_speed,’g’);
title(‘Speed’);
xlabel(‘Time [s]’)
ylabel(‘Speed [rpm]’);
legend(‘Sensor’,’SEW’,’Kistler’);
grid on
ax3=subplot(5,1,3);
plot(x,vibr,’.- r’, x,vibr_filt,’b’);
title(‘Vibrations’);
xlabel(‘Time [s]’)
ylabel(‘Acceleration [g]’);
grid on
ax4=subplot(5,1,4);
plot(x,Id_act_1,’r’);
title(‘Id act vs ref’);
xlabel(‘Time [s]’)
ylabel(‘Id [A]’);
grid on
ax5=subplot(5,1,5);
plot(x,Iq_act_1,’r’);
title(‘Iq act vs ref’);
xlabel(‘Time [s]’)
ylabel(‘Iq [A]’);
grid on
linkaxes([ax1,ax2,ax3,ax4,ax5],’x’);
%%% Making window to analyze FFT
speed_main_rpm_aver = abs(mean(speed_main(point_start_fft:(point_start_fft + speed_aver_wind_points))));
speed_rps = abs(speed_main_rpm_aver / 60);
period = 1 / speed_rps;
window_sec = fft_cycles * period;
window_poins = round(window_sec * Fs);
%%% Select and build desired window
switch lower(winType)
case ‘a’ % No window (raw data)
win = ones(window_poins,1);
case ‘b’ % Rectangular
win = rectwin(window_poins);
case ‘c’ % Hann
win = hann(window_poins);
case ‘d’ % Hamming
win = hamming(window_poins);
case ‘e’ % Flattop
win = flattopwin(window_poins);
case ‘f’ % Blackman-Harris
win = blackmanharris(window_poins);
case ‘g’ % Nuttall
win = nuttallwin(window_poins);
case ‘h’ % Chebyshev
win = chebwin(window_poins, 60); % 60 dB sidelobe suppression
otherwise
error(‘Unknown winType "%s"’, winType);
end
% Compute mean values for DC removal
torque_main_aver = mean(torque_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
vibr_main_aver = mean(vibr_main(point_start_fft:(point_start_fft + speed_aver_wind_points)));
% Preallocate arrays
time_wind = zeros(window_poins, 1);
torque_wind = zeros(window_poins, 1);
vibr_wind = zeros(window_poins, 1);
speed_wind = zeros(window_poins, 1);
% Create windowed signals
j = 1;
for i = point_start_fft:(point_start_fft + window_poins – 1)
time_wind(j) = x(i);
torque_wind(j) = (torque_main(i) – torque_main_aver) * win(j);
vibr_wind(j) = (vibr_main(i) – vibr_main_aver) * win(j);
speed_wind(j) = speed_main(i);
j = j + 1;
end
figure
set(gcf,’color’,’white’)
bx1=subplot(3,1,1);
plot(time_wind,torque_wind,’r’);
title(‘Torque window’);
xlabel(‘Time [s]’)
ylabel(‘Torque [Nm]’);
grid on
bx2=subplot(3,1,2);
plot(time_wind,vibr_wind,’r’);
title(‘Acceleretion window’);
xlabel(‘Time [s]’)
ylabel(‘Acceleration [g]’);
grid on
bx3=subplot(3,1,3);
plot(time_wind,speed_wind,’r’);
title(‘Speed window’);
xlabel(‘Time [s]’)
ylabel(‘Speed [rpm]’);
grid on
linkaxes([bx1,bx2,bx3],’x’);
%%% FFT of Vibrations
L=window_poins;
t = (0:L-1)*T; % Time vector
Y1=fft(vibr_wind,L);
P2 = abs(Y1/L).^2;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
if multipleORHz == 1
freq_ref = speed_main_rpm_aver/60;
freq_plot_lim = 100;
elseif multipleORHz == 2
freq_ref = 1;
freq_plot_lim = 4000;
end
f1 = (Fs*(0:(L/2))/L)/(freq_ref);
figure
set(gcf,’color’,’white’)
bar(f1,P1)
title(‘Single-Sided Amplitude Spectrum of X(t)’)
xlabel(‘f (multiple of mechanical frequency)’)
xlim([0 freq_plot_lim])
ylabel(‘Absolute value of Harmonic VIBRATIONS [g]’)
% — Define cutoff RPM for data selection —
rpm_min = 0;
rpm_max = 4200;
% Find indices corresponding to rpm_min and rpm_max in time vector x
ind_min = find(speed_main >= rpm_min, 1, ‘first’);
ind_max = find(speed_main <= rpm_max, 1, ‘last’);
% Cut signals based on these indices
time_cut = x(ind_min:ind_max);
rpm_cut = speed_main(ind_min:ind_max);
vibr_cut = vibr_main(ind_min:ind_max);
if any(rpm_cut <= 0)
warning(‘RPM data contains non-positive values. Fixing…’);
rpm_cut(rpm_cut <= 0) = NaN; % Or interpolate, or remove
% Option 1: interpolate missing values (if feasible)
rpm_cut = fillmissing(rpm_cut, ‘linear’);
% Option 2: truncate all rows with NaNs (if that’s acceptable)
validIdx = ~isnan(rpm_cut);
vibr_cut = vibr_cut(validIdx);
rpm_cut = rpm_cut(validIdx);
end
% Sampling frequency from original code
Fs = 4000; % Hz
% % —- Generate order map and waterfall plot —-
% Assuming you have the rpmordermap function in your path.
% If not, I can help to replace it with an equivalent.
% Frequency-based RPM map
[map, freq, rpm_axis, time_map, res] = rpmordermap(vibr_cut, Fs, rpm_cut, 2, …
‘Scale’, ‘dB’, ‘Window’, ‘hann’, ‘Amplitude’, ‘rms’);
[fr, rp] = meshgrid(freq, rpm_axis);
% Waterfall plot in frequency
figure;
waterfall(fr, rp, map’);
view(6, -60);
grid on;
xlabel(‘Frequency [Hz]’);
ylabel(‘RPM’);
zlabel(‘Amplitude [dB]’);
title(‘Waterfall Plot (Frequency Map)’);
% — FFT over speed steps (similar to your reference code) —
% Define speed steps (adjust based on your rpm range)
speed_steps = rpm_min:100:rpm_max;
% Define FFT window length in seconds (e.g., 0.5 s)
fft_window_sec = 0.5;
fft_window_points = round(fft_window_sec * Fs);
% Initialize figure for FFT results
for i = 1:length(speed_steps)
% Find start time index closest to current speed step
idx_start = find(rpm_cut >= speed_steps(i), 1, ‘first’);
if isempty(idx_start) || (idx_start + fft_window_points – 1) > length(vibr_cut)
continue; % Skip if index invalid or window exceeds data length
end
% Extract segment
segment = vibr_cut(idx_start : idx_start + fft_window_points – 1);
time_segment = time_cut(idx_start : idx_start + fft_window_points – 1);
% Remove DC
segment = segment – mean(segment);
% Apply Hann window
winvec = hann(length(segment));
% FFT
L = length(segment);
Y = fft(segment .* winvec);
P2 = abs(Y / L);
P1 = P2(1 : floor(L/2) + 1);
P1(2:end-1) = 2 * P1(2:end-1);
f = Fs * (0:(L/2)) / L;
% Plot time-domain and FFT for this step
subplot(2,1,1);
plot(time_segment, segment);
title(sprintf(‘Vibration at Speed Step: %d RPM’, speed_steps(i)));
xlabel(‘Time [s]’);
ylabel(‘Acceleration [g]’);
grid on;
hold on;
subplot(2,1,2);
bar(f, P1);
title(sprintf(‘FFT Spectrum at Speed Step: %d RPM’, speed_steps(i)));
xlabel(‘Frequency [Hz]’);
ylabel(‘Amplitude’);
grid on;
hold on;
end

code should be correct until the part of the water fall.

also the index of the data is no longer 1 to 9 but 1 to 3 being speed-raw vibration- filtered vibration

I would appreciete if somebody could help me make a code i can use with different data sets, maybe some have 1 to 6 and with the option of ploting the sigbals at the beginning and the fft how in the code.

Thanks in advanced matlab, fft, waterfalldiagramm, campbelldiagramm, electric_motor_control, pmsm MATLAB Answers — New Questions

​

Why do I receive License Manager Error -9?
Matlab News

Why do I receive License Manager Error -9?

PuTI / 2025-06-21

Why do I receive License Manager Error -9?Why do I receive License Manager Error -9? Why do I receive License Manager Error -9? error -9 MATLAB Answers — New Questions

​

Synchronizing the data of 2 subdevices within 1 device.
Matlab News

Synchronizing the data of 2 subdevices within 1 device.

PuTI / 2025-06-21

I am currently trying to figure out how to Sync up my IMU data. Within this IMU there is an Accelerometer and Gyroscope that get data with 200Hz and an GNSS that gets data with 10Hz. The ultimate goal is using the script(when its done) for event detection. Within the GNSS data I get is a POSIX time and a ‘timestamp[us]’ in the IMU data is only the ‘timestamp[us]. Does anyone have an idea on how I could manage to sync these two time variables so I can view real time events within matlab eventhough they have different frequencies.I am currently trying to figure out how to Sync up my IMU data. Within this IMU there is an Accelerometer and Gyroscope that get data with 200Hz and an GNSS that gets data with 10Hz. The ultimate goal is using the script(when its done) for event detection. Within the GNSS data I get is a POSIX time and a ‘timestamp[us]’ in the IMU data is only the ‘timestamp[us]. Does anyone have an idea on how I could manage to sync these two time variables so I can view real time events within matlab eventhough they have different frequencies. I am currently trying to figure out how to Sync up my IMU data. Within this IMU there is an Accelerometer and Gyroscope that get data with 200Hz and an GNSS that gets data with 10Hz. The ultimate goal is using the script(when its done) for event detection. Within the GNSS data I get is a POSIX time and a ‘timestamp[us]’ in the IMU data is only the ‘timestamp[us]. Does anyone have an idea on how I could manage to sync these two time variables so I can view real time events within matlab eventhough they have different frequencies. data synchronization MATLAB Answers — New Questions

​

importNetworkFromONNX did not recognize softmax layer
Matlab News

importNetworkFromONNX did not recognize softmax layer

PuTI / 2025-06-21

Hello,
I am using importNetworkFromONNX to import a neural network exported from pyTorch.
The network includes a softmax layer. (onnx_model.onnx file is attached as a zip file.)
Below picture shows the model’s netron view:

However, when I imported the onnx model, MATLAB did not recognize the softmax layer.
I know that I can relace the layer with MATLAB’s Softmax layer.
But, I want to know how to import the onnx model without replacing the layer.
Below is the code I used to import the onnx model.
clear
modelfile = "onnx_model.onnx";
localNet = importNetworkFromONNX(modelfile, InputDataFormats=’BC’);
InputSize = [4 1];
X = dlarray(rand(InputSize),"CB");
localNet = initialize(localNet, X);
localNet = expandLayers(localNet);
localNet.Layers

The results was:
>> test_import_onnx
ans =
9×1 Layer array with layers:

1 ‘onnx__Gemm_0’ Feature Input 4 features
2 ‘onnx__Gemm_0_BatchSizeVerifier’ Verify the fixed batch size Verify the fixed batch size of 1
3 ‘x_fc1_pi_Gemm’ Fully Connected 5 fully connected layer
4 ‘x_Relu’ ReLU ReLU
5 ‘x_fc2_pi_Gemm’ Fully Connected 2 fully connected layer
6 ‘SoftmaxLayer1003’ onnx_model.SoftmaxLayer1003 onnx_model.SoftmaxLayer1003
7 ‘x_fc1_v_Gemm’ Fully Connected 1 fully connected layer
8 ‘x11Output’ Custom output (‘CB’) See the OutputInformation property to find the output dimension ordering that is produced by this layer.
9 ‘x12Output’ Custom output (‘CB’) See the OutputInformation property to find the output dimension ordering that is produced by this layer.
Because onnx_model.SoftmaxLayer1003 did not work as a softmax layer, the outputs of SoftmaxLayer1003 were always [1; 1].Hello,
I am using importNetworkFromONNX to import a neural network exported from pyTorch.
The network includes a softmax layer. (onnx_model.onnx file is attached as a zip file.)
Below picture shows the model’s netron view:

However, when I imported the onnx model, MATLAB did not recognize the softmax layer.
I know that I can relace the layer with MATLAB’s Softmax layer.
But, I want to know how to import the onnx model without replacing the layer.
Below is the code I used to import the onnx model.
clear
modelfile = "onnx_model.onnx";
localNet = importNetworkFromONNX(modelfile, InputDataFormats=’BC’);
InputSize = [4 1];
X = dlarray(rand(InputSize),"CB");
localNet = initialize(localNet, X);
localNet = expandLayers(localNet);
localNet.Layers

The results was:
>> test_import_onnx
ans =
9×1 Layer array with layers:

1 ‘onnx__Gemm_0’ Feature Input 4 features
2 ‘onnx__Gemm_0_BatchSizeVerifier’ Verify the fixed batch size Verify the fixed batch size of 1
3 ‘x_fc1_pi_Gemm’ Fully Connected 5 fully connected layer
4 ‘x_Relu’ ReLU ReLU
5 ‘x_fc2_pi_Gemm’ Fully Connected 2 fully connected layer
6 ‘SoftmaxLayer1003’ onnx_model.SoftmaxLayer1003 onnx_model.SoftmaxLayer1003
7 ‘x_fc1_v_Gemm’ Fully Connected 1 fully connected layer
8 ‘x11Output’ Custom output (‘CB’) See the OutputInformation property to find the output dimension ordering that is produced by this layer.
9 ‘x12Output’ Custom output (‘CB’) See the OutputInformation property to find the output dimension ordering that is produced by this layer.
Because onnx_model.SoftmaxLayer1003 did not work as a softmax layer, the outputs of SoftmaxLayer1003 were always [1; 1]. Hello,
I am using importNetworkFromONNX to import a neural network exported from pyTorch.
The network includes a softmax layer. (onnx_model.onnx file is attached as a zip file.)
Below picture shows the model’s netron view:

However, when I imported the onnx model, MATLAB did not recognize the softmax layer.
I know that I can relace the layer with MATLAB’s Softmax layer.
But, I want to know how to import the onnx model without replacing the layer.
Below is the code I used to import the onnx model.
clear
modelfile = "onnx_model.onnx";
localNet = importNetworkFromONNX(modelfile, InputDataFormats=’BC’);
InputSize = [4 1];
X = dlarray(rand(InputSize),"CB");
localNet = initialize(localNet, X);
localNet = expandLayers(localNet);
localNet.Layers

The results was:
>> test_import_onnx
ans =
9×1 Layer array with layers:

1 ‘onnx__Gemm_0’ Feature Input 4 features
2 ‘onnx__Gemm_0_BatchSizeVerifier’ Verify the fixed batch size Verify the fixed batch size of 1
3 ‘x_fc1_pi_Gemm’ Fully Connected 5 fully connected layer
4 ‘x_Relu’ ReLU ReLU
5 ‘x_fc2_pi_Gemm’ Fully Connected 2 fully connected layer
6 ‘SoftmaxLayer1003’ onnx_model.SoftmaxLayer1003 onnx_model.SoftmaxLayer1003
7 ‘x_fc1_v_Gemm’ Fully Connected 1 fully connected layer
8 ‘x11Output’ Custom output (‘CB’) See the OutputInformation property to find the output dimension ordering that is produced by this layer.
9 ‘x12Output’ Custom output (‘CB’) See the OutputInformation property to find the output dimension ordering that is produced by this layer.
Because onnx_model.SoftmaxLayer1003 did not work as a softmax layer, the outputs of SoftmaxLayer1003 were always [1; 1]. importnetworkfromonnx, deep learning toolbox, onnx MATLAB Answers — New Questions

​

Issues in fitting a simple tumor growth model
Matlab News

Issues in fitting a simple tumor growth model

PuTI / 2025-06-21

Dear staff,
I’m learning how to implement a TGI model in Simbiology using a scripting approach. For that purpose I started to implement my model (based on the example tumor_growth_fitPKPD_completed.sbproj) just with a modified growth component and I’m trying to fit a dataset. Attached you can find the dataset (synth_data.mat) and the TGI script (test_code.m). As you can see the initial tumor mass of different patients varies quite a lot so I can’t use a unique initialization for the variable w0. I would like to understand if (and how) it is possible to pass to the model the initial tumor mass for each ID, in order to avoid to fit this parameter.
Thank you very much for your help.Dear staff,
I’m learning how to implement a TGI model in Simbiology using a scripting approach. For that purpose I started to implement my model (based on the example tumor_growth_fitPKPD_completed.sbproj) just with a modified growth component and I’m trying to fit a dataset. Attached you can find the dataset (synth_data.mat) and the TGI script (test_code.m). As you can see the initial tumor mass of different patients varies quite a lot so I can’t use a unique initialization for the variable w0. I would like to understand if (and how) it is possible to pass to the model the initial tumor mass for each ID, in order to avoid to fit this parameter.
Thank you very much for your help. Dear staff,
I’m learning how to implement a TGI model in Simbiology using a scripting approach. For that purpose I started to implement my model (based on the example tumor_growth_fitPKPD_completed.sbproj) just with a modified growth component and I’m trying to fit a dataset. Attached you can find the dataset (synth_data.mat) and the TGI script (test_code.m). As you can see the initial tumor mass of different patients varies quite a lot so I can’t use a unique initialization for the variable w0. I would like to understand if (and how) it is possible to pass to the model the initial tumor mass for each ID, in order to avoid to fit this parameter.
Thank you very much for your help. simbiology MATLAB Answers — New Questions

​

cooperative logistics – profit maximization
Matlab News

cooperative logistics – profit maximization

PuTI / 2025-06-20

I need help about this code dealing with a problem of cooperative logistics. I have to maximize profits, that is, demonstrate that with the three carriers together you get a profit greater than the sum of the profits given by the three carriers taken individually. MatLab takes over 30 minutes and gives infeasible. Thank you so much!
clc
clear all

%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%

Nr=20; % number of trips
Nc=3; % number of carriers
Vr=22; % number of vehicles

alpha=1; % extra charge on revenue;
cE= 0.2;
beta=0.8;
mu=80;

% dummy quantity
M=1000000;

% data about trips 20
d= [85, 80, 123, 111, 153, 112, 150, 107, 181, 135, 109, 106, 113, 148, 164, 127, 148, 141, 108, 118];
delta=[200, 200, 200, 200, 200, 200, 200, 200, 300, 300, 300, 300, 300, 300, 300, 400, 400, 400, 400, 400];
p=[25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25];

% data about carriers, INVENTATI
I=[3, 3, 3];
cc=[0.8, 1.0, 0.9];

% data about vehicles, 22
T= [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8];
CAP= [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30];
Eempty= [0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.2, 0.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.3, 0.2, 0.2, 0.3, 0.3, 0.2];
Efull= [0.3, 0.4, 0.3, 0.3, 0.3, 0.4, 0.3, 0.4, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.4, 0.4, 0.3, 0.4, 0.4]; % it is already Efull-Eempty
IsE=[0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0]; %1 if the vehicle is electric
Owner= [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3];

% repositioing distances between trips 20×20
eps= [14 16 6 19 26 27 21 10 16 18 29 7 30 27 28 10 8 11 17 10;
27 9 20 11 9 27 22 18 27 14 29 22 5 25 9 21 17 18 11 16;
25 11 10 30 19 14 7 9 19 26 27 22 29 7 12 29 24 7 9 30;
13 16 6 19 7 7 24 5 27 28 6 18 27 9 17 30 24 14 18 10;
22 23 21 13 16 16 25 21 14 8 14 13 18 17 11 7 16 29 6 6;
5 21 9 17 26 30 27 11 24 28 16 6 25 22 27 11 18 30 13 7;
13 11 15 24 20 30 26 11 7 21 7 17 29 5 11 12 22 22 8 5;
28 7 14 20 7 25 20 12 15 5 27 10 5 27 20 26 22 25 14 28;
16 19 11 12 8 25 24 9 26 26 17 15 8 5 20 19 7 21 28 24;
29 8 13 18 17 13 15 12 21 11 10 25 17 29 22 11 28 27 14 13;
10 22 26 16 22 19 7 28 10 14 12 28 27 10 12 6 23 30 28 26;
18 11 15 25 20 7 21 17 21 9 16 5 11 5 18 16 17 26 21 19;
14 15 25 19 11 16 8 25 21 9 11 8 13 5 10 7 5 28 30 8;
26 15 11 19 21 23 9 24 9 26 16 17 9 27 25 27 20 21 15 13;
27 27 14 19 28 9 22 24 18 11 22 15 15 26 27 27 27 19 17 19;
25 11 17 24 6 30 15 23 21 30 16 5 18 25 12 20 9 11 25 27;
29 9 17 11 24 8 29 28 21 16 20 15 6 28 8 22 7 13 26 22;
23 30 28 15 11 20 13 10 25 18 29 13 12 11 16 20 24 11 6 30;
28 23 12 29 16 29 26 22 28 17 20 8 30 22 8 25 12 7 13 9;
7 7 22 11 23 9 17 25 27 26 21 29 23 14 8 6 6 26 24 30];

% 1 if the pair of trips can be combined INVENTATO 20×20
O=[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0];
%O=[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0];

% distances between the depot of vehicle v (column) and
% the origin of trip n (row) 20×22 dati mancanti metto 1000 da 1 a 10
dO= [3 8 4 2 2 4 3 9 3 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
9 3 5 7 8 2 5 9 1 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
8 9 5 1 3 9 5 6 8 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 8 2 2 9 2 8 2 7 7 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 6 9 5 6 5 6 3 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 5 7 7 4 8 1 5 2 3 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 3 1 5 8 2 6 1 5 2 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 8 7 6 3 7 7 9 1 6 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 8 8 2 7 7 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7 7 8 2 8 3 9 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 2 3 3 6 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 6 8 5 4 2 8 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 2 6 5 6 7 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 1 2 7 8 8 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 9 8 9 8 9 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 6 7 9 8;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 5 7 9;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 1 8 6 4;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 4 1 3 8;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 8 9 7 3];

% distances between the destination of trip n (row)
% and the depot of vehicle v (column) 20×22 dati mancanti metto 1000 da 1 a 10

dD=[7 5 1 3 9 1 5 2 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
7 2 7 3 3 1 5 8 1 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 7 2 5 3 5 8 4 7 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
7 7 4 8 3 7 3 1 7 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 1 1 6 7 2 4 9 9 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
6 7 1 1 5 2 4 3 9 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
8 7 2 9 5 6 6 5 8 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
6 8 9 7 6 8 2 5 8 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 7 7 6 3 9 1 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 7 3 6 9 6 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 4 6 6 5 9 9 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7 7 1 7 2 1 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 9 8 6 2 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 1 5 4 1 3 3 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 6 6 6 4 6 9 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 6 7 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 9 6 8 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 4 6 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 5 4 9 1;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 8 6 9 9];

% operative cost 3×22
c=[0.7, 0.8, 0.8, 0.6, 0.9, 0.8, 0.6, 0.9, 0.6, 0.8, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100;
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.8, 0.8, 1, 0.7, 0.7, 0.8, 100, 100, 100, 100, 100;
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.9, 0.6, 0.7, 0.7];

% 1 if the trip (row) is owned by the carrier (column) INVENTATO 20×3
z=[1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
0 1 0;
0 1 0;
0 1 0;
0 1 0;
0 1 0
0 1 0;
0 1 0;
0 0 1;
0 0 1;
0 0 1;
0 0 1;
0 0 1];

% initial profit of carriers INVENTATO dati ottenuti dal prob1 per ciascun
% carrier
S0= [3380 2382 1561];
%S0 = [1000 1000 1000];

TOA=zeros(Nr,Vr);

% computation of TOA(n,v) for single trips;
for n=1:Nr
for v=1:Vr
TOA(n,v) = dO(n,v)+d(n)+dD(n,v);
end
end

E=zeros(Nr,Vr);
% computation of cost E;
for n=1:Nr
for v=1:Vr
E(n,v)= TOA(n,v)*Eempty(v) + (Efull(v)*d(n)*p(n)/CAP(v));
end
end

TOAC=zeros(Nr,Nr,Vr);
% computation of TOAC(n,k,v) for pairs of trips;
for n=1:Nr
for k=1:Nr
for v=1:Vr
TOAC(n,k,v) = dO(n,v)+d(n)+d(k)+eps(n,k)+dD(k,v);
end
end
end

EC=zeros(Nr,Nr,Vr);
% computation of cost E for combined trips;
for n=1:Nr
for k=1:Nr
for v=1:Vr
EC(n,k,v)= TOAC(n,k,v)*Eempty(v) + (Efull(v)*(d(n)*p(n)+d(k)*p(k))/CAP(v));
end
end
end

% % %%%%%%%%%%%%%%%% MAX %%%%%%%%%%%%%%%%%%%%%
%
%TO BE NOTED: I have not included the carrier in decision variables
% because here we are not includeing vehicles in the data structrue related
% to carriers but we are considering the whole list of vehicles
% than we associate vehicles with carriers by using the parameters Owner(v)

% decision variables with lower and upper bounds
x = optimvar(‘x’,Nr,Vr,’Type’,’integer’,’LowerBound’,0, ‘UpperBound’,1);
y = optimvar(‘y’,Nr,Nr,Vr,’Type’,’integer’,’LowerBound’,0, ‘UpperBound’,1);

% problema definition
prob = optimproblem(‘ObjectiveSense’,’max’);

% objective function computation

obj=0;
% %profit S1 of carriers;

for r=1:Nc
for n=1:Nr
for v=1:Vr
if (Owner(v) == r)
obj=obj+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
end
end
end
end
%
%profit S2 of carriers;

for r=1:Nc
for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (Owner(v) == r)
obj= obj+ (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
end
end
end
end
end
end

%profit S3 of carriers;

for r=1:Nc
%compensation profits and costs for single trips
for n=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
obj=obj+ delta(n)*z(n,r)*x(n,v);
else
if (Owner(v) == r)
obj=obj- delta(n)*z(n,rr)*x(n,v);
end
end
end
end
end
end
%compensation profits and costs for combined trips
for n=1:Nr
for k=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
obj = obj+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
else
if (Owner(v) == r)
obj = obj- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
end
end
end
end
end
end
end
end

%

prob.Objective = obj;

%%% CONSTRAINTS

% profits higher than initial profits

constr_profit = optimineq(Nc,1);
%
for r=1:Nc
S=0;
%S1 computation
for n=1:Nr
for v=1:Vr
if (Owner(v) == r)
S=S+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
end
end
end

%S2 computation
for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (Owner(v) == r)
S= S + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
end
end
end
end
end

% S3 computation;

%compensation profits and costs for single trips
for n=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
S=S+ delta(n)*z(n,r)*x(n,v);
else
if (Owner(v) == r)
S=S- delta(n)*z(n,rr)*x(n,v);
end
end
end
end
end
end
%compensation profits and costs for combined trips
for n=1:Nr
for k=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
S = S+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
else
if (Owner(v) == r)
S = S- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
end
end
end
end
end
end
end
constr_profit(r)=S>=S0(r);
end

prob.Constraints.constr_profit = constr_profit;

% all demand is executed (14);

constr_demand_sat=optimeq(1,1);
%
for n=1:Nr
temp=0;
for v=1:Vr
temp=temp+x(n,v);
for k=1:Nr
if (k ~= n)
temp=temp+y(n,k,v)+y(k,n,v);
end
end
end
constr_demand_sat(n)=temp==1;
end

prob.Constraints.constr_demand_sat = constr_demand_sat;

% Every vehicle must be assigned to at most one trip;

constr_assign=optimineq(Vr,1);

for v=1:Vr
temp=0;
for n=1:Nr
temp=temp+x(n,v);
for k=1:Nr
if (k ~= n)
temp=temp+y(n,k,v);
end
end

end
constr_assign(v)=temp<=1;
end

prob.Constraints.constr_assign = constr_assign;

% duration of the batteries of electric trucks;

constr_batteries_single=optimineq(Nr,Vr);

for n=1:Nr
for v=1:Vr
if (IsE(v) == 1)
constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=T(v)*mu;
else
constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=M;
end
end
end

prob.Constraints.constr_batteries_single = constr_batteries_single;

constr_batteries_combined=optimineq(Nr,Nr,Vr);

for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (IsE(v) == 1)
constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=T(v)*mu;
else
constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=M;
end
end
end
end
end

prob.Constraints.constr_batteries_combined = constr_batteries_combined;

% constraints on trips combination;

constr_combined=optimineq(Nr,Nr,Vr,1);

for n=1:Nr
for k=1:Nr
for v=1:Vr
constr_combined(n,k,v)=y(n,k,v)<=O(n,k);
end
end
end

prob.Constraints.constr_combined = constr_combined;

%
% % problem solution
show(prob)
[sol,val] = solve(prob);

% retrieve of optimal values
opt_variables_x=sol.x;
opt_variables_y=sol.y;

opt_obj=val;

% profits higher than initial profits

% %SS=zeros(Nc,1);
% %for r=1:Nc
% SS(r)=0;
% %S1 computation
% for n=1:Nr
% for v=1:Vr
% if (Owner(v) == r)
% SS(r)=SS(r)+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*sol.x(n,v);
% end
% end
% end
%
%
% %S2 computation
% for n=1:Nr
% for k=1:Nr
% if (k ~= n)
% for v=1:Vr
% if (Owner(v) == r)
% SS(r)= SS(r) + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*sol.y(n,k,v);
% end
% end
% end
% end
% end
%
%
% % S3 computation;
%
% %compensation profits and costs for single trips
% for n=1:Nr
% for rr=1:Nc
% if (rr ~= r)
% for v=1:Vr
% if (Owner(v) == rr)
% SS(r)=SS(r)+ delta(n)*z(n,r)*sol.x(n,v);
% else
% if (Owner(v) == r)
% SS(r)=SS(r)- delta(n)*z(n,rr)*sol.x(n,v);
% end
% end
% end
% end
% end
% end
% %compensation profits and costs for combined trips
% for n=1:Nr
% for k=1:Nr
% for rr=1:Nc
% if (rr ~= r)
% for v=1:Vr
% if (Owner(v) == rr)
% SS(r) = SS(r)+ sol.y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
% else
% if (Owner(v) == r)
% SS(r) = SS(r)- sol.y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
% end
% end
% end
% end
% end
% end
% end
%
% y2=0;
% for i=1:Nr
% for v=1:Vr
% y2=y2+opt_variables_y(i,2,v)+opt_variables_y(2,i,v);
% end;
% end;
%
% end
—-I need help about this code dealing with a problem of cooperative logistics. I have to maximize profits, that is, demonstrate that with the three carriers together you get a profit greater than the sum of the profits given by the three carriers taken individually. MatLab takes over 30 minutes and gives infeasible. Thank you so much!
clc
clear all

%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%

Nr=20; % number of trips
Nc=3; % number of carriers
Vr=22; % number of vehicles

alpha=1; % extra charge on revenue;
cE= 0.2;
beta=0.8;
mu=80;

% dummy quantity
M=1000000;

% data about trips 20
d= [85, 80, 123, 111, 153, 112, 150, 107, 181, 135, 109, 106, 113, 148, 164, 127, 148, 141, 108, 118];
delta=[200, 200, 200, 200, 200, 200, 200, 200, 300, 300, 300, 300, 300, 300, 300, 400, 400, 400, 400, 400];
p=[25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25];

% data about carriers, INVENTATI
I=[3, 3, 3];
cc=[0.8, 1.0, 0.9];

% data about vehicles, 22
T= [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8];
CAP= [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30];
Eempty= [0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.2, 0.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.3, 0.2, 0.2, 0.3, 0.3, 0.2];
Efull= [0.3, 0.4, 0.3, 0.3, 0.3, 0.4, 0.3, 0.4, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.4, 0.4, 0.3, 0.4, 0.4]; % it is already Efull-Eempty
IsE=[0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0]; %1 if the vehicle is electric
Owner= [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3];

% repositioing distances between trips 20×20
eps= [14 16 6 19 26 27 21 10 16 18 29 7 30 27 28 10 8 11 17 10;
27 9 20 11 9 27 22 18 27 14 29 22 5 25 9 21 17 18 11 16;
25 11 10 30 19 14 7 9 19 26 27 22 29 7 12 29 24 7 9 30;
13 16 6 19 7 7 24 5 27 28 6 18 27 9 17 30 24 14 18 10;
22 23 21 13 16 16 25 21 14 8 14 13 18 17 11 7 16 29 6 6;
5 21 9 17 26 30 27 11 24 28 16 6 25 22 27 11 18 30 13 7;
13 11 15 24 20 30 26 11 7 21 7 17 29 5 11 12 22 22 8 5;
28 7 14 20 7 25 20 12 15 5 27 10 5 27 20 26 22 25 14 28;
16 19 11 12 8 25 24 9 26 26 17 15 8 5 20 19 7 21 28 24;
29 8 13 18 17 13 15 12 21 11 10 25 17 29 22 11 28 27 14 13;
10 22 26 16 22 19 7 28 10 14 12 28 27 10 12 6 23 30 28 26;
18 11 15 25 20 7 21 17 21 9 16 5 11 5 18 16 17 26 21 19;
14 15 25 19 11 16 8 25 21 9 11 8 13 5 10 7 5 28 30 8;
26 15 11 19 21 23 9 24 9 26 16 17 9 27 25 27 20 21 15 13;
27 27 14 19 28 9 22 24 18 11 22 15 15 26 27 27 27 19 17 19;
25 11 17 24 6 30 15 23 21 30 16 5 18 25 12 20 9 11 25 27;
29 9 17 11 24 8 29 28 21 16 20 15 6 28 8 22 7 13 26 22;
23 30 28 15 11 20 13 10 25 18 29 13 12 11 16 20 24 11 6 30;
28 23 12 29 16 29 26 22 28 17 20 8 30 22 8 25 12 7 13 9;
7 7 22 11 23 9 17 25 27 26 21 29 23 14 8 6 6 26 24 30];

% 1 if the pair of trips can be combined INVENTATO 20×20
O=[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0];
%O=[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0];

% distances between the depot of vehicle v (column) and
% the origin of trip n (row) 20×22 dati mancanti metto 1000 da 1 a 10
dO= [3 8 4 2 2 4 3 9 3 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
9 3 5 7 8 2 5 9 1 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
8 9 5 1 3 9 5 6 8 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 8 2 2 9 2 8 2 7 7 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 6 9 5 6 5 6 3 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 5 7 7 4 8 1 5 2 3 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 3 1 5 8 2 6 1 5 2 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 8 7 6 3 7 7 9 1 6 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 8 8 2 7 7 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7 7 8 2 8 3 9 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 2 3 3 6 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 6 8 5 4 2 8 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 2 6 5 6 7 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 1 2 7 8 8 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 9 8 9 8 9 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 6 7 9 8;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 5 7 9;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 1 8 6 4;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 4 1 3 8;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 8 9 7 3];

% distances between the destination of trip n (row)
% and the depot of vehicle v (column) 20×22 dati mancanti metto 1000 da 1 a 10

dD=[7 5 1 3 9 1 5 2 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
7 2 7 3 3 1 5 8 1 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 7 2 5 3 5 8 4 7 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
7 7 4 8 3 7 3 1 7 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 1 1 6 7 2 4 9 9 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
6 7 1 1 5 2 4 3 9 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
8 7 2 9 5 6 6 5 8 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
6 8 9 7 6 8 2 5 8 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 7 7 6 3 9 1 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 7 3 6 9 6 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 4 6 6 5 9 9 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7 7 1 7 2 1 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 9 8 6 2 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 1 5 4 1 3 3 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 6 6 6 4 6 9 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 6 7 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 9 6 8 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 4 6 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 5 4 9 1;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 8 6 9 9];

% operative cost 3×22
c=[0.7, 0.8, 0.8, 0.6, 0.9, 0.8, 0.6, 0.9, 0.6, 0.8, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100;
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.8, 0.8, 1, 0.7, 0.7, 0.8, 100, 100, 100, 100, 100;
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.9, 0.6, 0.7, 0.7];

% 1 if the trip (row) is owned by the carrier (column) INVENTATO 20×3
z=[1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
0 1 0;
0 1 0;
0 1 0;
0 1 0;
0 1 0
0 1 0;
0 1 0;
0 0 1;
0 0 1;
0 0 1;
0 0 1;
0 0 1];

% initial profit of carriers INVENTATO dati ottenuti dal prob1 per ciascun
% carrier
S0= [3380 2382 1561];
%S0 = [1000 1000 1000];

TOA=zeros(Nr,Vr);

% computation of TOA(n,v) for single trips;
for n=1:Nr
for v=1:Vr
TOA(n,v) = dO(n,v)+d(n)+dD(n,v);
end
end

E=zeros(Nr,Vr);
% computation of cost E;
for n=1:Nr
for v=1:Vr
E(n,v)= TOA(n,v)*Eempty(v) + (Efull(v)*d(n)*p(n)/CAP(v));
end
end

TOAC=zeros(Nr,Nr,Vr);
% computation of TOAC(n,k,v) for pairs of trips;
for n=1:Nr
for k=1:Nr
for v=1:Vr
TOAC(n,k,v) = dO(n,v)+d(n)+d(k)+eps(n,k)+dD(k,v);
end
end
end

EC=zeros(Nr,Nr,Vr);
% computation of cost E for combined trips;
for n=1:Nr
for k=1:Nr
for v=1:Vr
EC(n,k,v)= TOAC(n,k,v)*Eempty(v) + (Efull(v)*(d(n)*p(n)+d(k)*p(k))/CAP(v));
end
end
end

% % %%%%%%%%%%%%%%%% MAX %%%%%%%%%%%%%%%%%%%%%
%
%TO BE NOTED: I have not included the carrier in decision variables
% because here we are not includeing vehicles in the data structrue related
% to carriers but we are considering the whole list of vehicles
% than we associate vehicles with carriers by using the parameters Owner(v)

% decision variables with lower and upper bounds
x = optimvar(‘x’,Nr,Vr,’Type’,’integer’,’LowerBound’,0, ‘UpperBound’,1);
y = optimvar(‘y’,Nr,Nr,Vr,’Type’,’integer’,’LowerBound’,0, ‘UpperBound’,1);

% problema definition
prob = optimproblem(‘ObjectiveSense’,’max’);

% objective function computation

obj=0;
% %profit S1 of carriers;

for r=1:Nc
for n=1:Nr
for v=1:Vr
if (Owner(v) == r)
obj=obj+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
end
end
end
end
%
%profit S2 of carriers;

for r=1:Nc
for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (Owner(v) == r)
obj= obj+ (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
end
end
end
end
end
end

%profit S3 of carriers;

for r=1:Nc
%compensation profits and costs for single trips
for n=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
obj=obj+ delta(n)*z(n,r)*x(n,v);
else
if (Owner(v) == r)
obj=obj- delta(n)*z(n,rr)*x(n,v);
end
end
end
end
end
end
%compensation profits and costs for combined trips
for n=1:Nr
for k=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
obj = obj+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
else
if (Owner(v) == r)
obj = obj- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
end
end
end
end
end
end
end
end

%

prob.Objective = obj;

%%% CONSTRAINTS

% profits higher than initial profits

constr_profit = optimineq(Nc,1);
%
for r=1:Nc
S=0;
%S1 computation
for n=1:Nr
for v=1:Vr
if (Owner(v) == r)
S=S+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
end
end
end

%S2 computation
for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (Owner(v) == r)
S= S + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
end
end
end
end
end

% S3 computation;

%compensation profits and costs for single trips
for n=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
S=S+ delta(n)*z(n,r)*x(n,v);
else
if (Owner(v) == r)
S=S- delta(n)*z(n,rr)*x(n,v);
end
end
end
end
end
end
%compensation profits and costs for combined trips
for n=1:Nr
for k=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
S = S+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
else
if (Owner(v) == r)
S = S- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
end
end
end
end
end
end
end
constr_profit(r)=S>=S0(r);
end

prob.Constraints.constr_profit = constr_profit;

% all demand is executed (14);

constr_demand_sat=optimeq(1,1);
%
for n=1:Nr
temp=0;
for v=1:Vr
temp=temp+x(n,v);
for k=1:Nr
if (k ~= n)
temp=temp+y(n,k,v)+y(k,n,v);
end
end
end
constr_demand_sat(n)=temp==1;
end

prob.Constraints.constr_demand_sat = constr_demand_sat;

% Every vehicle must be assigned to at most one trip;

constr_assign=optimineq(Vr,1);

for v=1:Vr
temp=0;
for n=1:Nr
temp=temp+x(n,v);
for k=1:Nr
if (k ~= n)
temp=temp+y(n,k,v);
end
end

end
constr_assign(v)=temp<=1;
end

prob.Constraints.constr_assign = constr_assign;

% duration of the batteries of electric trucks;

constr_batteries_single=optimineq(Nr,Vr);

for n=1:Nr
for v=1:Vr
if (IsE(v) == 1)
constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=T(v)*mu;
else
constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=M;
end
end
end

prob.Constraints.constr_batteries_single = constr_batteries_single;

constr_batteries_combined=optimineq(Nr,Nr,Vr);

for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (IsE(v) == 1)
constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=T(v)*mu;
else
constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=M;
end
end
end
end
end

prob.Constraints.constr_batteries_combined = constr_batteries_combined;

% constraints on trips combination;

constr_combined=optimineq(Nr,Nr,Vr,1);

for n=1:Nr
for k=1:Nr
for v=1:Vr
constr_combined(n,k,v)=y(n,k,v)<=O(n,k);
end
end
end

prob.Constraints.constr_combined = constr_combined;

%
% % problem solution
show(prob)
[sol,val] = solve(prob);

% retrieve of optimal values
opt_variables_x=sol.x;
opt_variables_y=sol.y;

opt_obj=val;

% profits higher than initial profits

% %SS=zeros(Nc,1);
% %for r=1:Nc
% SS(r)=0;
% %S1 computation
% for n=1:Nr
% for v=1:Vr
% if (Owner(v) == r)
% SS(r)=SS(r)+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*sol.x(n,v);
% end
% end
% end
%
%
% %S2 computation
% for n=1:Nr
% for k=1:Nr
% if (k ~= n)
% for v=1:Vr
% if (Owner(v) == r)
% SS(r)= SS(r) + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*sol.y(n,k,v);
% end
% end
% end
% end
% end
%
%
% % S3 computation;
%
% %compensation profits and costs for single trips
% for n=1:Nr
% for rr=1:Nc
% if (rr ~= r)
% for v=1:Vr
% if (Owner(v) == rr)
% SS(r)=SS(r)+ delta(n)*z(n,r)*sol.x(n,v);
% else
% if (Owner(v) == r)
% SS(r)=SS(r)- delta(n)*z(n,rr)*sol.x(n,v);
% end
% end
% end
% end
% end
% end
% %compensation profits and costs for combined trips
% for n=1:Nr
% for k=1:Nr
% for rr=1:Nc
% if (rr ~= r)
% for v=1:Vr
% if (Owner(v) == rr)
% SS(r) = SS(r)+ sol.y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
% else
% if (Owner(v) == r)
% SS(r) = SS(r)- sol.y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
% end
% end
% end
% end
% end
% end
% end
%
% y2=0;
% for i=1:Nr
% for v=1:Vr
% y2=y2+opt_variables_y(i,2,v)+opt_variables_y(2,i,v);
% end;
% end;
%
% end
—- I need help about this code dealing with a problem of cooperative logistics. I have to maximize profits, that is, demonstrate that with the three carriers together you get a profit greater than the sum of the profits given by the three carriers taken individually. MatLab takes over 30 minutes and gives infeasible. Thank you so much!
clc
clear all

%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%

Nr=20; % number of trips
Nc=3; % number of carriers
Vr=22; % number of vehicles

alpha=1; % extra charge on revenue;
cE= 0.2;
beta=0.8;
mu=80;

% dummy quantity
M=1000000;

% data about trips 20
d= [85, 80, 123, 111, 153, 112, 150, 107, 181, 135, 109, 106, 113, 148, 164, 127, 148, 141, 108, 118];
delta=[200, 200, 200, 200, 200, 200, 200, 200, 300, 300, 300, 300, 300, 300, 300, 400, 400, 400, 400, 400];
p=[25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25];

% data about carriers, INVENTATI
I=[3, 3, 3];
cc=[0.8, 1.0, 0.9];

% data about vehicles, 22
T= [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8];
CAP= [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30];
Eempty= [0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.2, 0.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.3, 0.2, 0.2, 0.3, 0.3, 0.2];
Efull= [0.3, 0.4, 0.3, 0.3, 0.3, 0.4, 0.3, 0.4, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.4, 0.4, 0.3, 0.4, 0.4]; % it is already Efull-Eempty
IsE=[0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0]; %1 if the vehicle is electric
Owner= [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3];

% repositioing distances between trips 20×20
eps= [14 16 6 19 26 27 21 10 16 18 29 7 30 27 28 10 8 11 17 10;
27 9 20 11 9 27 22 18 27 14 29 22 5 25 9 21 17 18 11 16;
25 11 10 30 19 14 7 9 19 26 27 22 29 7 12 29 24 7 9 30;
13 16 6 19 7 7 24 5 27 28 6 18 27 9 17 30 24 14 18 10;
22 23 21 13 16 16 25 21 14 8 14 13 18 17 11 7 16 29 6 6;
5 21 9 17 26 30 27 11 24 28 16 6 25 22 27 11 18 30 13 7;
13 11 15 24 20 30 26 11 7 21 7 17 29 5 11 12 22 22 8 5;
28 7 14 20 7 25 20 12 15 5 27 10 5 27 20 26 22 25 14 28;
16 19 11 12 8 25 24 9 26 26 17 15 8 5 20 19 7 21 28 24;
29 8 13 18 17 13 15 12 21 11 10 25 17 29 22 11 28 27 14 13;
10 22 26 16 22 19 7 28 10 14 12 28 27 10 12 6 23 30 28 26;
18 11 15 25 20 7 21 17 21 9 16 5 11 5 18 16 17 26 21 19;
14 15 25 19 11 16 8 25 21 9 11 8 13 5 10 7 5 28 30 8;
26 15 11 19 21 23 9 24 9 26 16 17 9 27 25 27 20 21 15 13;
27 27 14 19 28 9 22 24 18 11 22 15 15 26 27 27 27 19 17 19;
25 11 17 24 6 30 15 23 21 30 16 5 18 25 12 20 9 11 25 27;
29 9 17 11 24 8 29 28 21 16 20 15 6 28 8 22 7 13 26 22;
23 30 28 15 11 20 13 10 25 18 29 13 12 11 16 20 24 11 6 30;
28 23 12 29 16 29 26 22 28 17 20 8 30 22 8 25 12 7 13 9;
7 7 22 11 23 9 17 25 27 26 21 29 23 14 8 6 6 26 24 30];

% 1 if the pair of trips can be combined INVENTATO 20×20
O=[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0];
%O=[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0];

% distances between the depot of vehicle v (column) and
% the origin of trip n (row) 20×22 dati mancanti metto 1000 da 1 a 10
dO= [3 8 4 2 2 4 3 9 3 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
9 3 5 7 8 2 5 9 1 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
8 9 5 1 3 9 5 6 8 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 8 2 2 9 2 8 2 7 7 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 6 9 5 6 5 6 3 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 5 7 7 4 8 1 5 2 3 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 3 1 5 8 2 6 1 5 2 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 8 7 6 3 7 7 9 1 6 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 8 8 2 7 7 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7 7 8 2 8 3 9 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 2 3 3 6 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 6 8 5 4 2 8 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 2 6 5 6 7 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 1 2 7 8 8 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 9 8 9 8 9 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 6 7 9 8;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 5 7 9;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 1 8 6 4;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 4 1 3 8;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 8 9 7 3];

% distances between the destination of trip n (row)
% and the depot of vehicle v (column) 20×22 dati mancanti metto 1000 da 1 a 10

dD=[7 5 1 3 9 1 5 2 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
7 2 7 3 3 1 5 8 1 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1 7 2 5 3 5 8 4 7 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
7 7 4 8 3 7 3 1 7 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
5 1 1 6 7 2 4 9 9 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
6 7 1 1 5 2 4 3 9 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
8 7 2 9 5 6 6 5 8 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
6 8 9 7 6 8 2 5 8 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 7 7 6 3 9 1 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 7 3 6 9 6 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 4 6 6 5 9 9 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7 7 1 7 2 1 7 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 9 8 6 2 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 1 5 4 1 3 3 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 6 6 6 4 6 9 5 1000 1000 1000 1000 1000;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 6 7 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 9 6 8 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 4 6 2;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 5 4 9 1;
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 8 6 9 9];

% operative cost 3×22
c=[0.7, 0.8, 0.8, 0.6, 0.9, 0.8, 0.6, 0.9, 0.6, 0.8, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100;
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.8, 0.8, 1, 0.7, 0.7, 0.8, 100, 100, 100, 100, 100;
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.9, 0.6, 0.7, 0.7];

% 1 if the trip (row) is owned by the carrier (column) INVENTATO 20×3
z=[1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
1 0 0;
0 1 0;
0 1 0;
0 1 0;
0 1 0;
0 1 0
0 1 0;
0 1 0;
0 0 1;
0 0 1;
0 0 1;
0 0 1;
0 0 1];

% initial profit of carriers INVENTATO dati ottenuti dal prob1 per ciascun
% carrier
S0= [3380 2382 1561];
%S0 = [1000 1000 1000];

TOA=zeros(Nr,Vr);

% computation of TOA(n,v) for single trips;
for n=1:Nr
for v=1:Vr
TOA(n,v) = dO(n,v)+d(n)+dD(n,v);
end
end

E=zeros(Nr,Vr);
% computation of cost E;
for n=1:Nr
for v=1:Vr
E(n,v)= TOA(n,v)*Eempty(v) + (Efull(v)*d(n)*p(n)/CAP(v));
end
end

TOAC=zeros(Nr,Nr,Vr);
% computation of TOAC(n,k,v) for pairs of trips;
for n=1:Nr
for k=1:Nr
for v=1:Vr
TOAC(n,k,v) = dO(n,v)+d(n)+d(k)+eps(n,k)+dD(k,v);
end
end
end

EC=zeros(Nr,Nr,Vr);
% computation of cost E for combined trips;
for n=1:Nr
for k=1:Nr
for v=1:Vr
EC(n,k,v)= TOAC(n,k,v)*Eempty(v) + (Efull(v)*(d(n)*p(n)+d(k)*p(k))/CAP(v));
end
end
end

% % %%%%%%%%%%%%%%%% MAX %%%%%%%%%%%%%%%%%%%%%
%
%TO BE NOTED: I have not included the carrier in decision variables
% because here we are not includeing vehicles in the data structrue related
% to carriers but we are considering the whole list of vehicles
% than we associate vehicles with carriers by using the parameters Owner(v)

% decision variables with lower and upper bounds
x = optimvar(‘x’,Nr,Vr,’Type’,’integer’,’LowerBound’,0, ‘UpperBound’,1);
y = optimvar(‘y’,Nr,Nr,Vr,’Type’,’integer’,’LowerBound’,0, ‘UpperBound’,1);

% problema definition
prob = optimproblem(‘ObjectiveSense’,’max’);

% objective function computation

obj=0;
% %profit S1 of carriers;

for r=1:Nc
for n=1:Nr
for v=1:Vr
if (Owner(v) == r)
obj=obj+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
end
end
end
end
%
%profit S2 of carriers;

for r=1:Nc
for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (Owner(v) == r)
obj= obj+ (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
end
end
end
end
end
end

%profit S3 of carriers;

for r=1:Nc
%compensation profits and costs for single trips
for n=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
obj=obj+ delta(n)*z(n,r)*x(n,v);
else
if (Owner(v) == r)
obj=obj- delta(n)*z(n,rr)*x(n,v);
end
end
end
end
end
end
%compensation profits and costs for combined trips
for n=1:Nr
for k=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
obj = obj+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
else
if (Owner(v) == r)
obj = obj- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
end
end
end
end
end
end
end
end

%

prob.Objective = obj;

%%% CONSTRAINTS

% profits higher than initial profits

constr_profit = optimineq(Nc,1);
%
for r=1:Nc
S=0;
%S1 computation
for n=1:Nr
for v=1:Vr
if (Owner(v) == r)
S=S+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
end
end
end

%S2 computation
for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (Owner(v) == r)
S= S + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
end
end
end
end
end

% S3 computation;

%compensation profits and costs for single trips
for n=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
S=S+ delta(n)*z(n,r)*x(n,v);
else
if (Owner(v) == r)
S=S- delta(n)*z(n,rr)*x(n,v);
end
end
end
end
end
end
%compensation profits and costs for combined trips
for n=1:Nr
for k=1:Nr
for rr=1:Nc
if (rr ~= r)
for v=1:Vr
if (Owner(v) == rr)
S = S+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
else
if (Owner(v) == r)
S = S- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
end
end
end
end
end
end
end
constr_profit(r)=S>=S0(r);
end

prob.Constraints.constr_profit = constr_profit;

% all demand is executed (14);

constr_demand_sat=optimeq(1,1);
%
for n=1:Nr
temp=0;
for v=1:Vr
temp=temp+x(n,v);
for k=1:Nr
if (k ~= n)
temp=temp+y(n,k,v)+y(k,n,v);
end
end
end
constr_demand_sat(n)=temp==1;
end

prob.Constraints.constr_demand_sat = constr_demand_sat;

% Every vehicle must be assigned to at most one trip;

constr_assign=optimineq(Vr,1);

for v=1:Vr
temp=0;
for n=1:Nr
temp=temp+x(n,v);
for k=1:Nr
if (k ~= n)
temp=temp+y(n,k,v);
end
end

end
constr_assign(v)=temp<=1;
end

prob.Constraints.constr_assign = constr_assign;

% duration of the batteries of electric trucks;

constr_batteries_single=optimineq(Nr,Vr);

for n=1:Nr
for v=1:Vr
if (IsE(v) == 1)
constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=T(v)*mu;
else
constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=M;
end
end
end

prob.Constraints.constr_batteries_single = constr_batteries_single;

constr_batteries_combined=optimineq(Nr,Nr,Vr);

for n=1:Nr
for k=1:Nr
if (k ~= n)
for v=1:Vr
if (IsE(v) == 1)
constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=T(v)*mu;
else
constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=M;
end
end
end
end
end

prob.Constraints.constr_batteries_combined = constr_batteries_combined;

% constraints on trips combination;

constr_combined=optimineq(Nr,Nr,Vr,1);

for n=1:Nr
for k=1:Nr
for v=1:Vr
constr_combined(n,k,v)=y(n,k,v)<=O(n,k);
end
end
end

prob.Constraints.constr_combined = constr_combined;

%
% % problem solution
show(prob)
[sol,val] = solve(prob);

% retrieve of optimal values
opt_variables_x=sol.x;
opt_variables_y=sol.y;

opt_obj=val;

% profits higher than initial profits

% %SS=zeros(Nc,1);
% %for r=1:Nc
% SS(r)=0;
% %S1 computation
% for n=1:Nr
% for v=1:Vr
% if (Owner(v) == r)
% SS(r)=SS(r)+(alpha* d(n)*I(r) – TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*sol.x(n,v);
% end
% end
% end
%
%
% %S2 computation
% for n=1:Nr
% for k=1:Nr
% if (k ~= n)
% for v=1:Vr
% if (Owner(v) == r)
% SS(r)= SS(r) + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*sol.y(n,k,v);
% end
% end
% end
% end
% end
%
%
% % S3 computation;
%
% %compensation profits and costs for single trips
% for n=1:Nr
% for rr=1:Nc
% if (rr ~= r)
% for v=1:Vr
% if (Owner(v) == rr)
% SS(r)=SS(r)+ delta(n)*z(n,r)*sol.x(n,v);
% else
% if (Owner(v) == r)
% SS(r)=SS(r)- delta(n)*z(n,rr)*sol.x(n,v);
% end
% end
% end
% end
% end
% end
% %compensation profits and costs for combined trips
% for n=1:Nr
% for k=1:Nr
% for rr=1:Nc
% if (rr ~= r)
% for v=1:Vr
% if (Owner(v) == rr)
% SS(r) = SS(r)+ sol.y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
% else
% if (Owner(v) == r)
% SS(r) = SS(r)- sol.y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
% end
% end
% end
% end
% end
% end
% end
%
% y2=0;
% for i=1:Nr
% for v=1:Vr
% y2=y2+opt_variables_y(i,2,v)+opt_variables_y(2,i,v);
% end;
% end;
%
% end
—- #logistics #cooperation MATLAB Answers — New Questions

​

I am trying to use the PCBwriter function to output a Gerber file for PCB stack object with greater than 2 layers.
Matlab News

I am trying to use the PCBwriter function to output a Gerber file for PCB stack object with greater than 2 layers.

PuTI / 2025-06-20

I am trying to use the PCBwriter function to output a gerber file for pcb stack object with greater than 2 metal layers and I am getting the error:
[A, g] = gerberWrite(PcbWriteObj); Error using PCBWriter More than 2 layers of metal is not supported for generating gerber files.
Is there a way around this for PCBs with inner copper layers?
Does MATLAB support export of mulitlayer PCB designs (>2 Layers) as Gerber files?I am trying to use the PCBwriter function to output a gerber file for pcb stack object with greater than 2 metal layers and I am getting the error:
[A, g] = gerberWrite(PcbWriteObj); Error using PCBWriter More than 2 layers of metal is not supported for generating gerber files.
Is there a way around this for PCBs with inner copper layers?
Does MATLAB support export of mulitlayer PCB designs (>2 Layers) as Gerber files? I am trying to use the PCBwriter function to output a gerber file for pcb stack object with greater than 2 metal layers and I am getting the error:
[A, g] = gerberWrite(PcbWriteObj); Error using PCBWriter More than 2 layers of metal is not supported for generating gerber files.
Is there a way around this for PCBs with inner copper layers?
Does MATLAB support export of mulitlayer PCB designs (>2 Layers) as Gerber files? pcbwriter, gerber file MATLAB Answers — New Questions

​

Compare same frequency components of two signals
Matlab News

Compare same frequency components of two signals

PuTI / 2025-06-20

Hello,
I have two signals and I want to split each one of them in same number of components. After that, I want to find the components with the (almost) same frequency and compare their plots.
Any idea how could I find same frequency components of the signals?

Thanks.Hello,
I have two signals and I want to split each one of them in same number of components. After that, I want to find the components with the (almost) same frequency and compare their plots.
Any idea how could I find same frequency components of the signals?

Thanks. Hello,
I have two signals and I want to split each one of them in same number of components. After that, I want to find the components with the (almost) same frequency and compare their plots.
Any idea how could I find same frequency components of the signals?

Thanks. signal processing, frequency MATLAB Answers — New Questions

​

Interp1 is not working after applying unique because of rounding off
Matlab News

Interp1 is not working after applying unique because of rounding off

PuTI / 2025-06-20

I have a double array with many values. The first value is 3.569990000000001, the next multiple values is 3.570000000000001, then comes 3.570010000000001 and it works like this. I can see only these long values which i click on that cell, otherwise every value is being displayed as 3.57. Now I applied unique, so the multiple 3.570000000000001 is removed and there is now only one 3.570000000000001 and again and so on, but still all values are seen as 3.57. So if compared the first two values to see if it equal it is showing that they are not equal even though me as the user sees both values as 3.57. Now here comes the problem. I need to use interp1 on these values, but here interp1 is not reading them as 3.569990000000001 and 3.570000000000001. Interp1 is reading both of them as 3.57 and is giving me an error saying values need to be unique. I dont know what to do. I tried googling like converting to format long or something like that, i couldnt find anything. What can I do. Thank you.I have a double array with many values. The first value is 3.569990000000001, the next multiple values is 3.570000000000001, then comes 3.570010000000001 and it works like this. I can see only these long values which i click on that cell, otherwise every value is being displayed as 3.57. Now I applied unique, so the multiple 3.570000000000001 is removed and there is now only one 3.570000000000001 and again and so on, but still all values are seen as 3.57. So if compared the first two values to see if it equal it is showing that they are not equal even though me as the user sees both values as 3.57. Now here comes the problem. I need to use interp1 on these values, but here interp1 is not reading them as 3.569990000000001 and 3.570000000000001. Interp1 is reading both of them as 3.57 and is giving me an error saying values need to be unique. I dont know what to do. I tried googling like converting to format long or something like that, i couldnt find anything. What can I do. Thank you. I have a double array with many values. The first value is 3.569990000000001, the next multiple values is 3.570000000000001, then comes 3.570010000000001 and it works like this. I can see only these long values which i click on that cell, otherwise every value is being displayed as 3.57. Now I applied unique, so the multiple 3.570000000000001 is removed and there is now only one 3.570000000000001 and again and so on, but still all values are seen as 3.57. So if compared the first two values to see if it equal it is showing that they are not equal even though me as the user sees both values as 3.57. Now here comes the problem. I need to use interp1 on these values, but here interp1 is not reading them as 3.569990000000001 and 3.570000000000001. Interp1 is reading both of them as 3.57 and is giving me an error saying values need to be unique. I dont know what to do. I tried googling like converting to format long or something like that, i couldnt find anything. What can I do. Thank you. matlab MATLAB Answers — New Questions

​

Readtable and Readmatrix Ignore Specified Range and Produce Extra Variables
Matlab News

Readtable and Readmatrix Ignore Specified Range and Produce Extra Variables

PuTI / 2025-06-20

I have this very simple code:
stupid_matlab = readtable(‘PATH’, ‘ReadVariableNames’, false, ‘Range’, strcat(‘C2:D’, string(ri_length)));
However, it completely ignores the column portion of the range I give it. Changing the number part of the specified range changes the number of rows, but it always produces four extra columns. If I specify C2:C3, I get five columns. If I specify C2:D3, I get six columns. I have not idea where it is getting the extra columns from. I have also tried "readmatrix" and it does the exact same thing. I’ve attatched the .csv, but I’ve opened it directly and in excel and don’t see anything wrong.I have this very simple code:
stupid_matlab = readtable(‘PATH’, ‘ReadVariableNames’, false, ‘Range’, strcat(‘C2:D’, string(ri_length)));
However, it completely ignores the column portion of the range I give it. Changing the number part of the specified range changes the number of rows, but it always produces four extra columns. If I specify C2:C3, I get five columns. If I specify C2:D3, I get six columns. I have not idea where it is getting the extra columns from. I have also tried "readmatrix" and it does the exact same thing. I’ve attatched the .csv, but I’ve opened it directly and in excel and don’t see anything wrong. I have this very simple code:
stupid_matlab = readtable(‘PATH’, ‘ReadVariableNames’, false, ‘Range’, strcat(‘C2:D’, string(ri_length)));
However, it completely ignores the column portion of the range I give it. Changing the number part of the specified range changes the number of rows, but it always produces four extra columns. If I specify C2:C3, I get five columns. If I specify C2:D3, I get six columns. I have not idea where it is getting the extra columns from. I have also tried "readmatrix" and it does the exact same thing. I’ve attatched the .csv, but I’ve opened it directly and in excel and don’t see anything wrong. readtable, readmatrix, range MATLAB Answers — New Questions

​

Clean and automate indents in scripts an live scripts
Matlab News

Clean and automate indents in scripts an live scripts

PuTI / 2025-06-19

Hi everyone,
Please, do you have any tips to automate indentation in scripts?
Example by the picture below:Hi everyone,
Please, do you have any tips to automate indentation in scripts?
Example by the picture below: Hi everyone,
Please, do you have any tips to automate indentation in scripts?
Example by the picture below: automatic, indentation, alignment, text MATLAB Answers — New Questions

​

Varying device latency for ASIO Device
Matlab News

Varying device latency for ASIO Device

PuTI / 2025-06-19

Hi,
my goal is to measure impulse responses using Matlab and an ASIO device. I encountered the problem, that for different runs of my code, the peak of my impulse response is at a different location on the time axis. I was able to reproduce this behavior outside of my code in the "Impulse Respones Measurer". Looping back the output channel to an input channel and capturing a fiew imulse reponses while leaving everything else untouched gives me shifted impulse responses.

I also followed this example: Measure Audio Latency, which gave me the same results.
When measuring the same impulse response (i.e. same device, same configuration of channels etc.) with another software (REW), the delay remained constant between different runs.
I was also able to reproduce this behavior vor R2020a and R2024b
Any suggestions how to keep a constant latency between several runs of audioPlayerRecorder are appreciated.
Kind Regards,
SebastianHi,
my goal is to measure impulse responses using Matlab and an ASIO device. I encountered the problem, that for different runs of my code, the peak of my impulse response is at a different location on the time axis. I was able to reproduce this behavior outside of my code in the "Impulse Respones Measurer". Looping back the output channel to an input channel and capturing a fiew imulse reponses while leaving everything else untouched gives me shifted impulse responses.

I also followed this example: Measure Audio Latency, which gave me the same results.
When measuring the same impulse response (i.e. same device, same configuration of channels etc.) with another software (REW), the delay remained constant between different runs.
I was also able to reproduce this behavior vor R2020a and R2024b
Any suggestions how to keep a constant latency between several runs of audioPlayerRecorder are appreciated.
Kind Regards,
Sebastian Hi,
my goal is to measure impulse responses using Matlab and an ASIO device. I encountered the problem, that for different runs of my code, the peak of my impulse response is at a different location on the time axis. I was able to reproduce this behavior outside of my code in the "Impulse Respones Measurer". Looping back the output channel to an input channel and capturing a fiew imulse reponses while leaving everything else untouched gives me shifted impulse responses.

I also followed this example: Measure Audio Latency, which gave me the same results.
When measuring the same impulse response (i.e. same device, same configuration of channels etc.) with another software (REW), the delay remained constant between different runs.
I was also able to reproduce this behavior vor R2020a and R2024b
Any suggestions how to keep a constant latency between several runs of audioPlayerRecorder are appreciated.
Kind Regards,
Sebastian asio, matlab, audio, latency, delay, impulse response, loopback, signal processing MATLAB Answers — New Questions

​

COnnecting points on opposite end of mask
Matlab News

COnnecting points on opposite end of mask

PuTI / 2025-06-19

Given an mask image i have certain points as seen in the image, I want to connect the point to the point at the opposite side of it given the list of points/array of point.Well the output would look something like.Given an mask image i have certain points as seen in the image, I want to connect the point to the point at the opposite side of it given the list of points/array of point.Well the output would look something like. Given an mask image i have certain points as seen in the image, I want to connect the point to the point at the opposite side of it given the list of points/array of point.Well the output would look something like. image analysis, image processing, digital image processing, computer vision, contour MATLAB Answers — New Questions

​

Mass importing .txt files with sequential names to plot multiple graphs on same figure for easy comparison
Matlab News

Mass importing .txt files with sequential names to plot multiple graphs on same figure for easy comparison

PuTI / 2025-06-19

I have 30 .txt files numbered like so:
FreeField1.txt
FreeField2.txt
…
FreeField30.txt

I need to find a way to import the data from all the files at once so I can then plot a graph with all of them overlaid to allow for quick comparison.

The current method I have found which inputs all the data is like so
numfiles = 30;
S = dir(‘*.txt’)
N = length(S)
for i = 1:numel(S)
S(i).data = readtable(S(i).name)
end

But this gives a structure which I do not how to progress from. I am looking for any possible methods which could be used for this mass importation to then build off for the plotting of graphs.

Any help would be appreciated as I am very new to MATLab and am wishing to understand it more. I have attached a sample text file for convenience.I have 30 .txt files numbered like so:
FreeField1.txt
FreeField2.txt
…
FreeField30.txt

I need to find a way to import the data from all the files at once so I can then plot a graph with all of them overlaid to allow for quick comparison.

The current method I have found which inputs all the data is like so
numfiles = 30;
S = dir(‘*.txt’)
N = length(S)
for i = 1:numel(S)
S(i).data = readtable(S(i).name)
end

But this gives a structure which I do not how to progress from. I am looking for any possible methods which could be used for this mass importation to then build off for the plotting of graphs.

Any help would be appreciated as I am very new to MATLab and am wishing to understand it more. I have attached a sample text file for convenience. I have 30 .txt files numbered like so:
FreeField1.txt
FreeField2.txt
…
FreeField30.txt

I need to find a way to import the data from all the files at once so I can then plot a graph with all of them overlaid to allow for quick comparison.

The current method I have found which inputs all the data is like so
numfiles = 30;
S = dir(‘*.txt’)
N = length(S)
for i = 1:numel(S)
S(i).data = readtable(S(i).name)
end

But this gives a structure which I do not how to progress from. I am looking for any possible methods which could be used for this mass importation to then build off for the plotting of graphs.

Any help would be appreciated as I am very new to MATLab and am wishing to understand it more. I have attached a sample text file for convenience. data import, graph, text file MATLAB Answers — New Questions

​

3D surf plot for more than two quantities
Matlab News

3D surf plot for more than two quantities

PuTI / 2025-06-19

I want to plot the value in the same 3D graph .

syms x t r b %alpha
% Parameter values
a=(pi)/3;
g=9.8;
U=2.5;
O=7.29*10^(-5);
f=2*O*sin(a);
H=-(f/g)*U;
alpha=0.75; % fractional order
%%%%%%%%%initalization of variable
u_l=sym(zeros(1));
v_l=zeros(1,’sym’);
h_l=zeros(1,’sym’);
A_l=zeros(1,2,’sym’);
B_l=zeros(1,2,’sym’);
C_l=zeros(1,2,’sym’);
D_l=zeros(1,2,’sym’);
series1_l(x,t)=sym(zeros(1,1));
series2_l(x,t)=sym(zeros(1,1));
series3_l(x,t)=sym(zeros(1,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_up=sym(zeros(1));
v_up=zeros(1,’sym’);
h_up=zeros(1,’sym’);
A_up=zeros(1,2,’sym’);
B_up=zeros(1,2,’sym’);
C_up=zeros(1,2,’sym’);
D_up=zeros(1,2,’sym’);
series1_up(x,t)=sym(zeros(1,1));
series2_up(x,t)=sym(zeros(1,1));
series3_up(x,t)=sym(zeros(1,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_cr=sym(zeros(1));
v_cr=zeros(1,’sym’);
h_cr=zeros(1,’sym’);
A_cr=zeros(1,2,’sym’);
B_cr=zeros(1,2,’sym’);
C_cr=zeros(1,2,’sym’);
D_cr=zeros(1,2,’sym’);
series1_cr(x,t)=sym(zeros(1,1));
series2_cr(x,t)=sym(zeros(1,1));
series3_cr(x,t)=sym(zeros(1,1));
%%%%%%%% Initial condition fuzzy condition
R=0.5;
b_l=0 % lower bound
b_cr=0.5 % middle value
b_u=1 % upper value
%%%%%%%%%%%%%%%%%% lOWER VALUE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_l(1)=(2*b_l*(1-R)+R)*exp(x)*(sech(x))^2;
v_l(1)=(2*b_l*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_l(1)=(2*b_l*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%UPPER VALUE
u_up(1)=(2*b_u*(1-R)+R)*exp(x)*(sech(x))^2;
v_up(1)=(2*b_u*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_up(1)=(2*b_u*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Middle valur %%%%%%%%%%%%%%%%%%%%
u_cr(1)=(2*b_cr*(1-R)+R)*exp(x)*(sech(x))^2;
v_cr(1)=(2*b_cr*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_cr(1)=(2*b_cr*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%lOWER VALUR %%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:2
A_l=0;
B_l=0;
C_l=0;
D_l=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_l=A_l+u_l(r)*diff(u_l(k-r+1),x,1);
B_l=B_l+u_l(r)*diff(v_l(k-r+1),x,1);
C_l=C_l+u_l(r)*diff(h_l(k-r+1),x,1);
D_l=D_l+h_l(r)*diff(u_l(k-r+1),x,1);
end
u_l(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_l+f*v_l(k)-g*diff(h_l(k),x,1)));
v_l(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_l-f*u_l(k)-g*H*E);
h_l(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_l-H*v_l(k)-D_l);
end
% t1=simplify(u(2))
% var2 = vpa(t1)
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_l(x,t)=simplify(series1_l(x,t)+u_l(k)*(power(t,k-1)));
series2_l(x,t)=simplify(series2_l(x,t)+v_l(k)*(power(t,k-1)));
series3_l(x,t)=simplify(series3_l(x,t)+h_l(k)*(power(t,k-1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% UPPER VALUE %%%%%%%%%%%%%%%%
for k=1:2
A_up=0;
B_up=0;
C_up=0;
D_up=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_up=A_up+u_up(r)*diff(u_up(k-r+1),x,1);
B_up=B_up+u_up(r)*diff(v_up(k-r+1),x,1);
C_up=C_up+u_up(r)*diff(h_up(k-r+1),x,1);
D_up=D_up+h_up(r)*diff(u_up(k-r+1),x,1);
end
u_up(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_up+f*v_up(k)-g*diff(h_up(k),x,1)));
v_up(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_up-f*u_up(k)-g*H*E);
h_up(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_up-H*v_up(k)-D_up);
end
% t1=simplify(u(2))
% var2 = vpa(t1)_u
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_up(x,t)=simplify(series1_up(x,t)+u_up(k)*(power(t,k-1)));
series2_up(x,t)=simplify(series2_up(x,t)+v_up(k)*(power(t,k-1)));
series3_up(x,t)=simplify(series3_up(x,t)+h_up(k)*(power(t,k-1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%middle value %%%%%%%%%%%%%%
for k=1:2
A_cr=0;
B_cr=0;
C_cr=0;
D_cr=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_cr=A_cr+u_cr(r)*diff(u_up(k-r+1),x,1);
B_cr=B_cr+u_cr(r)*diff(v_up(k-r+1),x,1);
C_cr=C_cr+u_cr(r)*diff(h_up(k-r+1),x,1);
D_cr=D_up+h_cr(r)*diff(u_up(k-r+1),x,1);
end
u_cr(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_cr+f*v_up(k)-g*diff(h_cr(k),x,1)));
v_cr(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_cr-f*u_cr(k)-g*H*E);
h_cr(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_cr-H*v_cr(k)-D_cr);
end
% t1=simplify(u(2))
% var2 = vpa(t1)_u
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_cr(x,t)=simplify(series1_cr(x,t)+u_cr(k)*(power(t,k-1)));
series2_cr(x,t)=simplify(series2_cr(x,t)+v_cr(k)*(power(t,k-1)));
series3_cr(x,t)=simplify(series3_cr(x,t)+h_cr(k)*(power(t,k-1)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C1_l=zeros(1)
C2_l=zeros(1)
C3_l=zeros(1)
C1_up=zeros(1)
C2_up=zeros(1)
C3_up=zeros(1)
C1_cr=zeros(1)
C2_cr=zeros(1)
C3_cr=zeros(1)
row=0;
x=0:0.04:2
t=0:0.002:0.1
for i=1:length(x)
row=row+1;
col=0;
for j=1:length(t)
col=col+1;
C1_l(row,col)=series1_l(x(i),t(j));
C2_l(row,col)=series2_l(x(i),t(j));
C3_l(row,col)=series3_l(x(i),t(j));
%———————————-
C1_up(row,col)=series1_up(x(i),t(j));
C2_up(row,col)=series2_up(x(i),t(j));
C3_up(row,col)=series3_up(x(i),t(j));
%————————————
C1_cr(row,col)=series1_cr(x(i),t(j));
C2_cr(row,col)=series2_cr(x(i),t(j));
C3_cr(row,col)=series3_cr(x(i),t(j));

end
end
%——————————————————–
surf(x,t,C1_l,C1_up,C1_cr)
surf(x,t,C2_l,C2_up,C2_cr)
surf(x,t,C3_l,C3_up,C3_cr)
I want to plot all the graph in the same 3D plot. I dont want different graph using subplot. and I am using MATLAB version 2023aI want to plot the value in the same 3D graph .

syms x t r b %alpha
% Parameter values
a=(pi)/3;
g=9.8;
U=2.5;
O=7.29*10^(-5);
f=2*O*sin(a);
H=-(f/g)*U;
alpha=0.75; % fractional order
%%%%%%%%%initalization of variable
u_l=sym(zeros(1));
v_l=zeros(1,’sym’);
h_l=zeros(1,’sym’);
A_l=zeros(1,2,’sym’);
B_l=zeros(1,2,’sym’);
C_l=zeros(1,2,’sym’);
D_l=zeros(1,2,’sym’);
series1_l(x,t)=sym(zeros(1,1));
series2_l(x,t)=sym(zeros(1,1));
series3_l(x,t)=sym(zeros(1,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_up=sym(zeros(1));
v_up=zeros(1,’sym’);
h_up=zeros(1,’sym’);
A_up=zeros(1,2,’sym’);
B_up=zeros(1,2,’sym’);
C_up=zeros(1,2,’sym’);
D_up=zeros(1,2,’sym’);
series1_up(x,t)=sym(zeros(1,1));
series2_up(x,t)=sym(zeros(1,1));
series3_up(x,t)=sym(zeros(1,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_cr=sym(zeros(1));
v_cr=zeros(1,’sym’);
h_cr=zeros(1,’sym’);
A_cr=zeros(1,2,’sym’);
B_cr=zeros(1,2,’sym’);
C_cr=zeros(1,2,’sym’);
D_cr=zeros(1,2,’sym’);
series1_cr(x,t)=sym(zeros(1,1));
series2_cr(x,t)=sym(zeros(1,1));
series3_cr(x,t)=sym(zeros(1,1));
%%%%%%%% Initial condition fuzzy condition
R=0.5;
b_l=0 % lower bound
b_cr=0.5 % middle value
b_u=1 % upper value
%%%%%%%%%%%%%%%%%% lOWER VALUE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_l(1)=(2*b_l*(1-R)+R)*exp(x)*(sech(x))^2;
v_l(1)=(2*b_l*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_l(1)=(2*b_l*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%UPPER VALUE
u_up(1)=(2*b_u*(1-R)+R)*exp(x)*(sech(x))^2;
v_up(1)=(2*b_u*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_up(1)=(2*b_u*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Middle valur %%%%%%%%%%%%%%%%%%%%
u_cr(1)=(2*b_cr*(1-R)+R)*exp(x)*(sech(x))^2;
v_cr(1)=(2*b_cr*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_cr(1)=(2*b_cr*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%lOWER VALUR %%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:2
A_l=0;
B_l=0;
C_l=0;
D_l=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_l=A_l+u_l(r)*diff(u_l(k-r+1),x,1);
B_l=B_l+u_l(r)*diff(v_l(k-r+1),x,1);
C_l=C_l+u_l(r)*diff(h_l(k-r+1),x,1);
D_l=D_l+h_l(r)*diff(u_l(k-r+1),x,1);
end
u_l(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_l+f*v_l(k)-g*diff(h_l(k),x,1)));
v_l(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_l-f*u_l(k)-g*H*E);
h_l(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_l-H*v_l(k)-D_l);
end
% t1=simplify(u(2))
% var2 = vpa(t1)
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_l(x,t)=simplify(series1_l(x,t)+u_l(k)*(power(t,k-1)));
series2_l(x,t)=simplify(series2_l(x,t)+v_l(k)*(power(t,k-1)));
series3_l(x,t)=simplify(series3_l(x,t)+h_l(k)*(power(t,k-1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% UPPER VALUE %%%%%%%%%%%%%%%%
for k=1:2
A_up=0;
B_up=0;
C_up=0;
D_up=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_up=A_up+u_up(r)*diff(u_up(k-r+1),x,1);
B_up=B_up+u_up(r)*diff(v_up(k-r+1),x,1);
C_up=C_up+u_up(r)*diff(h_up(k-r+1),x,1);
D_up=D_up+h_up(r)*diff(u_up(k-r+1),x,1);
end
u_up(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_up+f*v_up(k)-g*diff(h_up(k),x,1)));
v_up(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_up-f*u_up(k)-g*H*E);
h_up(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_up-H*v_up(k)-D_up);
end
% t1=simplify(u(2))
% var2 = vpa(t1)_u
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_up(x,t)=simplify(series1_up(x,t)+u_up(k)*(power(t,k-1)));
series2_up(x,t)=simplify(series2_up(x,t)+v_up(k)*(power(t,k-1)));
series3_up(x,t)=simplify(series3_up(x,t)+h_up(k)*(power(t,k-1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%middle value %%%%%%%%%%%%%%
for k=1:2
A_cr=0;
B_cr=0;
C_cr=0;
D_cr=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_cr=A_cr+u_cr(r)*diff(u_up(k-r+1),x,1);
B_cr=B_cr+u_cr(r)*diff(v_up(k-r+1),x,1);
C_cr=C_cr+u_cr(r)*diff(h_up(k-r+1),x,1);
D_cr=D_up+h_cr(r)*diff(u_up(k-r+1),x,1);
end
u_cr(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_cr+f*v_up(k)-g*diff(h_cr(k),x,1)));
v_cr(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_cr-f*u_cr(k)-g*H*E);
h_cr(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_cr-H*v_cr(k)-D_cr);
end
% t1=simplify(u(2))
% var2 = vpa(t1)_u
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_cr(x,t)=simplify(series1_cr(x,t)+u_cr(k)*(power(t,k-1)));
series2_cr(x,t)=simplify(series2_cr(x,t)+v_cr(k)*(power(t,k-1)));
series3_cr(x,t)=simplify(series3_cr(x,t)+h_cr(k)*(power(t,k-1)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C1_l=zeros(1)
C2_l=zeros(1)
C3_l=zeros(1)
C1_up=zeros(1)
C2_up=zeros(1)
C3_up=zeros(1)
C1_cr=zeros(1)
C2_cr=zeros(1)
C3_cr=zeros(1)
row=0;
x=0:0.04:2
t=0:0.002:0.1
for i=1:length(x)
row=row+1;
col=0;
for j=1:length(t)
col=col+1;
C1_l(row,col)=series1_l(x(i),t(j));
C2_l(row,col)=series2_l(x(i),t(j));
C3_l(row,col)=series3_l(x(i),t(j));
%———————————-
C1_up(row,col)=series1_up(x(i),t(j));
C2_up(row,col)=series2_up(x(i),t(j));
C3_up(row,col)=series3_up(x(i),t(j));
%————————————
C1_cr(row,col)=series1_cr(x(i),t(j));
C2_cr(row,col)=series2_cr(x(i),t(j));
C3_cr(row,col)=series3_cr(x(i),t(j));

end
end
%——————————————————–
surf(x,t,C1_l,C1_up,C1_cr)
surf(x,t,C2_l,C2_up,C2_cr)
surf(x,t,C3_l,C3_up,C3_cr)
I want to plot all the graph in the same 3D plot. I dont want different graph using subplot. and I am using MATLAB version 2023a I want to plot the value in the same 3D graph .

syms x t r b %alpha
% Parameter values
a=(pi)/3;
g=9.8;
U=2.5;
O=7.29*10^(-5);
f=2*O*sin(a);
H=-(f/g)*U;
alpha=0.75; % fractional order
%%%%%%%%%initalization of variable
u_l=sym(zeros(1));
v_l=zeros(1,’sym’);
h_l=zeros(1,’sym’);
A_l=zeros(1,2,’sym’);
B_l=zeros(1,2,’sym’);
C_l=zeros(1,2,’sym’);
D_l=zeros(1,2,’sym’);
series1_l(x,t)=sym(zeros(1,1));
series2_l(x,t)=sym(zeros(1,1));
series3_l(x,t)=sym(zeros(1,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_up=sym(zeros(1));
v_up=zeros(1,’sym’);
h_up=zeros(1,’sym’);
A_up=zeros(1,2,’sym’);
B_up=zeros(1,2,’sym’);
C_up=zeros(1,2,’sym’);
D_up=zeros(1,2,’sym’);
series1_up(x,t)=sym(zeros(1,1));
series2_up(x,t)=sym(zeros(1,1));
series3_up(x,t)=sym(zeros(1,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_cr=sym(zeros(1));
v_cr=zeros(1,’sym’);
h_cr=zeros(1,’sym’);
A_cr=zeros(1,2,’sym’);
B_cr=zeros(1,2,’sym’);
C_cr=zeros(1,2,’sym’);
D_cr=zeros(1,2,’sym’);
series1_cr(x,t)=sym(zeros(1,1));
series2_cr(x,t)=sym(zeros(1,1));
series3_cr(x,t)=sym(zeros(1,1));
%%%%%%%% Initial condition fuzzy condition
R=0.5;
b_l=0 % lower bound
b_cr=0.5 % middle value
b_u=1 % upper value
%%%%%%%%%%%%%%%%%% lOWER VALUE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_l(1)=(2*b_l*(1-R)+R)*exp(x)*(sech(x))^2;
v_l(1)=(2*b_l*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_l(1)=(2*b_l*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%UPPER VALUE
u_up(1)=(2*b_u*(1-R)+R)*exp(x)*(sech(x))^2;
v_up(1)=(2*b_u*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_up(1)=(2*b_u*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Middle valur %%%%%%%%%%%%%%%%%%%%
u_cr(1)=(2*b_cr*(1-R)+R)*exp(x)*(sech(x))^2;
v_cr(1)=(2*b_cr*(1-R)+R+1)*2*x*(sech(2*x))^2;
h_cr(1)=(2*b_cr*(1-R)+R)*x^2*(sech(2*x))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%lOWER VALUR %%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:2
A_l=0;
B_l=0;
C_l=0;
D_l=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_l=A_l+u_l(r)*diff(u_l(k-r+1),x,1);
B_l=B_l+u_l(r)*diff(v_l(k-r+1),x,1);
C_l=C_l+u_l(r)*diff(h_l(k-r+1),x,1);
D_l=D_l+h_l(r)*diff(u_l(k-r+1),x,1);
end
u_l(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_l+f*v_l(k)-g*diff(h_l(k),x,1)));
v_l(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_l-f*u_l(k)-g*H*E);
h_l(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_l-H*v_l(k)-D_l);
end
% t1=simplify(u(2))
% var2 = vpa(t1)
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_l(x,t)=simplify(series1_l(x,t)+u_l(k)*(power(t,k-1)));
series2_l(x,t)=simplify(series2_l(x,t)+v_l(k)*(power(t,k-1)));
series3_l(x,t)=simplify(series3_l(x,t)+h_l(k)*(power(t,k-1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% UPPER VALUE %%%%%%%%%%%%%%%%
for k=1:2
A_up=0;
B_up=0;
C_up=0;
D_up=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_up=A_up+u_up(r)*diff(u_up(k-r+1),x,1);
B_up=B_up+u_up(r)*diff(v_up(k-r+1),x,1);
C_up=C_up+u_up(r)*diff(h_up(k-r+1),x,1);
D_up=D_up+h_up(r)*diff(u_up(k-r+1),x,1);
end
u_up(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_up+f*v_up(k)-g*diff(h_up(k),x,1)));
v_up(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_up-f*u_up(k)-g*H*E);
h_up(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_up-H*v_up(k)-D_up);
end
% t1=simplify(u(2))
% var2 = vpa(t1)_u
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_up(x,t)=simplify(series1_up(x,t)+u_up(k)*(power(t,k-1)));
series2_up(x,t)=simplify(series2_up(x,t)+v_up(k)*(power(t,k-1)));
series3_up(x,t)=simplify(series3_up(x,t)+h_up(k)*(power(t,k-1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%middle value %%%%%%%%%%%%%%
for k=1:2
A_cr=0;
B_cr=0;
C_cr=0;
D_cr=0;
if k==1
E=1;
else
E=0;
end
for r=1:k
A_cr=A_cr+u_cr(r)*diff(u_up(k-r+1),x,1);
B_cr=B_cr+u_cr(r)*diff(v_up(k-r+1),x,1);
C_cr=C_cr+u_cr(r)*diff(h_up(k-r+1),x,1);
D_cr=D_up+h_cr(r)*diff(u_up(k-r+1),x,1);
end
u_cr(k+1)=simplify(gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-A_cr+f*v_up(k)-g*diff(h_cr(k),x,1)));
v_cr(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-B_cr-f*u_cr(k)-g*H*E);
h_cr(k+1)=gamma(((k-1)*alpha)+1)/gamma((alpha*(k+1-1))+1)*(-C_cr-H*v_cr(k)-D_cr);
end
% t1=simplify(u(2))
% var2 = vpa(t1)_u
% t2=simplify(v(2))
% var = vpa(t2)
% t3=simplify(h(2))
for k=1:3
series1_cr(x,t)=simplify(series1_cr(x,t)+u_cr(k)*(power(t,k-1)));
series2_cr(x,t)=simplify(series2_cr(x,t)+v_cr(k)*(power(t,k-1)));
series3_cr(x,t)=simplify(series3_cr(x,t)+h_cr(k)*(power(t,k-1)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C1_l=zeros(1)
C2_l=zeros(1)
C3_l=zeros(1)
C1_up=zeros(1)
C2_up=zeros(1)
C3_up=zeros(1)
C1_cr=zeros(1)
C2_cr=zeros(1)
C3_cr=zeros(1)
row=0;
x=0:0.04:2
t=0:0.002:0.1
for i=1:length(x)
row=row+1;
col=0;
for j=1:length(t)
col=col+1;
C1_l(row,col)=series1_l(x(i),t(j));
C2_l(row,col)=series2_l(x(i),t(j));
C3_l(row,col)=series3_l(x(i),t(j));
%———————————-
C1_up(row,col)=series1_up(x(i),t(j));
C2_up(row,col)=series2_up(x(i),t(j));
C3_up(row,col)=series3_up(x(i),t(j));
%————————————
C1_cr(row,col)=series1_cr(x(i),t(j));
C2_cr(row,col)=series2_cr(x(i),t(j));
C3_cr(row,col)=series3_cr(x(i),t(j));

end
end
%——————————————————–
surf(x,t,C1_l,C1_up,C1_cr)
surf(x,t,C2_l,C2_up,C2_cr)
surf(x,t,C3_l,C3_up,C3_cr)
I want to plot all the graph in the same 3D plot. I dont want different graph using subplot. and I am using MATLAB version 2023a 3 d plot MATLAB Answers — New Questions

​

Previous 1 … 10 11 12 13 14 … 59 Next

Search

Categories

  • Matlab
  • Microsoft
  • News
  • Other
Application Package Repository Telkom University

Tags

matlab microsoft opensources
Application Package Download License

Application Package Download License

Adobe
Google for Education
IBM
Matlab
Microsoft
Wordpress
Visual Paradigm
Opensource

Sign Up For Newsletters

Be the First to Know. Sign up for newsletter today

Application Package Repository Telkom University

Portal Application Package Repository Telkom University, for internal use only, empower civitas academica in study and research.

Information

  • Telkom University
  • About Us
  • Contact
  • Forum Discussion
  • FAQ
  • Helpdesk Ticket

Contact Us

  • Ask: Any question please read FAQ
  • Mail: helpdesk@telkomuniversity.ac.id
  • Call: +62 823-1994-9941
  • WA: +62 823-1994-9943
  • Site: Gedung Panambulai. Jl. Telekomunikasi

Copyright © Telkom University. All Rights Reserved. ch

  • FAQ
  • Privacy Policy
  • Term