12
\$\begingroup\$

I am working for an Aero-Design team and we make RC planes that are capable of humanitarian aid. We wished to make a tail dragger and for stability reasons the "Angle of Incidence" (AOI) of the plane has to be between 10 to 20 degrees when it is grounded. Now the AOI depends upon a number of factors such as LHT, LVT, height of landing-gears etc... Thus the number of combination are too large for a human to manually calculate and thus I wrote MATLAB code that is able to consider all the possible combination and permutation and display all the sets of combinations that make the AOI fall within the desired range.

Now the code is extremely long and is filled with too many for loops, this makes slow to process and difficult to interpret. So I wish to know whether or not I could optimize the code, shorten it and make it faster (maybe by pre-defining variables in vectors form)?

Here is the in-depth description of how I developed the logic for this code.

Here is a gist of the logic used:

2D sketch describing a tail-dragger

All the variables that influence the AOI have been marked in the sketch. By using simple trigonometry we can form 3 equations that relate some of these variables with AOI:

Equation 1 $$ \tan(\alpha)= \frac{x}{(V+q)} $$ Equation 2 $$ \tan(\alpha)= \frac{T}{y} $$ Equation 3 $$ \tan(\alpha)= \frac{U}{(L+H+q)} $$

On rearranging the above 3 equations in terms of alpha we get the following equations:

Equation 1 (constraining tail gear): $$ \sin(\alpha) = \frac{-2V\cdot c + \sqrt{4V^2c^2-4(V^2+x^2)(c^2-x^2)}} {2(V^2+x^2)} $$

Equation 2 (constraining main landing gear): $$ \tan(\alpha) = \frac{T}{y} $$

