フィルターのクリア

rlocfind returns error: Execution of script poly as a function is not supported

6 ビュー (過去 30 日間)
Carlos Torelli
Carlos Torelli 2023 年 11 月 8 日
移動済み: Dyuman Joshi 2023 年 11 月 8 日
rlocfind has never worked for me, code provided works on collegues computers. I have trid unistalling controls toolbox to no avail.
full code provided below:
% Script to generate the full blown and approximate transfer functions for
% Lateral Directional dynamics of an aircraft.
% Housekeeping
clear
close all
clc
% Flight conditions
Alt = 30000;
Mach = .459;
U1 = 456;
qbar = 92.7;
g = 32.2;
theta1 = 2/57.3;
alpha1 = theta1;
% Aircraft Mass, Geometric and Inertial Properties
S = 182;
cbar = 5.47;
b = 33.8;
m = 6360;
Ixx_b = 7985*ones(1,8);
Iyy = 3326*ones(1,8);
Izz_b = 11183*ones(1,8);
Ixz_b = zeros(1,8);
% MoI Transformation (Roskam p346); This is used to transform the Moments of
% Inertia of the aircraft from the body fixed reference frame to the
% Stability axes. Rememebr that the stability axis is also a body fixed
% frame.
I_s = [cos(alpha1(1))^2 sin(alpha1(1))^2 -sin(2*alpha1(1));...
sin(alpha1(1))^2 cos(alpha1(1))^2 sin(2*alpha1(1));...
sin(2*alpha1(1))/2 -sin(2*alpha1(1))/2 cos(2*alpha1(1)) ]*[Ixx_b(1,:);Izz_b(1,:);Ixz_b(1,:)];
Ixx = I_s(1,:);
Izz = I_s(2,:);
Ixz = I_s(3,:);
% Lateral Directional coefficients and stability derivatives
c_l_b = -0.0944;
c_l_p = -0.442;
c_l_r = 0.0926;
c_y_b = -0.346;
c_y_p = -0.0827;
c_y_r = 0.300;
c_n_b = 0.1106;
c_n_Tb = zeros(1,8);
c_n_p = -0.0243;
c_n_r = -0.1390;
c_l_da = 0.1810;
c_l_dr = 0.015;
c_y_da = zeros(1,8);
c_y_dr = 0.2;
c_n_da = -0.0254;
c_n_dr = -0.0365;
% Lateral- Directional Dynamics
% % Generate the Dimensional Derivatives (Equations are in Table 5.7, Roskam p348)
Y_b = qbar(1)*S*c_y_b(1)/m(1);
Y_p = qbar(1)*S*b*c_y_p(1)/(2*m(1)*U1(1));
Y_r = qbar(1)*S*b*c_y_r(1)/(2*m(1)*U1(1));
Y_da = qbar(1)*S*c_y_da(1)/m(1);
Y_dr = qbar(1)*S*c_y_dr(1)/m(1);
L_b = qbar(1)*S*b*c_l_b(1)/Ixx(1);
L_p = qbar(1)*S*b*b*c_l_p(1)/(2*U1(1)*Ixx(1));
L_r = qbar(1)*S*b*b*c_l_r(1)/(2*U1(1)*Ixx(1));
L_da = qbar(1)*S*b*c_l_da(1)/Ixx(1);
L_dr = qbar(1)*S*b*c_l_dr(1)/Ixx(1);
N_b = qbar(1)*S*b*c_n_b(1)/Izz(1);
N_Tb = qbar(1)*S*b*c_n_Tb(1)/Izz(1);
N_p = qbar(1)*S*b*b*c_n_p(1)/(2*U1(1)*Izz(1));
N_r = qbar(1)*S*b*b*c_n_r(1)/(2*U1(1)*Izz(1));
N_da = qbar(1)*S*b*c_n_da(1)/Izz(1);
N_dr = qbar(1)*S*b*c_n_dr(1)/Izz(1);
% Full blown dynamics and corresponding transfer functions
% Generate the 4th order denominator polynomial (Roskam p322)
% Let us define the ratio of the moment and product of inertia according to
% Table 5.8 (Roskan p349)
A1_bar = Ixz/Ixx;
B1_bar = Ixz/Izz;
% Generate the denominator polynomial.
A_lat = U1(1)*(1 - A1_bar*B1_bar);
B_lat = -Y_b*(1 - A1_bar*B1_bar) -U1(1)*(L_p + N_r + A1_bar*N_p +B1_bar*L_r);
C_lat = U1(1)*(L_p*N_r - L_r*N_p) + Y_b*(N_r + L_p + A1_bar*N_p + B1_bar*L_r) ...
-Y_p*(L_b + N_b*A1_bar + N_Tb*A1_bar) +U1(1)*(L_b*B1_bar + N_b + N_Tb) ...
-Y_r*(L_b*B1_bar + N_b + N_Tb);
D_lat = -Y_b*(L_p*N_r - L_r*N_p) + Y_p*(L_b*N_r - N_b*L_r -N_Tb*L_r) ...
-g*cos(theta1(1))*(L_b + N_b*A1_bar + N_Tb*A1_bar) + U1(1)*(L_b*N_p - N_b*L_p - N_Tb*L_p) ...
-Y_r*(L_b*N_p - N_b*L_p -N_Tb*L_p);
E_lat = g*cos(theta1(1))*(L_b*N_r - N_b*L_r -N_Tb*L_r);
Den_lat = [A_lat B_lat C_lat D_lat E_lat];
% We are generating a transfer function with "delta" being a placeholder
% for any lateral directional control surface deflections (aileron or
% rudder)
Y_d = Y_dr;
L_d = L_dr;
N_d = N_dr;
%% Numerator polynomial of beta(s)/delta(s) transfer function
A_b = Y_d*(1 - A1_bar*B1_bar); % Verified.
B_b = -Y_d*(N_r + L_p + A1_bar*N_p + B1_bar*L_r) + Y_p*(L_d + N_d*A1_bar) ...
+Y_r*(L_d*B1_bar + N_d) -U1(1)*(L_d*B1_bar + N_d); % Verified.
C_b = Y_d*(L_p*N_r - N_p*L_r) +Y_p*(N_d*L_r - L_d*N_r) + g*cos(theta1(1))*(L_d + N_d*A1_bar) ...
+ Y_r*(L_d*N_p - N_d*L_p) -U1(1)*(L_d*N_p -N_d*L_p); % Verified.
D_b = g*cos(theta1(1))*(N_d*L_r - L_d*N_r); % Verified
% Generate the numerator polynomial
Num_b = [A_b B_b C_b D_b];
% Generate the beta(s)/delta(s) transfer function
sys_beta = tf(Num_b,Den_lat);
%% Numerator polynomial of phi(s)/delta(s) transfer function
A_phi = U1(1)*(L_d + N_d*A1_bar); % Verified.
B_phi = U1(1)*(N_d*L_r -L_d*N_r) -Y_b*(L_d + N_d*A1_bar) + Y_d*(L_b + N_b*A1_bar + N_Tb*A1_bar); % Verified.
C_phi = -Y_b*(N_d*L_r - L_d*N_r) +Y_d*(L_r*N_b + L_r*N_Tb -N_r*L_b) ...
+ (U1(1) - Y_r)*(N_b*L_d + N_Tb*L_d -L_b*N_d); % Verified.
% Generate the numerator polynomial
Num_phi = [A_phi B_phi C_phi];
% Generate the phi(s)/delta(s) transfer function
sys_phi = tf(Num_phi,Den_lat);
%% Numerator polynomial of the psi(s)/delta(s) or r(s)/delta(s) transfer function
A_psi = U1(1)*(N_d + L_d*B1_bar); % Verified.
B_psi = U1(1)*(L_d*N_p -N_d*L_p) -Y_b*(N_d + L_d*B1_bar) + Y_d*(L_b*B1_bar + N_b + N_Tb); % Verified.
C_psi = -Y_b*(L_d*N_p - N_d*L_p) + Y_p*(N_b*L_d +N_Tb*L_d -L_b*N_d) +Y_d*(L_b*N_p -N_b*L_p -N_Tb*L_p); % Verified.
D_psi = g*cos(theta1(1))*(N_b*L_d + N_Tb*L_d - L_b*N_d);
% Generate the numerator polynomial
Num_psi = [A_psi B_psi C_psi D_psi];
% Generate the psi(s)/delta(s) transfer function
sys_psi = tf(Num_psi, Den_lat); %Note: this is actually the yaw rate transfer function
%% Generate the dutch roll approximation of the dynamics
% Generate the denominator of the Dutch roll dynamics
A_lat_dr = 1;
B_lat_dr = -(N_r + Y_b/U1(1));
C_lat_dr = N_b+ (Y_b*N_r - N_b*Y_r)/U1(1);
% Generate the numerator for the Dutch Roll dynamics (yaw rate t.f)
A_dr = U1(1)*N_d; % Verified.
B_dr = N_b*Y_d - N_d*Y_b; % Verified.
%Dutch roll approximation
Den_dr = U1(1)*[A_lat_dr B_lat_dr C_lat_dr];
Num_dr = [A_dr B_dr];
%Generate the dutch roll, yaw rate transfer function
sys_dr = tf(Num_dr, Den_dr);
%% Approximation of Roll dynamics
Num_p = L_da;
Den_p = [1 -L_p];
% Generate the transfer function
sys_p = tf(Num_p,Den_p);
%% Now that we have the lateral directional dynamics and the transfer functions defined, we can simulate the dynamics in Matlab.
% previously, we have used the inbuilt Matlab function "step" to generate
% the response of the aircraft dynamics. Now, we will use another method to
% simulate the dynamics, using the "lsim" function. I would recommend that
% you take a look at this function by typing "doc lsim" at the command
% prompt.
% Let us create a "doublet maneuver here. Start by generating a time
% vector. This time vector is for 50 seconds, at a 100 Hz. That is the time
% step is 0.01s.
t = 0:0.01:50; % This creates a row vector.
t = t'; % It is ust my preference, but i like the time vector to be a column vector, so the row vector is transposed.
% let us create an input, 'u'.
u(1:100, 1) = 0; % The first 10 seconds of the input is zero.
u(101:110, 1) = 1; % The input is 1, for a duration of 1 second.
u(111:120,1) = -1; % The input is -1, for a duration of 1 second.
u(121:size(t),1) = 0; % The input is 0, for the rest of time.
% let us visualize the input.
figure;
plot(t,u);
% Now that we have defined the input and the time vector, we can simulate
% the dynamics of the system using the "lsim" command. The output of the
% simulation is stored in an output variable. We can use appropriate names
% for the outputs. For example, we can use the terms "sideslip", "p" or "r" for
% the sideslip, roll rate or yaw rate outputs.
sideslip = lsim(sys_beta,u,t); % Simulates the response of the beta transfer function.
p = lsim(sys_phi,u,t);
r = lsim(sys_psi,u,t);
% Plot the responses.
% Note that I have certain preferences
figure;
plot(t, sideslip,'LineWidth',1.5);
grid on;
xlabel('time (s)','FontName','Palatino Linotype', 'FontSize', 14);
ylabel('sideslip, \beta, radians','FontName','Palatino Linotype', 'FontSize', 14)
title('Sideslip response to a rudder doublet','FontName','Palatino Linotype', 'FontSize', 14)
figure; plot(t,p,'LineWidth',1.5)
grid on;
xlabel('time (s)','FontName','Palatino Linotype', 'FontSize', 14);
ylabel('roll rate, p (radians/s)','FontName','Palatino Linotype', 'FontSize', 14)
title('roll rate response to a rudder doublet','FontName','Palatino Linotype', 'FontSize', 14)
figure;plot(t,r,'LineWidth',1.5)
grid on;
xlabel('time (s)','FontName','Palatino Linotype', 'FontSize', 14);
ylabel('yaw rate, r (radians/s)','FontName','Palatino Linotype', 'FontSize', 14)
title('yaw rate response to a rudder doublet','FontName','Palatino Linotype', 'FontSize', 14)
%% You can use the denominator of the lateral directional transfer functions to
% determine the time domain specifications of the dynamics (as you did with
% the longitudinal dynamics)
damp(Den_lat)
%% New Stuff From Class 10/30
Num_rudder = 10;
Den_rudder = [1 10];
G_rudder = tf(Num_rudder,Den_rudder);
%adding Washout filter Two = 4 hence 0.25
Num_WO = [1 0];
Den_WO = [1 0.25];
G_WO = tf(Num_WO, Den_WO)
%Because in cascade we can multiply these two "boxes together" (look at 1.1 in notes for block diagram)
G_forward = G_rudder*sys_psi;
%Gyro Gain
k_r = -1
%Plotting Root Locus for system
figure;rlocus(G_forward*k_r*G_WO);
hold on;
sgrid
%Pick Pole location and always zom before using rlocfind
k = rlocfind(G_forward*k_r)
%Close the loop
cltf = feedback(k*G_forward,k_r);
%Get charecteristics of closed loop system
damp(cltf)
%simulation of closed loop response
figure;step(cltf,50)
%Dutch
G_forward_dr = G_rudder*sys_dr
figure;rlocus(G_forward_dr*k_r)
hold on;
sgrid
k1 = rlocfind(G_forward_dr*k_r)
%close loop
cltf1 = feedback(k1*G_forward_dr,k_r);
%charecteristics of closed loop system
damp(cltf1)
%Simulation of closed loop response
figure;step(cltf1,50)
  1 件のコメント
Carlos Torelli
Carlos Torelli 2023 年 11 月 8 日
and also this:
Error in DynamicSystem/rlocfind (line 59)
num = Gain * [zeros(1,lp-lz) poly(Zero)];

サインインしてコメントする。

採用された回答

Stephen23
Stephen23 2023 年 11 月 8 日
編集済み: Stephen23 2023 年 11 月 8 日
"Execution of script poly as a function is not supported"
You have created a script named POLY. Rename that script.
To find things named POLY:
which poly -all
  2 件のコメント
Dyuman Joshi
Dyuman Joshi 2023 年 11 月 8 日
編集済み: Dyuman Joshi 2023 年 11 月 8 日
Also, @Carlos Torelli, it is best to not use built-in function names as names for variables or scripts.
Carlos Torelli
Carlos Torelli 2023 年 11 月 8 日
移動済み: Dyuman Joshi 2023 年 11 月 8 日
fixed it perfectly!

サインインしてコメントする。

その他の回答 (0 件)

製品


リリース

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by