/r/matlab

Photograph via //r/matlab

Official MATLAB subreddit


MATLAB news, code tips and tricks, questions, and discussion! We are here to help, but won't do your homework or help you pirate software.

The effort you put into asking a question is often matched by the quality of our answers.

r/matlab discord channel


Sort By Topic

    Homework     Technical

    Code Share   News

    Tips             Misc


Places to learn Matlab
Matlab Resources

Try saturnapi to share and run MATLAB code in a web browser!

If you want flair simply Message the mods


/r/matlab

58,828 Subscribers

0

Help Needed: MATLAB/Simulink Modeling and Simulation of Aircraft Fuel Management System

Hi everyone,

I’m currently working on a challenging assignment for my Aerospace Systems Design module and could really use some help from anyone with experience in MATLAB and Simulink, especially those familiar with aircraft fuel systems.

Assignment Overview: The task involves analyzing, simulating, verifying, and validating an aircraft fuel management system (ATA 28). The goal is to model the system using MATLAB and Simulink, focusing on several key areas such as system architecture, fuel gauging, control systems, and failure management.

Key Components of the System:

  • Single fuel tank in each wing.
  • Main and standby pumps for each engine.
  • Fuel shut-off valves.
  • Cross-feed capability between tanks.
  • Ground refueling system.
  • Fuel transfer between wings.
  • Fuel gauging using capacitance probes.

What I Need Help With:

  1. Fuel Gauging Calculations:
    • Accurate methods for volumetric and mass gauging of fuel using capacitance probe data.
    • Alternate methods (e.g., pressure sensors, ultrasonic sensors) for verifying fuel quantity.
  2. GUI Development:
    • Creating an effective GUI for the pilot’s flight deck to display system information.
    • Ensuring real-time updates and interactive controls.
  3. Control Systems Implementation:
    • Implementing control logic for pumps and valves.
    • Handling system failures and maintaining balance (Center of Gravity).

Example Calculation for Reference: Here’s an example of the calculations I’m working with:

  1. Fuel Height Calculation (using capacitance probes):
    • Probe P1 (bottom, fuselage end): Wetted ratio r1=0.8r_1 = 0.8r1​=0.8, length = 0.5 meters
    • Probe P2 (bottom, wing-tip end): Wetted ratio r2=0.6r_2 = 0.6r2​=0.6, length = 0.4 meters
    • Height h1=0.4h_1 = 0.4h1​=0.4 meters, h2=0.24h_2 = 0.24h2​=0.24 meters
    • Average height havg=0.32h_{avg} = 0.32havg​=0.32 meters
  2. Volume and Mass Calculation:
    • Volume V=5 m×1 m×0.32 m=1.6 m3V = 5 \, \text{m} \times 1 \, \text{m} \times 0.32 \, \text{m} = 1.6 \, \text{m}^3V=5m×1m×0.32m=1.6m3
    • Mass m=1.6 m3×800 kg/m3=1280 kgm = 1.6 \, \text{m}^3 \times 800 \, \text{kg/m}^3 = 1280 \, \text{kg}m=1.6m3×800kg/m3=1280kg

https://preview.redd.it/3cpqfklwdo9d1.png?width=1011&format=png&auto=webp&s=1ddc76ea979e46af96d656d51255f342f4a0a7bb

Any guidance, sample code, or references to similar projects would be incredibly helpful. Thanks in advance for your assistance!

TL;DR: Need help with MATLAB/Simulink modeling of an aircraft fuel management system for a coursework assignment. Looking for advice on system analysis, modeling strategy, fuel gauging calculations, GUI development, and control systems implementation.

1 Comment
2024/06/30
09:11 UTC

1

The question about dir function and sorting string/char arrays

Hello everyone,

I've got a question for you. While using dir function in Matlab, I needed to define my variable name which is given below:

 ListName = dir([FilePath, '/', 'VariableName*', 'Property*']);

My variable names in the file are a bit long. Thus, I am trying to search the names in two separations such as "VariableName..." and "Poperty..." using "*".

Is there a common name for such indexes to define sorting the string/character arrays?

Good day!

1 Comment
2024/06/29
18:44 UTC

2

having some problems with packages (?) that i dont know how to solve

1 Comment
2024/06/29
16:18 UTC

2

Lyapunov approach to full-state feedback control on LTV system

