%Group 17 Optimiser file
clc
clear
format compact
format short g
close all
%--------------------------------------
r_i=0.01; %Internal radius (m)
t=0.001; %Thickness (m)
r_o=r_i+t;
P=3; %Pressure (bars)


S_y=170; %Yield stress (MPa)
n_d=2; %Safety factor
stress_all=S_y/n_d %Allowable stress (MPa)
%--------------------------------------
%Constraints
C=1400*17/1000 % Capacity or Vi (m^3)
H=<3.5; %Height(m)
CF=10 %Cost factor from sections
%--------------------------------------
sol=[];
for i=1:length(section)
    r_i=section(i).r_i;
    for j=1:length(section(i).t
        t=section(i).t(j);
        for k=1:length(section(i).L
            L=section(i).L(k)
        if t/r_i<0.05 %Thin walled
            sigma_t=P*r_i/t;
            sigma_r=0;
            sigma_l=P*r_i/(2*t);
        else %Thick walled
            sigma_t=r_i^2*P/(r_o^2-r_i^2)*(1+r_o^2/r_i^2);
            sigma_r=r_i^2*P/(r_o^2-r_i^2)*(1-r_o^2/r_i^2);
            sigma_l=r_i^2*P/(r_o^2-r_i^2);
        end
        sigma_vM=sqrt(((sigma_t-sigma_r)^2+(sigma_t-sigma_l)^2+(sigma_l-sigma_r)^2)/2) %von Mises stress (MPa)
        if L<0 || H<0 || S_y<0 || r_i<0 || stress_all<0 || t<0 || P<0
            errorlg(' All inputs must be greater than zero')
        elseif sigma_vM<=stress_all %Check if von Mises stress is greater than allowable stress
            V_mat=2*pi*r_i*L+4*pi*r_i^2*t; %Volume of material (m^3)
            H=2*r_o %Height (m)
            L_T=L+H %Total length (m)
            v_cost=V_mat*CF %Vessel cost (£)
            sol=[sol;r_i,t,H,L_T,V_mat,v_cost]
        end
        end
    end
end


sol_sorted=sortrows(sol,4); %Sort feasible solutions based on vessel cost
op_sol=sol_sorted(1,:);
%--------------------------------------
%Display optimum values
disp(['optimum internal radius (m)=',num2str(sol_sorted(1,1))])
disp(['optimum thickness (m)=',num2str(sol_sorted(1,2))])
disp(['optimum height (m)=',num2str(sol_sorted(1,3))])
disp(['optimum total length (m)=',num2str(sol_sorted(1,4))])
disp(['optimum material volume (m^3)=',num2str(sol_sorted(1,5))])
disp(['optimum vessel cost (£)=',num2str(sol_sorted(1,6))])