Optimization of pump speed, impeller diameter and pipe diameter

I need assistance trying to find optimized values for pump speed, impeller diameter, and pipe diameter using a loop. The starting values for each are: 3000 RPM, 0.5m, 0.5m. Each value and be increased or decreased by 1 RPM or 1 mm respectively. Once the optimized values are found, i need them to be plugged into the script to find the minimum cost.

Here is the code:

clear all
clc

Q = input('Are you going to use English units? Y/N \n','s');
if Q == 'Y'
    Dp=input('Diameter of pipe in ft: ');
    Di=input('impeller diameter in ft: ');
elseif Q == 'N'
    Dp=input('Diameter of pipe in m: ') *3.28;
    Di=input('impeller diameter in m: ')*3.28;
else
    error('ya dun goofed, Beevers')
    
end
W=input('impeller rotational speed in rpm: ');
Wi=W*2*pi/60;

Temp= 50; %deg F
gamma_w = 62.4; %specific weight of the water from the aquifer
h=100*3.28; %height m to ft
l=5000*3.28; %length ft
h2=40*3.28; %second height ft
Vol=5000*(3.28^3); %cubic  volume ft^3
p=1.940; %slug/ft^3
mu=(27.4*10^-6); %lb*s/ft^2
L=5000*3.28; %ft
Zb=140*3.28; %ft
g=32.2; %ft/s^2
e_gi=0.0005; %ft

ne=2+17;%number of elbows
nv=4; %number of valves
Kl90elb=0.2;
Kl45bnd=0.2;
Klscv=2;
KlRee=.8;
Klext=1;
Kl=2*Kl90elb+15*Kl45bnd+4*Klscv+KlRee+Klext;

D_1 = 7.75/12; %impeller diameter in feet
omega_1 = 1750*2*pi/60; %rotational speed of pump one in rads/s
Q_1 = 0.002228.*[0;100;200;300;400;500;600;630]; %volume flow rate for pump one in ft^3/s
H_1 = [55;54;50.5;46;39;30;17;13]; %pump one head in feet
eta_1 = [NaN;NaN;.415;.48;.505;.48;.39;NaN]; %efficiency of pump one as a percentage
Pump_1 = [Q_1 H_1 eta_1]; %table for pump one
Power_1 = gamma_w.*Q_1.*H_1; %Pump power for pump one

D_2 = 7.75/12; %impeller diameter in feet
omega_2 = 2850*2*pi/60; %rotational speed of pump one in rads/s
Q_2 = 0.002228.*[0;100;200;300;400;500;600;680]; %volume flow rate for pump two in ft^3/s
H_2 = [151;149.5;155.5;140;130;119;102;87]; %pump two head in feet
eta_2 = [NaN;NaN;NaN;.42;.485;.525;.525;.48]; %efficiency of pump two as a percentage
Pump_2 = [Q_2 H_2 eta_2]; %table for pump two
Power_2 = gamma_w.*Q_2.*H_2; %Pump power for pump two

D_3 = 7.75/12; %impeller diameter in feet
omega_3 = 3500*2*pi/60; %rotational speed of pump one in rads/s
Q_3 = 0.002228.*[0;100;200;300;400;500;600;700]; %volume flow rate for pump three in ft^3/s
H_3 = [225;222;218;212;205;195;179;142]; %pump three head in feet
eta_3 = [NaN;NaN;NaN;.40;.455;.51;.53;.49]; %efficiency of pump three as a percentage
Pump_3 = [Q_3 H_3 eta_3]; %table for pump three
Power_3 = gamma_w.*Q_3.*H_3; %Pump power for pump three

D_4 = 6/12; %impeller diameter in feet
omega_4 = 3500*2*pi/60; %rotational speed of pump one in rads/s
Q_4 = 0.002228.* [0;100;200;300;400;500;580]; %volume flow rate for pump four in ft^3/s
H_4 = [113;110.5;105;95;81;68;39]; %pump four head in feet
eta_4 = [NaN;NaN;NaN;.47;.49;.40;.34]; %efficiency of pump four as a percentage
Pump_4 = [Q_4 H_4 eta_4]; %table for pump four
Power_4 = gamma_w.*Q_4.*H_4; %Pump power for pump four

