Optimizing for a function output with multiple variables having different ranges.

Hello! I have a function that takes multiple inputs of temperature, pressure and so on to calculate net profit of a steam power plant. I am trying to create a new function that lets me run my first function multiple times with different ranges of my inputs and store or present the results of each iteration.

For example:

function [Net_Profit] = SteamPlant(T1, T2, P1, P2);

% This function would take two temperatures and pressures and output a Net Profit value%

Since it relies on an external function to calculate variables, I can't do "T1 = 35:45:10;" and then run the function. I need to find a way to run my function multiple times with singular input values across a range for each input.

Let me know if you need any clarification or details. I really appreciate the help and or suggestions! (:

Answer 1

[aara]

Hi, a for loop would be very useful, you could use it to cycle through all the datapoints you want:

%define T1, T2, P1,P2 Ranges:
T1 = 20:45;
T2 = 30:60;
%place values for P1 and P2
P1 = ...;
P2 = ...;
%you would also require a matrix to place all your results in
%the results would be equal to length(T1)*length(T2)*length(P1)*length(P2).
%define a counter for to store you results:
counter = 1;

for i = 1:length(T1)
    for j = 1:length(T2)
        for k = 1:length(P1)
            for m = 1:length(P2)
                NetProfit(counter,:) = SteamPlant(i,j,k,m);
                Input_data(counter,:) = [counter, i, j, k, m];
                counter = counter+1; 
            end
        end
    end
end

This would be very time consuming for large values. It will cycle through all values of P2 at constant P1,T1,T2 then change P1 by one and cycle again through P2. It continues until it has cycled through all values T1.

[Josef Woloschek]

When I input this type of code into my script, it says that I have a duplicate. Attached is my full code, the first two blocks are my attempts and making some sort of looping function and the answer you gave. Below that is my main function.

When I tried putting the code in its own script file and under a function, it says it doesn't recognize the original SteamPlant function I'm trying to call.

I'm really new to MATLAB so sorry if I'm making a really obvious mistake. I greatly appreciate the help!

% T1 = [40,45];
% boiler_T_range = [400,450,500,550,600];
% boiler_P_range = [5,6,7,8,9,10,11,12,13,14,15];
% turbine_P_range = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15];
% turbine1_eff = [.80,.85,.90,.95];
% turbine2_eff = [.80,.85,.90,.95];
% 
% for a = 1:2
%     condenser_T = T1(a);
%     [Net_profit(a)] = SteamPlant(T1,600,15,8,1,13,.85,.85);
% end
% Total(a,:) = [Net_profit];

%define T1, T2, P1,P2 Ranges:
T1 = 40:45;
T3 = 550:600;
%place values for P1 and P2
P3 = 5:15;
P4 = 1:15;
%you would also require a matrix to place all your results in
%the results would be equal to length(T1)*length(T2)*length(P1)*length(P2).
%define a counter for to store you results:
counter = 1;
for i = 1:length(T1)
    for j = 1:length(T3)
        for k = 1:length(P3)
            for m = 1:length(P4)
                Net_profit(counter,:) = SteamPlant(i,j,k,m);
                Input_data(counter,:) = [counter, i, j, k, m];
                counter = counter+1; 
            end
        end
    end
end

function [Table,Specific_Work,Net_profit] = SteamPlant(T1,T3,P3,P4,Reheat1,Reheat2,Turbine1,Turbine2)
%Function using input parameters to give cycle properties

state3 = steam_properties('P', P3, 'T', T3+273.15, 'SI');
state4ideal = steam_properties('P', P4, 's', state3.s, 'SI');
h4a = 0.15*state3.h+Turbine1*state4ideal.h;
state1 = steam_properties('T', T1+273.15, '', [ ], 'SI');
state4actual = steam_properties('P', P4, 'h', h4a, 'SI');
state2ideal = steam_properties('P', P3, 's', state1.s_l, 'SI');
h2a = ((state2ideal.h-state1.h_l)/.75)+state1.h_l;
state2actual = steam_properties('P', P3, 'h', h2a, 'SI');
state5 = steam_properties('P', P4, 'T', T3+273.15, 'SI');
state6ideal = steam_properties('T', T1+273.15, 's', state5.s, 'SI');
h6a = 0.15*state5.h+Turbine2*state6ideal.h;
state6actual = steam_properties('P', state1.P, 'h', h6a, 'SI');
state2b = steam_properties('P', state2ideal.P, '', 0, 'SI');

%Table Vectors
States = {'State 1'; 'State 2 Ideal'; 'State 2 Actual'; 'State 3'; 'State 4 Ideal'; 'State 4 Actual'; 'State 5'; 'State 6 Ideal'; 'State 6 Actual'};
Temperature = [state1.T; state2ideal.T; state2actual.T; state3.T; state4ideal.T; state4actual.T; state5.T; state6ideal.T; state6actual.T];
Pressure = [state1.P; state2ideal.P; state2actual.P; state3.P; state4ideal.P; state4actual.P; state5.P; state6ideal.P; state6actual.P];
Enthalpy = [state1.h_l; state2ideal.h; state2actual.h; state3.h; state4ideal.h; state4actual.h; state5.h; state6ideal.h; state6actual.h];
Entropy = [state1.s_l; state2ideal.s; state2actual.s; state3.s; state4ideal.s; state4actual.s; state5.s; state6ideal.s; state6actual.s];
Quality = [0; state2ideal.x; state2actual.x; state3.x; state4ideal.x; state4actual.x; state5.x; state6ideal.x; state6actual.x];

