Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-79201.

2 次查看（过去 30 天）
Turgut Ataseven
Turgut Ataseven 2022-1-2
编辑： Torsten 2022-1-3
Hi. First of all sorry for asking the same question again and again but I am really close to the solution (maybe not). I have 2 questions.
1) This 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. Different flight maneuvers are performed for the specified time intervals.
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,den.m,den_cr.m,den_des.m,drag.m,drag_cr.m,drag_des.m,lift.m,lift_cr.m,lift_des.m are functions in seperate tabs. Code in the main file (Flight.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]
x6 = W;
% solving 1st order ODE using numerical methods
t0=0;
tend=3960;
h=0.05;
N=(tend-t0)/h;
t=t0:h:tend;
x = zeros(6,length(t));
% Initial conditions
x(:,1)=[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];
for i=2:length(t)
if and (x(3,i-1) > 3608.92,x(3,i-1)<=5249.3) % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
t =0:60;
x4 = Velocities(3608.9:5249.3); % Changing speed [m/s]
x5 = 0; % Changing head angle [deg]
f = fuel_cl(3608.9:5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cl(3608.9:5249.3); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 5; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag(3608.9:5249.3,x4); % Changing drag coefficient
Cl = lift(3608.9:5249.3,x4); % Changing lift coefficient
p = den(3608.9:5249.3); % Changing density [kg/m3]
elseif and (x(3,i-1) > 5249.3,x(3,i-1)<=5249.3) % Perform cruise flight for t=60 minutes.
for t=60:3660
x4 = Cruise_Vel(5249.3); % Changing speed [m/s]
x5 = 0; % Changing head angle [deg]
f = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cr(5249.3); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,x4); % Changing drag coefficient
Cl = lift_cr(5249.3,x4); % Changing lift coefficient
p = den_cr(5249.3); % Changing density [kg/m3]
end
% Turn with β=30 bank angle until heading is changed by η=270◦.
for t=3660:3720
x4 = Cruise_Vel(5249.3); % Changing speed [m/s]
x5 = 0:30:270; % Changing head angle [deg]
f = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cr(5249.3); % Changing thrust [N]
u2 = 30; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,x4); % Changing drag coefficient
Cl = lift_cr(5249.3,x4); % Changing lift coefficient
p = den_cr(5249.3); % Changing density [kg/m3]
end
elseif and (x(3,i-1) < 5249.3,x(3,i-1)>=3608.9) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
t=3720:3780;
x4 = Des_Vel(5249.3,3608.9); % Changing speed [m/s]
x5 = 270; % Changing head angle [deg]
f = fuel_des(5249.3,3608.9); % Changing fuel flow [kg/min]
u1 = Thr_des(5249.3,3608.9); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 4; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_des(5249.3,3608.9,x4); % Changing drag coefficient
Cl = lift_des(5249.3,3608.9,x4); % Changing lift coefficient
p = den_des(5249.3,3608.9); % Changing density [kg/m3]
elseif and (x(3,i-1) < 3608.9,x(3,i-1)>=3608.9) % Complete a 360◦ turn (loiter) at level flight.
t = 3780:3900;
x4 = Cruise_Vel(3608.9); % Changing speed [m/s]
lon = [270 300 360 60 120 180 240 270];
x5 = wrapTo360(lon); % Changing head angle [deg]
f = fuel_cr(3608.9); % Changing fuel flow [kg/min]
u1 = Thr_cr(3608.9); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(3608.9,x4); % Changing drag coefficient
Cl = lift_cr(3608.9,x4); % Changing lift coefficient
p = den_cr(3608.9); % Changing density [kg/m3]
elseif and (x(3,i-1) < 3608.9,x(3,i-1)>=2624.67) % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
t=3900:3960;
x4 = Des_Vel(3608.9,2624.67); % Changing speed [m/s]
x5 = 270; % Changing head angle [deg]
f = fuel_des(3608.9,2624.67); % Changing fuel flow [kg/min]
u1 = Thr_des(3608.9,2624.67); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 4.5; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_des(3608.9,2624.67,x4); % Changing drag coefficient
Cl = lift_des(3608.9,2624.67,x4); % Changing lift coefficient
p = den_des(3608.9,2624.67); % Changing density [kg/m3]
else
fprintf("A problem occured.");
end
dx1dt = x4 .* cos(x5) .* cos(u3);
dx2dt = x4 .* sin(x5) .* cos(u3);
dx3dt = x4 .* sin(u3);
dx4dt = -C_D.*S.*p.*(x4.^2)./(2.*x6)-g0.*sin(u3)+u1./x6;
dx5dt = -Cl.*S.*p.*x4./(2.*x6).*sin(u2);
dx6dt = -f;
x(1,i)= x(1,i-1) + h * dx1dt; %%% line 136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x(2,i)= x(2,i-1) + h * dx2dt;
x(3,i)= x(3,i-1) + h * dx3dt;
x(4,i)= x(4,i-1) + h * dx4dt;
x(5,i)= x(5,i-1) + h * dx5dt;
x(6,i)= x(6,i-1) + h * dx6dt;
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/60);
I got this error:
Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-79201.
Error in Flight (line 136)
x(1,i)= x(1,i-1) + h * dx1dt;
dx1dt is the 1-by-79201 one there. Is it because I have set the output lenghts of aforementioned functions in seperate tabs as linspace(H1,H2,79201)? If yes what should I hange, if no what should I do? 79201 is the length of t. Before doing that I have faced with array size problems. I am open to better suggestions.
2) With these lines I've tried to state that altitude [x(3,i-1)] remains constant:
elseif and (x(3,i-1) > 5249.3,x(3,i-1)<=5249.3) % Perform cruise flight for t=60 minutes.
and to state descent:
elseif and (x(3,i-1) < 3608.9,x(3,i-1)>=2624.67) % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
Are they correct? If no, can you provide the corrected versions?
Thanks.
  3 个评论
