Лабораторная работа: Исследование нелинейных систем
Проведите качественное исследование системы: оцените возможность возникновения периодических режимов и их устойчивость. Там, где это возможно, запишите выражения для критического коэффициента усиления ЛЧ, при котором в нелинейной системе возникают автоколебания.
Изучите текст программы (файл GB_prog.m) и структуру моделей (файлы GB_mod.mdl и R_Fourie.mdl).
Графический расчет параметров периодических режимов c использованием метода гармонического баланса производится по сценарию, записанному в файле GB_prog.m. Моделирование нелинейной системы осуществляется при помощи файлов GB_prog.m и GB_mod.mdl, а анализ спектрального состава периодического режима на выходе линейной части – при помощи файлов GB_prog.m и R_Fourie.mdl.
Cодержание файла GB_prog.m:
%Исследование нелинейных систем методом гармонического баланса
%Используемые файлы: GB_prog.m, GB_mod.mdl и R_Fourie.mdl.
%Используемые обозначениЯ: НЭ – нелинейный элемент, ЛЧ – линейнаЯ часть.
%Очистка всех переменных в памЯти
clear all
%закрытие всех предыдущих рисунков
set (0,'ShowHiddenHandles', 'on')
delete (get(0,'Children'))
%Задание значениЯ переменной, определЯющей положение переключателЯ
%конфигурации нелинейной части в файле GB_mod.mdl
%1 – НЭ-двухпозиционное реле с гистерезисом
%2 – НЭ-трехпозиционное реле без гистерезиса
%3 – НЭ-люфт
config_nlin = 2;
%Задание значениЯ переменной, определЯющей положение переключателЯ
%конфигурации линейной части в файле GB_mod.mdl
%1 – ЛЧ – Wл1 (p)=k/[(T1*p+1) (T2*p+1) (T3p+1)]
%2 – ЛЧ – Wл2 (p)=k/[(T1*p+1) (T2*p+1) p]
%3 – ЛЧ – Wл3 (p)=[k (T1*p+1)]/[(T2*p‑1)^2 (T3*p+1)^2]
%4 – ЛЧ – Wл4 (p)=[k (T1*p+1)]/[(T2*p‑1) p]
config_lin = 2;
k = 5;
T1 = 0.1;
T2 = 0.2;
T3 = 1;
b = 0.1;
c = 1;
%Графический расчет параметров периодических режимов c использованием
%метода гармонического баланса
%описание линейной части
switch config_lin
case 1,
%начальное значение, шаг и конечное значение частоты (в рад)
w = [0.02:0.01:100];
%комплексный передаточный коэффициент ЛЧ
W_lin = k./ ((T1*j*w+1).* (T2*j*w+1).*(T3*j*w+1));
%определение значениЯ строковой переменной lin
lin = 'Wл1 (p)';
case 2,
%начальное значение, шаг и конечное значение частоты (в рад)
w = [2:0.01:100];
%комплексный передаточный коэффициент ЛЧ
W_lin = k./ ((T1*j*w+1).* (T2*j*w+1).*(j*w));
%определение значениЯ строковой переменной lin
lin = 'Wл2 (p)';
case 3,
%начальное значение, шаг и конечное значение частоты (в рад)
w = [0.01:0.01:300];
%комплексный передаточный коэффициент ЛЧ
W_lin = (k*(T1*j*w+1))./ ((T2*j*w‑1).^2.*(T3*j*w+1).^2);
%определение значениЯ строковой переменной lin
lin = 'Wл3 (p)';
case 4,
%начальное значение, шаг и конечное значение частоты (в рад)
w = [1:0.01:100];
%комплексный передаточный коэффициент ЛЧ
W_lin = (k*(T1*j*w+1))./ ((T2*j*w‑1).*(j*w));
%определение значениЯ строковой переменной lin
lin = 'Wл4 (p)';
end
%описание ЛЧ как частотной характеристики через W_lin(jw)
SYSL = frd (W_lin, w);
%–
%Описание нелинейной части
switch config_nlin
case 1,
%начальное значение, шаг и конечное значение амплитуды на входе НЭ
A = [b:0.01:3.0];
%коэффициенты гармонической линеаризации
q = 2*c/(pi*b).* ((2*b./A).* sqrt (1 – (b./A).^2));
q1 = -2*c/(pi*b).* 2*(b./A).^2;
%определение значениЯ строковой переменной nlin
nlin = '2‑х поз. реле с гист.';
case 2,
%начальное значение, шаг и конечное значение амплитуды на входе НЭ
A = [b+0.0001:0.005:b*sqrt(2)];
A_ = [b*sqrt(2):0.005:3.0];
%коэффициенты гармонической линеаризации
q = (4*c./ (pi*A)).* sqrt (1 – (b./A).^2);
q1 = 0;
q_ = (4*c./ (pi*A_)).* sqrt (1 – (b./A_).^2);
q1_ = 0;
%эквивалентный комплексный передаточный коэффициент W_nlin(jA) НЭ
W_nlin_ = q_ + j*q1_;
%характеристика -1/W_nlin(jA)
S_ = -1./W_nlin_;
%описание НЭ как амплитудной характеристики через -1/W_nlin(jA)
SYSN_ = frd (S_, A_);
%определение значениЯ строковой переменной nlin
nlin = '3‑х поз. реле без гист.';
case 3,
%начальное значение, шаг и конечное значение амплитуды на входе НЭ
A = [b+0.001:0.01:100];
%коэффициенты гармонической линеаризации
alfa = asin (1–2*b./A);
q = (1/pi)*(pi/2+alfa+0.5*sin (2*alfa));
q1 = – (4*b).* (1‑b./A)./ (pi*A);
%определение значениЯ строковой переменной nlin
nlin = 'люфт';
end
%эквивалентный комплексный передаточный коэффициент W_nlin(jA) НЭ
W_nlin = q + j*q1;
%характеристика -1/W_nlin(jA)
S = -1./W_nlin;
%описание НЭ как амплитудной характеристики через -1/W_nlin(jA)
SYSN = frd (S, A);
%–
%ВизуализациЯ
%Построение W_lin(jw) и -1/W_nlin(jA) при помощи plot (общий вид)
figure(1)
gr_W_lin = plot (real(W_lin), imag (W_lin));
set (gr_W_lin, {'Color'}, {'r'});
hold on
gr_S = plot (real(S), imag(S));
set (gr_S, {'Color'}, {'b'});
title(['Расчет гарм. баланса. НЭ – ', nlin, ', b=', num2str(b), ', c=', num2str(c), '; ЛЧ – ',…
lin, ', k=', num2str(k), ', T1=', num2str(T1), ', T2=', num2str(T2),…
', T3=', num2str(T3)])
xlabel ('re(W lin), re(S)');
ylabel ('im(W lin), im(S)');
legend ('W lin(jw)', 'S(jA)', 0);
grid on
%–
%построение W_lin(jw) и -1/W_nlin(jA) при помощи LTI Viewer
if config_nlin == 2
ltiview({'nyquist'}, SYSL, '-b', SYSN, '-r')
ltiview({'nyquist'}, SYSL, '-b', SYSN_, '-r')
else
ltiview({'nyquist'}, SYSL, '-b', SYSN, '-r')
end
%–
%Завершение графического расчета параметров периодических режимов
%c использованием метода гармонического баланса
%====================================================
%^^^^^^^^^^^^^^Моделирование нелинейной системы (НС)^^^^^^^^^^^^^^^^^^^^^
% с определением параметров автоколебаний
%задание времени моделированиЯ
t_end = 15;
%ограничение шага моделированиЯ сверху
step_max = 0.005;
%задание амплитуды сигнала на выходе ЛЧ в начальный момент времени
y0 = 0.5;
%вызов модели
open_system ('GB_mod.mdl');
%запуск модели
sim ('GB_mod');
%считывание фактических параметров моделированиЯ
Max_Step_Size = get_param ('GB_mod', 'MaxStep');
Stop_Time = get_param ('GB_mod', 'StopTime');
%извлечение последнего элемента из векторов vec_period и vec_amp
%(определение установившихсЯ значений периода и амплитуды автоколебаний
%с использованием данных, записанных в рабочую область памЯти)
clc % очистка командного окна
period = vec_period (length(vec_period));
clear vec_period;
amp_kol = vec_amp (length(vec_amp))%вывод в командное окно
clear vec_amp;
%расчет частоты сигнала на выходе ЛЧ (используетсЯ длЯ последующего
%расчета амплитуды гармоник)
frequency = 1/period;
w_kol = 2*pi*frequency %вывод в командное окно
%Построение процесса во времени на выходе ЛЧ
figure(2)
gr = plot (t_and_y(:, 1), t_and_y(:, 2));
set (gr, {'Color'}, {'r'});
title(['Процесс y(t). НЭ – ', nlin, ', b=', num2str(b), ', c=', num2str(c), '; ЛЧ – ',…
lin, ', k=', num2str(k), ', T1=', num2str(T1), ', T2=', num2str(T2),…
', T3=', num2str(T3), ', y0=', num2str(y0)])
xlabel(['t, cек Амп.кол.=', num2str (amp_kol), ', Wкол=', num2str (w_kol), 'c^-1']);
ylabel('y');
grid on
%–
%закрытие модели
%close_system ('GB_mod', 1)
%=====================================================
%^^^^^^^^^^^^^^^^Анализ спектра автоколебательного процесса^^^^^^^^^^^^^^
%вызов программы, определЯющей амплитуды гармоник сигнала на выходе ЛЧ
open_system ('R_Fourie.mdl');
%транслЯциЯ параметров моделированиЯ НС 'GB_mod.mdl' в файл R_Fourie.mdl
set_param ('R_Fourie', 'MaxStep', Max_Step_Size);
set_param ('R_Fourie', 'StopTime', Stop_Time);
%запуск программы, вычислЯющей амплитуды гармоник
sim ('R_Fourie');
%закрытие программы
%close_system ('R_Fourie', 1)
%извлечение последнего элемента из каждого вектора, содержащего амплитуды
%гармоник (определение установившихсЯ значений амплитуд гармоник с
%использованием данных, записанных в рабочую область памЯти)
magn_0=vec_magn_0 (length(vec_magn_0));
clear vec_magn_0;
magn_1=vec_magn_1 (length(vec_magn_1));
clear vec_magn_1;
magn_2=vec_magn_2 (length(vec_magn_2));
clear vec_magn_2;
magn_3=vec_magn_3 (length(vec_magn_3));
clear vec_magn_3;
magn_4=vec_magn_4 (length(vec_magn_4));
clear vec_magn_4;
magn_5=vec_magn_5 (length(vec_magn_5));
clear vec_magn_5;
%Проверка гипотезы фильтра
filtration = magn_1/magn_3
%построение гистограммы
figure(3);
bar([0 1 2 3 4 5], [magn_0 magn_1 magn_2 magn_3 magn_4 magn_5]);
grid on
title(['Гарм. состав y(t). НЭ – ', nlin, ', b=', num2str(b),…
', c=', num2str(c), '; ЛЧ – ', lin, ', k=', num2str(k), ', T1=', num2str(T1),…
', T2=', num2str(T2), ', T3=', num2str(T3)])
xlabel(['Номер гармоники ГФ: A1/A3=', num2str(filtration)]);
ylabel ('Амплитуда гармоники');
Все команды Matlab, использованные при составлении данной программы, описаны в приложении.
В m‑файле задаются значения параметров линейной части и нелинейного элемента, а также начальные условия по выходу линейной части, указывается время моделирования и шаг моделирования.
Для проведения графического расчета параметров периодических
режимов m‑файл содержит описания ЛЧ (через комплексный передаточный
коэффициент) и нелинейности (через коэффициенты гармонической линеаризации и
эквивалентный комплексный передаточный коэффициент). Совместное воспроизведение
характеристик Wл(jw) и 
осуществляется
двумя способами: при помощи команды plot и команды ltiview. Первая позволяет
совместно вывести на рисунок указанные характеристики без вывода
комплексно-сопряженных характеристик Wл(-jw) и 
, но не позволяет по точке
их пересечения определить частоту и амплитуду периодических режимов. При помощи
второй одновременно выводятся все четыре характеристики, но имеется возможность
определения частоты и амплитуды по точке пересечения: команда zoom контекстного
меню, вызываемого правой клавишей мыши, позволяет увеличить нужную область
рисунка, после чего нужно поместить указатель мыши на требуемую характеристику
и нажать на ее левую клавишу – на графике появятся отметка и информационное
окно, которые можно переместить мышью в точку пересечения. Для случая системы с
идеальным трехпозиционным реле предусмотрено раздельное построение
характеристики 
до и после точки
экстремума при помощи двух команд ltiview.
После завершения графического расчета производится вызов и запуск модели нелинейной системы (GB_mod.mdl). По окончании моделирования строится временной процесс y(t) на выходе линейной части, определяются частота и амплитуда автоколебаний (их значения можно увидеть в командном окне и на созданном рисунке).
На этапе анализа спектрального состава периодического режима на выходе линейной части производится вызов и запуск модели (R_Fourie.mdl), определяющей амплитуды гармоник сигнала y(t), записанного в рабочую область памяти, и строится результирующая гистограмма «Номер гармоники – амплитуда гармоники». В командное окно и окно с гистограммой выводится отношение амплитуд первой и третьей гармоник, что позволяет сделать заключение о выполнении / невыполнении гипотезы фильтра.
При составлении моделей в Simulink (GB_mod.mdl и R_Fourie.mdl) используются элементы библиотек Simulink (Math, Linear, Nonlinear, Signals & Systems, Sinks и Sources), Simulink Extras (Additional Linear) и Power System Blockset (Extra Library\Measurements), доступные через Simulink Library Browser (рис. 3.3 – 3.6). В модели нелинейной системы GB_mod.mdl, показанной на рис. 3.3, предусмотрены переключатели конфигурации линейной и нелинейной частей, управление ими осуществляется соответственно через переменные config_lin и config_nlin, значения которых задаются в m‑файле GB_prog.m.
Параметры моделирования должны быть указаны в окне Simulation parameters, доступном через меню Simulation\Simulation parameters окна, в котором открыт mdl‑файл (рис. 3.7, 3.8).
Установку параметров различных функциональных блоков моделей GB_mod.mdl и R_Fourie.mdl поясняют рис. 3.9 – 3.15.
Рис. 3.4. Организация трехпозиционного реле без гистерезиса
(файл GB_mod.mdl)
Рис. 3.5. Подсистема определения периода и амплитуды автоколебаний (файл GB_mod.mdl)
Рис. 3.6. Модель в Simulink, вычисляющая амплитуду гармоник в составе периодического процесса на выходе линейной части (файл R_Fourie.mdl)
Рис. 3.7. Параметры моделирования для файла GB_mod.mdl
Рис. 3.8. Параметры моделирования для файла R_Fourie.mdl
а б
Рис. 3.9. Параметры блоков в составе трехпозиционного реле без гистерезиса (файл GB_mod.mdl): а – блока Relay1; б – блока Relay2
а б
Рис. 3.10. Параметры нелинейностей (файл GB_mod.mdl): а люфта; б – двухпозиционного реле с гистерезисом
а б
Рис. 3.11. Параметры блоков линейной части (файл GB_mod.mdl): а – инерционного звена; б – интегрирующего звена
Рис. 3.12. Параметры блока To Workspace (файл GB_mod.mdl)
а б
Рис. 3.13. Параметры блоков подсистемы определения периода и амплитуды автоколебаний (файл GB_mod.mdl): а – блока ограничения сигнала Saturation2; б – блока памяти Memory1
Рис. 3.14. Параметры блока Fourier (файл R_Fourie.mdl)
Рис. 3.15. Параметры блока From Workspace (файл R_Fourie.mdl)
Выполнение работы
Модернизируйте файл сценария (m‑файл) в соответствии со своим вариантом: установите значения переменных конфигурации, т.е. выберите передаточную функцию линейной части и первую из двух нелинейностей вашего варианта; установите заданные значения параметров ЛЧ и нелинейности.
Запустите m‑файл на исполнение. Зафиксируйте совместно воспроизведенные командой plot частотную характеристику линейной части Wл(-jw) и характеристику нелинейности –1/Wн(jA). По точкам их пересечения определите Ап и wп, используя LTI Viewer.
