I keep getting the solve stops prematurely for fsolve?
    4 visualizzazioni (ultimi 30 giorni)
  
       Mostra commenti meno recenti
    
Why do I keep getting this error? I think the equations are solvable though? 
%% ROC ESTIMATION FOR PS = 3 WITH DIFFERENT THRUST VECTORING SETTINGS
% from sea level
clc
clear all
% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;
csfc = 200/60/60/1000;
%% ------ Functions to solve for FSOLVE (steady climbing flight) ------ %%
% func1 = mass' = (csfc*T)
% func2 = gamma' = 1/mu * (T*sin(epsilon)+L-mgcos(gamma)) = 0
% func3 = Cm = 0 
%% ------ FSOLVE variables ------ %%
% x(1) = velocity
% x(2) = deltae
% x(3) = gamma
%% ------ Altitude range from flight envelope ------ %%
% Using the lowest max alt from level flight envelope so easier for comparison
hlist_3 = linspace(0,14000,300);
%% ------ Array for storing ------ %%
% 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. altitude   7. thrust   8.
% 1/ROC   9.c*T/roc  
% ORIGINAL NO THRUST VECTORING
final_list_3_ORIGINAL = zeros(length(hlist_3),9);
% PS = 2.5 WITH THRUST VECTORING
final_list_3_thrustvec_5 = zeros(length(hlist_3),9);
final_list_3_thrustvec_10 = zeros(length(hlist_3),9);
final_list_3_thrustvec_15 = zeros(length(hlist_3),9);
final_list_3_thrustvec_20 = zeros(length(hlist_3),9);
%% PS = 3 NO THRUST VEC ORIGINAL
for j = 1:length(hlist_3)
    PS = 3;
    h = hlist_3(j);
    [rho, speedsound] = atmos(h);
    epsilon = 0;    
    [aoa_start,aoa_end] = aoa_guess(PS,epsilon);
    aoalist = linspace(aoa_start,aoa_end,150);
    store_v_de_gam_aoa_roc = zeros(length(aoalist),6); 
    % 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. thrust
    x = [40 0 0];
    for i = 1:length(aoalist)
        T = @(x) PS *( (7+x(1)/speedsound)*200/3 + h/1000*(2*x(1)/speedsound -11));
        qbar = @(x) 0.5*rho*x(1)^2*S;
        CL = @(x) 6.44*aoalist(i) + 0.355*x(2);
        L = @(x) CL(x)*qbar(x);
        Cm = @(x) 0.05 - 0.683*aoalist(i) - 0.923*x(2);
        f1 = @(x) csfc*T(x);
        f2 = @(x) (1/(m*x(1)))*(T(x)*sin(epsilon)+L(x)-m*g*cos(x(3)));
        f3 = @(x) Cm(x);
        func = @(x) [f1(x);f2(x);f3(x)];
        x = fsolve(@(x) func(x), x);
    end   
end
%% ------ Atmos function ------ %%
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos(h)
h1 = 11000;    % Height of tropopause
h2 = 20000;    % End height of table
g = 9.81;
R = 287; 
c = 6.51e-3;       % temperature lapse dt/dh = - c = -6.51 degcelcius/km
T0 = 15+273.15;       % Temperature sea level
p0 = 101325;       % pressure sealevel
rho0 = 101325/R/T0;      % density sealevel = pressure / R*T, R=287, T = 15 degcelcius       
T1 = T0 - c*h1;    % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506;   % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506;   % Density at 11km
T2 = T1;    % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1));  % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1));  % Density at 20km
if h <= h1
%       disp('Troposphere');
    T = T0 - c*h;
    p = p0 * (T/T0)^5.2506;
    rho = rho0 * (T/T0)^4.2506;
    speedsound = (1.4*R*T)^0.5;
elseif h <= h2
%       disp('Tropopause');
    T = T1;
    p = p1 * exp(-g/(R*T)*(h-h1));
    rho = rho1 * exp(-g/(R*T)*(h-h1));
    speedsound = (1.4*R*T)^0.5;
end
return 
end
%% ------ AOA guess function ------ %%
% To guess the range of AOA at each altitude
function [aoa_start,aoa_end] = aoa_guess(PS,epsilon)
h = 0;
W = 1248.5*9.81;        % Weight of UAV
S = 17.1;       % Wing area        
%     CL = (6.44*aoa + 0.355*eledefl);
%     CD = (0.03 + 0.05*(6.44*aoa + 0.355*eledefl)^2);
%     Cm = (0.05 - 0.683*aoa - 0.923*eledefl);
%     thrust = (PS * ( (7 + V/speedsound )*200/3  + h*(2*(V/speedsound) - 11) ));
%%
% V = x(1);
% aoa = x(2);
% eledefl = x(3);
[rho, speedsound] = atmos(h);
qbar = @(x) 0.5*rho*x(1)^2*S;
CL = @(x) (6.44*x(2) + 0.355*x(3));
CD = @(x) (0.03 + 0.05*(CL(x))^2);
Cm = @(x) 0.05 - 0.683*x(2) - 0.923*x(3);
thrust = @(x) ( PS * ( 200/3*( 7 + x(1)/speedsound ) + h/1000*( 2*x(1)/speedsound - 11) ) );
func_1 = @(x) qbar(x)*CL(x) + thrust(x)*sin(x(2)+epsilon) - W;      % Vertical
func_2 = @(x) qbar(x)*CD(x) - thrust(x)*cos(x(2)+epsilon);      % Horizontal
func_3 = @(x) Cm(x);
FUNC = @(x) [func_1(x); func_2(x); func_3(x)];
% ------ Guess --> Range of AOA ------ %
Xlow = fsolve(@(x) FUNC(x), [10 0 -30/180*pi]); % high AOA
aoa_end = Xlow(1,2);
Xhigh = fsolve(@(x) FUNC(x), [100 25/180*pi -30/180*pi]); % low AOA
aoa_start = Xhigh(1,2);
end
2 Commenti
  Torsten
      
      
 il 17 Mar 2022
				I can't reproduce the error. Under Ocatve's "fsolve", your code finishes without error message.
Risposte (1)
  Star Strider
      
      
 il 17 Mar 2022
        Possibly — 
% ------ Guess --> Range of AOA ------ %
opts = optimoptions('fsolve', 'MaxIterations', 5E+3);                           % <— ADD THIS
Xlow = fsolve(@(x) FUNC(x), [10 0 -30/180*pi], opts); % high AOA
aoa_end = Xlow(1,2);
Xhigh = fsolve(@(x) FUNC(x), [100 25/180*pi -30/180*pi], opts); % low AOA
aoa_start = Xhigh(1,2);
It might be necessary to change other options as well.  See the options documentation section for those details.  
.
0 Commenti
Vedere anche
Categorie
				Scopri di più su Line Plots in Help Center e File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!