Hi everyone. I'm wondering if it's possible to solve the Lyapunov equation on Matlab for a full state feedback control in order to find the stabilizing K for a linear time-varying system, and in case how? I can't find any source about this.

Thanks for the attention.

0 Comments
2024/06/28
15:25 UTC

0

Do anyone know how to use the amdf?

5 Comments
2024/06/28
10:45 UTC

1

Curve Fitting toolbox limitations

Is there a way to modify the graph view in the Curve Fitting toolbox (e.g. plot as loglog)? If not, why not‽

3 Comments
2024/06/28
03:25 UTC

3

Need help with simulink!

Hey guys, I feel like I'm being stupid for having trouble with such a simple thing.

I'm making a smart lock system and testing out a PoC on simulink. As you can see in the screenshot, it will have an NFC embedded system into it. The thing is, My C block function won't read the variables as it should. I have 2 constant blocks with "100" and "101" but my diagnostic will always print "0" no matter the values I change the variable to. So yeah I'd appreciate if anyone could give me a hand <3 First time using Simulink don't judge me ;-;

https://preview.redd.it/ou4uo8o1i59d1.png?width=153&format=png&auto=webp&s=1ea5ed587221a23ff3fa44b965409f34a0cd603c

https://preview.redd.it/khirl4qyh59d1.png?width=519&format=png&auto=webp&s=b277a3b648342241a933d0c2dcb864aae8ebb3d1

5 Comments
2024/06/27
17:41 UTC

2

Index access problem with parfor loop

m = 1e4;
n = 1e3;

Mat = rand(m,n);
Changes = rand(m,n);
Mat_copy = Mat;

parfor idx = [1:n/2; 1+n/2:n]
  % changes some values of Mat_copy
  Mat_copy(:,idx) = Changes(:,idx);

  % some other calculations

  % changes values of Mat_copy back
  Mat_copy(:,idx) = Mat(:,idx);
end

This is my example code.

As a simple for loop it works, but if I switch to parfor it no longer works.

I know what the error is, but wanted to ask if anyone has an idea how I can get something like this to work in parallel anyway

2 Comments
2024/06/27
17:02 UTC

19

New Onramp courses available Core MATLAB Skills and Simulink Fundamentals

MATLAB Onramp and SImulink Onramp have been very popular and people wanted to continue learning using the same format. Unfortunately, we only had a longer-format training courses.

Now we have new "learning paths" that let you take those courses in small chunk at a time.

Core MATLAB Skills

https://matlabacademy.mathworks.com/details/core-matlab-skills/lpmlcms

Simulink Fundamentals

https://matlabacademy.mathworks.com/details/simulink-fundamentals/slbe

1 Comment
2024/06/27
13:13 UTC

0

Start a timer (from 0 to 10) after a variable is set to 1

13 Comments
2024/06/26
20:13 UTC

1

Intensive

Does anyone know of any in person intensives to learn MatLab? Ideally I'm imagining a 1 month intensive code every day for like 8 hours come out knowing MatLab.

Ps I live in the uk so something uk based would be ideal.

4 Comments
2024/06/26
06:47 UTC

0

how can i disable the sign in button and stop the popup window while this fuction is useless for me

1 Comment
2024/06/26
03:52 UTC

11

Any tips for starting matlab?

I have no coding experience whatsoever. My lab only uses matlab for analyzing its primary form of data collection. I’ve tried attempting the analyses following the GUI but it makes no sense to me. I tried a guided matlab workshop to help but it confused me as well because it wasn’t using data that is meaningful to me and the functions just overwhelmed me. Idk how to overcome this as I am expected to create a “script” for my project soon. Does anyone have tips to get more familiar with it? Perhaps an online resource or workshop. I’d need it to be dumbed down as much as possible. I’ve dabbled with this every couple weeks over the last year. I get frustrated from not knowing what code or functions means and then I set it aside. I’d really like to tackle this issue head on. I appreciate any help!!

9 Comments
2024/06/26
02:56 UTC

1

Help Making a Wire Length calculator

Imagine I have a wire with a set diameter in milimeters, and I want to make this wire into a multi layer multi column coil, with a set inner radius and a set outer radius. How can I make a calculator that calculates the length of the wire used. I made a code in Scilab but it turned out wrong so I came here.

4 Comments
2024/06/25
21:49 UTC

2

Guidance for GRBL CNC Measurement Device project

Good Morning Everyone,

First of all thank you to anyone that can provide me with guidance.

