I keep getting the solve stops prematurely for fsolve?
6 views (last 30 days)
Show older comments
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 Comments
Torsten
on 17 Mar 2022
I can't reproduce the error. Under Ocatve's "fsolve", your code finishes without error message.
Answers (1)
Star Strider
on 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 Comments
See Also
Categories
Find more on Nonlinear ARX Models in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!