how to code transfer function
Show older comments
Hi, I would like to know how to code and plot this tf:
This in my attempt;
d = (e^-(t*s)/t1*s + 1)*dc;
not sure about it, any help will be greatly appreciated. Thanks much.
10 Comments
Cesar Cardenas
on 19 Sep 2022
Paul
on 19 Sep 2022
Hi Cesar,
I assume that deltac represents the control actuator commands: de0, da0, and dr0. In the code, these are treated as global constants. In reality they are functions of time. How are de0, da0, and dr0 to be devloped as functions of time as the simulation model runs?
Cesar Cardenas
on 22 Sep 2022
Paul
on 22 Sep 2022
those are actuator commands as global constants, not sure how they develop,
Do you want a solution that works only if the actuator commands are constant? Or do you need a more general solution? If the latter, we'll need model equations for how the actuator commands evolve as a function of time.
Cesar Cardenas
on 22 Sep 2022
Cesar Cardenas
on 22 Sep 2022
Cesar Cardenas
on 22 Sep 2022
Star Strider
on 22 Sep 2022
Cesar Cardenas
on 22 Sep 2022
Star Strider
on 22 Sep 2022
My pleasure!
Accepted Answer
There seems to be two questions in this Question.
The first question is: How to plot the transfer function. For values of tau0 and tau1, the transfer function from delta_c to delta is
H = tf(1,[tau1 1],'InputDelay',tau0);
The concept of "plot the transfer function" is ambiguous. Typical plots of a transfer function might be
bode(H)
impulse(H)
What kind of plot is desired?
Or, if you really want to plot the actuator position in response to an actuator command (which isn't really "plot the transfer function"), then we need to define the actuator command of interest as a function of time and then use functions like step or lsim as apprropriate for the specific actuator command input of interest.
The second question is how to add those actuator dynamics to the simulation. For constant actuator command inputs, and assuming that the actuator commands are zero for t < 0 (i.e., we are injecting a step command into the actuator) the approach would be as follows. Referring to the code in this comment ...
Add three additional elements to the state vector, y. These will be the values of de, da, and dr. Make sure to initialize these values in y0
In function STOL_Eom:
Extract de, da, and dr from the state vector.
Compute:
dec = de0*(t >tau0). Note that this approach is a little loose. Technically, we'd have to integrate from t = 0 to t = tau0, stop the simulation, and then restart it. But we can try to make this a little simpler, assuming that tau0 is reasonably small as would normally be the case for an actuator.
dedot = (dec - de) / tau1
Repeat for da and dr.
Add dedot, dadot, and drdot to the end of the xdot vector
In the main code, call ode45 with with additional option of InitialStep = tau0. See ode45 to see how to use the odeset function to specify this option.
The idea here is that ode45 will hopefully take a small first step from t = 0 to t = tau0 to catch the step change in the actuator commands. This can be checked in the t vector after the simulation completes.This approach is a bit loose, but may be sufficient for your needs. If it doesn't work as well as you'd like, come back and we can address this issue a bit better.
Based on the structure of the code, I guess tau0 and tau1 would be added to the global variable list with values assigned in the main script.
18 Comments
Cesar Cardenas
on 22 Sep 2022
Paul
on 22 Sep 2022
Hi Cesar,
I suggest you try it yourself and then come back if having a problem. Is there something that's unclear about the suggested approachj? If so, ask away and we'll try to get it sorted out.
Cesar Cardenas
on 22 Sep 2022
Cesar Cardenas
on 22 Sep 2022
Paul
on 22 Sep 2022
In STOL_EOM:
at the top:
de = x(13);
da = x(14);
dr = x(15);
Comment out this, at least for now
% Control deflections
%de = deEq;
%da = 0;
%da = - 0.5*phi - 2*pr; % Add some roll stiffness and damping
%dr = 0;
%dT = dTEq;
In the middle:
dec = de0*(t >= tau0);
dac = da0*(t >= tau0);
drc = dr0*(t >= tau0);
dedot = (dec - de)/tau_1;
dadot = (dec - de)/tau_1;
drdot = (dec - de)/tau_1;
At the bottom:
xdot(13) = dedot;
xdot(14) = dadot;
xdot(15) = drdot;
In the main script:
y0 = [X0; Theta0; V0; omega0; de_initial; da_initial; dr_initial];
where *_initial are the initial values of the actuator deflections.
[t,y] = ode45('STOL_EOM',[0:0.1:t_final]',y0,odeset('InitialStep'),tau0);
Cesar Cardenas
on 22 Sep 2022
I had typos in my response. Corrected below. I also noticed some other issues, corrected below. Let's see if it at least runs now.
Nope, still an error. What is hat(e3), etc. in the computation of RIB? Is that just supposed to be e3?
Also, I don't think de0, da0, or dr0 are defined anywhere.
% GenericFixedWingScript.m solves the nonlinear equations of motion for a
% fixed-wing aircraft flying in ambient wind with turbulence.
global e1 e2 e3 rho m g S b c Inertia Power ...
Cx0 Cxu Cxw Cxw2 ...
Cz0 Czw Czw2 Czq Czde ...
Cm0 Cmw Cmq Cmde ...
Cy0 Cyv Cyp Cyr Cyda Cydr ...
Cl0 Clv Clp Clr Clda Cldr ...
Cn0 Cnv Cnv2 Cnp Cnr Cnda Cndr ...
dTEq deEq ...
VonKarmanFlag Amp Phase Omega Phase2 Omega2 ...
tau0 tau_1
tau0 and tau_1 added to the global list
tau0 = 1; % actuator delay
tau_1 = 2; % time constant
% Basis vectors
e1 = [1;0;0];
e2 = [0;1;0];
e3 = [0;0;1];
% Atmospheric and gravity parameters (Constant altitude: Sea level)
a = 340.3; % Speed of sound (m/s)
rho = 1.225; % Density (kg/m^3)
g = 9.80665; % Gravitational acceleration (m/s^2)
% Aircraft parameters (Bix3)
m = 1.202; % Mass (kg)
W = m*g; % Weight (Newtons)
Ix = 0.095; % Roll inertia (kg-m^2)
Iy = 215e3; % Pitch inertia (kg-m^2)
Iz = 447e3; % Yaw inertia (kg-m^2)
Ixz = 0; % Roll/Yaw product of inertia (kg-m^2)
Inertia = [Ix, 0, -Ixz; 0, Iy, 0; -Ixz, 0, Iz];
S = 0.285; % Wing area (m^2)
b = 1.54; % Wing span (m)
c = 0.188; % Wing chord (m)
AR = (b^2)/S; % Aspect ratio
% Longitudinal nondimensional stability and control derivatives
Cx0 = 0; % Define CT0 = 0.197 below. (Cx0 = 0.197 in [Simmons et al, 2019])
Cxu = -0.156;
Cxw = 0.297;
Cxw2 = 0.960;
Cz0 = -0.179;
Czw = -5.32;
Czw2 = 7.02;
Czq = -8.20;
Czde = -0.308;
Cm0 = 0.0134;
Cmw = -0.240;
Cmq = -4.49;
Cmde = -0.364;
% Lateral-directional nondimensional stability and control derivatives
Cy0 = 0;
Cyv = -0.251;
Cyp = 0.170;
Cyr = 0.350;
Cyda = 0.103;
Cydr = 0.0157;
Cl0 = 0;
Clv = -0.0756;
Clp = -0.319;
Clr = 0.183;
Clda = 0.170;
Cldr = -0.0117;
Cn0 = 0;
Cnv = 0.0408;
Cnv2 = 0.0126;
Cnp = -0.242;
Cnr = -0.166;
Cnda = 0.0416;
Cndr = -0.0618;
% Solve for wings-level, constant-altitude equilibrium flight condition
VEq = 12; % Initial guess for speed (m/s)
thetaEq = 2*(pi/180); % Initial guess for pitch (rad)
deEq = -(Cm0 + Cmw*sin(thetaEq))/Cmde; % Initial guess for elevator (rad)
Power = 200; % Maximum Power (W)
CT0 = 0.197; % CT0 = Cx0 in [Simmons et al, 2019]
dTEq = (CT0*rho*S*(VEq^3))/(2*Power); % Nominal throttle setting
% Define initial state
X0 = zeros(3,1); % (m)
Theta0 = [0;thetaEq;0]; % (rad)
V0 = [VEq*cos(thetaEq);0;VEq*sin(thetaEq)]; % (m/s)
omega0 = zeros(3,1); % (rad/s)
de, da, and dr are scalars and should be initialized as such
%de_initial = zeros(3,1);
%da_initial = zeros(3,1);
%dr_initial = zeros(3,1);
de_initial = 0;
da_initial = 0;
dr_initial = 0;
% y0 = [X0; Theta0; V0; omega0];
y0 = [X0; Theta0; V0; omega0; de_initial; da_initial; dr_initial];
t_final = 10;
% [t,y] = ode45('GenericFixedWingEOM',(0:0.1:t_final)',y0);
Fix the typo in the call to odeset. The odefun argument should be a function handle.
%[t,y] = ode45('GenericFixedWingEOM',[0:0.1:t_final]',y0,odeset('InitialStep'),tau0);
[t,y] = ode45(@GenericFixedWingEOM,[0:0.1:t_final]',y0,odeset('InitialStep',tau0));
%%%% CODE TO COMPUTE FLOW RELATIVE VELOCITY, ETC %%%%
figure(1)
subplot(2,1,1)
plot(t,y(:,1:3))
legend('X','Y','Z')
ylabel('Position (m)')
subplot(2,1,2)
plot(t,y(:,4:6)*(180/pi))
legend('\phi','\theta','\psi')
ylabel('Attitude (deg)')
figure(2)
subplot(2,1,1)
plot(t,y(:,7:9))
legend('u','v','w')
ylabel('Velocity (m/s)')
subplot(2,1,2)
plot(t,y(:,10:12)*(180/pi))
ylabel('Angular Velocity (deg/s)')
legend('p','q','r')
function xdot = GenericFixedWingEOM(t,x)
% GenericFixedWingEOM.m
%
% Nonlinear equations of motion for a rigid, fixed-wing airplane.
global e1 e2 e3 rho m g S b c Inertia Power ...
Cx0 Cxu Cxw Cxw2 ...
Cz0 Czw Czw2 Czq Czde ...
Cm0 Cmw Cmq Cmde ...
Cy0 Cyv Cyp Cyr Cyda Cydr ...
Cl0 Clv Clp Clr Clda Cldr ...
Cn0 Cnv Cnv2 Cnp Cnr Cnda Cndr ...
dTEq deEq ...
VonKarmanFlag Amp Phase Omega Phase2 Omega2 ...
tau0 tau_1 % added
% Parse out state vector components
X = x(1:3);
Theta = x(4:6);
phi = Theta(1);
theta = Theta(2);
psi = Theta(3);
V = x(7:9);
u = V(1);
v = V(2);
w = V(3);
omega = x(10:12);
p = omega(1);
q = omega(2);
r = omega(3);
RIB = expm(psi*hat(e3))*expm(theta*hat(e2))*expm(phi*hat(e1));
LIB = [1, sin(phi)*tan(theta), cos(phi)*tan(theta);
0, cos(phi), -sin(phi);
0, sin(phi)/cos(theta), cos(phi)/cos(theta)];
% if VonKarmanFlag == 1
%
% %%% WIND VELOCITY (1D VON KARMAN) %%%
% WX = Amp(:,1)'*cos(Omega*X(1)+Phase(:,1));
% WY = Amp(:,2)'*cos(Omega*X(1)+Phase(:,2));
% WZ = Amp(:,3)'*cos(Omega*X(1)+Phase(:,3));
% W = [WX; WY; WZ];
%
% %%% WIND GRADIENT (1D VON KARMAN) %%%
% WXX = -(Amp(:,1).*Omega)'*sin(Omega*X(1)+Phase(:,1));
% WYX = -(Amp(:,2).*Omega)'*sin(Omega*X(1)+Phase(:,2));
% WZX = -(Amp(:,3).*Omega)'*sin(Omega*X(1)+Phase(:,3));
% GradW = [WXX, 0, 0; WYX, 0, 0; WZX, 0, 0];
%
% elseif VonKarmanFlag == 2
%
% %%% WIND VELOCITY (2D VON KARMAN) %%%
% WX = Amp(:,1)'*(cos(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WY = Amp(:,2)'*(cos(Omega*X(1)+Phase(:,2)).*cos(Omega2*X(2)+Phase2(:,2)));
% WZ = Amp(:,3)'*(cos(Omega*X(1)+Phase(:,3)).*cos(Omega2*X(2)+Phase2(:,3)));;
% W = [WX; WY; WZ];
%
% %%% WIND GRADIENT (2D VON KARMAN) %%%
% WXX = -(Amp(:,1).*Omega)'*(sin(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WYX = -(Amp(:,2).*Omega)'*(sin(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WZX = -(Amp(:,3).*Omega)'*(sin(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WXY = -(Amp(:,1).*Omega2)'*(cos(Omega*X(1)+Phase(:,1)).*sin(Omega2*X(2)+Phase2(:,1)));
% WYY = -(Amp(:,2).*Omega2)'*(cos(Omega*X(1)+Phase(:,1)).*sin(Omega2*X(2)+Phase2(:,1)));
% WZY = -(Amp(:,3).*Omega2)'*(cos(Omega*X(1)+Phase(:,1)).*sin(Omega2*X(2)+Phase2(:,1)));
% GradW = [WXX, WXY, 0; WYX, WYY, 0; WZX, WZY, 0];
%
% end
% Wind velocity in body frame, air-relative velocity, and dynamic pressure
% W = RIB'*W; % Convert wind velocity to body frame
W = [0;0;0];
Vr = V - W;
ur = Vr(1);
vr = Vr(2);
wr = Vr(3);
PDyn = (1/2)*rho*norm(Vr)^2;
% Wind gradient in body frame
%Phi = RIB'*GradW'*RIB;
% Effective angular velocity of the fluid (Note: Etkin argues for
% a different formulation of the "angular velocity of the wind")
%temp = -(Phi-Phi');
%omegaf = [-temp(2,3); temp(1,3); -temp(1,2)];
omegaf =[0;0;0];
omegar = omega - omegaf;
pr = omegar(1);
qr = omegar(2);
rr = omegar(3);
%%% AERODYNAMIC FORCES AND MOMENTS %%%
% Aerodynamic angles
%beta = asin(Vr(2)/norm(Vr));
%alpha = atan2(Vr(3),Vr(1));
% Control deflections
% %de = deEq;
% da = 0;
% %da = - 0.5*phi - 2*pr; % Add some roll stiffness and damping
% dr = 0;
% dT = dTEq;
Why not pull these out of x at the top of the function with the other state variables?
de = x(13);
da = x(14);
dr = x(15);
Defined in the global script
%tau0 = 1;
%tau_1 = 2;
dec = de0*(t >= tau0);
dac = da0*(t >= tau0);
drc = dr0*(t >= tau0);
dadot and drdot equations are incorrect. Looks like I had these incorrect in my comment above (copy/paste error, obviously).
dedot = (dec - de)/tau_1;
%dadot = (dec - de)/tau_1;
%drdot = (dec - de)/tau_1;
dadot = (dac - da)/tau_1;
drdot = (drc - dr)/tau_1;
The xdot vector has to be constructed below.
%xdot(13) = dedot;
%xdot(14) = dadot;
%xdot(15) = drdot;
% Aerodynamic force coefficients
Cx = Cx0 + Cxu*(ur/norm(Vr)) + Cxw*(wr/norm(Vr)) + Cxw2*(wr/norm(Vr))^2;
Cy = Cy0 + Cyv*(vr/norm(Vr)) + Cyp*((b*pr)/(2*norm(V))) + Cyr*((b*rr)/(2*norm(Vr))) ...
+ Cyda*da + Cydr*dr;
Cz = Cz0 + Czw*(wr/norm(Vr)) + Czw2*(wr/norm(Vr))^2 + Czq*((c*qr)/(2*norm(Vr))) ...
+ Czde*de;
Thrust = dT*Power/norm(Vr); % Constant power
X = PDyn*S*Cx + Thrust;
Y = PDyn*S*Cy;
Z = PDyn*S*Cz;
Force_Aero = [X; Y; Z];
% Components of aerodynamic moment (modulo "unsteady" terms)
Cl = Cl0 + Clv*(vr/norm(Vr)) + Clp*((b*pr)/(2*norm(Vr))) + Clr*((b*rr)/(2*norm(Vr))) + ...
Clda*da + Cldr*dr;
Cm = Cm0 + Cmw*(wr/norm(Vr)) + Cmq*((c*qr)/(2*norm(Vr))) + Cmde*de;
Cn = Cn0 + Cnv*(vr/norm(Vr)) + Cnv2*(vr/norm(Vr))^2 + Cnp*((b*pr)/(2*norm(Vr))) + Cnr*((b*rr)/(2*norm(Vr))) + ...
Cnda*da + Cndr*dr;
L = PDyn*S*b*Cl;
M = PDyn*S*c*Cm;
N = PDyn*S*b*Cn;
Moment_Aero = [L; M; N];
% Sum of forces and moments
Force = RIB'*(m*g*e3) + Force_Aero;
Moment = Moment_Aero;
%%% EQUATIONS OF MOTION %%%
%%% KINEMATIC EQUATIONS %%%
XDot = RIB*V;
ThetaDot = LIB*omega;
%%% DYNAMIC EQUATIONS %%%
VDot = (1/m)*(cross(m*V,omega) + Force);
omegaDot = inv(Inertia)*(cross(Inertia*omega,omega) + Moment);
Add the actuator rates to the xdot vector
xdot = [XDot; ThetaDot; VDot; omegaDot; dedot; dadot; drdot];
end
Cesar Cardenas
on 23 Sep 2022
Cesar Cardenas
on 23 Sep 2022
Cesar Cardenas
on 23 Sep 2022
Cesar Cardenas
on 23 Sep 2022
It looks like you didn't actually use all of code changes I provided above. They have been reimplemented below, along with adding the actuation commands to the script. The code below runs, as demonstrated, for constant actuator commands.
You can try to implement a doublet actuator command as follows:
a. Remove de0, da0, dr0 from the script and global variable list
b. Create a function with the signature
dc = doublet(t, ...) that computes the acutator command as a function of time.
c. In GenericFixedWingEOM, prior the computation of the actuator derviatives:
de0 = doublet(t-tau0, ....). Do the same for da0 and dr0.
This approach might work ok for your needs if the width of the doublet steps are fairly wide. Again, plot the achieved actuator deflections to see if they are what they should be.
Good luck with your project.
global e1 e2 e3 rho m g S b c Inertia Power ...
Cx0 Cxu Cxw Cxw2 ...
Cz0 Czw Czw2 Czq Czde ...
Cm0 Cmw Cmq Cmde ...
Cy0 Cyv Cyp Cyr Cyda Cydr ...
Cl0 Clv Clp Clr Clda Cldr ...
Cn0 Cnv Cnv2 Cnp Cnr Cnda Cndr ...
dTEq deEq ...
VonKarmanFlag Amp Phase Omega Phase2 Omega2 tau0 tau_1 ...
de0 da0 dr0 % added
Define small constant actuator deflection commands for testing
de0 = 0.001;
da0 = 0.0005;
dr0 = 0.0002;
tau0 = 1; % actuator delay
tau_1 = 2; % time constant
% Basis vectors
e1 = [1;0;0];
e2 = [0;1;0];
e3 = [0;0;1];
% Atmospheric and gravity parameters (Constant altitude: Sea level)
a = 340.3; % Speed of sound (m/s)
rho = 1.225; % Density (kg/m^3)
g = 9.80665; % Gravitational acceleration (m/s^2)
% Aircraft parameters (Bix3)
m = 1.202; % Mass (kg)
W = m*g; % Weight (Newtons)
Ix = 0.095; % Roll inertia (kg-m^2)
Iy = 215e3; % Pitch inertia (kg-m^2)
Iz = 447e3; % Yaw inertia (kg-m^2)
Ixz = 0; % Roll/Yaw product of inertia (kg-m^2)
Inertia = [Ix, 0, -Ixz; 0, Iy, 0; -Ixz, 0, Iz];
S = 0.285; % Wing area (m^2)
b = 1.54; % Wing span (m)
c = 0.188; % Wing chord (m)
AR = (b^2)/S; % Aspect ratio
% Longitudinal nondimensional stability and control derivatives
Cx0 = 0; % Define CT0 = 0.197 below. (Cx0 = 0.197 in [Simmons et al, 2019])
Cxu = -0.156;
Cxw = 0.297;
Cxw2 = 0.960;
Cz0 = -0.179;
Czw = -5.32;
Czw2 = 7.02;
Czq = -8.20;
Czde = -0.308;
Cm0 = 0.0134;
Cmw = -0.240;
Cmq = -4.49;
Cmde = -0.364;
% Lateral-directional nondimensional stability and control derivatives
Cy0 = 0;
Cyv = -0.251;
Cyp = 0.170;
Cyr = 0.350;
Cyda = 0.103;
Cydr = 0.0157;
Cl0 = 0;
Clv = -0.0756;
Clp = -0.319;
Clr = 0.183;
Clda = 0.170;
Cldr = -0.0117;
Cn0 = 0;
Cnv = 0.0408;
Cnv2 = 0.0126;
Cnp = -0.242;
Cnr = -0.166;
Cnda = 0.0416;
Cndr = -0.0618;
% Solve for wings-level, constant-altitude equilibrium flight condition
VEq = 13.7; % Initial guess for speed (m/s)
thetaEq = 2*(pi/180); % Initial guess for pitch (rad)
deEq = -(Cm0 + Cmw*sin(thetaEq))/Cmde; % Initial guess for elevator (rad)
Power = 200; % Maximum Power (W)
CT0 = 0.197; % CT0 = Cx0 in [Simmons et al, 2019]
dTEq = (CT0*rho*S*(VEq^3))/(2*Power); % Nominal throttle setting
% Define initial state
X0 = zeros(3,1); % (m)
Theta0 = [0;thetaEq;0]; % (rad)
V0 = [VEq*cos(thetaEq);0;VEq*sin(thetaEq)]; % (m/s)
omega0 = zeros(3,1); % (rad/s)
de_initial = 0;
da_initial = 0;
dr_initial = 0;
%y0 = [X0; Theta0; V0; omega0];
y0 = [X0; Theta0; V0; omega0; de_initial; da_initial; dr_initial];
t_final = 10;
[t,y] = ode45(@GenericFixedWingEOM,[0:0.1:t_final]',y0,odeset('InitialStep',tau0));
%t_final = 30;
%[t,y] = ode45('GenericFixedWingEOM',[0:0.1:t_final]',y0);
%%%% CODE TO COMPUTE FLOW RELATIVE VELOCITY, ETC %%%%
figure(1)
subplot(2,1,1)
plot(t,y(:,1:3))
legend('X','Y','Z')
ylabel('Position (m)')
subplot(2,1,2)
plot(t,y(:,4:6)*(180/pi))
legend('\phi','\theta','\psi')
ylabel('Attitude (deg)')
figure(2)
subplot(2,1,1)
plot(t,y(:,7:9))
legend('u','v','w')
ylabel('Velocity (m/s)')
subplot(2,1,2)
plot(t,y(:,10:12)*(180/pi))
ylabel('Angular Velocity (deg/s)')
legend('p','q','r')
See if actuator reponse to constant input is delaybe by tau0 = 1. It is.
figure
plot(t,y(:,13:15))
legend('de','da','dr')
function xdot = GenericFixedWingEOM(t,x)
% GenericFixedWingEOM.m
%
% Nonlinear equations of motion for a rigid, fixed-wing airplane.
global e1 e2 e3 rho m g S b c Inertia Power ...
Cx0 Cxu Cxw Cxw2 ...
Cz0 Czw Czw2 Czq Czde ...
Cm0 Cmw Cmq Cmde ...
Cy0 Cyv Cyp Cyr Cyda Cydr ...
Cl0 Clv Clp Clr Clda Cldr ...
Cn0 Cnv Cnv2 Cnp Cnr Cnda Cndr ...
dTEq deEq ...
VonKarmanFlag Amp Phase Omega Phase2 Omega2 ...
tau0 tau_1 ...
de0 da0 dr0 % added
% Parse out state vector components
X = x(1:3);
Theta = x(4:6);
phi = Theta(1);
theta = Theta(2);
psi = Theta(3);
V = x(7:9);
u = V(1);
v = V(2);
w = V(3);
omega = x(10:12);
p = omega(1);
q = omega(2);
r = omega(3);
Extract actuator deflections from the state vector.
de = x(13);
da = x(14);
dr = x(15);
RIB = expm(psi*hat(e3))*expm(theta*hat(e2))*expm(phi*hat(e1));
LIB = [1, sin(phi)*tan(theta), cos(phi)*tan(theta);
0, cos(phi), -sin(phi);
0, sin(phi)/cos(theta), cos(phi)/cos(theta)];
W = [0;0;0];
Vr = V - W;
ur = Vr(1);
vr = Vr(2);
wr = Vr(3);
PDyn = (1/2)*rho*norm(Vr)^2;
omegaf =[0;0;0];
omegar = omega - omegaf;
pr = omegar(1);
qr = omegar(2);
rr = omegar(3);
%%% AERODYNAMIC FORCES AND MOMENTS %%%
% Aerodynamic angles
%beta = asin(Vr(2)/norm(Vr));
%alpha = atan2(Vr(3),Vr(1));
% Control deflections
%de = deEq;
%da = 0;
%da = - 0.5*phi - 2*pr; % Add some roll stiffness and damping
%dr = 0;
Need to assign to dT
dT = dTEq;
dec = de0*(t >= tau0);
dac = da0*(t >= tau0);
drc = dr0*(t >= tau0);
dedot = (dec - de)/tau_1;
%dadot = (dec - de)/tau_1;
%drdot = (dec - de)/tau_1;
dadot = (dac - da)/tau_1;
drdot = (drc - dr)/tau_1;
% Aerodynamic force coefficients
Cx = Cx0 + Cxu*(ur/norm(Vr)) + Cxw*(wr/norm(Vr)) + Cxw2*(wr/norm(Vr))^2;
Cy = Cy0 + Cyv*(vr/norm(Vr)) + Cyp*((b*pr)/(2*norm(V))) + Cyr*((b*rr)/(2*norm(Vr))) ...
+ Cyda*da + Cydr*dr;
Cz = Cz0 + Czw*(wr/norm(Vr)) + Czw2*(wr/norm(Vr))^2 + Czq*((c*qr)/(2*norm(Vr))) ...
+ Czde*de;
Thrust = dT*Power/norm(Vr); % Constant power
X = PDyn*S*Cx + Thrust;
Y = PDyn*S*Cy;
Z = PDyn*S*Cz;
Force_Aero = [X; Y; Z];
% Components of aerodynamic moment (modulo "unsteady" terms)
Cl = Cl0 + Clv*(vr/norm(Vr)) + Clp*((b*pr)/(2*norm(Vr))) + Clr*((b*rr)/(2*norm(Vr))) + ...
Clda*da + Cldr*dr;
Cm = Cm0 + Cmw*(wr/norm(Vr)) + Cmq*((c*qr)/(2*norm(Vr))) + Cmde*de;
Cn = Cn0 + Cnv*(vr/norm(Vr)) + Cnv2*(vr/norm(Vr))^2 + Cnp*((b*pr)/(2*norm(Vr))) + Cnr*((b*rr)/(2*norm(Vr))) + ...
Cnda*da + Cndr*dr;
L = PDyn*S*b*Cl;
M = PDyn*S*c*Cm;
N = PDyn*S*b*Cn;
Moment_Aero = [L; M; N];
% Sum of forces and moments
Force = RIB'*(m*g*e3) + Force_Aero;
Moment = Moment_Aero;
%%% EQUATIONS OF MOTION %%%
%%% KINEMATIC EQUATIONS %%%
XDot = RIB*V;
ThetaDot = LIB*omega;
%%% DYNAMIC EQUATIONS %%%
VDot = (1/m)*(cross(m*V,omega) + Force);
omegaDot = inv(Inertia)*(cross(Inertia*omega,omega) + Moment);
%xdot = [XDot; ThetaDot; VDot; omegaDot];
Include actuator rates in state derivative vector
xdot = [XDot; ThetaDot; VDot; omegaDot; dedot; dadot; drdot];
end
function xhat = hat(x)
xhat = [0,-x(3),x(2);x(3),0,-x(1);-x(2),x(1),0];
end
Cesar Cardenas
on 24 Sep 2022
Paul
on 24 Sep 2022
Why not just copy and paste the code I posted? It runs without error.
I don't understand your questiion about de, da, and dr. They are now state variables in the state vector x and are pulled out of x to use inc ompuations in the exact same way you were already extracting omega, for example.
The error message sounds like it results because you didn't add the variables dedot, dadot, and drdot to the xdot vector at the bottom of GenericFixedWingEOM, as I've now shown twice, and for some reason you commented out in the code you posted in this comment.
Cesar Cardenas
on 24 Sep 2022
Edited: Cesar Cardenas
on 24 Sep 2022
Cesar Cardenas
on 24 Sep 2022
Edited: Cesar Cardenas
on 24 Sep 2022
Cesar Cardenas
on 24 Sep 2022
The call in GenericFixedWingEOM should be to "doublet", not "double." I realize you copied exactly what I typed in that comment (since corrected).
What is your definition of a doublet?
I assumed that it's some sort of actuator command defined as a function of time. For example, step the actuator command to 0.001 for 1 second, down to -0.001 for 1 second, and back to zero.
t = 0:.001:3;
figure
plot(t,doublet(t))
ylim([-1.1 1.1]*1e-3)
xlabel('Time (sec)');
ylabel('Actuator Command')
Perhaps my assumption was incorrect.
function dc = doublet(t)
dc = 0.001*(t < 1) - 0.001*(t >=1 & t < 2) + 0*(t >= 2);
end
More Answers (1)
Try something like this:
tau0 = 2;
tau1 = 3;
Gp = tf(1, [tau1 1],'InputDelay', tau0)
For more info, please check out:
Categories
Find more on Unit Conversions 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!