Equation 3 (constraining the whole system): $$ \sin(\alpha) = \frac{-(2Hc+2Lc) + \sqrt{(2Hc+2Lc)^2-4(L^2+H^2+2LH+U^2)(C^2- U^2)}}{2(L^2+H^2+2LH+U^2} $$ On solving these 3 equations simultaneously, we can find the AOI of the plane for given value of the variables.

AoA of a Tail-Dragger

clc   
clear all    
L=100;           %Distance between LEw to LEt   
H=20;            %Chord of H-stab  
B=15;            %Height of the fuselage   
w=(-5);          %Height of center of wing integration above the fuselage   
c=5;             %Clearance of the trailing edge of H-stab from the ground   
k=0;             %Distance betweem MLG from LEw (+ve if MLG is after LEw)  
for x=7:0.5:10;     %Height of TG  
    for y=60:0.5:130;       %Distance btw MLG and TG      
        for z=20:0.5:25;        %Height of MLG  
           for E=5:0.5:10;      %Height difference of LEw to LEt  
                T=z+B+w-E-x;    %Height of TG from the ground  
                U=z+B+w-E;      %Height of tail boom above the ground  
                V=L+H-y-k;      %Distance btw TG and trailing edge of H-stab  
                dis1=4*(V*c)^2-4*(V^2+x^2)*(c^2-x^2);  
                dis2=(2*H*c+2*L*c)^2-4*(L^2+H^2+2*L*H+U^2)*(c^2-U^2);  
                if (dis1>0 && dis2>0)  
                    alpha1=asin((-2*V*c+sqrt(dis1))/(2*(V^2+x^2)));%for eqn1  
                    alpha11=round(RD(alpha1));                     %for eqn1  
                    alpha2=atan(T/y);               %for eqn_2  
                    alpha22=round(RD(alpha2));      %for eqn_2  
                    alpha3=asin((-(2*H*c+2*L*c)+sqrt(dis2))...
                        /(2*(L^2+H^2+2*L*H+U^2)));       %for eqn_3  
                    alpha33=round(RD(alpha3));           %for eqn_3  
                    if(alpha11==alpha22 && alpha22==alpha33)  
                        alpha=alpha11;     
                    if(alpha>10 && alpha<20)    
                    disp(['The value of theta1 is:    ',num2str(alpha)]);   
                    disp(['The value of x is:    ',num2str(x)]);   
                    disp(['The value of y is:    ',num2str(y)]);     
                    disp(['The value of z is:    ',num2str(z)]);     
                    disp(['The value of E is:    ',num2str(E)]);     
                    disp('####################');   
                    end    
                end    
            end   
       end   
    end   
end   
end         
\$\endgroup\$
  • \$\begingroup\$ What is RD? Can you include it in the code? \$\endgroup\$ – Stewie Griffin Jan 13 '17 at 10:38
  • \$\begingroup\$ RD is a simple function that concerts radians into degree. I have included it in the link. \$\endgroup\$ – Neelay Doshi Jan 16 '17 at 8:29
3
\$\begingroup\$

In this answer I'll address only one specific improvement you can make: use vectorized operators instead of for loops. This gave me a speedup of about x70. Here is how:

L = 100;         % Distance between LEw to LEt   
H = 20;          % Chord of H-stab  
B = 15;          % Height of the fuselage   
w = -5;          % Height of center of wing integration above the fuselage   
c = 5;           % Clearance of the trailing edge of H-stab from the ground   
k = 0;           % Distance betweem MLG from LEw (+ve if MLG is after LEw)  
x = 7:0.5:10;    % Height of TG  
y = 60:0.5:130;  % Distance btw MLG and TG      
z = 20:0.5:25;   % Height of MLG  
E = 5:0.5:10;    % Height difference of LEw to LEt  

% span each of the vectors along a separate dimension
x_ = x(:);              %  size(x_) = [numel(x),1]
y_ = shiftdim(y(:),-1); %  size(y_) = [1,numel(y)]
z_ = shiftdim(z(:),-2); %  size(z_) = [1,1,numel(z)]
E_ = shiftdim(E(:),-3); %  size(E_) = [1,1,1,numel(E)]

U = bsxfun(@minus,z_,E_)+B+w;  % Height of tail boom above the ground. size(U) = [1,1,numel(z),numel(E)]
T = bsxfun(@minus,U,x_); % Height of TG from the ground.  size(T) = [numel(x),1,numel(z),numel(E)]
V = L+H-y_-k; % Distance btw TG and trailing edge of H-stab.   size(V) = [1,numel(y)]
dis1 = bsxfun(@minus,4*(V*c).^2,4*bsxfun(@times,bsxfun(@plus,V.^2,x_.^2),(c^2-x_.^2))); % size(dis1) = [numel(x),numel(y)]
dis2=(2*H*c+2*L*c)^2-4*(L^2+H^2+2*L*H+U.^2).*(c^2-U.^2);  % size(dis2) = [1,1,numel(z),numel(E)]
inds = bsxfun(@and,dis1>0,dis2>0); % logical indices of where dis1>0 and dis2>0.  size(inds) = [numel(x),numel(y),numel(z),numel(E)]
alpha1 = asin(bsxfun(@plus,-2*V*c,sqrt(dis1))./(2*bsxfun(@plus,V.^2,x_.^2))); % for eqn1.  size(alpha1) = [numel(x),numel(y)]
alpha11 = round(rad2deg(alpha1)); % for eqn1.  size(alpha11) = [numel(x),numel(y)]
alpha2 = atan(bsxfun(@rdivide,T,y_)); % for eqn_2. size(alpha2) = [numel(x),numel(y),numel(z),numel(E)]
alpha22 = round(rad2deg(alpha2)); %for eqn2. size(alpha2) = [numel(x),numel(y),numel(z),numel(E)]  
alpha3 = asin((-(2*H*c+2*L*c)+sqrt(dis2))./(2*(L^2+H^2+2*L*H+U.^2))); %for eqn_3. size(alpha3) = [1,1,numel(z),numel(E)]
alpha3 = round(rad2deg(alpha3)); %for eqn_3. size(alpha3) = [1,1,numel(z),numel(E)]  
alpha = repmat(alpha11,[1,1,numel(z),numel(E)]); % size(alpha) = [numel(x),numel(y),numel(z),numel(E)]
alpha(bsxfun(@ne,alpha11,alpha22) | bsxfun(@ne,alpha22,alpha3) | (inds==0)) = nan;  % this replaces the if statements
inds2 = find(~isnan(alpha) & (alpha>10) & (alpha<20)); % inds2 is a column vector
[inds_x,inds_y,inds_z,inds_E] = ind2sub(size(inds),inds2);

% ther results will be stored in the following matrix, where each row is a result, and the columns are alpha,x,y,z,E
results = [alpha(inds2),x(inds_x)',y(inds_y)',z(inds_z)',E(inds_E)'];

Here is a short explanation: The key idea is to span each of the vectors which you run over in the for loop, along a different direction. We span x along the first dimension, y along the second, z along the third, and E along the fourth dimension. This lets you use the bsxfun function to operate on these vectors giving you all the combinations you need. So for example, in the code above, T contains all combinations of the expression "z+B+w-E-x". Since it depends on the vectors x,z and E, it's dimensions will be [numel(x),1,numel(z),numel(E)] (where numel is the number of elements in the vector). The 1 is there since T does not depend on y.

The if statements in your code are replaced by vectorized conditioning, in the line the starts with alpha(bsxfun(. A drawback of this approach is that I calculated alpha for all combinations of the input, and then I replace with nan the values which do not meet the conditions of your if statements. This means that I waste time on unnecessary calculations. But as I said, I get a big speed improvement on your example so I think we're okay.

Another drawback of the vectorized approach is that it consumes more memory, since all the calculation results are stored at the same time. This may be somewhat mitigated by making sure to clear variables that are not needed any more. But in the case where there is not enough memory to even start the vectorized calculation, a hybrid approach can be taken, where you vectorize (say) x,y and z, but run over E with a for loop.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.