%Table
Table = table(States,Temperature,Pressure,Enthalpy,Entropy,Quality);

%Work Calculations
Net_work = 20000
Specific_Work_1to2 = h2a - state1.h_l
Specific_HT_2to3 = state3.h-h2a
Specific_Work_3to4 = state3.h - h4a
Specific_HT_4to5 = state5.h - h4a
Specific_Work_5to6 = state5.h - h6a
MassFlow = Net_work/(Specific_Work_3to4+Specific_Work_5to6-Specific_Work_1to2)
Total_Work_1to2 = MassFlow*Specific_Work_1to2
Total_HT_2to3 = MassFlow*Specific_HT_2to3
Total_Work_3to4 = MassFlow*Specific_Work_3to4
Total_HT_4to5 = MassFlow*Specific_HT_4to5
Total_Work_5to6 = MassFlow*Specific_Work_5to6

%T-s vapour dome
Ts = linspace (20+273.15, 373+273.15, 500);
for i = 1:1:length(Ts)
    satts = steam_properties('T', Ts(i), '', [], 'SI');
    x1t(i) = satts.s_l;
    x2t(i) = satts.s_v;
    yts(i) = Ts(i);
end
Ts_pointx = [x1t(end), x2t(end)];
Ts_pointy = [yts(end), yts(end)];

%P-v vapour dome
Pv = linspace (5e-3, 22, 10000);

for i = 1:1:length(Pv)
    satpv = steam_properties('P', Pv(i), '', [], 'SI');
    x1p(i) = satpv.v_l;
    x2p(i) = satpv.v_v;
    ypv(i) = Pv(i);
end
Pv_pointx = [x1p(end), x2p(end)];
Pv_pointy = [ypv(end), ypv(end)];

%Graphs
figure(1)
xT = [state1.s_l, state2ideal.s, state2b.s_l, state2b.s_v, state3.s, state4ideal.s, state5.s, state6ideal.s, state1.s_l];
yT = [state1.T, state2ideal.T, state2b.T, state2b.T, state3.T, state4ideal.T, state5.T, state6ideal.T, state1.T]
xTa = [state1.s_l, state2actual.s, state2b.s_l, state2b.s_v, state3.s, state4actual.s, state5.s, state6actual.s, state1.s_l];
yTa = [state1.T, state2actual.T, state2b.T, state2b.T, state3.T, state4actual.T, state5.T, state6actual.T, state1.T];
TsPlot = plot(xT,yT,'b--o')
hold on
plot(x1t,yts,'r')
hold on
plot(x2t,yts,'r')
hold on
plot(xTa,yTa,'m-.x')
hold on
plot(Ts_pointx,Ts_pointy,'r')
xlabel('Entropy (kJ/kg*k)')
ylabel('Temperature (k)')

figure(2)
xP = [state1.v_l, state2ideal.v, state2b.v_l, state2b.v_v, state3.v, state4ideal.v, state5.v, state6ideal.v, state1.v_l];
yP = [state1.P, state2ideal.P,state2b.P, state2b.P, state3.P, state4ideal.P, state5.P, state6ideal.P, state1.P];
xPa = [state1.v_l, state2actual.v, state2b.v_l, state2b.v_v, state3.v, state4actual.v, state5.v, state6actual.v, state1.v_l];
yPa = [state1.P, state2actual.P,state2b.P, state2b.P, state3.P, state4actual.P, state5.P, state6actual.P, state1.P];
PvPlot = plot(xP,yP,'b--o')
hold on 
plot(xPa,yPa,'m-.x')
hold on
plot(x1p,ypv,'r')
hold on
plot(x2p,ypv,'r')
hold on
plot(Pv_pointx,Pv_pointy,'r')
set(gca,'Xscale','log');
set(gca,'Yscale','log');
xlabel('Specific Volume (m^3/kg)')
ylabel('Pressure (MPa)')

%Cycle Efficieny
%%Ideal_efficieny = 1 - ((state3.h - state2ideal.h_l)+(state5)/(state3.h-state2ideal.h)
Actual_efficiency = ((state3.h-state2actual.h)+(state5.h-state4actual.h)-(state6actual.h-state1.h_l))/((state3.h-state2actual.h)+(state5.h-state4actual.h))

%Gas Cost
BTU = (Total_HT_2to3 + Total_HT_4to5)/1.06
mmBTU=BTU/(10^6)
NG_cost_year = mmBTU*3*86400*365

%Power Cost
MW = ((Total_Work_5to6 + Total_Work_3to4) - Total_Work_1to2)/1000
Power_revenue_year = MW * 25 * 24 * 365
Net_profit = Power_revenue_year - NG_cost_year

%Trade off curves
Trade_off = linspace (Reheat1, Reheat2, 500);
for i = 1:1:length(Trade_off)
    Reheat_prop = steam_properties('P', Trade_off(i), '', [], 'SI');
    %x1p(i) = satpv.s_l;
    %x2p(i) = satpv.s_v;
    %ypv(i) = Pv(i);
end

end

[Walter Roberson]

I predict you tried to store all of this in SteamPlant.m . When you define a function in a script then the function cannot have the same name as the script .