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

Question

Carlos Torelli 2023-11-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2044407-rlocfind-returns-error-execution-of-script-poly-as-a-function-is-not-supported

移动：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)