I have two variables that I need to loop over calculations to find the best combination

2 views (last 30 days)
Basically, I have written a code that does some calculations using two variables. It must then iterate through these variables to find the best combination of these variables. I would then like to plot these values in a surface plot. I will add the full code that I have have. I am terrible at coding in general and have tried my best to get this to work.
This might be a long problem with a lot of issues, so I apologise in advance. Even small tips or fixes will be of great help. Here is my code, I have added comments to show my general idea to approach the problem.
%Constants
p1= 97000; %Pressure 1
t1=285.15; %Temp 1
t5=1493.15; %Temp 5
y=1.4; %Gamma
cp=1005;
n=0.9; %Small Stage Efficiency
W_out= (3.6*(10^7))/0.96; %Work out
a=50; %Incrrement Size
optimal_cost=1000000000000000000; %Very large value for money
best_p_ratio=0;
best_p_int_ratio=0;
optimal_efficiency=0;
p_ratio_array= linspace(1.5,10,a); %Creation of array for Pressure Ratio
p_ratio_int_array= linspace(0.1,0.9,a)'; %Creation of array for Intercooler Pressure ratio as a percent
%Intercooler pressure ratio is set as a percentage as the intercooler ratio
%should not be higher or equal to the Pressure Ratio
t_cost= linspace(0,0,a); %Array creation for Total Cost
n_th_count= linspace(0,0,a); %array creation for efficiency counting
for kk=1:a
p_ratio=p_ratio_array(kk);
for kk=1:a
p_ratio_int=p_ratio_int_array(kk);
%Compressor
t2=t1*((p_ratio*p_ratio_int)^((y-1)/(n*y)));
w1_2=cp*(t1-t2);
t3=t1;
q2_3=cp*(t3-t2);
t4=t3*((p_ratio/(p_ratio*p_ratio_int))^((y-1)/(n*y)));
w3_4=cp*(t3-t4);
%Combustion and Turbine
q4_5=cp*(t5-t4);
t6=t5*((1/p_ratio)^((n*(y-1))/y));
w5_6=cp*(t5-t6);
%Efficiency and Mass Flow Rate
n_th_count(kk)=((w5_6+(w1_2+w3_4))./q4_5).*100;
n_th=n_th_count(kk);
m_dot= W_out./ (w5_6+(w1_2+w3_4));
%Costs
t_cost(kk)= ( 15000.*((m_dot).^0.85).*((p_ratio).^2.5) )+( (25000.*m_dot).*25 );
total_cost=t_cost(kk);
if total_cost<0 %I was getting negative values for total cost, so I added this statement to stop that from occuring
break
end
if total_cost<optimal_cost %My if statement to find the best combination by finding the cheapest price
best_p_ratio=p_ratio;
best_p_int_ratio=p_ratio_int;
optimal_cost=total_cost;
optimal_efficiency=n_th;
end
end
end
hold on
%Plots
figure(1)
AX=plot(p_ratio_array,n_th_count,'bo-', ...
'LineWidth', 1, 'MarkerSize', 5, ...;
'MarkerFaceColor', 'r', 'MarkerEdgeColor', [.5, .9, .7]);
xlabel('Pressure Ratio');
ylabel('Efficiency (%)');
title('Pressure Ratio Vs Efficiency')
figure(2)
plot(p_ratio_int_array,n_th_count,'bo-', ...
'LineWidth', 1, 'MarkerSize', 5, ...;
'MarkerFaceColor', 'r', 'MarkerEdgeColor', [.5, .9, .7]);
xlabel('Pressure Ratio Intercooler');
ylabel('Efficiency (%)');
title('Pressure Ratio Intercooler Vs Efficiency')
figure(3)
plot(n_th_count,t_cost,'bo-', ...
'LineWidth', 1, 'MarkerSize', 5, ...;
'MarkerFaceColor', 'r', 'MarkerEdgeColor', [.5, .9, .7]);
xlabel('Efficiency (%)');
ylabel('Total Cost');
title('Efficiency Vs Total Cost');
figure(4)
surf(p_ratio_array,t_cost,p_ratio_int_array)
xlabel('Pressure Ratio');
ylabel('Total Cost');
zlabel('Intercooler Pressure Ratio');

Accepted Answer

