如何用正弦波输入模拟系统输出?
How to simulate a system output with a sine wave input?
我想模拟我拥有的某个齿轮系统的输出。齿轮系统的外观对问题来说并不是特别重要,我设法从机械系统中获得了所需的微分方程。
这是我的代码
% parameters
N2 = 90;
N1 = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2 = 2;
D = 8;
K = 5;
J = (N2/N1)^2 * Jn1 + Jn2 + J2;
% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);
% initial state: (position, velocity) [rad; rad/s]
x0 = [0; 0];
% define the time span
t = linspace(0, 15, 10000)';
% define the input step
T1 = zeros(length(t), 1);
T1(t>=0) = 1;
% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);
% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);
% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
ylabel('[rad]');
grid minor;
legend({'$T_1$', '$\theta_1$', '$\theta_2$'}, 'Interpreter', 'latex',...
'location', 'southeast');
hold off;
这应该可以生成一个图表,显示 Heaviside/step 输入的位置,我的输出。我的问题是,我将如何为正弦波输入执行此操作。我认为我应该使用 sin(w*t)
而不是 (t>=0)
,其中 w 是我的脉冲频率。不过,我似乎无法完成这项工作。任何帮助将非常感激! :)
测试运行
这些是我目前在MATLAB R2019b上得到的输出结果。正如 Luis' 的评论所建议的那样,我还声明了一个正弦曲线作为 T1
作为输入。目前不确定这个结果是否是预期的输出。
代码片段:
t = linspace(0, 15, 10000)';
f = 0.1;
phi = 0;
T1 = sin(2*pi*f*t + phi);
f
→ 正弦输入频率(本例中为 0.1Hz)。
phi
→ 正弦波 input/initial 相位的相位偏移(本例中为 0)。
t
→ 时间向量决定正弦曲线的样本。
0
→ 开始时间(本例中为 0 秒)。
15
→ 结束时间(本例中为 15 秒)。
10000
→ 开始时间 (0s) 和结束时间 (15s) 之间的样本数。
脚本中的实现:
% parameters
N2 = 90;
N1 = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2 = 2;
D = 8;
K = 5;
J = (N2/N1)^2 * Jn1 + Jn2 + J2;
% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);
% initial state: (position, velocity) [rad; rad/s]
x0 = [0; 0];
% define the time span
t = linspace(0, 15, 10000)';
% define the input step
T1 = zeros(length(t), 1);
T1(t>=0) = 1;
f = 0.1; %Sinusoid frequency = 0.1Hz%
phi = 0; %Phase = 0%
T1 = sin(2*pi*f*t + phi);
% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);
% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);
% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
ylabel('[rad]');
grid minor;
legend({'$T_1$', '$\theta_1$', '$\theta_2$'}, 'Interpreter', 'latex',...
'location', 'southeast');
hold off;
这是我的问题的解决方案:)
function x = integrate(t, u, x0)
% parameters
N2 = 90;
N1 = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2 = 2;
D = 8;
K = 5;
J = (N2/N1)^2 * Jn1 + Jn2 + J2;
% integrate the differential equation
[t, x] = ode23(@fun, t, x0);
% plot results
figure('color', 'white');
% plot position
yyaxis left;
plot(t, x(:, 1));
ylabel('$x$ [rad]', 'Interpreter', 'latex');
% plot velocity
yyaxis right;
plot(t, x(:, 2));
ylabel('$\dot{x}$ [rad/s]', 'Interpreter', 'latex');
grid minor;
xlabel('$t$ [s]', 'Interpreter', 'latex');
function g = fun(t, x)
g = zeros(2, 1);
g(1) = x(2);
g(2) = (-K/J)*x(1) + (-D/J)*x(2) + (N2/(N1*J)*u(t));
end
end
现在我们可以使用匿名函数了,例如:
t = linspace(0, 120, 10000)';
x0 = [0.1; 0];
x = integrate(t, @(t)(sin(1.5*t)), x0);
我想模拟我拥有的某个齿轮系统的输出。齿轮系统的外观对问题来说并不是特别重要,我设法从机械系统中获得了所需的微分方程。 这是我的代码
% parameters
N2 = 90;
N1 = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2 = 2;
D = 8;
K = 5;
J = (N2/N1)^2 * Jn1 + Jn2 + J2;
% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);
% initial state: (position, velocity) [rad; rad/s]
x0 = [0; 0];
% define the time span
t = linspace(0, 15, 10000)';
% define the input step
T1 = zeros(length(t), 1);
T1(t>=0) = 1;
% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);
% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);
% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
ylabel('[rad]');
grid minor;
legend({'$T_1$', '$\theta_1$', '$\theta_2$'}, 'Interpreter', 'latex',...
'location', 'southeast');
hold off;
这应该可以生成一个图表,显示 Heaviside/step 输入的位置,我的输出。我的问题是,我将如何为正弦波输入执行此操作。我认为我应该使用 sin(w*t)
而不是 (t>=0)
,其中 w 是我的脉冲频率。不过,我似乎无法完成这项工作。任何帮助将非常感激! :)
测试运行
这些是我目前在MATLAB R2019b上得到的输出结果。正如 Luis' 的评论所建议的那样,我还声明了一个正弦曲线作为 T1
作为输入。目前不确定这个结果是否是预期的输出。
代码片段:
t = linspace(0, 15, 10000)';
f = 0.1;
phi = 0;
T1 = sin(2*pi*f*t + phi);
f
→ 正弦输入频率(本例中为 0.1Hz)。
phi
→ 正弦波 input/initial 相位的相位偏移(本例中为 0)。
t
→ 时间向量决定正弦曲线的样本。
0
→ 开始时间(本例中为 0 秒)。
15
→ 结束时间(本例中为 15 秒)。
10000
→ 开始时间 (0s) 和结束时间 (15s) 之间的样本数。
脚本中的实现:
% parameters
N2 = 90;
N1 = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2 = 2;
D = 8;
K = 5;
J = (N2/N1)^2 * Jn1 + Jn2 + J2;
% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);
% initial state: (position, velocity) [rad; rad/s]
x0 = [0; 0];
% define the time span
t = linspace(0, 15, 10000)';
% define the input step
T1 = zeros(length(t), 1);
T1(t>=0) = 1;
f = 0.1; %Sinusoid frequency = 0.1Hz%
phi = 0; %Phase = 0%
T1 = sin(2*pi*f*t + phi);
% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);
% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);
% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
ylabel('[rad]');
grid minor;
legend({'$T_1$', '$\theta_1$', '$\theta_2$'}, 'Interpreter', 'latex',...
'location', 'southeast');
hold off;
这是我的问题的解决方案:)
function x = integrate(t, u, x0)
% parameters
N2 = 90;
N1 = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2 = 2;
D = 8;
K = 5;
J = (N2/N1)^2 * Jn1 + Jn2 + J2;
% integrate the differential equation
[t, x] = ode23(@fun, t, x0);
% plot results
figure('color', 'white');
% plot position
yyaxis left;
plot(t, x(:, 1));
ylabel('$x$ [rad]', 'Interpreter', 'latex');
% plot velocity
yyaxis right;
plot(t, x(:, 2));
ylabel('$\dot{x}$ [rad/s]', 'Interpreter', 'latex');
grid minor;
xlabel('$t$ [s]', 'Interpreter', 'latex');
function g = fun(t, x)
g = zeros(2, 1);
g(1) = x(2);
g(2) = (-K/J)*x(1) + (-D/J)*x(2) + (N2/(N1*J)*u(t));
end
end
现在我们可以使用匿名函数了,例如:
t = linspace(0, 120, 10000)';
x0 = [0.1; 0];
x = integrate(t, @(t)(sin(1.5*t)), x0);