Category: Matlab
Category Archives: Matlab
Can I use ”PLACE” matlab function to find feedback gain for discrete state space controller?
I am using PLACE function in MATLAB to find feedback gain for continuous system and I am getting desired output. But when I am using PLACE function for discrete system then I am not getting correct feedback gain. Can I use PLACE function for discrete?
If not then suggest me some another way to find feedback gain matrix for discrete system.I am using PLACE function in MATLAB to find feedback gain for continuous system and I am getting desired output. But when I am using PLACE function for discrete system then I am not getting correct feedback gain. Can I use PLACE function for discrete?
If not then suggest me some another way to find feedback gain matrix for discrete system. I am using PLACE function in MATLAB to find feedback gain for continuous system and I am getting desired output. But when I am using PLACE function for discrete system then I am not getting correct feedback gain. Can I use PLACE function for discrete?
If not then suggest me some another way to find feedback gain matrix for discrete system. control, state space MATLAB Answers — New Questions
I don’t know why my Monte Carlo Localization simulation doesn’t work, I’m using the Monte Carlo Localization toolbox
I think it’s a relatively simple tool and I can’t find my error. Here’s the code:
timeStep = 0.1; % Time step (s)
lin_vel = 1; % Linear velocity (m/s)
resolution = 1; % Map resolution (m)
initial_position = [2 10 90]; % [m m degrees]
robot_pose = initial_position;
lidar_res = 0.1; % Lidar resolution (m)
lidar_range = 5.5;
num_beams = 173;
angles = linspace(-pi/2,pi/2,num_beams);
map_matrix = [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 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 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];
% Binary occupancy map matrix
map = binaryOccupancyMap(map_matrix,resolution);
% In order to follow the red path, the robot has to perform the following
% movements:
% – Up for 5 seconds
% – Right for 8 seconds
% – Down for 6 seconds
% – Right for 6 seconds
% – Down for 7 seconds
% – Left for 7 seconds
% – Up for 5 seconds
% – Left for 7 seconds
% – Up for 2 seconds
% Total simulation time: 53 seconds, we’ll simulate for 1 whole minute
totalTime = 60; % Total simulation time (s)
numParticles = 5000; % Initial number of particles
mcl = monteCarloLocalization; % Initialize MCL object
mcl.UseLidarScan = true;
mcl.GlobalLocalization = true;
mcl.ParticleLimits = [500 numParticles];
mcl.UpdateThresholds = [0.2 0.2 0.2];
mcl.ResamplingInterval = 1;
motionModel = odometryMotionModel; % Initialise motion model
motionModel.Noise = [0.05 0.05 0.05 0.05]; % Add motion error
mcl.MotionModel = motionModel;
sensorModel = likelihoodFieldSensorModel;
sensorModel.Map = map;
sensorModel.SensorLimits = [0 lidar_range];
sensorModel.NumBeams = num_beams;
mcl.SensorModel = sensorModel;
% Simulation loop
for t=0:timeStep:totalTime
if t < 5
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
elseif (t>=5)&&(t<13)
robot_pose(1) = robot_pose(1) + lin_vel * timeStep;
robot_pose(3) = 0;
elseif (t>=13)&&(t<19)
robot_pose(2) = robot_pose(2) – lin_vel * timeStep;
robot_pose(3) = 270;
elseif (t>=19)&&(t<25)
robot_pose(1) = robot_pose(1) + lin_vel * timeStep;
robot_pose(3) = 0;
elseif (t>=25)&&(t<32)
robot_pose(2) = robot_pose(2) – lin_vel * timeStep;
robot_pose(3) = 270;
elseif (t>=32)&&(t<39)
robot_pose(1) = robot_pose(1) – lin_vel * timeStep;
robot_pose(3) = 180;
elseif (t>=39)&&(t<44)
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
elseif (t>=44)&&(t<51)
robot_pose(1) = robot_pose(1) – lin_vel * timeStep;
robot_pose(3) = 180;
elseif (t>=51)&&(t<53)
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
else
% Robot stopped
end
% rayIntersection func. works with radians
robot_pose_radians = robot_pose;
robot_pose_radians(1,3) = robot_pose_radians(1,3)*(pi/180);
intersectionPts = rayIntersection(map,robot_pose_radians,angles,lidar_range);
ranges = sqrt((intersectionPts(:,1)-robot_pose(1)).^2+(intersectionPts(:,2)-robot_pose(2)).^2);
scan = lidarScan(ranges,angles);
[isUpdated,estimatedPose,covariance] = mcl(robot_pose,scan);
particles = getParticles(mcl);
% Update particles plot every 2 "seconds"
if mod(t,2) == 0
particles_plot = particles;
end
figure(1)
show(map)
hold on
scatter(particles_plot(:,1),particles_plot(:,2),’b.’)
hold on
scatter(estimatedPose(:,1),estimatedPose(:,2),’ro’)
hold on
scatter(intersectionPts(:,1),intersectionPts(:,2),’r.’)
hold on
scatter(robot_pose(:,1),robot_pose(:,2),’bo’)
hold off
xlim([0 19])
ylim([0 17])
endI think it’s a relatively simple tool and I can’t find my error. Here’s the code:
timeStep = 0.1; % Time step (s)
lin_vel = 1; % Linear velocity (m/s)
resolution = 1; % Map resolution (m)
initial_position = [2 10 90]; % [m m degrees]
robot_pose = initial_position;
lidar_res = 0.1; % Lidar resolution (m)
lidar_range = 5.5;
num_beams = 173;
angles = linspace(-pi/2,pi/2,num_beams);
map_matrix = [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 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 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];
% Binary occupancy map matrix
map = binaryOccupancyMap(map_matrix,resolution);
% In order to follow the red path, the robot has to perform the following
% movements:
% – Up for 5 seconds
% – Right for 8 seconds
% – Down for 6 seconds
% – Right for 6 seconds
% – Down for 7 seconds
% – Left for 7 seconds
% – Up for 5 seconds
% – Left for 7 seconds
% – Up for 2 seconds
% Total simulation time: 53 seconds, we’ll simulate for 1 whole minute
totalTime = 60; % Total simulation time (s)
numParticles = 5000; % Initial number of particles
mcl = monteCarloLocalization; % Initialize MCL object
mcl.UseLidarScan = true;
mcl.GlobalLocalization = true;
mcl.ParticleLimits = [500 numParticles];
mcl.UpdateThresholds = [0.2 0.2 0.2];
mcl.ResamplingInterval = 1;
motionModel = odometryMotionModel; % Initialise motion model
motionModel.Noise = [0.05 0.05 0.05 0.05]; % Add motion error
mcl.MotionModel = motionModel;
sensorModel = likelihoodFieldSensorModel;
sensorModel.Map = map;
sensorModel.SensorLimits = [0 lidar_range];
sensorModel.NumBeams = num_beams;
mcl.SensorModel = sensorModel;
% Simulation loop
for t=0:timeStep:totalTime
if t < 5
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
elseif (t>=5)&&(t<13)
robot_pose(1) = robot_pose(1) + lin_vel * timeStep;
robot_pose(3) = 0;
elseif (t>=13)&&(t<19)
robot_pose(2) = robot_pose(2) – lin_vel * timeStep;
robot_pose(3) = 270;
elseif (t>=19)&&(t<25)
robot_pose(1) = robot_pose(1) + lin_vel * timeStep;
robot_pose(3) = 0;
elseif (t>=25)&&(t<32)
robot_pose(2) = robot_pose(2) – lin_vel * timeStep;
robot_pose(3) = 270;
elseif (t>=32)&&(t<39)
robot_pose(1) = robot_pose(1) – lin_vel * timeStep;
robot_pose(3) = 180;
elseif (t>=39)&&(t<44)
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
elseif (t>=44)&&(t<51)
robot_pose(1) = robot_pose(1) – lin_vel * timeStep;
robot_pose(3) = 180;
elseif (t>=51)&&(t<53)
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
else
% Robot stopped
end
% rayIntersection func. works with radians
robot_pose_radians = robot_pose;
robot_pose_radians(1,3) = robot_pose_radians(1,3)*(pi/180);
intersectionPts = rayIntersection(map,robot_pose_radians,angles,lidar_range);
ranges = sqrt((intersectionPts(:,1)-robot_pose(1)).^2+(intersectionPts(:,2)-robot_pose(2)).^2);
scan = lidarScan(ranges,angles);
[isUpdated,estimatedPose,covariance] = mcl(robot_pose,scan);
particles = getParticles(mcl);
% Update particles plot every 2 "seconds"
if mod(t,2) == 0
particles_plot = particles;
end
figure(1)
show(map)
hold on
scatter(particles_plot(:,1),particles_plot(:,2),’b.’)
hold on
scatter(estimatedPose(:,1),estimatedPose(:,2),’ro’)
hold on
scatter(intersectionPts(:,1),intersectionPts(:,2),’r.’)
hold on
scatter(robot_pose(:,1),robot_pose(:,2),’bo’)
hold off
xlim([0 19])
ylim([0 17])
end I think it’s a relatively simple tool and I can’t find my error. Here’s the code:
timeStep = 0.1; % Time step (s)
lin_vel = 1; % Linear velocity (m/s)
resolution = 1; % Map resolution (m)
initial_position = [2 10 90]; % [m m degrees]
robot_pose = initial_position;
lidar_res = 0.1; % Lidar resolution (m)
lidar_range = 5.5;
num_beams = 173;
angles = linspace(-pi/2,pi/2,num_beams);
map_matrix = [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 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1;
1 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 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];
% Binary occupancy map matrix
map = binaryOccupancyMap(map_matrix,resolution);
% In order to follow the red path, the robot has to perform the following
% movements:
% – Up for 5 seconds
% – Right for 8 seconds
% – Down for 6 seconds
% – Right for 6 seconds
% – Down for 7 seconds
% – Left for 7 seconds
% – Up for 5 seconds
% – Left for 7 seconds
% – Up for 2 seconds
% Total simulation time: 53 seconds, we’ll simulate for 1 whole minute
totalTime = 60; % Total simulation time (s)
numParticles = 5000; % Initial number of particles
mcl = monteCarloLocalization; % Initialize MCL object
mcl.UseLidarScan = true;
mcl.GlobalLocalization = true;
mcl.ParticleLimits = [500 numParticles];
mcl.UpdateThresholds = [0.2 0.2 0.2];
mcl.ResamplingInterval = 1;
motionModel = odometryMotionModel; % Initialise motion model
motionModel.Noise = [0.05 0.05 0.05 0.05]; % Add motion error
mcl.MotionModel = motionModel;
sensorModel = likelihoodFieldSensorModel;
sensorModel.Map = map;
sensorModel.SensorLimits = [0 lidar_range];
sensorModel.NumBeams = num_beams;
mcl.SensorModel = sensorModel;
% Simulation loop
for t=0:timeStep:totalTime
if t < 5
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
elseif (t>=5)&&(t<13)
robot_pose(1) = robot_pose(1) + lin_vel * timeStep;
robot_pose(3) = 0;
elseif (t>=13)&&(t<19)
robot_pose(2) = robot_pose(2) – lin_vel * timeStep;
robot_pose(3) = 270;
elseif (t>=19)&&(t<25)
robot_pose(1) = robot_pose(1) + lin_vel * timeStep;
robot_pose(3) = 0;
elseif (t>=25)&&(t<32)
robot_pose(2) = robot_pose(2) – lin_vel * timeStep;
robot_pose(3) = 270;
elseif (t>=32)&&(t<39)
robot_pose(1) = robot_pose(1) – lin_vel * timeStep;
robot_pose(3) = 180;
elseif (t>=39)&&(t<44)
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
elseif (t>=44)&&(t<51)
robot_pose(1) = robot_pose(1) – lin_vel * timeStep;
robot_pose(3) = 180;
elseif (t>=51)&&(t<53)
robot_pose(2) = robot_pose(2) + lin_vel * timeStep;
robot_pose(3) = 90;
else
% Robot stopped
end
% rayIntersection func. works with radians
robot_pose_radians = robot_pose;
robot_pose_radians(1,3) = robot_pose_radians(1,3)*(pi/180);
intersectionPts = rayIntersection(map,robot_pose_radians,angles,lidar_range);
ranges = sqrt((intersectionPts(:,1)-robot_pose(1)).^2+(intersectionPts(:,2)-robot_pose(2)).^2);
scan = lidarScan(ranges,angles);
[isUpdated,estimatedPose,covariance] = mcl(robot_pose,scan);
particles = getParticles(mcl);
% Update particles plot every 2 "seconds"
if mod(t,2) == 0
particles_plot = particles;
end
figure(1)
show(map)
hold on
scatter(particles_plot(:,1),particles_plot(:,2),’b.’)
hold on
scatter(estimatedPose(:,1),estimatedPose(:,2),’ro’)
hold on
scatter(intersectionPts(:,1),intersectionPts(:,2),’r.’)
hold on
scatter(robot_pose(:,1),robot_pose(:,2),’bo’)
hold off
xlim([0 19])
ylim([0 17])
end monte carlo localization, particle filter, lidar MATLAB Answers — New Questions
Read excel data in Simulink and convert it into an array
I am trying to import data from an excel sheet as an array into Simulink and turn it into an array. I doing it using a matlab function
the function is written as
function data = dataserial(file,range)
%#codegen
data = readmatrix(file,’Sheet1′,range);
I am getting the following error:
Function ‘readmatrix’ not supported for code generation.
I tried looking this up but I am still a bit uncertain on the excact solution.
I just wanted to enquire if anyone has done the same: Importing excel data as an array
Also, I have tried the "From Spreadsheet" block in simulink but it was not working as wellI am trying to import data from an excel sheet as an array into Simulink and turn it into an array. I doing it using a matlab function
the function is written as
function data = dataserial(file,range)
%#codegen
data = readmatrix(file,’Sheet1′,range);
I am getting the following error:
Function ‘readmatrix’ not supported for code generation.
I tried looking this up but I am still a bit uncertain on the excact solution.
I just wanted to enquire if anyone has done the same: Importing excel data as an array
Also, I have tried the "From Spreadsheet" block in simulink but it was not working as well I am trying to import data from an excel sheet as an array into Simulink and turn it into an array. I doing it using a matlab function
the function is written as
function data = dataserial(file,range)
%#codegen
data = readmatrix(file,’Sheet1′,range);
I am getting the following error:
Function ‘readmatrix’ not supported for code generation.
I tried looking this up but I am still a bit uncertain on the excact solution.
I just wanted to enquire if anyone has done the same: Importing excel data as an array
Also, I have tried the "From Spreadsheet" block in simulink but it was not working as well excel, readmatrix MATLAB Answers — New Questions
How to ploting input signal by “Digilent Analog Discovery Pro 3450”
I am currently trying to use a library to interface Analog Discovery with MATLAB. However, I am encountering an issue. I am practicing the Getting Started code from the Digilent toolbox, and there are no problems when using the Analog Discovery 2 (AD2_0) for Blocking Input Stream. However, when using the "Digilent Analog Discovery Pro 3450" (ADP3450_0), the Blocking Input Stream does not work. I have attached the plotted image below.
The plotting keeps getting stuck like the image below. How can I fix this? The code is as follows:
daqreset;
devices = daqlist %list devices
devices.DeviceInfo(1)
%Single Scan
Read data from 2 analog input channels the display data
daqreset;
dq = daq("digilent");
addinput(dq,"ADP3450_0", "ai0" ,’Voltage’);
addinput(dq,"ADP3450_0", "ai1" ,’Voltage’);
ss = read(dq)
%Blocking Input Stream
Read data for 2 analog input channels then plot
daqreset;
dq = daq("digilent");
addinput(dq, "ADP3450_0", "ai0" ,’Voltage’);
addinput(dq, "ADP3450_0", "ai1" ,’Voltage’);
systemRate = 1e8; %100Mhz max frequency for AD2_0
rDiv = uint32(1) % Rate divisor is expected to be integer
rate = systemRate/rDiv;
dq.Rate = rate; % 100Mhz
data = read(dq, 100);
figure(1);
title(‘Blocking Input Stream’);
xlabel(‘Samples’);
ylabel(‘Volts’);
plot([data.ADP3450_0_ai0, data.ADP3450_0_ai1]);I am currently trying to use a library to interface Analog Discovery with MATLAB. However, I am encountering an issue. I am practicing the Getting Started code from the Digilent toolbox, and there are no problems when using the Analog Discovery 2 (AD2_0) for Blocking Input Stream. However, when using the "Digilent Analog Discovery Pro 3450" (ADP3450_0), the Blocking Input Stream does not work. I have attached the plotted image below.
The plotting keeps getting stuck like the image below. How can I fix this? The code is as follows:
daqreset;
devices = daqlist %list devices
devices.DeviceInfo(1)
%Single Scan
Read data from 2 analog input channels the display data
daqreset;
dq = daq("digilent");
addinput(dq,"ADP3450_0", "ai0" ,’Voltage’);
addinput(dq,"ADP3450_0", "ai1" ,’Voltage’);
ss = read(dq)
%Blocking Input Stream
Read data for 2 analog input channels then plot
daqreset;
dq = daq("digilent");
addinput(dq, "ADP3450_0", "ai0" ,’Voltage’);
addinput(dq, "ADP3450_0", "ai1" ,’Voltage’);
systemRate = 1e8; %100Mhz max frequency for AD2_0
rDiv = uint32(1) % Rate divisor is expected to be integer
rate = systemRate/rDiv;
dq.Rate = rate; % 100Mhz
data = read(dq, 100);
figure(1);
title(‘Blocking Input Stream’);
xlabel(‘Samples’);
ylabel(‘Volts’);
plot([data.ADP3450_0_ai0, data.ADP3450_0_ai1]); I am currently trying to use a library to interface Analog Discovery with MATLAB. However, I am encountering an issue. I am practicing the Getting Started code from the Digilent toolbox, and there are no problems when using the Analog Discovery 2 (AD2_0) for Blocking Input Stream. However, when using the "Digilent Analog Discovery Pro 3450" (ADP3450_0), the Blocking Input Stream does not work. I have attached the plotted image below.
The plotting keeps getting stuck like the image below. How can I fix this? The code is as follows:
daqreset;
devices = daqlist %list devices
devices.DeviceInfo(1)
%Single Scan
Read data from 2 analog input channels the display data
daqreset;
dq = daq("digilent");
addinput(dq,"ADP3450_0", "ai0" ,’Voltage’);
addinput(dq,"ADP3450_0", "ai1" ,’Voltage’);
ss = read(dq)
%Blocking Input Stream
Read data for 2 analog input channels then plot
daqreset;
dq = daq("digilent");
addinput(dq, "ADP3450_0", "ai0" ,’Voltage’);
addinput(dq, "ADP3450_0", "ai1" ,’Voltage’);
systemRate = 1e8; %100Mhz max frequency for AD2_0
rDiv = uint32(1) % Rate divisor is expected to be integer
rate = systemRate/rDiv;
dq.Rate = rate; % 100Mhz
data = read(dq, 100);
figure(1);
title(‘Blocking Input Stream’);
xlabel(‘Samples’);
ylabel(‘Volts’);
plot([data.ADP3450_0_ai0, data.ADP3450_0_ai1]); digilent, data acquisition MATLAB Answers — New Questions
Optimization of EPANET network using GA toolkit
A pipe network distribution layout for irrigation purpose is prepared in EPANET. This network have 80 pipes and 80 nodes plus one reservoir Optimisation of pipe network for minimum cost is required to be done using GA tool box of MATLAB.
The pipe network (EPANET file i.e .inp file) needs to be coupled with MATLAB. And then the same network needs to be optimised. That means diameters of all 80 pipes need to be determined using GA tool box of MATLAB, such that it that will suffice the limiting conditions of design viz: minimum pressure head at outlet nodes must be at least 1 m and velocity in all pipes must be in between 0.6 and 1.8 m/s.
the set of commercially available diameters and there cost is available, need to use those diameters only.A pipe network distribution layout for irrigation purpose is prepared in EPANET. This network have 80 pipes and 80 nodes plus one reservoir Optimisation of pipe network for minimum cost is required to be done using GA tool box of MATLAB.
The pipe network (EPANET file i.e .inp file) needs to be coupled with MATLAB. And then the same network needs to be optimised. That means diameters of all 80 pipes need to be determined using GA tool box of MATLAB, such that it that will suffice the limiting conditions of design viz: minimum pressure head at outlet nodes must be at least 1 m and velocity in all pipes must be in between 0.6 and 1.8 m/s.
the set of commercially available diameters and there cost is available, need to use those diameters only. A pipe network distribution layout for irrigation purpose is prepared in EPANET. This network have 80 pipes and 80 nodes plus one reservoir Optimisation of pipe network for minimum cost is required to be done using GA tool box of MATLAB.
The pipe network (EPANET file i.e .inp file) needs to be coupled with MATLAB. And then the same network needs to be optimised. That means diameters of all 80 pipes need to be determined using GA tool box of MATLAB, such that it that will suffice the limiting conditions of design viz: minimum pressure head at outlet nodes must be at least 1 m and velocity in all pipes must be in between 0.6 and 1.8 m/s.
the set of commercially available diameters and there cost is available, need to use those diameters only. epanet-matlab toolkit, optimization using ga toolkit MATLAB Answers — New Questions
Difference between hist and imhist
Hi,
I have a fundamental question about the hist and imhist and an explanation will be great as I am a bit lost.
I am running this script:
I = imread(‘pout.tif’);
figure;
x=imhist(I)
imhist(I)
figure(2);
hist(x)
and when I compare the two histograms they do not look alike, why?
How can I save the imhist data so that the hist will show the same plot?
ThanksHi,
I have a fundamental question about the hist and imhist and an explanation will be great as I am a bit lost.
I am running this script:
I = imread(‘pout.tif’);
figure;
x=imhist(I)
imhist(I)
figure(2);
hist(x)
and when I compare the two histograms they do not look alike, why?
How can I save the imhist data so that the hist will show the same plot?
Thanks Hi,
I have a fundamental question about the hist and imhist and an explanation will be great as I am a bit lost.
I am running this script:
I = imread(‘pout.tif’);
figure;
x=imhist(I)
imhist(I)
figure(2);
hist(x)
and when I compare the two histograms they do not look alike, why?
How can I save the imhist data so that the hist will show the same plot?
Thanks hist, imhist MATLAB Answers — New Questions
Paramter optimization by curve fitting
I tried paramter optimization by curve fitting. I have five paramters to optimize. If possible, I would kindly like to know the below things.
I kindly would like to know, how opts could be chnaged to have good fit? I have this code. It says options are not excecuted but doesn’t give an error.
How to decide lower and upper bounds?Is it soley based on literature?(In my case, I have only one literature similar to my material, but it didnt give a good fit)
Main file
%rng default % For reproducibility
% Constants
AD = 1.2; % Air density
p0 = 101320; % Air pressure
SH = 1.4; % Specific heat of air
PN = 0.7; % Prandtl constant
c0 = 343; % Air velocity
d = 0.05; % Assuming a constant value for d
% Define the model equation
SAC_fun = @(x, f) your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d);
% Load experimental data from Excel file
exp_data = readmatrix(‘/Users/pramodya/Desktop/PhD/Modelling/MATLAB Final codes/My research/Exp.results/Done/C-1.xlsx’); % Assuming first column contains frequencies and second column contains corresponding SAC values
f = exp_data(:, 1); % Frequency data
ydata = exp_data(:, 2); % SAC experimental data
% Plot the experimental data
figure
plot(f, ydata, ‘r*’)
xlabel ‘Frequency’
ylabel ‘SAC’
title(‘Experimental Data’)
% Define the objective function for GA
objectiveFunction = @(x) sum((SAC_fun(x, f) – ydata).^2);
% Define bounds
lb = [0, 0, 0, 0.1]; % Lower bounds for parameters
ub = [2,20, 1e8, 1]; % Upper bounds for parameters
% Define GA options with a specified number of generations
opts = optimoptions(@ga,’SelectionFcn’,@selectiontournament, …
‘FitnessScalingFcn’,@fitscalingprop);
% Run GA to fit the model to the data
nvars = 4; % Number of variables to optimize
[sol, fval] = ga(objectiveFunction, nvars, [], [], [], [], lb, ub, [], options);
% Plot the fitted curve
figure
responsedata = SAC_fun(sol, f);
plot(f, ydata, ‘y*’, f, responsedata, ‘g-‘)
legend(‘Experimental Data’, ‘Fitted Curve’)
xlabel ‘Frequency’
ylabel ‘SAC’
title(‘Fitted Response’)
% Display optimized parameter values
disp(‘Optimized Parameter Values:’);
disp([‘b: ‘, num2str(sol(1))]);
disp([‘T: ‘, num2str(sol(2))]);
disp([‘FR: ‘, num2str(sol(3))]);
disp([‘P: ‘, num2str(sol(4))]);
% Calculate R-squared value
SS_total = sum((ydata – mean(ydata)).^2);
SS_residual = sum((ydata – responsedata).^2);
R_squared = 1 – (SS_residual / SS_total);
% Display R-squared value
disp([‘R-squared Value: ‘, num2str(R_squared)]);
Function file:
function SAC = your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d)
% Parameters
%VL = x(1);
b = x(1);%shape factor
T = x(2);
FR = x(3);
P = x(4);
% SAC calculation based on provided equations
s = (1/2*b)*sqrt((8*T*AD*2*pi*f)/(FR*P));
T_bes = besselj(1, s*sqrt(-1i)) ./ besselj(0, s*sqrt(-1i));
p = (T*AD/P) * (1 – (2./(s*sqrt(-1i))) .* T_bes).^-1;
T_bes2 = besselj(1, PN^0.5*s*sqrt(-1i)) ./ besselj(0, PN^0.5*s*sqrt(-1i));
K_w = (SH*p0/P) * (SH – ((2*(SH-1))./(PN^0.5*s*sqrt(-1i))) .* T_bes).^-1;
%Characteristic Impedance
Zc = sqrt(K_w.* p);
%Complex wave number
k_c = 2 * pi * f .* sqrt(p./K_w);
%Surface impedance of sample
Zs = -1i * Zc .* cot(k_c * d);
z0 = AD * c0;
%Normalised impedance
Zn = Zs / z0;
%Reflection coefficient
R = (Zs – z0) ./ (Zs + z0);
%Absorption coefficient
SAC = 1 – abs(R).^2;
end
Thank you in advance.I tried paramter optimization by curve fitting. I have five paramters to optimize. If possible, I would kindly like to know the below things.
I kindly would like to know, how opts could be chnaged to have good fit? I have this code. It says options are not excecuted but doesn’t give an error.
How to decide lower and upper bounds?Is it soley based on literature?(In my case, I have only one literature similar to my material, but it didnt give a good fit)
Main file
%rng default % For reproducibility
% Constants
AD = 1.2; % Air density
p0 = 101320; % Air pressure
SH = 1.4; % Specific heat of air
PN = 0.7; % Prandtl constant
c0 = 343; % Air velocity
d = 0.05; % Assuming a constant value for d
% Define the model equation
SAC_fun = @(x, f) your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d);
% Load experimental data from Excel file
exp_data = readmatrix(‘/Users/pramodya/Desktop/PhD/Modelling/MATLAB Final codes/My research/Exp.results/Done/C-1.xlsx’); % Assuming first column contains frequencies and second column contains corresponding SAC values
f = exp_data(:, 1); % Frequency data
ydata = exp_data(:, 2); % SAC experimental data
% Plot the experimental data
figure
plot(f, ydata, ‘r*’)
xlabel ‘Frequency’
ylabel ‘SAC’
title(‘Experimental Data’)
% Define the objective function for GA
objectiveFunction = @(x) sum((SAC_fun(x, f) – ydata).^2);
% Define bounds
lb = [0, 0, 0, 0.1]; % Lower bounds for parameters
ub = [2,20, 1e8, 1]; % Upper bounds for parameters
% Define GA options with a specified number of generations
opts = optimoptions(@ga,’SelectionFcn’,@selectiontournament, …
‘FitnessScalingFcn’,@fitscalingprop);
% Run GA to fit the model to the data
nvars = 4; % Number of variables to optimize
[sol, fval] = ga(objectiveFunction, nvars, [], [], [], [], lb, ub, [], options);
% Plot the fitted curve
figure
responsedata = SAC_fun(sol, f);
plot(f, ydata, ‘y*’, f, responsedata, ‘g-‘)
legend(‘Experimental Data’, ‘Fitted Curve’)
xlabel ‘Frequency’
ylabel ‘SAC’
title(‘Fitted Response’)
% Display optimized parameter values
disp(‘Optimized Parameter Values:’);
disp([‘b: ‘, num2str(sol(1))]);
disp([‘T: ‘, num2str(sol(2))]);
disp([‘FR: ‘, num2str(sol(3))]);
disp([‘P: ‘, num2str(sol(4))]);
% Calculate R-squared value
SS_total = sum((ydata – mean(ydata)).^2);
SS_residual = sum((ydata – responsedata).^2);
R_squared = 1 – (SS_residual / SS_total);
% Display R-squared value
disp([‘R-squared Value: ‘, num2str(R_squared)]);
Function file:
function SAC = your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d)
% Parameters
%VL = x(1);
b = x(1);%shape factor
T = x(2);
FR = x(3);
P = x(4);
% SAC calculation based on provided equations
s = (1/2*b)*sqrt((8*T*AD*2*pi*f)/(FR*P));
T_bes = besselj(1, s*sqrt(-1i)) ./ besselj(0, s*sqrt(-1i));
p = (T*AD/P) * (1 – (2./(s*sqrt(-1i))) .* T_bes).^-1;
T_bes2 = besselj(1, PN^0.5*s*sqrt(-1i)) ./ besselj(0, PN^0.5*s*sqrt(-1i));
K_w = (SH*p0/P) * (SH – ((2*(SH-1))./(PN^0.5*s*sqrt(-1i))) .* T_bes).^-1;
%Characteristic Impedance
Zc = sqrt(K_w.* p);
%Complex wave number
k_c = 2 * pi * f .* sqrt(p./K_w);
%Surface impedance of sample
Zs = -1i * Zc .* cot(k_c * d);
z0 = AD * c0;
%Normalised impedance
Zn = Zs / z0;
%Reflection coefficient
R = (Zs – z0) ./ (Zs + z0);
%Absorption coefficient
SAC = 1 – abs(R).^2;
end
Thank you in advance. I tried paramter optimization by curve fitting. I have five paramters to optimize. If possible, I would kindly like to know the below things.
I kindly would like to know, how opts could be chnaged to have good fit? I have this code. It says options are not excecuted but doesn’t give an error.
How to decide lower and upper bounds?Is it soley based on literature?(In my case, I have only one literature similar to my material, but it didnt give a good fit)
Main file
%rng default % For reproducibility
% Constants
AD = 1.2; % Air density
p0 = 101320; % Air pressure
SH = 1.4; % Specific heat of air
PN = 0.7; % Prandtl constant
c0 = 343; % Air velocity
d = 0.05; % Assuming a constant value for d
% Define the model equation
SAC_fun = @(x, f) your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d);
% Load experimental data from Excel file
exp_data = readmatrix(‘/Users/pramodya/Desktop/PhD/Modelling/MATLAB Final codes/My research/Exp.results/Done/C-1.xlsx’); % Assuming first column contains frequencies and second column contains corresponding SAC values
f = exp_data(:, 1); % Frequency data
ydata = exp_data(:, 2); % SAC experimental data
% Plot the experimental data
figure
plot(f, ydata, ‘r*’)
xlabel ‘Frequency’
ylabel ‘SAC’
title(‘Experimental Data’)
% Define the objective function for GA
objectiveFunction = @(x) sum((SAC_fun(x, f) – ydata).^2);
% Define bounds
lb = [0, 0, 0, 0.1]; % Lower bounds for parameters
ub = [2,20, 1e8, 1]; % Upper bounds for parameters
% Define GA options with a specified number of generations
opts = optimoptions(@ga,’SelectionFcn’,@selectiontournament, …
‘FitnessScalingFcn’,@fitscalingprop);
% Run GA to fit the model to the data
nvars = 4; % Number of variables to optimize
[sol, fval] = ga(objectiveFunction, nvars, [], [], [], [], lb, ub, [], options);
% Plot the fitted curve
figure
responsedata = SAC_fun(sol, f);
plot(f, ydata, ‘y*’, f, responsedata, ‘g-‘)
legend(‘Experimental Data’, ‘Fitted Curve’)
xlabel ‘Frequency’
ylabel ‘SAC’
title(‘Fitted Response’)
% Display optimized parameter values
disp(‘Optimized Parameter Values:’);
disp([‘b: ‘, num2str(sol(1))]);
disp([‘T: ‘, num2str(sol(2))]);
disp([‘FR: ‘, num2str(sol(3))]);
disp([‘P: ‘, num2str(sol(4))]);
% Calculate R-squared value
SS_total = sum((ydata – mean(ydata)).^2);
SS_residual = sum((ydata – responsedata).^2);
R_squared = 1 – (SS_residual / SS_total);
% Display R-squared value
disp([‘R-squared Value: ‘, num2str(R_squared)]);
Function file:
function SAC = your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d)
% Parameters
%VL = x(1);
b = x(1);%shape factor
T = x(2);
FR = x(3);
P = x(4);
% SAC calculation based on provided equations
s = (1/2*b)*sqrt((8*T*AD*2*pi*f)/(FR*P));
T_bes = besselj(1, s*sqrt(-1i)) ./ besselj(0, s*sqrt(-1i));
p = (T*AD/P) * (1 – (2./(s*sqrt(-1i))) .* T_bes).^-1;
T_bes2 = besselj(1, PN^0.5*s*sqrt(-1i)) ./ besselj(0, PN^0.5*s*sqrt(-1i));
K_w = (SH*p0/P) * (SH – ((2*(SH-1))./(PN^0.5*s*sqrt(-1i))) .* T_bes).^-1;
%Characteristic Impedance
Zc = sqrt(K_w.* p);
%Complex wave number
k_c = 2 * pi * f .* sqrt(p./K_w);
%Surface impedance of sample
Zs = -1i * Zc .* cot(k_c * d);
z0 = AD * c0;
%Normalised impedance
Zn = Zs / z0;
%Reflection coefficient
R = (Zs – z0) ./ (Zs + z0);
%Absorption coefficient
SAC = 1 – abs(R).^2;
end
Thank you in advance. paramter optimizaiton, ga, curve fitting MATLAB Answers — New Questions
How can I download MATLAB & Simulink R2013a student version without cd?
How can I download MATLAB & Simulink R2013a student version without cd? I bought a new computer without an optical drive, so I cannot download MATLAB using the installation disc that came when I purchased the software. Just wondering how I could download the software for a mac (I already have a key code for authorization)?How can I download MATLAB & Simulink R2013a student version without cd? I bought a new computer without an optical drive, so I cannot download MATLAB using the installation disc that came when I purchased the software. Just wondering how I could download the software for a mac (I already have a key code for authorization)? How can I download MATLAB & Simulink R2013a student version without cd? I bought a new computer without an optical drive, so I cannot download MATLAB using the installation disc that came when I purchased the software. Just wondering how I could download the software for a mac (I already have a key code for authorization)? software, download, version upgrade MATLAB Answers — New Questions
Can I use splitapply with two grouping variables and how can it be done in matlab app designer?
Hi ! 🙂
I am referring to this example:
load carsmall
whos
Name Size Bytes Class Attributes
Acceleration 100×1 800 double
Cylinders 100×1 800 double
Displacement 100×1 800 double
Horsepower 100×1 800 double
MPG 100×1 800 double
Mfg 100×13 2600 char
Model 100×33 6600 char
Model_Year 100×1 800 double
Origin 100×7 1400 char
Weight 100×1 800 double
cmdout 1×33 66 char
[g,id]=findgroups(Model_Year);
hAx=axes; hold(hAx,’on’)
splitapply(@(x,y)scatter(x,y,’filled’),Displacement,Horsepower,g)
legend(hAx,"Model Year "+string(1900+id),’location’,’northwest’)
xlabel(‘Displacement’), ylabel(‘Horsepower’)
——————————–
I have a table which is heterogenous. And I want to show the names of the first preferred salt user has chosen and second preferred salt (these are two of the column in my table) along with their properties on the x and y axis. Their properties are also the different columns in my tableHi ! 🙂
I am referring to this example:
load carsmall
whos
Name Size Bytes Class Attributes
Acceleration 100×1 800 double
Cylinders 100×1 800 double
Displacement 100×1 800 double
Horsepower 100×1 800 double
MPG 100×1 800 double
Mfg 100×13 2600 char
Model 100×33 6600 char
Model_Year 100×1 800 double
Origin 100×7 1400 char
Weight 100×1 800 double
cmdout 1×33 66 char
[g,id]=findgroups(Model_Year);
hAx=axes; hold(hAx,’on’)
splitapply(@(x,y)scatter(x,y,’filled’),Displacement,Horsepower,g)
legend(hAx,"Model Year "+string(1900+id),’location’,’northwest’)
xlabel(‘Displacement’), ylabel(‘Horsepower’)
——————————–
I have a table which is heterogenous. And I want to show the names of the first preferred salt user has chosen and second preferred salt (these are two of the column in my table) along with their properties on the x and y axis. Their properties are also the different columns in my table Hi ! 🙂
I am referring to this example:
load carsmall
whos
Name Size Bytes Class Attributes
Acceleration 100×1 800 double
Cylinders 100×1 800 double
Displacement 100×1 800 double
Horsepower 100×1 800 double
MPG 100×1 800 double
Mfg 100×13 2600 char
Model 100×33 6600 char
Model_Year 100×1 800 double
Origin 100×7 1400 char
Weight 100×1 800 double
cmdout 1×33 66 char
[g,id]=findgroups(Model_Year);
hAx=axes; hold(hAx,’on’)
splitapply(@(x,y)scatter(x,y,’filled’),Displacement,Horsepower,g)
legend(hAx,"Model Year "+string(1900+id),’location’,’northwest’)
xlabel(‘Displacement’), ylabel(‘Horsepower’)
——————————–
I have a table which is heterogenous. And I want to show the names of the first preferred salt user has chosen and second preferred salt (these are two of the column in my table) along with their properties on the x and y axis. Their properties are also the different columns in my table matlab app designer, split apply, two grouping variables MATLAB Answers — New Questions
simulink -sinx integral
As I know, integral -sinx dx is cosx
but, in sinulink there is offset !As I know, integral -sinx dx is cosx
but, in sinulink there is offset ! As I know, integral -sinx dx is cosx
but, in sinulink there is offset ! simulink MATLAB Answers — New Questions
I’m trying to solve two separate systems using ODE23 while controlling the number of indexes so that the indexes are equal for both systems, but I couldn’t. Can you help? Than
function [sol1,sol2] = two_system
clc
close all
clear all
opt = odeset(‘RelTol’,1e-3);
sol1 = ode23(@system1,[0 5],[0],opt);
sol2 = ode23(@system2,[0 5],[0],opt);
figure(1)
plot(sol1.x(1,:),sol1.y(1,:),’b’,’LineWidth’,1.5)
hold on
plot(sol2.x(1,:),sol2.y(1,:),’–r’,’LineWidth’,1.5)
L1=length(sol1.x(1,:))
L2=length(sol2.x(1,:))
end
function dydt = system1(t,y)
M1 = y(1);
F_M1= 2;
theta1=0*pi/180;
A1 = 0.01;
%%%%%%% equation_system1 %%%%%%%
dMdt=(A1*(2*pi*(F_M1))*cos(2*pi*(F_M1)*(t)+theta1));
dydt = [dMdt];
end
function dydt = system2(t,y)
M2 = y(1);
F_M2= 2;
theta2=90*pi/180;
A2 = 0.01;
%%%%%%% equation_system2 %%%%%%%
dMdt=(A2*(2*pi*(F_M2))*cos(2*pi*(F_M2)*(t)+theta2));
dydt = [dMdt];
endfunction [sol1,sol2] = two_system
clc
close all
clear all
opt = odeset(‘RelTol’,1e-3);
sol1 = ode23(@system1,[0 5],[0],opt);
sol2 = ode23(@system2,[0 5],[0],opt);
figure(1)
plot(sol1.x(1,:),sol1.y(1,:),’b’,’LineWidth’,1.5)
hold on
plot(sol2.x(1,:),sol2.y(1,:),’–r’,’LineWidth’,1.5)
L1=length(sol1.x(1,:))
L2=length(sol2.x(1,:))
end
function dydt = system1(t,y)
M1 = y(1);
F_M1= 2;
theta1=0*pi/180;
A1 = 0.01;
%%%%%%% equation_system1 %%%%%%%
dMdt=(A1*(2*pi*(F_M1))*cos(2*pi*(F_M1)*(t)+theta1));
dydt = [dMdt];
end
function dydt = system2(t,y)
M2 = y(1);
F_M2= 2;
theta2=90*pi/180;
A2 = 0.01;
%%%%%%% equation_system2 %%%%%%%
dMdt=(A2*(2*pi*(F_M2))*cos(2*pi*(F_M2)*(t)+theta2));
dydt = [dMdt];
end function [sol1,sol2] = two_system
clc
close all
clear all
opt = odeset(‘RelTol’,1e-3);
sol1 = ode23(@system1,[0 5],[0],opt);
sol2 = ode23(@system2,[0 5],[0],opt);
figure(1)
plot(sol1.x(1,:),sol1.y(1,:),’b’,’LineWidth’,1.5)
hold on
plot(sol2.x(1,:),sol2.y(1,:),’–r’,’LineWidth’,1.5)
L1=length(sol1.x(1,:))
L2=length(sol2.x(1,:))
end
function dydt = system1(t,y)
M1 = y(1);
F_M1= 2;
theta1=0*pi/180;
A1 = 0.01;
%%%%%%% equation_system1 %%%%%%%
dMdt=(A1*(2*pi*(F_M1))*cos(2*pi*(F_M1)*(t)+theta1));
dydt = [dMdt];
end
function dydt = system2(t,y)
M2 = y(1);
F_M2= 2;
theta2=90*pi/180;
A2 = 0.01;
%%%%%%% equation_system2 %%%%%%%
dMdt=(A2*(2*pi*(F_M2))*cos(2*pi*(F_M2)*(t)+theta2));
dydt = [dMdt];
end mischa kim MATLAB Answers — New Questions
Simulate Fick’s 1st and 2nd laws of mass diffusion.
Dear a person who is looking this.
I am trying to simulate the mass diffusion based on Fick’s laws on 2D plane.
The initial condition is that the concentration on the left wall is 1.
But the result seems to be odd as below.
I guess the problem is from the code for the 2nd law. So I have tried to modify that like
C(i,j) = C(i,j) …
– dt * ((Jx(i, j) – Jx(i-1,j))/(dx) …
+ (Jy(i, j) – Jy(i,j-1))/(dy));
But the result was the same.
How can I solve this problem. I really want some help..
Thank you for reading this.
======Here is my code (you can operate by just copy and paste)=======
% Parameters
Nx = 20; % Number of grid points in x-direction
Ny = 20; % Number of grid points in y-direction
D = 0.1; % Diffusion coefficient
dx = 1.0; % Spatial step in x-direction
dy = 1.0; % Spatial step in y-direction
dt = 0.1; % Time step
T = 200.0; % Total simulation time
C = zeros(Nx, Ny); % Initial concentration
Jx = zeros(Nx, Ny);
Jy = zeros(Nx, Ny);
C(:, 1) = 1.0; % Initial condition: concentrated mass at the center
% Time loop
for t = 0:dt:T
% Create a copy of the concentration array to hold the new values
% C_new = C;
% Iterate over each interior point
% flux
% Fick’s 1st law: J = -D * del c
for i = 2:Nx-1
for j = 2:Ny-1
Jx(i,j) = – D * (C(i,j) – C(i-1,j))/dx – D * (C(i+1,j) – C(i,j))/dx ;
Jy(i,j) = – D * (C(i,j) – C(i,j-1))/dy – D * (C(i,j+1) – C(i,j))/dy ;
end
end
% concentration
% Fick’s 2nd law: partial dc/dt = – div J
for i = 2:Nx-1
for j = 2:Ny-1
C(i,j) = C(i,j) …
– dt * ((Jx(i+1, j) – Jx(i-1,j))/(2*dx) …
+ (Jy(i, j+1) – Jy(i,j-1))/(2*dy));
end
end
% Optionally plot the concentration at each time step
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
colormap(‘jet’);
title(sprintf(‘Time = %.2f’, t));
drawnow;
end
% Final plot
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
title(‘Final Concentration’);
xlabel(‘X’);
ylabel(‘Y’);Dear a person who is looking this.
I am trying to simulate the mass diffusion based on Fick’s laws on 2D plane.
The initial condition is that the concentration on the left wall is 1.
But the result seems to be odd as below.
I guess the problem is from the code for the 2nd law. So I have tried to modify that like
C(i,j) = C(i,j) …
– dt * ((Jx(i, j) – Jx(i-1,j))/(dx) …
+ (Jy(i, j) – Jy(i,j-1))/(dy));
But the result was the same.
How can I solve this problem. I really want some help..
Thank you for reading this.
======Here is my code (you can operate by just copy and paste)=======
% Parameters
Nx = 20; % Number of grid points in x-direction
Ny = 20; % Number of grid points in y-direction
D = 0.1; % Diffusion coefficient
dx = 1.0; % Spatial step in x-direction
dy = 1.0; % Spatial step in y-direction
dt = 0.1; % Time step
T = 200.0; % Total simulation time
C = zeros(Nx, Ny); % Initial concentration
Jx = zeros(Nx, Ny);
Jy = zeros(Nx, Ny);
C(:, 1) = 1.0; % Initial condition: concentrated mass at the center
% Time loop
for t = 0:dt:T
% Create a copy of the concentration array to hold the new values
% C_new = C;
% Iterate over each interior point
% flux
% Fick’s 1st law: J = -D * del c
for i = 2:Nx-1
for j = 2:Ny-1
Jx(i,j) = – D * (C(i,j) – C(i-1,j))/dx – D * (C(i+1,j) – C(i,j))/dx ;
Jy(i,j) = – D * (C(i,j) – C(i,j-1))/dy – D * (C(i,j+1) – C(i,j))/dy ;
end
end
% concentration
% Fick’s 2nd law: partial dc/dt = – div J
for i = 2:Nx-1
for j = 2:Ny-1
C(i,j) = C(i,j) …
– dt * ((Jx(i+1, j) – Jx(i-1,j))/(2*dx) …
+ (Jy(i, j+1) – Jy(i,j-1))/(2*dy));
end
end
% Optionally plot the concentration at each time step
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
colormap(‘jet’);
title(sprintf(‘Time = %.2f’, t));
drawnow;
end
% Final plot
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
title(‘Final Concentration’);
xlabel(‘X’);
ylabel(‘Y’); Dear a person who is looking this.
I am trying to simulate the mass diffusion based on Fick’s laws on 2D plane.
The initial condition is that the concentration on the left wall is 1.
But the result seems to be odd as below.
I guess the problem is from the code for the 2nd law. So I have tried to modify that like
C(i,j) = C(i,j) …
– dt * ((Jx(i, j) – Jx(i-1,j))/(dx) …
+ (Jy(i, j) – Jy(i,j-1))/(dy));
But the result was the same.
How can I solve this problem. I really want some help..
Thank you for reading this.
======Here is my code (you can operate by just copy and paste)=======
% Parameters
Nx = 20; % Number of grid points in x-direction
Ny = 20; % Number of grid points in y-direction
D = 0.1; % Diffusion coefficient
dx = 1.0; % Spatial step in x-direction
dy = 1.0; % Spatial step in y-direction
dt = 0.1; % Time step
T = 200.0; % Total simulation time
C = zeros(Nx, Ny); % Initial concentration
Jx = zeros(Nx, Ny);
Jy = zeros(Nx, Ny);
C(:, 1) = 1.0; % Initial condition: concentrated mass at the center
% Time loop
for t = 0:dt:T
% Create a copy of the concentration array to hold the new values
% C_new = C;
% Iterate over each interior point
% flux
% Fick’s 1st law: J = -D * del c
for i = 2:Nx-1
for j = 2:Ny-1
Jx(i,j) = – D * (C(i,j) – C(i-1,j))/dx – D * (C(i+1,j) – C(i,j))/dx ;
Jy(i,j) = – D * (C(i,j) – C(i,j-1))/dy – D * (C(i,j+1) – C(i,j))/dy ;
end
end
% concentration
% Fick’s 2nd law: partial dc/dt = – div J
for i = 2:Nx-1
for j = 2:Ny-1
C(i,j) = C(i,j) …
– dt * ((Jx(i+1, j) – Jx(i-1,j))/(2*dx) …
+ (Jy(i, j+1) – Jy(i,j-1))/(2*dy));
end
end
% Optionally plot the concentration at each time step
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
colormap(‘jet’);
title(sprintf(‘Time = %.2f’, t));
drawnow;
end
% Final plot
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
title(‘Final Concentration’);
xlabel(‘X’);
ylabel(‘Y’); physics, diffusion, mass diffusion, differential equations MATLAB Answers — New Questions
How to set Matlab do display all calculation results?
Hi,
How can I set Matlab to display all the calculation results?
I am solving the system of ODEs with large spatial variable and time step distribution. If I set total time to 100 with time step of 1 only the results starting from t=50 are displayed.
How to fix it?
ThanksHi,
How can I set Matlab to display all the calculation results?
I am solving the system of ODEs with large spatial variable and time step distribution. If I set total time to 100 with time step of 1 only the results starting from t=50 are displayed.
How to fix it?
Thanks Hi,
How can I set Matlab to display all the calculation results?
I am solving the system of ODEs with large spatial variable and time step distribution. If I set total time to 100 with time step of 1 only the results starting from t=50 are displayed.
How to fix it?
Thanks odes, results display MATLAB Answers — New Questions
histogram
how to calculate threshold value using histogram ,that is average of two hightest velles in in histogramhow to calculate threshold value using histogram ,that is average of two hightest velles in in histogram how to calculate threshold value using histogram ,that is average of two hightest velles in in histogram histogram MATLAB Answers — New Questions
plot in designated subfigures in loops
I have the following code. I want to plotsubfigures as required in the comments in main.m. How to plot in designated subfigures in such loops below?
% main.m
tiledlayout(2,3)
for i=1:2
for j=1:3
function myfun;
end
end
% myfun.m
plot % in subplot(1,j)
plot % in subplot(2,j)I have the following code. I want to plotsubfigures as required in the comments in main.m. How to plot in designated subfigures in such loops below?
% main.m
tiledlayout(2,3)
for i=1:2
for j=1:3
function myfun;
end
end
% myfun.m
plot % in subplot(1,j)
plot % in subplot(2,j) I have the following code. I want to plotsubfigures as required in the comments in main.m. How to plot in designated subfigures in such loops below?
% main.m
tiledlayout(2,3)
for i=1:2
for j=1:3
function myfun;
end
end
% myfun.m
plot % in subplot(1,j)
plot % in subplot(2,j) plot, subplot, tilelayout, axes MATLAB Answers — New Questions
How do I install on Ubuntu?
I am a novice Ubuntu user, I know how to use Terminal. I am trying to install MATLAB 2020a on my Ubuntu machine, I have downloaded the Linux Installer which is a zip folder. What do I do from here? Thanks!I am a novice Ubuntu user, I know how to use Terminal. I am trying to install MATLAB 2020a on my Ubuntu machine, I have downloaded the Linux Installer which is a zip folder. What do I do from here? Thanks! I am a novice Ubuntu user, I know how to use Terminal. I am trying to install MATLAB 2020a on my Ubuntu machine, I have downloaded the Linux Installer which is a zip folder. What do I do from here? Thanks! 2020a, linux, ubuntu, install, installation MATLAB Answers — New Questions
Solve an ODE with the parameters defined in the function changing in a for loop in the main script
I have to cycle trough different values of stifness kp in the main script with a for loop and solve an Ode. How can I modify the kp value in the function defined used as imput to an ode45 in the main script?
Thanks in advance
function xdot = Time_domain_system(t,x)
theta_ref = 1;
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
kp = 50;
%equivalent parameters
mx = J+m*L^2/4;
cx = c*L^2/4;
kx = k*L^2/4 – m*g*L/2;
xdot(1) = x(2);
xdot(2) = kp*theta_ref/mx -(cx/mx)*x(2) – ((kx + kp)/mx)*x(1);
xdot = xdot’;
endI have to cycle trough different values of stifness kp in the main script with a for loop and solve an Ode. How can I modify the kp value in the function defined used as imput to an ode45 in the main script?
Thanks in advance
function xdot = Time_domain_system(t,x)
theta_ref = 1;
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
kp = 50;
%equivalent parameters
mx = J+m*L^2/4;
cx = c*L^2/4;
kx = k*L^2/4 – m*g*L/2;
xdot(1) = x(2);
xdot(2) = kp*theta_ref/mx -(cx/mx)*x(2) – ((kx + kp)/mx)*x(1);
xdot = xdot’;
end I have to cycle trough different values of stifness kp in the main script with a for loop and solve an Ode. How can I modify the kp value in the function defined used as imput to an ode45 in the main script?
Thanks in advance
function xdot = Time_domain_system(t,x)
theta_ref = 1;
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
kp = 50;
%equivalent parameters
mx = J+m*L^2/4;
cx = c*L^2/4;
kx = k*L^2/4 – m*g*L/2;
xdot(1) = x(2);
xdot(2) = kp*theta_ref/mx -(cx/mx)*x(2) – ((kx + kp)/mx)*x(1);
xdot = xdot’;
end ode45, ordinary differential equation controls, control MATLAB Answers — New Questions
About While Loops problems
Is there anyway to end the program by entering "quit"? As the program keeps going despite displays "Exiting…"
The following is the syntax.
function Individual_Project_Final()
disp(‘Welcome to Unit Converter interface! ‘);
disp(‘ ‘);
while true
% Accept user input
input_str = input(‘Enter a number to be converted (or type "quit" to exit program): ‘, ‘s’);
% Check if the user wants to quit
if strcmp(input_str, ‘quit’)
disp(‘Exiting…’)
break;
end
% Convert input string to a number
numin = str2double(input_str);
% Check if the input is a valid number
if isnan(numin)
disp(‘Input invalid. Please enter a whole or decimal number (or type "quit to exit program): ‘);
else
disp(‘Input valid. Now please select the unit conversion type (or type "quit to exit program): ‘);
break;
end
end
% Display the menu for selecting conversion type
disp(‘ ‘)
disp(‘Select conversion type:’)
disp(‘1) Celsius to Fahrenheit’)
disp(‘2) Fahrenheit to Celsius’)
disp(‘3) Centimeters to Inches’)
disp(‘4) Inches to Centimeters’)
disp(‘5) Meters to Foot’)
disp(‘6) Foot to Meters’)
disp(‘7) Kilometers to Miles’)
disp(‘8) Miles to Kilometers’)
disp(‘9) Grams to Ounces’)
disp(’10) Ounces to Grams’)
disp(’11) Kilograms to Pounds’)
disp(’12) Pounds to Kilograms’)
disp(’13) Tonnes to Tons’)
disp(’14) Tons to Tonnes’)
while true
% Accept user input
input_str = input(‘Please select conversion type (1-14) (or type "quit" to exit): ‘, ‘s’);
% Check if the user wants to quit
if strcmp(input_str, ‘quit’)
disp(‘Exiting…’);
break;
end
% Convert input string to a number
convtype = str2double(input_str);
% Check if input is a valid conversion type
if ~ismember(convtype, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
disp(‘Conversion type invalid. Please enter a number from 1 to 14 (or type "quit to exit program): ‘);
else
disp(‘Input valid. We are now finalizing the resluts based on your input. ‘);
break;
end
end
% Call the user-defined function from Module 2
[numout,unit] = MyUnitConverter(numin,convtype);
% Output generation
disp(‘ ‘);
fprintf(‘The converted number is: %.2f %sn’, numout, unit);
endIs there anyway to end the program by entering "quit"? As the program keeps going despite displays "Exiting…"
The following is the syntax.
function Individual_Project_Final()
disp(‘Welcome to Unit Converter interface! ‘);
disp(‘ ‘);
while true
% Accept user input
input_str = input(‘Enter a number to be converted (or type "quit" to exit program): ‘, ‘s’);
% Check if the user wants to quit
if strcmp(input_str, ‘quit’)
disp(‘Exiting…’)
break;
end
% Convert input string to a number
numin = str2double(input_str);
% Check if the input is a valid number
if isnan(numin)
disp(‘Input invalid. Please enter a whole or decimal number (or type "quit to exit program): ‘);
else
disp(‘Input valid. Now please select the unit conversion type (or type "quit to exit program): ‘);
break;
end
end
% Display the menu for selecting conversion type
disp(‘ ‘)
disp(‘Select conversion type:’)
disp(‘1) Celsius to Fahrenheit’)
disp(‘2) Fahrenheit to Celsius’)
disp(‘3) Centimeters to Inches’)
disp(‘4) Inches to Centimeters’)
disp(‘5) Meters to Foot’)
disp(‘6) Foot to Meters’)
disp(‘7) Kilometers to Miles’)
disp(‘8) Miles to Kilometers’)
disp(‘9) Grams to Ounces’)
disp(’10) Ounces to Grams’)
disp(’11) Kilograms to Pounds’)
disp(’12) Pounds to Kilograms’)
disp(’13) Tonnes to Tons’)
disp(’14) Tons to Tonnes’)
while true
% Accept user input
input_str = input(‘Please select conversion type (1-14) (or type "quit" to exit): ‘, ‘s’);
% Check if the user wants to quit
if strcmp(input_str, ‘quit’)
disp(‘Exiting…’);
break;
end
% Convert input string to a number
convtype = str2double(input_str);
% Check if input is a valid conversion type
if ~ismember(convtype, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
disp(‘Conversion type invalid. Please enter a number from 1 to 14 (or type "quit to exit program): ‘);
else
disp(‘Input valid. We are now finalizing the resluts based on your input. ‘);
break;
end
end
% Call the user-defined function from Module 2
[numout,unit] = MyUnitConverter(numin,convtype);
% Output generation
disp(‘ ‘);
fprintf(‘The converted number is: %.2f %sn’, numout, unit);
end Is there anyway to end the program by entering "quit"? As the program keeps going despite displays "Exiting…"
The following is the syntax.
function Individual_Project_Final()
disp(‘Welcome to Unit Converter interface! ‘);
disp(‘ ‘);
while true
% Accept user input
input_str = input(‘Enter a number to be converted (or type "quit" to exit program): ‘, ‘s’);
% Check if the user wants to quit
if strcmp(input_str, ‘quit’)
disp(‘Exiting…’)
break;
end
% Convert input string to a number
numin = str2double(input_str);
% Check if the input is a valid number
if isnan(numin)
disp(‘Input invalid. Please enter a whole or decimal number (or type "quit to exit program): ‘);
else
disp(‘Input valid. Now please select the unit conversion type (or type "quit to exit program): ‘);
break;
end
end
% Display the menu for selecting conversion type
disp(‘ ‘)
disp(‘Select conversion type:’)
disp(‘1) Celsius to Fahrenheit’)
disp(‘2) Fahrenheit to Celsius’)
disp(‘3) Centimeters to Inches’)
disp(‘4) Inches to Centimeters’)
disp(‘5) Meters to Foot’)
disp(‘6) Foot to Meters’)
disp(‘7) Kilometers to Miles’)
disp(‘8) Miles to Kilometers’)
disp(‘9) Grams to Ounces’)
disp(’10) Ounces to Grams’)
disp(’11) Kilograms to Pounds’)
disp(’12) Pounds to Kilograms’)
disp(’13) Tonnes to Tons’)
disp(’14) Tons to Tonnes’)
while true
% Accept user input
input_str = input(‘Please select conversion type (1-14) (or type "quit" to exit): ‘, ‘s’);
% Check if the user wants to quit
if strcmp(input_str, ‘quit’)
disp(‘Exiting…’);
break;
end
% Convert input string to a number
convtype = str2double(input_str);
% Check if input is a valid conversion type
if ~ismember(convtype, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
disp(‘Conversion type invalid. Please enter a number from 1 to 14 (or type "quit to exit program): ‘);
else
disp(‘Input valid. We are now finalizing the resluts based on your input. ‘);
break;
end
end
% Call the user-defined function from Module 2
[numout,unit] = MyUnitConverter(numin,convtype);
% Output generation
disp(‘ ‘);
fprintf(‘The converted number is: %.2f %sn’, numout, unit);
end while loop MATLAB Answers — New Questions
Unable to outerjoin two tables, Error message of ‘Left and right key variables have incompatible types’
Hi, it seems that no matter what I try, my left and right key variables will always be of incompatible types.
It’s gotten so frustrating that I literally renamed Table A to have the exact same category names as Table B.
For example, each Table now has the same categories, One, Two, Three, Four, Five, and Six
But for some reason, matlab.mathworks.com keeps saying that ‘Ordinal arrays must have the same categories, including their order’ – which I believe I already made sure the tables were.
Dear community, what exactly am I doing wrong? Is MATLAB Online not compatible with joining tables?
My version is 2024a – Thanks.Hi, it seems that no matter what I try, my left and right key variables will always be of incompatible types.
It’s gotten so frustrating that I literally renamed Table A to have the exact same category names as Table B.
For example, each Table now has the same categories, One, Two, Three, Four, Five, and Six
But for some reason, matlab.mathworks.com keeps saying that ‘Ordinal arrays must have the same categories, including their order’ – which I believe I already made sure the tables were.
Dear community, what exactly am I doing wrong? Is MATLAB Online not compatible with joining tables?
My version is 2024a – Thanks. Hi, it seems that no matter what I try, my left and right key variables will always be of incompatible types.
It’s gotten so frustrating that I literally renamed Table A to have the exact same category names as Table B.
For example, each Table now has the same categories, One, Two, Three, Four, Five, and Six
But for some reason, matlab.mathworks.com keeps saying that ‘Ordinal arrays must have the same categories, including their order’ – which I believe I already made sure the tables were.
Dear community, what exactly am I doing wrong? Is MATLAB Online not compatible with joining tables?
My version is 2024a – Thanks. outerjoin, table, ordinal arrays, incompatible types MATLAB Answers — New Questions
problem in function inputs
Hi guys. I have this function that gives me the equations for my DAE system.
function S = Sec_model_fun(t,x,Q0,T1)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%FGout = x(16); %Gas outlet[mol/s]
%% Constants
Kc_Total_gases = 1;
tauI_Total_gases =1;
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
%T1=230+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;…
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
kL=1; kH=1;
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
S9= VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
S10= VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
S11= VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
S12= VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
S13= VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
I want to solve the S1,S2,..,S15 equations for different Q0 and T1. I have tried this code:
clear variables
clc
%Initial Condition
Q0=350; T1=230;
%S = @Sec_model_fun_for_optimization;
%S1 = S(Q0,T1);
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0 Q0 T1];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@Sec_model_fun_for_optimization,[0 18000],y0);
but I have got this error:
Not enough input arguments.
Error in Sec_model_fun_for_optimization (line 27)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
As you can see I have done it in a wrong way and MATLAB can’t recognize Q0 and T1 as my function inputs. What should I do to solve this problem??Hi guys. I have this function that gives me the equations for my DAE system.
function S = Sec_model_fun(t,x,Q0,T1)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%FGout = x(16); %Gas outlet[mol/s]
%% Constants
Kc_Total_gases = 1;
tauI_Total_gases =1;
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
%T1=230+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;…
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
kL=1; kH=1;
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
S9= VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
S10= VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
S11= VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
S12= VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
S13= VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
I want to solve the S1,S2,..,S15 equations for different Q0 and T1. I have tried this code:
clear variables
clc
%Initial Condition
Q0=350; T1=230;
%S = @Sec_model_fun_for_optimization;
%S1 = S(Q0,T1);
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0 Q0 T1];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@Sec_model_fun_for_optimization,[0 18000],y0);
but I have got this error:
Not enough input arguments.
Error in Sec_model_fun_for_optimization (line 27)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
As you can see I have done it in a wrong way and MATLAB can’t recognize Q0 and T1 as my function inputs. What should I do to solve this problem?? Hi guys. I have this function that gives me the equations for my DAE system.
function S = Sec_model_fun(t,x,Q0,T1)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%FGout = x(16); %Gas outlet[mol/s]
%% Constants
Kc_Total_gases = 1;
tauI_Total_gases =1;
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
%T1=230+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;…
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
kL=1; kH=1;
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
S9= VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
S10= VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
S11= VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
S12= VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
S13= VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
I want to solve the S1,S2,..,S15 equations for different Q0 and T1. I have tried this code:
clear variables
clc
%Initial Condition
Q0=350; T1=230;
%S = @Sec_model_fun_for_optimization;
%S1 = S(Q0,T1);
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0 Q0 T1];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@Sec_model_fun_for_optimization,[0 18000],y0);
but I have got this error:
Not enough input arguments.
Error in Sec_model_fun_for_optimization (line 27)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
As you can see I have done it in a wrong way and MATLAB can’t recognize Q0 and T1 as my function inputs. What should I do to solve this problem?? function, ode MATLAB Answers — New Questions