Index exceeds the number of array elements (1).

1 次查看（过去 30 天）
Turgut Ataseven
Turgut Ataseven 2021-12-26
回答： Walter Roberson 2021-12-26
Hi. My code tries to solve 6 ODEs with 6 state variables [horizontal position (x1 and x2), altitude (x3), the true airspeed (x4), the heading angle (x5) and the mass of the aircraft (x6)] and 3 control inputs [engine thrust (u1), the bank angle (u2) and the flight path angle (u3)] by using Euler's method.
Velocities.m, Cruise_Vel.m, Des_Vel.m, Thr_cl.m, Thr_cr.m, Thr_des.m, fuel_cl.m, fuel_cr.m, fuel_des.m,drag.m,lift.m,den.m are functions in seperate files. Code in the main file (FlightPlan.m) is:
% Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
% Perform cruise flight for t=60 minutes.
% Turn with β=30 bank angle until heading is changed by η=270◦.
% Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
% Complete a 360◦ turn (loiter) at level flight.
% Descent to h3=800 [m] with κ=4.5◦ flight path angle.
% Aircraft Properties
W = .44225E+06; % .44225E+03 tons = .44225E+06 kg
S = .51097E+03; % Surface Area [m^2]
g0 = 9.80665; % Gravitational acceleration [m/s2]
% Equations
% dx1 = x4*cos(x5/180*pi)*cos(u3/180*pi);
% dx2 = x4*sin(x5/180*pi)*cos(u3/180*pi);
% dx3 = x4*sin(u3/180*pi);
% dx4 = -C_D*S*p*x4.*x4/(2*x6)-g0*sin(u3/180*pi)+u1/x6;
% dx5 = -Cl*S*p*x4/(2*x6)/sin(u2/180*pi);
% dx6 = -f;
% Initial conditions
x1(1) = 0; % Initial position [m]
x2(1) = 0; % Initial position [m]
x3(1) = 3608.92; % Initial altitude [ft]
x4(1) = 1.0e+02 * 1.161544478045788; % Initial speed [m/s]
x5(1) = 0; % Assuming aircraft headed to North initially.
x6(1) = W; % Initial mass [kg]
u1(1) = 1.0e+05 * 7.078900012026454; % Initial thrust [N]
u2(1) = 0; % Initial bank angle [deg]
u3(1) = 5; % Initial flight path angle [deg]
C_D(1) = 0.026278212364444; % Initial drag coefficient
Cl(1) = 0.124852132431623; % Initial lift coefficient
p(1) = 1.027625283894403; % Initial density [kg/m3]
f(1) = 1.0e+02 * 4.497203403070063; % Initial fuel flow [kg/min]
dx1dt = @(x4,x5,u3) x4.*cos(x5*pi/180).*cos(u3*pi/180);
dx2dt = @(x4,x5,u3) x4*sin(x5*pi/180)*cos(u3*pi/180);
dx3dt = @(x4,u3) x4*sin(u3*pi/180);
dx4dt = @(C_D,p,x6,x4,u3,u1) -C_D*S*p*x4.*x4/(2*x6)-g0*sin(u3)+u1/x6;
dx5dt = @(Cl,p,x6,x4,u2) -Cl*S*p*x4/(2*x6)/sin(u2);
dx6dt = @(f) -f;
% solving 1st order ODE using numerical methods
t0=0;
tend=66;
y0=0;
h=0.2;
N=(tend-t0)/h;
t=t0:h:tend;
if and (t >= 0,t<=1) % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
x3(t) = linspace(3608.9,5249.3,2); % Changing altitude [m] -> [ft]
x4(t) = Velocities(3608.9,5249.3); % Changing speed [m/s]
x5(t) = 0; % Changing head angle [deg]
f(t) = fuel_cl(3608.9,5249.3); % Changing fuel flow [kg/min]
u1(t) = Thr_cl(3608.9,5249.3); % Changing thrust [N]
u2(t) = 0; % Changing bank angle [deg]
u3(t) = 5; % Changing flight path angle [deg]
V_ver(t) = x4(t)*sin(u3(t)*pi/180); % Changing vertical speed [m/s]
C_D = drag(3608.9,5249.3,x4(t)); % Changing drag coefficient
Cl = lift(3608.9,5249.3,x4(t)); % Changing lift coefficient
p = den(3608.9,5249.3); % Changing density [kg/m3]
end
if and (t >1,t<=61) % Perform cruise flight for t=60 minutes.
x3(t) = 5249.3; % Changing altitude [m] -> [ft]
x4(t) = Cruise_Vel(5249.3); % Changing speed [m/s]
x5(t) = 0; % Changing head angle [deg]
f(t) = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1(t) = Thr_cr(5249.3); % Changing thrust [N]
u2(t) = 0; % Changing bank angle [deg]
u3(t) = 0; % Changing flight path angle [deg]
V_ver(t) = x4(t)*sin(u3(t)*pi/180); % Changing vertical speed [m/s]
C_D = drag(5249.3,5249.3,x4(t)); % Changing drag coefficient
Cl = lift(5249.3,5249.3,x4(t)); % Changing lift coefficient
p = den(5249.3,5249.3); % Changing density [kg/m3]
end
if and (t >61,t<=62) % Turn with β=30 bank angle until heading is changed by η=270◦.
x3(t) = 5249.3; % Changing altitude [m] -> [ft]
x4(t) = Cruise_Vel(5249.3); % Changing speed [m/s]
x5(t) = 0:30:270; % Changing head angle [deg]
f(t) = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1(t) = Thr_cr(5249.3); % Changing thrust [N]
u2(t) = 30; % Changing bank angle [deg]
u3(t) = 0; % Changing flight path angle [deg]
V_ver(t) = x4(t)*sin(u3(t)*pi/180); % Changing vertical speed [m/s]
C_D = drag(5249.3,5249.3,x4(t)); % Changing drag coefficient
Cl = lift(5249.3,5249.3,x4(t)); % Changing lift coefficient
p = den(5249.3,5249.3); % Changing density [kg/m3]
end
if and (t >62,t<=63) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
x3(t) = linspace(5249.3,3608.9,2); % Changing altitude [m] -> [ft]
x4(t) = Des_Vel(5249.3,3608.9); % Changing speed [m/s]
x5(t) = 270; % Changing head angle [deg]
f(t) = fuel_des(5249.3,3608.9); % Changing fuel flow [kg/min]
u1(t) = Thr_des(5249.3,3608.9); % Changing thrust [N]
u2(t) = 0; % Changing bank angle [deg]
u3(t) = 4; % Changing flight path angle [deg]
V_ver(t) = x4(t)*sin(u3(t)*pi/180); % Changing vertical speed [m/s]
C_D = drag(5249.3,3608.9,x4(t)); % Changing drag coefficient
Cl = lift(5249.3,3608.9,x4(t)); % Changing lift coefficient
p = den(5249.3,3608.9); % Changing density [kg/m3]
end
if and (t >63,t<=65) % Complete a 360◦ turn (loiter) at level flight.
x3(t) = 3608.9; % Changing altitude [m] -> [ft]
x4(t) = Cruise_Vel(3608.9); % Changing speed [m/s]
lon = [270 300 360 60 120 180 240 270];
x5(t) = wrapTo360(lon); % Changing head angle [deg]
f(t) = fuel_cr(3608.9); % Changing fuel flow [kg/min]
u1(t) = Thr_cr(3608.9); % Changing thrust [N]
u2(t) = 0; % Changing bank angle [deg]
u3(t) = 0; % Changing flight path angle [deg]
V_ver(t) = x4(t)*sin(u3(t)*pi/180); % Changing vertical speed [m/s]
C_D = drag(3608.9,3608.9,x4(t)); % Changing drag coefficient
Cl = lift(3608.9,3608.9,x4(t)); % Changing lift coefficient
p = den(3608.9,3608.9); % Changing density [kg/m3]
end
if and (t >65,t<=66) % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
x3(t) = linspace(3608.9,2624.67,2); % Changing altitude [m] -> [ft]
x4(t) = Des_Vel(3608.9,2624.67); % Changing speed [m/s]
x5(t) = 270; % Changing head angle [deg]
f(t) = fuel_des(3608.9,2624.67); % Changing fuel flow [kg/min]
u1(t) = Thr_des(3608.9,2624.67); % Changing thrust [N]
u2(t) = 0; % Changing bank angle [deg]
u3(t) = 4.5; % Changing flight path angle [deg]
V_ver(t) = x4(t)*sin(u3(t)*pi/180); % Changing vertical speed [m/s]
C_D = drag(3608.9,2624.67,x4(t)); % Changing drag coefficient
Cl = lift(3608.9,2624.67,x4(t)); % Changing lift coefficient
p = den(3608.9,2624.67); % Changing density [kg/m3]
end
for i=1:N
x1(i+1)=x1(i)+h*dx1dt(x4(i),x5(i),u3(i)); % line 152 %%%%%%%%%%%%%%%%%%%%%
x2(i+1)=x2(i)+h*dx2dt(x4(i),x5(i),u3(i));
x3(i+1)=x3(i)+h*dx3dt(x4(i),u3(i));
x4(i+1)=x4(i)+h*dx4dt(C_D(i),p(i),x6(i),x4(i),u3(i),u1(i));
x5(i+1)=x5(i)+h*dx5dt(Cl(i),p(i),x6(i),x4(i),u2(i));
x6(i+1)=x6(i)+h*dx6dt(f(i));
end
tot=cell2mat(f); % Total fuel consumption during mission [kg/min]
Tot_fuel=sum(tot);
figure(1)
plot3(x1(:),x2(:),x3(:)); % 3D position graph
figure(2)
plot(t,x4(:)); % Vtas − Time graph
figure(3)
plot(t,V_ver(:)); % V_vertical − Time graph
figure(4)
plot(t,x5(:)); % Heading − Time graph
figure(5)
plot(t,x6(:)); % Mass − Time graph
figure(6)
plot(t,u1(:)); % Thrust − Time graph
figure(7)
plot(t,u2(:)); % Bank Angle − Time graph
figure(8)
plot(t,u3(:)); % Flight Path Angle − Time graph
fprintf('Total fuel consumption during mission is %.2f [kg]',Tot_fuel*tend);
legend('Euler');
When I run the code, this error displays:
Index exceeds the number of array elements (1).
Error in FlightPlan (line 152)
x1(i+1)=x1(i)+h*dx1dt(x4(i),x5(i),u3(i));
After checking the Workspace, I saw that state variables (x1,x2,x3,x4,x5,x6) were updated as 1x2 matrices (first columns are initial conditions).
I want your help to utilize the Euler's Method on my code. How should I update the loops to plot Variable vs Time graphs? Is using (t) in loops is problematic? I need your help.
Thanks.

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2021-12-26
t=t0:h:tend;
if and (t >= 0,t<=1) % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
t is a vector. and (t >= 0,t<=1) is a vector. if of a vector is considered true only if all of the values in the vector are non-zero; in the context of a logical test, that means that the condition had to be met for all of the values.
You need to switch from if to using logical indexing.
https://blogs.mathworks.com/loren/2013/02/20/logical-indexing-multiple-conditions/

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by