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

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');

Answer 1

[Prachi Kulkarni]

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');