Effacer les filtres
Effacer les filtres

I keep getting the solve stops prematurely for fsolve?

4 vues (au cours des 30 derniers jours)
Rachel Ong
Rachel Ong le 17 Mar 2022
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 commentaires
Torsten
Torsten le 17 Mar 2022
I can't reproduce the error. Under Ocatve's "fsolve", your code finishes without error message.
Rachel Ong
Rachel Ong le 17 Mar 2022
I dont know why I keep getting the error
Solver stopped prematurely.
fsolve stopped because it exceeded the iteration limit,
options.MaxIterations = 4.000000e+02.

Connectez-vous pour commenter.

Réponses (1)

Star Strider
Star Strider le 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.
.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by