While loop that checks if difference is greater than zero not working

Hello everyone,

Here is the basic algorithm:

D is incremented.

L is incremented until the difference between variables x_cp - x_cg - D*SM = 0

Once x_cp - x_cg - D*SM = 0, the value lambda is compared to the previous value of lambda, if its larger it is set as the current lambda value.

D is incremented.

My problem is the while loop that checks the difference allows the difference to beome negative. I've tried a bunch of different approaches and none seem to work.

The full code is below, its hard to create a reproducible example just because of the problem's many nuances.

clear all
close all
clc

M_l = 5; %payload mass in [kg]
h_max = 3048; %max height in [meters]
a_max = 5; %normalized max acceleration
SM = 1; %Static margin for this test
rho_s = 1820; %density of structural material [kg/m^3]
rho_p = 1772; %density of propellant [kg/m^3]
sig_s = 8E7; %Shell working stress [Pa]
g = 9.81; %gravitational acceleration [m/s^2]
R = 287; %gas constant [J/kg-K]
T = 273.15; %standard temperature [K]
P = 100000; %standard pressure [pa]
gamma = 1.4; %coefficient ratio

R_opt = 1 + a_max; %Mass ratio: Pg. 182 typeset notes
W_eq = sqrt(h_max*g*(1/((log(R_opt)/2)*(log(R_opt)-2)+(R_opt-1)/R_opt))); %Equivalent V [m/s]: Pg. 183 typeset notes
t_b = ((R_opt-1)*W_eq)/(g*R_opt); %Burn time [s]: Pg. 183 typeset notes
a = sqrt(gamma*T*R); %Speed of sound [m/s]
M = W_eq/a; %Mach number
P_o = P*(1 + (gamma-1)/2*M^2)^(gamma/(gamma-1)); %Stagnation pressure [pa]: Pg. 

i = 0;
j = 0;
k = 0;

D = 0.1; %Starting D value [m]

condition = false; %While loop condition variable
lambda = 0; %Initial lambda value
lambdaf = 0.001;

while condition == false
    %Loop increment for D
    
    if lambdaf > lambda
        lambda = lambdaf;
    else
        break
    end
    i = i+1;
    D = D + 0.05;
    check = 1;
    L = 0.1; %Starting L value [m]
    L_p = 0;
    
    while check > 0
        
        j = j+1;
        %Loop increment for L
        L = L + 0.05;
        
        if L + D > L_p
            
            k = k+1;
            
            del = (D*P_o)/(2*sig_s); %Thickness of shell [m]
            
            M_fin = D^2/2*del*rho_s; %Mass of fin [kg]
            
            M_fus = pi*D*del*(L+D)*rho_s; %Mass of fuselage [kg], total length: L + D
            
            M_cone = (sqrt(D^2+(D/2)^2) + D)*pi*D*del*rho_s; %Mass of cone [kg]
            
            M_s = M_fus + M_cone + 4*M_fin; %Total mass of structure [kg]
            
            M_p = (R_opt - 1)*(M_s + M_l); %Mass of propellant [kg]
            
            L_p = M_p/(pi*D^2*(rho_p/4)); %Length of propellant section [m]
            
            %Ensures current increment is valid
            %The following values are derived from rocket dimensions based on
            %Centuri noes
            
            S = D;
            d = D;
            R = D/2;
            l = sqrt((D/2)^2+D^2);
            N = 4;
            a = D;
            b = 0;
            m = D;
            
            Cnan = 2; %nose coefficient
            x_nbar = (2/3)*D; %nose pressure moment arm for cone [m]
            
            Cnaf = (4*N*(S/d)^2)/(1+sqrt(1+((2*l)/(a+b))^2)); %Fin coefficient
            k_fb = 1 + R/(S + R); %Fin correction factor
            Cnafb = k_fb*Cnaf; %Final fin coefficient
            
            x_f = L+D; %Fin pressure moment arm [m]
            delx_f = (m*a)/(3*a) + (1/6)*a; %Calculation of delta X [m]
            x_fbar = x_f + delx_f; %Total pressure fin moment arm [m]
            
            Cna = Cnan + Cnafb; %Total pressure coefficient
            
            x_cp = (Cnan*x_nbar + Cnafb*x_fbar)/(Cna); %Center of pressure [m]
            
            L1 = (2/3)*D; %CG moment arm for nose cone and payload
            L2 = D+(L+D)/2; %CG moment arm for fuselage
            L3 = 5/3*D + L; %CG moment arm for fins
            L4 = 2*D+L-(L_p)/2; %CG moment arm for propellant
            
            %Center of gravity calculation
            x_cg = (M_l*L1 + M_cone*L1 + M_fus*L2 + N*M_fin*L3 + M_p*L4)/(M_l + M_fus + N*M_fin + M_p + M_cone);
            
            %Prints out check value for loop
            
            lambdaf = M_l/(M_p + M_s);
            
            check = x_cp - x_cg - D*SM
            
        else
            continue
        end
        
    end
    
end

m_dot =  M_p/t_b;
Thrust = m_dot*W_eq;

[Stephen23]

"My problem is the while loop that checks the difference allows the difference to beome negative."

The value will be negative in order for the WHILE loop to exit, so the final CHECK value will be negative. Are you expecting to get the penultimate CHECK value, the last (positive) one which did not fail that logical comparison ?

[raymond bryant]

Yes penultimate.

I've tried several different methods for it to be penultimate as well and I still run into similar problems.

Ideally, I'd have it check directly against zero. I really thought I could just round the difference but apparently it's not that simple.

[Stephen23]

You could store the penultimate values that you need, e.g.:

while val>0
    ..
    pen_val = val;
    val = .. new value;
end
pen_val % will be the last positive value

[raymond bryant]

I looked at this problem again and I'm going to try using fmincon. I think I have the necessary information for two nonlinear constraints.

Sometimes the while loop allows multiple negative values...

[Stephen23]

"Sometimes the while loop allows multiple negative values..."

No, the inner WHILE loop stops as soon as the first negative value occurs, just as it should.

However your code logic certainly allows multiple negative values:

while ..
    ..
    check = 1;
    ..
    while check > 0
        ..
    end
end

After getting one negative value you jump back to the outer WHILE loop, but in doing so you reset CHECK to a positive value. Thus the cycle continues...