Back to index.
% 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