how to code transfer function

35 views (last 30 days)
Cesar Cardenas
Cesar Cardenas on 19 Sep 2022
Commented: Paul on 26 Sep 2022
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
Show 8 older comments
Cesar Cardenas
Cesar Cardenas on 22 Sep 2022
@Star Strider it is running and plotting now. thanks much.

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Paul
Paul on 22 Sep 2022
Edited: Paul on 22 Sep 2022
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
Show 16 older comments
Cesar Cardenas
Cesar Cardenas on 24 Sep 2022
Hi Paul, sorry for bothering so much, This is what I'm trying to do for the doublet based on your comment
Not really sure if it is the right approach. Thanks for your feeddback in advance.
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 ...
%removed de0 da0 dr0
% 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);
%Actuator deflections taken from 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;
%added
de0 = double(t:tau0:tau_1);
da0 = double(t:tau0:tau_1);
dr0 = double(t:tau0:tau_1);
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];
xdot = [XDot; ThetaDot; VDot; omegaDot; dedot; dadot; drdot];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dc = doublet(t,u)
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 ...
%removed de0 da0 dr0
%Actuator deflections are defined
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 = 0.211; % Mass (slug)
W = m*g; % Weight (Newtons)
Ix = 0.2163; % Roll inertia (slug-ft^2)
Iy = 0.1823; % Pitch inertia (slug-ft^2)
Iz = 0.3396; % Yaw inertia (slug-ft^2)
Ixz = 0.0364; % Roll/Yaw product of inertia (slug-ft^2)
Inertia = [Ix, 0, -Ixz; 0, Iy, 0; -Ixz, 0, Iz];
S = 4.92; % Wing area (ft^2)
b = 5.91; % Wing span (ft)
c = 0.833; % Wing chord (ft)
AR = (b^2)/S; % Aspect ratio
% Longitudinal nondimensional stability and control derivatives
Cx0 = 0;
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) 45ft/s converted to (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;
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 = 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,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')
figure
plot(t,y(:,13:15))
legend('de','da','dr')
Paul
Paul on 26 Sep 2022
Ran in:
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

Sign in to comment.

More Answers (1)

Sam Chak
Sam Chak on 19 Sep 2022
Ran in:
@Cesar Cardenas,
Try something like this:
tau0 = 2;
tau1 = 3;
Gp = tf(1, [tau1 1],'InputDelay', tau0)
Gp = 1 exp(-2*s) * ------- 3 s + 1 Continuous-time transfer function.
For more info, please check out:
https://www.mathworks.com/help/control/ug/time-delays-in-linear-systems.html

Categories

Find more on Environmental Models in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!