I keep getting the solve stops prematurely for fsolve?
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
le 17 Mar 2022
I can't reproduce the error. Under Ocatve's "fsolve", your code finishes without error message.
Réponses (1)
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.
.
0 commentaires
Voir également
Catégories
En savoir plus sur Nonlinear ARX Models dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!