Currently working on a home project that involves the use of a customized CNC machine with a gantry that has a measurement device attached.

The basis of the MATLAB script is too:

  • Start and initialize devices (GRBL v1.1F CNC Shield V3 and four channel oscilloscope)
  • MATLAB to 'setup' both devices where GRBL homes and oscilloscope is setup
  • Begin experiment after the above tasks are completed and the "experiment" parameters are set.
  • GCODE is generated based on the algorithm, for each line of GCODE positional instructions that are sent, halt the system until the measurements from the oscilloscope are complete and saved to a csv with gcode XY coordinates.
  • Move to next position and proceed to take measurements, continue until finished.
  • Gantry moves for example for every 10 positions on X, Y moves one times. Grid plane is around 40cm x 40cm.

I have a rough working script, but I cannot fathom on how to do develop the generate GCODE, send gcode, wait for measurements, proceed section. I'm using it's a nested for loop with gcode generator being the top most loop with the measurements for loop nested?

Any guidance would be much appreciated, might be my ADHD but I cannot see the wood in the trees at the moment!

Here is my current script;

(Disclosure, the script is roughly 40% ChatGPT and the rest myself with GPT fixing my errors)

function WIP_VLC() 
[deviceObj, serialObj] = startDevices(); % Start device 
[gcodeCommands] = grblSerial(serialObj); % GRBL serial monitor 
%%positionSetup(deviceObj, serialObj, gcodeCommands); % Positioning Measurement Setup closeDevices(deviceObj, serialObj); % Close devices and delete end 

%% Initialise GRLB v1.1f & Visa device % Configure grbl and instrument settings 

function [deviceObj, serialObj] = startDevices() 
disp(['Initialising VLC Positioning Measurement \n']); 

% VISA-TCIP Settings 
% %%%%%%%%%%%%%%%%%%%%%%%% 

fprintf('Initialising DS1054Z...'); 
instrreset; 
deviceObj = visa('ni', 'TCPIP::169.254.237.158::INSTR'); 
pause(1); 

% Set the buffer size deviceObj.InputBufferSize = 60e6; 
deviceObj.OutputBufferSize = 20e6; 

% Set the timeout value 
deviceObj.Timeout = 10; 

% Set the Byte order 
deviceObj.ByteOrder = 'littleEndian'; 

% Open the connection 
fopen(deviceObj); 

% Identify the device 
fprintf(deviceObj, '*IDN?'); 
idn = fscanf(deviceObj); 
disp(['Oscilloscope Identification: ', idn]); 

% Reset & clear the oscilloscope 

fprintf(deviceObj, '*RST'); 
pause(0.1); 
fprintf(deviceObj, '*CLS'); 
pause(0.1); 

% Switch four channels ON 
channels = [1, 2, 3, 4]; 
ch = channels 
fprintf(deviceObj, [':CHANNEL', num2str(ch), ':DISPLAY ON']); 
pause(0.1); 
end 

% device format setup 
fprintf(deviceObj, ':SYST:HEAD OFF'); 
% Disable the header 
fprintf(deviceObj, ':ACQ:TYPE NORM'); 
% Set the acquisition type to normal 
fprintf(deviceObj, ':ACQ:MEMD LONG'); 
% Set the number of acquisition points 
fprintf(deviceObj, ':WAV:FORM WORD'); 
% Set the waveform format to WORD 
fprintf(deviceObj, ':WAV:BYTE LSB'); 

% Set the byte order 
% Preamble 
% Get the preamble block which contains information about the waveform 

preambleBlock = query(deviceObj, ':WAV:PRE?');
% Parse the preamble block 
preambleBlock = regexp(preambleBlock, ',', 'split'); 
points = str2double(preambleBlock{3}); 
% Number of data points 
XIncrement = str2double(preambleBlock{5}); 
% Time interval between points 
XOrigin = str2double(preambleBlock{6}); 
% The start time of the waveform 
YIncrement = str2double(preambleBlock{8}); 
% Voltage increment per bit 
YOrigin = str2double(preambleBlock{9}); 
% The start voltage of the waveform 
YReference = str2double(preambleBlock{10}); 
% Reference voltage disp('DS1054Z is ready...\n'); 
pause(0.2); 

% GRBL Settings 
% %%%%%%%%%%%%%%%%%%%%% 

