Back to [[kyung_ha_thesis_code|index]].

<file matlab Figs3_4.m>
% Optimize the mass distribution for a quadruple pendulum for longitudinal
% seismic isolation. The script optimizes equations equation 6 in 
% T1300786-v7 based on equations 2 and 3. Note, this does not take into 
% account other important factors such as bounce and roll modes; and 
% dilution factors; and wire lengths. 
%
% This script optimizes the mass distribution in two different ways.
% Agreement between the two enhances the confidence that they are correct.
% The first is a brute force approach that simply scans many values of the
% mass distribution and picks the best result. The second iterates between
% equations 2 and 3. 
% 
% The script generates Figures 3 and 4 in T1300786-v7.
%
% Brett Shapiro
% Updated 8 Dec 2014

%% Given parameters
list_constants;

totalmass = total_mass; % the assumed permitted total suspended mass (I recently heard, Dec 2014, from Brian Lantz that we may be able to push this up based on a better understanding of the ISI pyaload limit, which would improve overal performance)
testmass = M4;  % assumed test mass mass


% totalmass = 248; % the assumed permitted total suspended mass (I recently heard, Dec 2014, from Brian Lantz that we may be able to push this up based on a better understanding of the ISI pyaload limit, which would improve overal performance)
% testmass = 80;  % assumed test mass mass

%% Brute force optimization
% This section scans many mass values of the UIM and PUM and simply picks
% the best result. The top mass value is calculated from the remaining
% permitted mass in the sus.

pointnum = 1000; % number of mass values to try for each stage
uimvec = logspace(log10(30),log10(50),pointnum); % uim mass values to try
pumvec = logspace(log10(30),log10(50),pointnum); % pum mass values to try

isolationfactor = NaN(pointnum,pointnum); % initialize the isolation factor before entering the loop. This is the figure of merit. Smaller numbers are better.

disp(' ')
disp('Looping over UIM and PUM mass values')

for uimind = 1:pointnum % looping over UIM mass values
    
    % display progress on command line for slow loops
%     if mod(uimind,pointnum/100) == 0
%         
% %         disp([num2str(100*uimind/pointnum),'%'])
%         
%     end
    
    for pumind = 1:pointnum     % loop over PUM mass values
        
        topmass = totalmass - uimvec(uimind) - pumvec(pumind) - testmass;
        
        if topmass < 0
            
            continue
            
        else
            
            isolationfactor(uimind,pumind) = (uimvec(uimind) + pumvec(pumind) + testmass) * (pumvec(pumind) + testmass) / ( topmass*uimvec(uimind)*pumvec(pumind) ); % the mass term of equation 6 in T1300786-v7
            
        end
        
    end
    
end

% find minimum in each column, i.e. minimum isolation for a given pum mass
[min1,min1ind] = min(isolationfactor);

% find minimum isolation overall
[min2,min2ind] = min(min1);

% Best mass values
uimmass = uimvec(min1ind(min2ind));
pummass = pumvec(min2ind);
topmass = totalmass - uimmass - pummass - testmass; % top mass is just the remainder from the total permitted mass

% display results to the command line
disp(' ')
disp(['Best TOP is ',num2str(topmass),' kg'])
disp(['Best UIM is ',num2str(uimmass),' kg'])
disp(['Best PUM is ',num2str(pummass),' kg'])


%% This is a much more efficient optimization. It iterates between the solutions m2 and m3 in equations 2 and 3 of T1300786-v7.
% Equation 2 gives the best value of the UIM given the other mass values.
% Equation 3 gives the best value of the PUM given the other mass values.
% Iterating between the two should give the best value overall.

% Initializations before loop
iter = 10; % number of loop iterations
pumitervec = zeros(iter,1); % PUM mass values
uimitervec = zeros(iter,1); % UIM mass values
topitervec = zeros(iter,1); % top mass mass values
pumitervec(1) = 80; % pick/guess a PUM mass value to start the loop
isolationiterfactor = zeros(iter,1); % the isolation factor, this is the figure of merit, smaller is better

% aLIGO10HzLisolation = 52.267252286217634; % (m1+m2+m3+m4)*(m2+m3+m4)*(m3+m4)/(m1*m2*m3), production fiber model

% iterate between UIM and PUM optimal equations to converge the optimal
for ii = 1:iter-1;
    
    % optimal UIM mass given the PUM mass: Eq 2
    uimitervec(ii) = -(pumitervec(ii) + testmass) + sqrt(totalmass*(pumitervec(ii) + testmass));
    
    topitervec(ii) = totalmass - testmass - uimitervec(ii) - pumitervec(ii);
    
    isolationiterfactor(ii) = (uimitervec(ii) + pumitervec(ii) + testmass) * (pumitervec(ii) + testmass) / ( topitervec(ii)*uimitervec(ii)*pumitervec(ii) );
    
    % Optimal PUM mass given the UIM mass: Eq 3
    A = testmass*(uimitervec(ii) + testmass)/(totalmass + testmass);
    pumitervec(ii+1) = -A + sqrt( A^2 + A*(totalmass-uimitervec(ii)-testmass) ); 
    
end

 % final iteration for some terms
uimitervec(iter) = -(pumitervec(iter) + testmass) + sqrt(totalmass*(pumitervec(iter) + testmass));
topitervec(iter) = totalmass - testmass - uimitervec(iter) - pumitervec(iter);
isolationiterfactor(iter) = (uimitervec(iter) + pumitervec(iter) + testmass) * (pumitervec(iter) + testmass) / ( topitervec(iter)*uimitervec(iter)*pumitervec(iter) );
g = 9.81; % gravity
freqval = 10; % frequency at which to evaluate isolation
%L4 = 2.14/4; L1=(2.14-L4)/3; L2=L1; L3=L1; % lengths of wires for isolation calculation
totalL=total_length;
L4=F_l; 
L1=(totalL-L4)/3; 
L2=L1; 
L3=L1;
isolationiterfactor = ((g^4)/((2*pi*freqval)^8))*(1/(L1*L2*L3*L4))*totalmass*isolationiterfactor; % isolation factor, equation 1

% plotting evolution of mass values
figure,plot(1:iter,topitervec,'-*',1:iter,uimitervec,'-*',1:iter,pumitervec,'-*','LineWidth',5,'MarkerSize',10),grid on  %maxfig
set(gca,'FontSize',30)
%title(['m_4 = ',num2str(testmass),'kg, Final m_1 = ',num2str(round(100*topitervec(iter))/100),' kg, Final m_2 = ',num2str(round(100*uimitervec(iter))/100),' kg, final m_3 = ',num2str(round(100*pumitervec(iter))/100),' kg'])
title('m_4 = 160 kg')
ylabel('Mass (kg)')
xlabel('Loop iterations')
%axis(1,10,35,80)
legend('m_1 (TOP)','m_2 (UIM)','m_3 (PUM)')
%FillPage('w')
idplot

% plotting evolution of isolation factor
figure,plot(1:iter,isolationiterfactor,'-*','LineWidth',5,'MarkerSize',10), grid on  %maxfig
set(gca,'FontSize',30)
xlabel('Loop iterations')
title(['L4 = ',num2str(L4),'m, Minimizing C_1 by iterating between optimal m_2 and m_3 values'])
ylabel('C_1 (m/m)')
%FillPage('w')
idplot
</file>