Prachi Kulkarni
Prachi Kulkarni on 20 Oct 2021
Hi,
In the code you have provided, the same variable kk is being used to loop over both the variable arrays p_ratio_array and p_ratio_int_array. Further, to iterate through all combinations of the 2 input variables, t_cost and n_th_count need to be 2D matrices.
Here is a modified version of your code which will give the surface plot. Comments have been added only where changes are made. You may need to make some changes to get the other 3 line plots as per your requirement.
p1= 97000;
t1=285.15;
t5=1493.15;
y=1.4;
cp=1005;
n=0.9;
W_out= (3.6*(10^7))/0.96;
a=50;
optimal_cost=1000000000000000000;
best_p_ratio=0;
best_p_int_ratio=0;
optimal_efficiency=0;
p_ratio_array= linspace(1.5,10,a);
p_ratio_int_array= linspace(0.1,0.9,a)';
t_cost= zeros(a,a); %2D Array creation for Total Cost
n_th_count= zeros(a,a); %2D array creation for efficiency counting
for kk1=1:a % Using variable kk1 instead of kk
p_ratio=p_ratio_array(kk1); % Using variable kk1 instead of kk
for kk2=1:a % Using variable kk2 instead of kk
p_ratio_int=p_ratio_int_array(kk2); % Using variable kk2 instead of kk
t2=t1*((p_ratio*p_ratio_int)^((y-1)/(n*y)));
w1_2=cp*(t1-t2);
t3=t1;
q2_3=cp*(t3-t2);
t4=t3*((p_ratio/(p_ratio*p_ratio_int))^((y-1)/(n*y)));
w3_4=cp*(t3-t4);
q4_5=cp*(t5-t4);
t6=t5*((1/p_ratio)^((n*(y-1))/y));
w5_6=cp*(t5-t6);
n_th_count(kk1,kk2)=((w5_6+(w1_2+w3_4))./q4_5).*100; % Using 2D indexing (kk1,kk2) instead of kk
n_th=n_th_count(kk1,kk2); % Using 2D indexing (kk1,kk2) instead of kk
m_dot= W_out./ (w5_6+(w1_2+w3_4));
t_cost(kk1,kk2)= ( 15000.*((m_dot).^0.85).*((p_ratio).^2.5) )+( (25000.*m_dot).*25 ); % Using 2D indexing (kk1,kk2) instead of kk
total_cost=t_cost(kk1,kk2); % Using 2D indexing (kk1,kk2) instead of kk
if total_cost<0
break
end
iftotal_cost<optimal_cost
best_p_ratio=p_ratio;
best_p_int_ratio=p_ratio_int;
optimal_cost=total_cost;
optimal_efficiency=n_th;
end
end
end
hold on
figure(1)
AX=plot(p_ratio_array,n_th_count,'bo-', ...
'LineWidth', 1, 'MarkerSize', 5, ...;
'MarkerFaceColor', 'r', 'MarkerEdgeColor', [.5, .9, .7]);
xlabel('Pressure Ratio');
ylabel('Efficiency (%)');
title('Pressure Ratio Vs Efficiency')
figure(2)
plot(p_ratio_int_array,n_th_count,'bo-', ...
'LineWidth', 1, 'MarkerSize', 5, ...;
'MarkerFaceColor', 'r', 'MarkerEdgeColor', [.5, .9, .7]);
xlabel('Pressure Ratio Intercooler');
ylabel('Efficiency (%)');
title('Pressure Ratio Intercooler Vs Efficiency')
figure(3)
plot(n_th_count,t_cost,'bo-', ...
'LineWidth', 1, 'MarkerSize', 5, ...;
'MarkerFaceColor', 'r', 'MarkerEdgeColor', [.5, .9, .7]);
xlabel('Efficiency (%)');
ylabel('Total Cost');
title('Efficiency Vs Total Cost');
figure(4)
[X,Y]=meshgrid(p_ratio_array,p_ratio_int_array); % Using meshgrid to get X,Y coordinates for surface plots
surf(X,Y,abs(t_cost)) % Surface plot of magnitude of t_cost (magnitude taken because t_cost is of type complex double)
xlabel('Pressure Ratio');
ylabel('Total Cost');
zlabel('Intercooler Pressure Ratio');

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by