fprintf('Initialising GRBL...\n'); 
%%% Checks if port is already opened %%% 
serialPort = instrfind('Port', 'COM3'); 
% Close serial port 
if ~isempty(serialPort) 
if strcmp(serialPort.Status, 'open') 
fclose(serialPort); 
disp('Serial port COM closed.'); 
end 
end 

% Delete serial port object 
if ~isempty(serialPort) 
delete(serialPort); 
end 

% Visa Device Parameters 
serialPort = 'COM3'; 
baudRate = 115200; 
serialObj = serial(serialPort, 'BaudRate', baudRate, 'Terminator', 'LF'); 
fopen(serialObj); 
pause(2); 
% Initialization 
fprintf(serialObj, '\r\n'); 
pause(2); 

% Wait for GRBL to initialize
% Clear GRBL's buffer by reading the output 

while serialObj.BytesAvailable > 0 fscanf(serialObj); 
end 
fprintf('GRBL is ready... \n'); 
fprintf('System initialised... \n'); 
end 

%% GRBL Function function [gcodeCommands] = grblSerial(serialObj) 
% List of example G-code commands for testing 

gcodeCommands = { 
'G28', % Move to home position 'G92 X0 Y0 Z0',
% Set work coordinates to zero 'G0 X10 Y10 Z0', 
% Rapid move to X=10, Y=10, Z=0 'G1 X20 Y20 Z-1 F100', 
% Linear move to X=20, Y=20, Z=-1 at feed rate 100 'G0 X0', 
% Move to X=0 'G0 Y0', 
% Move to Y=0 'G0 Z0', 
% Move to Z=0 'F200', 
% Set feed rate to 200 mm/min '?', 
% Status report '$$', 
% List GRBL settings 'M3 S1000', 
% Turn on spindle clockwise at 1000 RPM 'M5', 
% Turn off spindle 'M8', 
% Turn on coolant 'M9', 
% Turn off coolant 'G1 X10 Y0', 
% Move to X=10 'G1 X10 Y10', 
% Move to Y=10 'G1 X0 Y10', 
% Move to X=0 'G1 X0 Y0' 
% Move to starting position }; 
% Function to check limit status function checkLimits(s) fprintf(s, '$$'); 
% Send command to list GRBL settings 

pause(1); 
if s.BytesAvailable > 0 
response = fscanf(s); 
disp('Current GRBL settings including limit settings:'); 
disp(response); 
end 
end 

% Main interactive loop for sending G-code commands 
while true 
% Prompt the user for a G-code command gcodeCommand = input(['Enter G-code command (' ... '"test" to run test commands, ' ... '"limits" to check limits, ' ... '"exit" to quit)' ... 'Command: '], 's'); 

% Exit the loop if the user types "exit" 
if strcmpi(gcodeCommand, 'exit') 
break; 
end 

% Run predefined test commands if the user types "test"
if strcmpi(gcodeCommand, 'test') 
for i = 1:length(gcodeCommands) 
fprintf(serialObj, '%s\n', gcodeCommands{i}); 
pause(1); 

% Give some time for GRBL to process the command 
% Read and display the response from GRBL 

while serialObj.BytesAvailable > 0 
response = fscanf(serialObj); 
disp(['GRBL: ', response]);
end 
end 

% Check limit status if the user types "limits" 
elseif strcmpi(gcodeCommand, 'limits') checkLimits(serialObj); 
else 

% Send the user-entered G-code command to GRBL 
fprintf(serialObj, '%s\n', gcodeCommand); 
pause(1); 

% Give some time for GRBL to process the command 
% Read and display the response from GRBL 

while serialObj.BytesAvailable > 0 
response = fscanf(serialObj); 
disp(['GRBL: ', response]); 
end 
end 
end 

%% Experiment/Measurement Parameters 
% function positionSetup(deviceObj, serialObj, gcodeCommands) 
% end

%% Close, delete and clear Objects 
function closeDevices(deviceObj, serialObj) 
fclose(deviceObj); 
delete(deviceObj); 
clear deviceObj; 
fprintf('deviceObj deleted... \n'); 

fclose(serialObj); 
delete(serialObj); 
clear serialObj; 
fprintf('serialObj deleted... \n'); 
end 

%% Useful Commands % fprintf(deviceObj, ':CHANNEL2:DISPLAY OFF');
0 Comments
2024/06/25
10:02 UTC

1

Help installing program

I am new to MatLab can someone please help me install https://github.com/ZurichNCH/Automatic-High-Frequency-Oscillation-Detector I’ve been trying to understand how for an extremely long time. Could you please give detailed instructions?

