3. Исследование автоколебаний в нелинейной системе
Для существования периодических колебаний в замкнутой системе характеристическое уравнение системы должно иметь мнимые корни:
,
где
– амплитуда колебаний.
Коэффициенты гармонической линеаризации для данной нелинейности:
Эту систему можно решать в Matlab с помощью функции fsolve, но этой функции надо задавать начальные значения переменных. При разных начальных значениях получается разный результат. Поэтому будем перебирать начальные значения до тех пор, пока не получим вещественные и положительные корни, которые будут соответствовать устойчивым колебаниям. Для устойчивости найденных периодических решений должно быть выполнено следующее неравенство:
Для предварительной оценки диапазона, в котором следует перебирать начальные значения амплитуды и частоты можно воспользоваться частотным методом и получить значения амплитуды и частоты при заданных параметрах системы.
Частотный метод:
Построим графики этих функций с помощью Matlab:
K=20; D=2; T1=0.5; T2=0.1; T3=2; T4=1; b=1; c=2;
w = 3.5:0.0001:100;
sys1 = tf( [ K*D*T3, K*D ], [ T1*T2*T4, T1*T2+T2*T4+T1*T4, T1+T2+T4, 1, 0 ] );
figure
nyquist( sys1, w );
hold on
a = 0.5:0.5:25;
WN = 4*c/pi./a.*sqrt(1-(b/2)^2./a.^2)-j*2*c*b/pi./a.^2;
plot(-1./WN)
Обозначим
вещественную часть функции
буквой U:
Зависимость амплитуды и частоты колебаний от K2:
Зависимость амплитуды и частоты колебаний от Ti:
Код для расчёта и построения зависимостей в Matlab:
K=20; D=2; T1=0.5; T2=0.1; T3=2; T4=1; b=1; c=2;
%% зависимость от К
K = 1:40;
x0 = [1;1]; % начальные значения
res = zeros(length(K),2); % для записи результатов
for i=1:length(K) % цикл по значениям К
found = false;
for a=1:30 % перебираем начальные значения амплитуды
for w=1:10 % перебираем начальные значения частоты
x0(1)=a; x0(2)=w;
x = fsolve( @(X) [ T1*T2*T4*X(2)^4 - (T1+T2+T4)*X(2)^2 + 2*c*K(i)*D/pi/X(1)*(2*sqrt(1-(b/2)^2/X(1)^2) + b*T3*X(2)/X(1));
2*c*K(i)*D/pi/X(1)*(2*T3*X(2)*sqrt(1-(b/2)^2/X(1)^2)-b/X(1))-(T1*T2+T1*T4+T2*T4)*X(2)^3+X(2) ], x0, optimset('Display','off') );
if ( imag(x(1))==0 && imag(x(2))==0 && x(1)>0 && x(2)>0 ) % если полученные а и w вещественны и > 0
sq = sqrt(x(1)^2-(b/2)^2); % то проверим их на устойчивость
dx_da = 4*c*K(i)*D/pi/x(1)^3*(x(1)^2/sq-2*sq-b*T3*x(2)); % производная уравнения вещественной части по а
dx_dw = 4*T1*T2*T4*x(2)^3-2*(T1+T2+T4)*x(2)+2*K(i)*D*T3*b*c/pi/x(1)^2;
dy_da = 4*c*K(i)*D/pi/x(1)^3*(T3*x(2)*x(1)^2/sq-2*T3*x(2)*sq+b); % производная уравнения мнимой части по а
dy_dw = -3*(T1*T2+T1*T4+T2*T4)*x(2)^2+1+4*c*K(i)*D*T3/pi/x(1)^2*sq;
if ( dx_da*dy_dw > dx_dw*dy_da ) % это условие должно быть выполнено для устойчивости
found = true; % если оно выполнено, значит мы нашли подходящее решение
break % прекращаем перебирать начальные значения
end
end
end
if found
break
end
end
if ~found
disp([ 'Устойчивое решение для К=', num2str(K(i)), ' не найдено.' ])
continue
end
res(i,1) = x(1);
res(i,2) = x(2);
end
figure
subplot(2,1,1)
plot(K,res(:,1))
grid on
title('a(K_2)')
subplot(2,1,2)
plot(K,res(:,2))
grid on
title('w(K_2)')
%% зависимость от Т
T = [ 0.01:0.01:0.5, 0.6:0.1:5];
x0 = [1;1]; % начальные значения
res = zeros(length(T),2); % для записи результатов
for i=1:length(T) % цикл по значениям T
found = false;
for a=1:20 % перебираем начальные значения амплитуды
for w=1:50 % перебираем начальные значения частоты
x0(1)=a; x0(2)=w;
x = fsolve( @(X) [ T1*T2*T(i)*X(2)^4 - (T1+T2+T(i))*X(2)^2 + 2*c*K*D/pi/X(1)*(2*sqrt(1-(b/2)^2/X(1)^2) + b*T3*X(2)/X(1));
2*c*K*D/pi/X(1)*(2*T3*X(2)*sqrt(1-(b/2)^2/X(1)^2)-b/X(1))-(T1*T2+T1*T(i)+T2*T(i))*X(2)^3+X(2) ], x0, optimset('Display','off') );
if ( imag(x(1))==0 && imag(x(2))==0 && x(1)>0 && x(2)>0 )
sq = sqrt(x(1)^2-(b/2)^2);
dx_da = 4*c*K*D/pi/x(1)^3*(x(1)^2/sq-2*sq-b*T3*x(2));
dx_dw = 4*T1*T2*T(i)*x(2)^3-2*(T1+T2+T(i))*x(2)+2*K*D*T3*b*c/pi/x(1)^2;
dy_da = 4*c*K*D/pi/x(1)^3*(T3*x(2)*x(1)^2/sq-2*T3*x(2)*sq+b);
dy_dw = -3*(T1*T2+T1*T(i)+T2*T(i))*x(2)^2+1+4*c*K*D*T3/pi/x(1)^2*sq;
if ( dx_da*dy_dw > dx_dw*dy_da ) % это условие должно быть выполнено для устойчивости
found = true; % если оно выполнено, значит мы нашли подходящее решение
break % прекращаем перебирать начальные значения
end
end
end
if found
break
end
end
if ~found
disp([ 'Устойчивое решение для Ti=', num2str(T(i)), ' не найдено.' ])
continue
end
res(i,1) = x(1);
res(i,2) = x(2);
end
figure
subplot(2,1,1)
plot(T,res(:,1))
grid on
title('a(T_i)')
subplot(2,1,2)
plot(T,res(:,2))
grid on
title('w(T_i)')