D_5 = 7/12; %impeller diameter in feet
omega_5 = 3500*2*pi/60; %rotational speed of pump one in rads/s
Q_5 = 0.002228.*[0;100;200;300;400;500;600;640]; %volume flow rate for pump five in ft^3/s
H_5 = [170;165;160;151.5;141;128;111;102]; %pump five head in feet
eta_5 = [NaN;NaN;NaN;.405;.47;.525;.525;.52]; %efficiency of pump five as a percentage
Pump_5 = [Q_5 H_5 eta_5]; %table for pump five
Power_5 = gamma_w.*Q_5.*H_5; %Pump power for pump five

C_Q1 = Q_1.*(1/(omega_1*D_1^3)); %Flow coefficient for Pump 1
C_Q2 = Q_2.*(1/(omega_2*D_2^3)); %Flow coefficient for Pump 2
C_Q3 = Q_3.*(1/(omega_3*D_3^3)); %Flow coefficient for Pump 3
C_Q4 = Q_4.*(1/(omega_4*D_4^3)); %Flow coefficient for Pump 4
C_Q5 = Q_5.*(1/(omega_5*D_5^3)); %Flow coefficient for Pump 5

C_H1 = (32.2.*H_1)/(omega_1^2*D_1^2); %Head coefficient for Pump 1
C_H2 = (32.2.*H_2)/(omega_2^2*D_2^2); %Head coefficient for Pump 2
C_H3 = (32.2.*H_3)/(omega_3^2*D_3^2); %Head coefficient for Pump 3
C_H4 = (32.2.*H_4)/(omega_4^2*D_4^2); %Head coefficient for Pump 4
C_H5 = (32.2.*H_5)/(omega_5^2*D_5^2); %Head coefficient for Pump 5

C_P1 = Power_1/(1.94*omega_1^3*D_1^5); %Power coefficient for Pump 1
C_P2 = Power_2/(1.94*omega_2^3*D_2^5); %Power coefficient for Pump 2
C_P3 = Power_3/(1.94*omega_3^3*D_3^5); %Power coefficient for Pump 3
C_P4 = Power_4/(1.94*omega_4^3*D_4^5); %Power coefficient for Pump 4
C_P5 = Power_5/(1.94*omega_5^3*D_5^5); %Power coefficient for Pump 5