2 Comments
2024/06/25
01:56 UTC

1

Building a Participant Study GUI in MATLAB

Hello everyone,

I'm fairly new to MATLAB and haven't really used it for participant study before. I'm working on an audio perception experiment and wanted to build a GUI where participants would be able to listen to different audio clips and give them a specific rating. I would ideally like the responses by the participants to be anonymous.

In terms of execution, I'm planning on building a web app that could allow all participants to easily access the study. Could you suggest me some resources or guides I should take a look at for that?

9 Comments
2024/06/24
21:05 UTC

1

Matlab Code

We are currently trying to find code that would help with communication between MATLAB and the skin conductance software, BIOPAC or AcqKnowledge. More specifically, we want to create code on MATLAB that will let our computer handling BIOPAC know when certain events happen. That way, we can look back on data and pay attention to skin conductance levels when important events happened in the MATLAB-coded tasks. Is anyone aware of what this code might look like?

2 Comments
2024/06/24
18:51 UTC

0

How can i make this?

2 Comments
2024/06/24
17:10 UTC

1

Parameter optimisation with globalsearch (fmincon)

Hello everybody,

Im trying to optimise some parameters with the globalsearch fmincon algorithm. The equation is basicly x1**p1+x2**p2+... +x25*p25 = y

x and y are 31842x1 matrixes.

The overall goal is to eventually use less than 25 data sets to get the same y. The sum of the parameters is 1. Y was calculated with previously determined factors p, to proof the algorithm works i want the optimization to reach the same p's previously determined to show the algorithm works. With the constraint that the sum of the parameters is 1, the individual parameters differ greatly from my desired value. If i lift the constrain that the sum of parameters is 1, the parameters match, but the overall sum is 1.002. Ive tried to see what happens if i set sum p <1.001 but then the parameters still greatly differ.

This is my code:

num_parameters = 25;
A = ones(1,num_parameters);
b = 1; x0 = ones(num_parameters, 1)/num_parameters;
%Optimierung
rng default % For reproducibility
gs = GlobalSearch; Gleichung = @(p) rmse(zwva7 * p,BSW_Probandenreihe);
opts = optimoptions('fmincon','Display','iter','StepTolerance',1e-75,'MaxFunctionEvaluations',5e+06,MaxIterations=40e+02,ConstraintTolerance=0.001,Algorithm='active-set');
problem = createOptimProblem('fmincon','x0',x0,'objective',Gleichung,'lb',lb,'ub',ub,'Aeq',A,'beq',b,'options',opts);
[p_opt,fval] = run(gs,problem)

Any ideas how i can come closer to my desired outcome? The stop code is fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance.,

2 Comments
2024/06/24
16:31 UTC

1

How can I plot a constant value that divides at certain point?

It's an image from old University manual and I have to resolve a certain application based on it.

I have to write an Matlab script to plot this graph. It's purpose is to simulate the graph's behavior logically, not just drawing the lines through xline/yline.

As example I am aware that plotting this will show incorrect results:

D1 = 0.04;   D2 = 0.05;   D3 = 0.06; 

L_Al = 0.1095;
xNod = 0.0547;

x = 0 : 0.003: L_Al;
y1 = D1 * x;  y2 = D2 * x; y3 = D3 * x; 
plot(x, y1, x, y2, x, y3), axis([0 0.12 0.0 0.65])

I need to plot every value (D1, D2, D3) going constant from point ( X: 0 Y: 0.04/0.05/0.06 ) until they reach xNod value ( X: 0.0547 Y: 0.06) then /2 * 0.867 (it divides in half then it's being reduced to 86.7% of it's size) and then going constant until reaching L_Al = 0.1095 on X axis, so the result should be similar to what image above shows. Any suggestions on how it can be done?

Thank you in advance.

PS : I am new here, sorry if I am posting unrelevant question, I just really need help with that.

3 Comments
2024/06/24
12:14 UTC

2

Unable to Plot X axis according to required conditions

According to attached image from a research paper, I wish to plot the graph where the x axis is displayed as such. I have a graph where the x axis is in the linear scale and I am confused as to how to adjust the scale from interval to interval. Guidance would be appreciated.

https://preview.redd.it/3m0y5t74zh8d1.jpg?width=481&format=pjpg&auto=webp&s=0d57d1c53c594a6acdbcd660cb0db3ace169e58f

