Как имитировать выход системы с синусоидальным входом?

Я хочу смоделировать выход определенной зубчатой ​​​​системы, которая у меня есть. То, как выглядит зубчатая передача, не имеет особого значения для задачи, мне удалось получить необходимое дифференциальное уравнение из механической системы. Вот код, который у меня есть

% 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;

Это должно работать при создании графика, показывающего позиции, мои результаты, для ввода Хевисайда/шага. Мой вопрос в том, как бы я сделал это для синусоидального входа. Думаю, у меня должно быть sin(w*t) вместо (t>=0), где w — частота пульса. Тем не менее, я не могу заставить эту работу работать. Любая помощь могла бы быть полезна! :)


person Marx    schedule 02.02.2021    source источник
comment
Я никогда не использовал lsim, но разве недостаточно определить T1 как синусоидальную волну вместо ступенчатой ​​функции? Например, T1 = sin(2*pi*f*t + phi), где f — нужная частота, а phi — начальная фаза.   -  person Luis Mendo    schedule 03.02.2021
comment
lsim на самом деле простая в использовании команда, в основном вы моделируете систему в MATLAB с произвольным количеством входных данных и временем. Вот полезная ссылка об этом mathworks.com/help/control/ref /lti.lsim.html Также да; достаточно определить T1 как синусоидальную волну вместо ступенчатой ​​функции, но по какой-то причине у меня это не работает. Я рекомендую вам использовать мой код и реализовать свое решение, дайте мне отзыв. Я запускаю свой код на R2018b   -  person Marx    schedule 03.02.2021
comment
Что вы имеете в виду, говоря, что это не работает для меня? Вы получаете ошибку? Вывод неверный?   -  person Matteo V    schedule 03.02.2021


Ответы (2)


Тестовый забег

Это выходные результаты, которые я сейчас получаю в MATLAB R2019b. Как предложил комментарий Луиса, я также объявил синусоиду как T1, чтобы она служила входом. В настоящее время не уверен, является ли этот результат ожидаемым результатом. Вывод результатов

Фрагмент кода:

t = linspace(0, 15, 10000)';
f = 0.1; 
phi = 0; 
T1 = sin(2*pi*f*t + phi);

f → Частота синусоидального входа (0,1 Гц в этом примере).
phi → Фазовый сдвиг синусоидального входа/начальная фаза (0 в этом примере).
t → Временной вектор, определяющий выборки синусоиды.
0 → Время начала (0 секунд в этом примере).
15 → Время окончания (в этом примере 15 секунд).
10000 → Количество выборок между временем начала (0 секунд) и временем окончания (15 секунд).


Реализация в скрипте:

% 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;
person MichaelTr7    schedule 03.02.2021
comment
Абсолютно верно! Спасибо! :) - person Marx; 03.02.2021
comment
Нет проблем, рад, что это сработало. - person MichaelTr7; 03.02.2021

Вот решение моей проблемы :)

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);
person Marx    schedule 03.02.2021