显示 1更早的评论
Turgut Ataseven
Turgut Ataseven 2022-1-2
@Torsten Thanks for the response and patience.
Sorry but I have thought about your suggestions at this and previous questions but could not come up with a reasonable solution. Do you want me to use preallocations like:
x1 = zeros(1);
x2 = zeros(1);
...
Cl = zeros(1);
p = zeros(1);
If not (probably), can you give just one example on that? [x4; x5; f;u1; C_D; Cl; p] varies with altitude (x3) and I don't know how to state that without using arrays which gets the value from functions like this (79201 is the length(t)):
function [Vtas_cl] = Velocities(H1,H2)
%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%
Vcl_1 = 335; % Standard calibrated airspeed [kt]
Vcl_2 = 172.3; % Standard calibrated airspeed [kt] -> [m/s] (To find the Mach transition altitude)
Vcl_2_in_knots = 335; % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cl = 0.86; % Standard calibrated airspeed [kt]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
g0 = 9.80665; % Gravitational acceleration [m/s2]
a0= 340.294; % Speed of Sound [m/s]
Vd_CL1 = 5; % Climb speed increment below 1500 ft (jet)
Vd_CL2 = 10; % Climb speed increment below 3000 ft (jet)
Vd_CL3 = 30; % Climb speed increment below 4000 ft (jet)
Vd_CL4 = 60; % Climb speed increment below 5000 ft (jet)
Vd_CL5 = 80; % Climb speed increment below 6000 ft (jet)
CV_min = 1.3; % Minimum speed coefficient
Vstall_TO = .14200E+03; % Stall speed at take-off [KCAS]
CAS_climb = Vcl_2;
Mach_climb = M_cl;
delta_trans = (((1+((K-1)/2)*(CAS_climb/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_climb^2)^(K/(K-1)))-1); % Pressure ratio at the transition altitude
teta_trans = delta_trans ^ (-Bt*R/g0); % Temperature ratio at the transition altitude
H_p_trans_climb = (1000/0.348/6.5)*(T0*(1-teta_trans)); % Transition altitude for climb [ft]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H_climb = linspace(H1,H2,79201); %%%%% Input %%%%%%%%%%%
Vnom_climb_jet = zeros(1, length(H_climb));
for k1 = 1:length(H_climb)
if (0<=H_climb(k1)&&H_climb(k1)<1500)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;
elseif (1500<=H_climb(k1)&&H_climb(k1)<3000)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
elseif (3000<=H_climb(k1)&&H_climb(k1)<4000)
Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL3;
elseif (4000<=H_climb(k1)&&H_climb(k1)<5000)
Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL4;
elseif (5000<=H_climb(k1)&&H_climb(k1)<6000)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL5;
elseif (6000<=H_climb(k1)&&H_climb(k1)<10000)
Vnom_climb_jet (k1)= min(Vcl_1,250);
elseif (10000<=H_climb(k1)&&H_climb(k1)<=H_p_trans_climb)
Vnom_climb_jet(k1) = Vcl_2_in_knots;
elseif (H_p_trans_climb<H_climb(k1))
Vnom_climb_jet(k1) = M_cl;
end
Vcas_cl(k1) = Vnom_climb_jet(k1)* 0.514; % [kn] -> [m/s]
H_climb (k1)= H_climb(k1) * 0.3048; % [feet] -> [m]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0; % Value of the real temperature T in ISA conditions [K]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
P0 = 101325; % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665; % Gravitational acceleration [m/s2]
p0 = 1.225; % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;
T(k1) = T0 + deltaT + Bt * H_climb(k1); % Temperature [K]
P (k1)= P0*((T(k1)-deltaT)/T0).^((-g0)/(Bt*R)); % Pressure [Pa]
p (k1)= P(k1) ./ (R*T(k1)); % Density [kg/m^3]
Vtas_cl(k1) = (2*P(k1)/visc/p(k1)*((1 + P0/P(k1)*((1 + visc*p0*Vcas_cl(k1)*Vcas_cl(k1)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
% Output %%%%%%%%%%%%
end
end
Torsten
Torsten 2022-1-3
编辑：Torsten 2022-1-3
I mean:
Instead of setting
H_climb = linspace(H1,H2,79201);
in "Velocities", simply set
H_climb = x3
where x3 is your solution for the altitude at time t(i-1) (thus one single value).
This will return one single value for x4.
Same for the other deduced model variables x5, f, u1, u2, u3, C_D, Cl and p.

请先登录，再进行评论。

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

回答（0 个）

类别

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

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by