3 Comments
2024/06/24
10:32 UTC

2

Vector plot generation

I'm having trouble with the quiver function since there is a size problem. Issue is, all my four variables (x,y, u_at_t1 and v_at_t1) have the same size 381x101 and type for all is double. Still in getting an error saying the size of Y must match the size of U

1 Comment
2024/06/24
07:29 UTC

3

A* Code Review Request: function is slow

Just posted this on Code Reviewer down here (contains entirety of the function and more details):

performance - Working on Bi-Directional A* function in Matlab, want to speed it up - Code Review Stack Exchange

Currently it takes a significant amount of time (5+ minutes) to compute a path that includes 1000 nodes, as my environment gets more complex and more nodes are added to the environment the slower the function becomes. Since my last post asking a similar question, I have changed to a bi-directional approach, and changed to 2 MiniHeaps (1 for each direction). Wanted to see if anyone had any ideas on how to improve the speed of the function or if there were any glaring issues.

function [path, totalCost, totalDistance, totalTime, totalRE, nodeId] = AStarPathTD(nodes, adjacencyMatrix3D, heuristicMatrix, start, goal, Kd, Kt, Ke, cost_calc, buildingPositions, buildingSizes, r, smooth)
    % Find index of start and goal nodes
    [~, startIndex] = min(pdist2(nodes, start));
    [~, goalIndex] = min(pdist2(nodes, goal));

    if ~smooth
        connectedToStart = find(adjacencyMatrix3D(startIndex,:,1) < inf & adjacencyMatrix3D(startIndex,:,1) > 0); %getConnectedNodes(startIndex, nodes, adjacencyMatrix3D, r, buildingPositions, buildingSizes);
        connectedToEnd = find(adjacencyMatrix3D(goalIndex,:,1) < inf & adjacencyMatrix3D(goalIndex,:,1) > 0); %getConnectedNodes(goalIndex, nodes, adjacencyMatrix3D, r, buildingPositions, buildingSizes);
        if isempty(connectedToStart) || isempty(connectedToEnd)
            if isempty(connectedToEnd) && isempty(connectedToStart)
                nodeId = [startIndex; goalIndex];
            elseif isempty(connectedToEnd) && ~isempty(connectedToStart)
                nodeId = goalIndex;
            elseif isempty(connectedToStart) && ~isempty(connectedToEnd)
                nodeId = startIndex;
            end
            path = [];
            totalCost = [];
            totalDistance = [];
            totalTime = [];
            totalRE = [];
            return;
        end
    end

    % Bidirectional search setup
    openSetF = MinHeap(); % From start
    openSetB = MinHeap(); % From goal
    openSetF = insert(openSetF, startIndex, 0);
    openSetB = insert(openSetB, goalIndex, 0);

    numNodes = size(nodes, 1);
    LARGENUMBER = 10e10;
    gScoreF = LARGENUMBER * ones(numNodes, 1); % Future cost from start
    gScoreB = LARGENUMBER * ones(numNodes, 1); % Future cost from goal
    fScoreF = LARGENUMBER * ones(numNodes, 1); % Total cost from start
    fScoreB = LARGENUMBER * ones(numNodes, 1); % Total cost from goal
    gScoreF(startIndex) = 0;
    gScoreB(goalIndex) = 0;
    
    cameFromF = cell(numNodes, 1); % Path tracking from start
    cameFromB = cell(numNodes, 1); % Path tracking from goal

    % Early exit flag
    isPathFound = false;
    meetingPoint = -1;
    
    %pre pre computing costs
    heuristicCosts = arrayfun(@(row) calculateCost(heuristicMatrix(row,1), heuristicMatrix(row,2), heuristicMatrix(row,3), Kd, Kt, Ke, cost_calc), 1:size(heuristicMatrix,1));
    costMatrix = inf(numNodes, numNodes);
    for i = 1:numNodes
        for j = i +1: numNodes
            if adjacencyMatrix3D(i,j,1) < inf
                costMatrix(i,j) = calculateCost(adjacencyMatrix3D(i,j,1), adjacencyMatrix3D(i,j,2), adjacencyMatrix3D(i,j,3), Kd, Kt, Ke, cost_calc);
                costMatrix(j,i) = costMatrix(i,j);
            end
        end
    end
    costMatrix = sparse(costMatrix);

    %initial costs
    fScoreF(startIndex) = heuristicCosts(startIndex);
    fScoreB(goalIndex) = heuristicCosts(goalIndex);
    
    %KD Tree
    kdtree = KDTreeSearcher(nodes);

    % Main loop
    while ~isEmpty(openSetF) && ~isEmpty(openSetB)
        % Forward search
        [openSetF, currentF] = extractMin(openSetF);
        if isfinite(fScoreF(currentF)) && isfinite(fScoreB(currentF))
            if fScoreF(currentF) + fScoreB(currentF) < LARGENUMBER % Possible meeting point
                isPathFound = true;
                meetingPoint = currentF;
                break;
            end
        end
        % Process neighbors in parallel
        neighborsF = find(adjacencyMatrix3D(currentF, :, 1) < inf & adjacencyMatrix3D(currentF, :, 1) > 0);
        tentative_gScoresF = inf(1, numel(neighborsF));
        tentativeFScoreF = inf(1, numel(neighborsF));
        validNeighborsF = false(1, numel(neighborsF));
        gScoreFCurrent = gScoreF(currentF);
        parfor i = 1:numel(neighborsF)
            neighbor = neighborsF(i);
            tentative_gScoresF(i) = gScoreFCurrent +  costMatrix(currentF, neighbor);
            if  ~isinf(tentative_gScoresF(i))
                validNeighborsF(i) = true;   
                tentativeFScoreF(i) = tentative_gScoresF(i) +  heuristicCosts(neighbor);
            end
        end

        for i = find(validNeighborsF)
            neighbor = neighborsF(i);
            tentative_gScore = tentative_gScoresF(i);
            if tentative_gScore < gScoreF(neighbor)
                cameFromF{neighbor} = currentF;
                gScoreF(neighbor) = tentative_gScore;
                fScoreF(neighbor) = tentativeFScoreF(i);
                openSetF = insert(openSetF, neighbor, fScoreF(neighbor));
            end
        end