Q = [C_Q1' C_Q2' C_Q3' C_Q4' C_Q5']; %All flow data in a cell array
H = [C_H1' C_H2' C_H3' C_H4' C_H5']; %All head data in a cell array
P = [C_P1' C_P2' C_P3' C_P4' C_P5']; %All power data in a cell array
eta = [eta_1' eta_2' eta_3' eta_4' eta_5']; %All efficiency data in a cell array

%Plotting all datapoints in single graph
x1 = Q;
y1 = H;
y2 = eta;
ind = ~isnan(y2);
y2 = y2(ind);
x2 = x1;
x2 = x2(ind);

%now we do the curve fits:
%QHE
%we found them in matlab using the polyfit and printing the coefficients to
%find the curve fit lines, then plotted them with 100 intervals between 0
%and 0.3
QH = polyfit(x1,H,2);
format compact;
QE = polyfit(x2,y2,4);
format compact;
figure('Name','Fluid Flow of My Tears')
plot(x1,y1,'^');
title('Head Rise & Efficiency vs. Flow Coefficient')
hold on
yyaxis right
xlabel('Flow Coefficient');
ylabel('Efficiency')
plot(x2,y2,'x');
hold on
QHcap= zeros(1,100); %empty vec
Qspace= 0:0.0003:0.03; %100 points for line
QHline= polyval(QH,Qspace);
yyaxis left;
ylabel('Head Rise Coefficient')
plot1 = plot(Qspace,QHline,'b');
hold on
QEline= polyval(QE,Qspace);
yyaxis right;
plot2 = plot(Qspace,QEline,'r');
legend({'HeadC data','HeadC fit','Efficiency data','Efficiency fit'},'Location','south')

%now we do fig 2 with Cp instead of eta

%now we do the curve fits:
%QHP
%we found them in matlab using the polyfit and polyval the coefficients to
%find the curve fit lines, then plotted them with 100 intervals between 0
%and 0.3
QH2 = polyfit(x1,H,2);
QHline= polyval(QH2,Qspace);
QP = polyfit(x1,P,2);
format compact;
figure(2);
plot(x1,y1,'^');
title('Head Rise & Power Coefficient vs. Flow Coefficient');
xlabel('Flow Coefficient')
hold on;
yyaxis right;
ylabel('Power Coefficient')
plot(x1,P,'o');
Qspace= 0:0.0003:0.03; %100 points for line
yyaxis left;
ylabel('Head Rise Coefficient')
plot3 = plot(Qspace,QHline,'b');
hold on;
QPline= polyval(QP,Qspace);
yyaxis right;
plot4 = plot(Qspace,QPline,'r');
legend({'HeadC data','HeadC fit','PowerC data','PowerC fit'},'Location','south')

%uses discriminant and quad formula to find cq intersection
disco=QP-QH2;
dis=disco(2)^2-4*disco(1)*disco(3);
if dis>0
    x_coord1=(-(disco(2)^2)+(disco(2)^2-4*disco(1)*disco(3))^(1/2))/(2*disco(1));
    x_coord2=(-(disco(2)^2)-(disco(2)^2-4*disco(1)*disco(3))^(1/2))/(2*disco(1));
       %if function below finds real root closest to 0
    if x_coord1 > x_coord2 && x_coord1 > 0
        cq=x_coord1;
    elseif x_coord2 > x_coord1 && x_coord2 > 0
        cq=x_coord2;
    else
        fprintf('Something''s wrong')
    end
else
    fprintf('Something''s wrong')
end
%This will relate the Reynolds Number and relative roughness to
%a corresponding darcy friction factor
Re=(p * cq * Wi * Di^3)/(pi/4) /Dp/mu;
%^ Cq found found by the intersection point between two curves
%in 2nd plot
%fprintf('Reynold''s number is: %d\n',Re);
rel_rough=e_gi/Dp;
df=1;

if Re < 2300
    f = 64/Re;
    
%elseif (2300<Re) && (Re<4000)
    %disp('flow is transitional and cannot be predicted');
    
else Re > 4000;
    i=0;
    f=0.023; %we'll use this as a first guess bc it's like in the middle of the moody graph
    while(abs(df)>10^-6) && (i<100)
    f1=f;
    f = ((-2.0*log10(((rel_rough/(3.7))+(2.51/(Re*f^(1/2))))))^(-1))^2;
    %fprintf('friction fact is: %d\n', f); %just to check to see if loops are working
    df=f1-f;
    i=i+1;
    end
end

v=(mu*Re)/(p*Dp); 
h= (v^2/(2*g))*(1 + f*(L/Dp)+Kl)+Zb;
Qop= v*pi*(Dp/2)^2;
fprintf('The operating volume flow rate is: %d ft^3/s \n',Qop);
eop=2.635*10^6*(cq^4)-1.734*10^5*(cq^3)+2721*(cq^2)+4.328*(cq)+0.2851;%efficiency at optimal point
fprintf('The operating efficiency is: %d\n',eop)
Power_p=(h*Qop*g*p/550*745.7/1000); %pump power in kwatts
t=(Vol/Qop)/3600; %time to reach daily volume needed to supply town
fprintf('This is how long the pump will run to fill the town water tower: %d hours \n',t)

%let's go shopping:
pmc=3500+600*Di*30.48; %pump and motor cost
vlv=nv*(300+80*Dp*30.48); %valve cost
elbow=ne*(50+20*Dp*30.48); %elbow cost
pipc=(l/3.28)*(1.50*Dp*30.48);%pipe cost

elec=0.115*Power_p*t; %electricity cost per day
elec_total = elec*365*5; %cost for 5 years
fprintf('total cost of electricity over 5 years is: $ %.2f\n',elec_total);

total_cost=elec_total + pmc +vlv + elbow + pipc;

capital_cost=pmc+vlv+elbow+pipc;
%results:
fprintf('Flow coefficient for the system is: %d\n',cq)
ch=-94.3*(cq)^2-0.5624*(cq)+0.1251;
fprintf('Head rise coefficient for the system is: %d\n',ch)
fprintf('Pump head is: %d ft \n',h)
fprintf('Total 5 year power consumption is: %d kWh \n',Power_p*t*5*365)
fprintf('Capital cost of system with input parameters is: $ %.2f\n',capital_cost)
fprintf('Annual operating cost for a pump with the input parameters is:$ %.2f\n',(total_cost-capital_cost)/5)
fprintf('Total cost for a pump with the input parameters is: $ %.2f\n',total_cost)
%%
RPM_range = {1800:1:3600};
Dp_range = {0:0.01:1};
Di_range = {0:0.01:1};

Answer 1

[Kholood]

Optimization linear programming MATLAB for pipe line of water