% Backward search

% Backward search
        [openSetB, currentB] = extractMin(openSetB);
        if isfinite(fScoreF(currentB)) && isfinite(fScoreB(currentB))
            if fScoreF(currentB) + fScoreB(currentB) < LARGENUMBER % Possible meeting point
                isPathFound = true;
                meetingPoint = currentB;
                break;
            end
        end
        % Process neighbors in parallel
        neighborsB = find(adjacencyMatrix3D(currentB, :, 1) < inf & adjacencyMatrix3D(currentB, :, 1) > 0);
        tentative_gScoresB = inf(1, numel(neighborsB));
        tentativeFScoreB = inf(1, numel(neighborsB));
        validNeighborsB = false(1, numel(neighborsB));
        gScoreBCurrent = gScoreB(currentB);
        parfor i = 1:numel(neighborsB)
            neighbor = neighborsB(i);
            tentative_gScoresB(i) = gScoreBCurrent + costMatrix(currentB, neighbor);
            if ~isinf(tentative_gScoresB(i))
                validNeighborsB(i) = true;
                tentativeFScoreB(i) = tentative_gScoresB(i) + heuristicCosts(neighbor)
            end
        end
        
        for i = find(validNeighborsB)
            neighbor = neighborsB(i);
            tentative_gScore = tentative_gScoresB(i);
            if tentative_gScore < gScoreB(neighbor)
                cameFromB{neighbor} = currentB;
                gScoreB(neighbor) = tentative_gScore;
                fScoreB(neighbor) = tentativeFScoreB(i);
                openSetB = insert(openSetB, neighbor, fScoreB(neighbor));
            end
        end

    end

    if isPathFound
        pathF = reconstructPath(cameFromF, meetingPoint, nodes);
        pathB = reconstructPath(cameFromB, meetingPoint, nodes);
        pathB = flipud(pathB);
        path = [pathF; pathB(2:end, :)]; % Concatenate paths
        totalCost = fScoreF(meetingPoint) + fScoreB(meetingPoint);
        
        pathIndices = knnsearch(kdtree, path, 'K', 1);
        totalDistance = 0;
        totalTime = 0;
        totalRE = 0;
        for i = 1:(numel(pathIndices) - 1)
            idx1 = pathIndices(i);
            idx2 = pathIndices(i+1);
            totalDistance = totalDistance + adjacencyMatrix3D(idx1, idx2, 1);
            totalTime = totalTime + adjacencyMatrix3D(idx1, idx2, 2);
            totalRE = totalRE + adjacencyMatrix3D(idx1, idx2, 3);

        end
       
        nodeId = [];
    else
        path = [];
        totalCost = [];
        totalDistance = [];
        totalTime = [];
        totalRE = [];
        nodeId = [currentF; currentB];
    end
end

function path = reconstructPath(cameFrom, current, nodes)
    path = current;
    while ~isempty(cameFrom{current})
        current = cameFrom{current};
        path = [current; path];
    end
    path = nodes(path, :);
end

function [cost] = calculateCost(RD,RT,RE, Kt,Kd,Ke,cost_calc)       
    % Time distance and energy cost equation constants can be modified on needs
            if cost_calc == 1
            cost = RD/Kd; % weighted cost function
            elseif cost_calc == 2
                cost = RT/Kt;
            elseif cost_calc == 3
                cost = RE/Ke;
            elseif cost_calc == 4
                cost = RD/Kd + RT/Kt;
            elseif cost_calc == 5
                cost = RD/Kd +  RE/Ke;
            elseif cost_calc == 6
                cost =  RT/Kt + RE/Ke;
            elseif cost_calc == 7
                cost = RD/Kd + RT/Kt + RE/Ke;
            elseif cost_calc == 8
                cost =  3*(RD/Kd) + 4*(RT/Kt) ;
            elseif cost_calc == 9
                cost =  3*(RD/Kd) + 5*(RE/Ke);
            elseif cost_calc == 10
                cost =  4*(RT/Kt) + 5*(RE/Ke);
            elseif cost_calc == 11
                cost =  3*(RD/Kd) + 4*(RT/Kt) + 5*(RE/Ke);
            elseif cost_calc == 12
                cost =  4*(RD/Kd) + 5*(RT/Kt) ;
            elseif cost_calc == 13
                cost =  4*(RD/Kd) + 3*(RE/Ke);
            elseif cost_calc == 14
                cost =  5*(RT/Kt) + 3*(RE/Ke);
            elseif cost_calc == 15
                cost =  4*(RD/Kd) + 5*(RT/Kt) + 3*(RE/Ke);
            elseif cost_calc == 16
                cost =  5*(RD/Kd) + 3*(RT/Kt) ;
            elseif cost_calc == 17
                cost =  5*(RD/Kd) + 4*(RE/Ke);
            elseif cost_calc == 18
                cost =  3*(RT/Kt) + 4*(RE/Ke);
            elseif cost_calc == 19
                cost =  5*(RD/Kd) + 3*(RT/Kt) + 4*(RE/Ke);
            end  
end

Update 1:  main bottleneck is the parfor loop for neighborsF and neighborsB, I have updated the code form the original post; for a basic I idea of how the code works is that the A* function is inside of a for loop to record the cost, distance, time, RE, and path of various cost function weights.
20 Comments
2024/06/24
00:35 UTC

1

matlab R2007a version

All,

Is it ok to install matlab R2007a version on the latest IOS on mac pro? I checked the system support and requirements page on matlab website but I didn't get the information I was looking for. It says Mac OS X 10.5 so I am assuming that this is not supported on the latest IOS. Does that mean this file is pretty much useless at this point since it's so old and outdated? Is there any way that I can use this version of matlab or should I just throw it away? I am not sure if this version is no longer supported since it's too old. Any recommendation and inputs are appreciated.

7 Comments
2024/06/23
20:19 UTC

1

Wind vector plot generation

I have the 30 day data for U component, V component and mean sea level pressure of a certain location with given lat and long. How do I plot the data of the horizontal and vertical component of wind with respect to msl pressure in vector form over an map of the location

1 Comment
2024/06/23
05:44 UTC

1

Simulink, Library doesn’t show in Library browser

Hello, I just installed PLC coder on simulink for Count up block, but the library doesn’t show in library browser. How I can fix it?

1 Comment
2024/06/22
23:36 UTC

2

Plotting data over pixels which have logical values = 1 and leaving the pixels with logical = 0 untouched

I am working on a code which displays the output of pressure sensitive paint over a 1024x1024 image. I need to plot the pressure at each pixel where the logical = 1 only. Any tips on how to do this efficiently? I also want to display this with imagesc() or maybe imshow(). Is there a way to make the pixels that have logical = 0 black instead of complying with the upper/lower limit of the colormap?

2 Comments
2024/06/22
14:38 UTC

Back To Top