Матлаб. К. Ю. Петрова введение в matlab учебное пособие
 Скачать 2.57 Mb.
|
|
Моделирование в SIMULINK Необходимые начальные сведения о работе в SIMULINK были приведены в разд. 3. В этом параграфе описываются некоторые дополнительные возможности SIMULINK по составлению и анализу схем моделирования, а также его взаимодействию с MATLAB. Редактор дифференциальных уравнений DEE Эффективным средством для решения нелинейных дифференциальных уравнений с известными начальными условиями является моделирование в SIMULINK. Оно требует построения эквивалентной структурной схемы и ее реализации на стандартных линейных и нелинейных блоках. Определенную помощь в этом процессе может оказать редактор дифференциальных уравнений DEE (Differential Equation Editor). Этот редактор сам строит схему моделирования для SIMULINK по системе дифференциальных уравнений, записанной в форме 97 Коши. Это немного эффективней, чем "ручное" построение схемы на интеграторах, и несколько проще, чем написание блока «с нуля». Редактор вызывается командой dee. Появляется окно с несколькими блоками, один из них называется Differential Equation Editor, его надо скопировать (перетащить мышкой) на рабочую страницу SIMULINK, войти в него двойным щелчком мыши и ввести параметры дифференциального уравнения (левые части, начальные условия и др.). Приведем пример заполнения параметров для моделирования дифференциального уравнения математического маятника 0, (0) (0) 0, sin y y y y записанного в форме Коши: 1 2 2 1 sin x x x x Name: pendulum # of inputs 0 First order equations, f(x,u): x0 dx/dt= x(2) -sin(x(1)) 1 0 Output equations, f(x,u) y= x(1) x(2) В дальнейшем входить в блок для изменения параметров можно с помощью команды diffeqed. Схема моделирования приведена слева на рис. 5.6. Чтобы одновременно наблюдать графики двух сигналов в блоке Scope в параметрах осциллографа устанавливаем поле ”Number of axes” равным двум. Параметры осциллографа можно редактировать, нажав кнопку . Если есть желание взглянуть непосредственно на схему моделирования, соответствующую заданному уравнению, надо, выделив блок, нажать правую кнопку мыши и выбрать пункт меню look under mask, при этом появится схема, показанная справа на рис. 5.6. Рис. 5.6 Scope Differential Equation Editor DEE 2 Port2 1 Port1 u(2) y2 u(1) y1 f(u) x2 u(2) x1 Mux SysMux 1/s Integ2 1/s Integ1 98 Анализ Simulink-моделей Для получения графиков временных или частотных характеристик схемы, нарисованной в SIMULINK, и построения диаграммы нулей и полюсов ее передаточной функции, используется команда Linear Analysis из пункта Tools меню, расположенного вверху рабочей страницы SIMULINK. Ответ выдается в виде диаграмм средства LTIview. Предварительно на схеме надо пометить вход и выход, пользуясь метками входных и выходных точек Input Point, Output Point из библиотечки Control_System_Toolbox (в ранних версиях – Model_Inputs_and_Outputs) (см. рис. 5.7). Для модели (в общем случае – нелинейной), реализованной в SIMULINK, можно получить описание в пространстве состояний. Это делается при помощи команды linmod. Поясним ее применение на примере схемы следящей системы, показанной на рис. 5.8, которой мы присвоим имя test_sys.mdl. Предварительно на схеме надо пометить вход и выход, используя блоки In и Out. Найдем описание схемы в пространстве состояний при k=1. Набирая в командной строке код [A,B,C,D]=linmod('test_ sys'), получим ответ: A = -1 1 -1 0 B = 0 1 C = 1 0 D = 0 Число выходных параметров команды linmod можно брать любым от 0 до 4. Применяя эту команду для нелинейных систем, получим описание в пространстве состояний линеаризованной системы. В теории управления часто встает задача определения стационарных точек системы, т.е. положений равновесия, при которых все производные равны нулю 0 X . Для Simulink-моделей эту задачу можно решить при помощи команды trim (от trim – приводить в порядок; уравновешивать). Входным параметром функции является имя mdl-файла, а выходными – координаты точки равновесия в фазовом пространстве X, а также установившиеся значения входного и выходного сигналов u и y. Возвращаемое значение dX в случае успеха должно быть равно нулю. Применительно к нашей схеме результатом команды [X,u,y,dX]=trim('test_ sys') будет: X = 0.5 0.5 u = –0.5 y = 0.5000 dX = 0 0 Output Point Input Point Рис . 5.7 1 Out1 1 s+1 Transfer Fcn Step Scope s 1 Integrator k Gain 1 In1 x 1 x 2 y v u Рис. 5.8 99 Действительно, если входной сигнал взять равным u= –0,5, то при у=0,5 и v=1 на выходе сумматора получим нуль, и выходной сигнал интегратора х 2 =0,5 не будет изменяться. Производные от переменных х 1 , х 2 будут равны нулю. Изменим значение k, взяв k=2, и вновь выполним команду [X,u,y,dX]=trim('test_cx'). Результат будет иным: X = 0.3333 0.3333 u = –0.3333 y = 0.3333 dX = 1.0e-015 * 0 0.1110 Теперь стационарное состояние схемы достигается при u= –1/3, х 1 = х 2 =1/3. Для проверки получим те же результаты теоретически. При k = 1 схема описывается уравнением 1 , , , 1 1 2 2 1 1 v x y v u x x x x x Полагая 0 2 x , получаем уравнение x 1 – u = 1. Это одно уравнение с двумя неизвестными x 1 и u, оно имеет бесчисленное множество решений. Команда trim выдает решение с минимальной нормой x 1 = 0,5; u = 0,5 (оно соответствует псевдообращению матрицы [1 –1]). Остальные переменные получаются из равенств –x 1 + x 2 = 0, y = x 1 , что дает x 2 = y = 0,5. При k = 2 схема описывается уравнениями , 1 2 , 1 1 2 2 1 1 x y u x x x x x Полагая 0 2 x , получаем уравнение 2x 1 – u = 1. По команде trim получим одно из его решений x 1 = 1/3, u = –1/3 (заметим, что оно не совпадает с решением x 1 = 2/5, u = –1/5, полученным псевдообращением матрицы [2 1]). Маскирование подсистем в SIMULINK При моделировании сложных систем, когда Simulink-схема становится слишком большой, бывает полезно перейти к иерархической модели, разбив систему на подсистемы. Выделив мышью группу блоков, можно изобразить их в виде одного блока подсистемы, используя пункт меню Edit/Create subsystem. В дальнейшем пользователь может, щелкнув мышью по блоку подсистемы, войти внутрь и редактировать эту подсистему. Иногда такую возможность следует запретить, произведя так называемое «маскирование» системы (Edit/Mask subsystem). Содержимое маскированной системы можно просмотреть по команде меню Edit/Look under mask. По команде Edit/Edit mask открывается редактор масок. Он содержит кнопку «Unmask», позволяющую удалить маску, а также управлять параметрами маскированной подсистемы и процедурой отрисовки блока в редакторе SIMULINK. Пример. Создадим маскированную подсистему, состоящую из двух последовательно соединенных усилителей. В качестве коэффициентов усиления зададим не конкретные числа, а переменные g1 и g2. Отредактируем маску так, чтобы, во-первых, изображение блока содержало окружность и надпись «2 Gains» (рис. 5.9, слева) и, во-вторых, чтобы при щелчке мышью по маскированной подсистеме открывалось диалоговое окно, запрашивающее оба коэффициента усиления (рис. 5.9, справа). 100 2 gains Subsystem Scope 1 Constant Рис. 5.9 Для решения первой задачи в редакторе масок на панели Icon в поле Draw commands следует ввести plot( (1+sin(0:0.2:6.2))/2,(1+cos(0:0.2:6.2))/2) text(0.1,0.5,'2 gains') В списке Units выберем Normalized, это означает, что при отрисовке размер блока считается равным 1x1. Для решения второй задачи на панели Parameters добавим два параметра – g1 и g2 (это имена параметров усилителей) с приглашениями “Gain 1” и “Gain 2”, которые будут отображены на форме. Пример . Моделирование двойного маятника. Консервативная колебательная система из двух взаимосвязанных маятников описывается уравнениями 0 , 0 2 1 2 2 1 1 y y y y ky y Схема моделирования этих уравнений в SIMULINK приведена на рис. 5.10, а. Ее левая часть отвечает первому маятнику, правая – второму. Выделяя поочередно эти части и выполняя их маскирование, получаем схему, изображенную на рис. 5.10, б. Она занимает меньше места и наглядно показывает структуру взаимосвязи маятников. Чтобы войти внутрь подсистемы и изменить ее параметры, нужен двойной щелчок мышью. 101 y1 y2 In1 Out1 Subsystem2 In1 Out1 Subsystem1 Scope2 Рис. 5.10 На рис. 5.11 приведен результат моделирования для начальных условий 1 ) 0 ( , 0 ) 0 ( 2 1 y y и значения k=1.2. 0 5 10 15 20 25 30 35 40 -1 -0.5 0 0.5 1 t y 2 Рис. 5.11 Управление Simulink-моделью из MATLAB При организации численных экспериментов в SIMULINK может возникнуть необходимость многократно запускать одну и ту же модель с различными параметрами и входными сигналами. Это можно обеспечить, управляя процессом моделирования из командного окна. Gain2 Gain1 y1 y2 Scope1 1 s Integrator3 1 s Integrator2 1 s Integrator1 1 s Integrator -1 -K- a) б) 102 Прежде всего, необходимо открыть Simulink-модель вручную или при помощи команды open. Так, команда open('vdp') открывает модель vdp.mdl (схема моделирования уравнения Вандерполя). Чтобы запустить процесс моделирования из командной строки используется команда sim, например, sim('vdp'). В качестве параметров этой команде можно передать время моделирования, входной сигнал (при наличии входных портов у модели), начальные условия и опции решателя. В этом случае синтаксис ее вызова будет иметь вид y =sim('model',timespan,options,ut), где timespan – интервал моделирования, ut – массив, состоящий из отсчетов времени и значений входного сигнала. Пример. Создадим модель, показанную на рис. 5.12 и сохраним ее под именем simple_model. Найдем ее реакцию на синусоиду на интервале 30с. Опций моделирования использовать не будем – в качестве аргумента возьмем пустой массив. 1 Out1 Scope 1 s Integrator 3 Gain 1 In1 Рис. 5.12 >>open('simple_model'); >>t=0:0.1:30; >>sim('simple_model',[0 30],[],[t' sin(t')]) >> [t1,y1]=sim('simple_model',[0 30],[],[t' sin(t')]); >>plot(t1,y1) После выполнения этих команд можно открыть блок Scope в схеме и увидеть результаты моделирования. Поместить результаты моделирования в рабочее пространство и затем посмотреть график с помощью команды plot можно, используя вызов sim с выходными параметрами. При вызове команды sim с опциями решателя структура опций задается командой simset. Для примера изменим тип решателя на ode23: >> [t2,y2]=sim('simple_model',[0 30],simset ('Solver','ode23'),[t' sin(t')]); С полным списком опций моделирования можно ознакомиться, используя команду simget, например, o=simget ('simple_model'). Помимо задания опций моделирования, существует возможность программно изменять параметры блоков. Для их чтения и записи используют команды get_param и set_param. Поскольку названия параметров блоков не всегда совпадают с текстами пояснений в окне параметров, бывает полезна конструкция get_param(gcb, 'objectparameters'). Мышью в окне модели выбирается нужный блок. После этого команда gcb возвращает имя текущего блока. Так, если будет выделен блок Gain , команда gcb вернет строку 'simple_model/Gain'. При помощи той же функции get_param можно получить любой конкретный параметр, например, >> p=get_param('simple_model/Gain', 'Gain'), class(p) p =3 ans = сhar Обратите внимание на то, что тип параметра – строка. Это значит, что если, например, захочется сменить коэффициент усиления с 3 на 4 при помощи функции set_param , то придется набрать не set_param('simple_model/Gain', 'Gain',4), а set_param('simple_model/Gain', 'Gain',num2str(4)) ! Использование этих команд позволяет, в частности, организовать циклическое изменение параметров схемы моделирования (коэффициентов усиления, начальных условий) с помощью оператора MATLAB for. Пример (задача об оптимальных начальных условиях). Движение маятника описывается нелинейным дифференциальным уравнением 103 0 sin 3 , 0 y y y Его начальные условия удовлетворяют условию 2 2 0 2 0 R y y (лежат на окружности радиуса R в фазовой плоскости). Требуется найти значения , , 0 0 y y при которых энергия выходного сигнала будет максимальной либо минимальной 0 2 extremum dt y E Решение. Составим схему моделирования в SIMULINK и назовем ее sys.mdl (рис. 5.13). y 1 Out1 1 s y 1 s x sin Trigonometric Function Scope 0.3 Gain Рис. 5.13 Рассмотрим два способа поиска оптимальных начальных условий. Способ 1 (одномерный линейный поиск). Будем запускать схему из различных начальных условий, вычисляя каждый раз энергию Е и запоминая результат, а затем выберем экстремальные значения. Начальные условия будем формировать в MATLAB по формулам 0 , cos , sin 0 R y R y Программу для управления процессом моделирования оформим в виде отдельного m-файла optic.m. Program optic %Optimal Initial Condition for sys.mdl open('sys'); % вызов модели sys.mdl E=zeros(1,200); R=1; % массив для записи энергий for i=1:200 x0=R*cosd(i); y0=R*sind(i); % очередные начальные условия set_param('sys/x','initial',num2str(x0)); % установка начальных условий в модель set_param('sys/y','initial',num2str(y0)); sim('sys',[0,40],[],[]); % запуск моделирования E(i)=trapz(tout,yout.^2); % вычисление энергии disp(i); % вывод текущего i на экран end plot(E,'LineWidth',2);grid % график зависимости энергии title('ENERGY vs. IC angle'), % от угла на единичной окружности Результат моделирования для R=1 представлен слева на рис. 5.14. Из графика видно, что максимальная энергия выходного сигнала примерно равна 2,4. 104 0 5 10 15 20 25 -1 -0.5 0 0.5 1 t y y min y max Рис. 5.14 Более точные данные получаем, добавляя в программу следующие команды: [EM,iM]=max(E);x0M=cosd(iM);y0M=sind(iM); [Em,im]=min(E);x0m=cosd(im);y0m=sind(im); disp('iM,x0M y0M EM='),[iM,x0M y0M EM], disp('im,x0m y0m Em='),[im,x0m y0m Em], Результат их выполнения: iM,x0M y0M EM = 47 0.6820 0.7314 2.4118 im,x0m y0m Em = 136 -0.7193 0.6947 1.2637 Таким образом, максимальная энергия выходного сигнала, равная Е=2,412, достигается при φ=47 о , т.е. при начальных условиях 731 , 0 , 682 , 0 0 0 y y (маятник отклоняем вправо и толкаем от положения равновесия). Минимальная энергия выходного сигнала, равная Е=1,264, достигается при φ=136 о , т.е. при начальных условиях 695 , 0 , 719 , 0 0 0 y y (маятник отклоняем влево и толкаем к положению равновесия). Соответствующие графики приведены на рис. 5.14, справа. Изменим радиус R окружности начальных условий, взяв его близким к критическому значению, при котором угол у отклонения маятника приближается к . В этом случае маятник начинает «зависать» в верхнем положении, что приводит к резкому увеличению энергии выходного сигнала (при у → энергия Е→ ). На рис. 5.15 показаны графики энергии Е в зависимости от угла φ (слева) и максимального и минимального выходного сигнала у в зависимости от времени (справа) для R=2,29, там же приведены численные значения экстремальных параметров. 0 50 100 150 200 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 ENERGY vs. IC angle 105 0 50 100 150 200 0 10 20 30 40 50 60 ENERGY vs. IC angle 0 5 10 15 20 25 -2 -1 0 1 2 3 4 R=2.29 R=2,29 iM,x0M y0M EM= 41 1.7283 1.5024 51.8763 im,x0m y0m Em= 122 -1.2135 1.9420 6.5953 Рис. 5.15 Способ 2 (поиск экстремума с помощью функции fmincon). В разд. 4.2.1 были описаны функции MATLAB для решения задач на условный и безусловный экстремум. В данном случае надо найти экстремумы функции Е (энергии выходного сигнала), которая зависит от двух аргументов (начальных условий маятника) при наличии ограничения (сумма квадратов начальных условий фиксирована): 0 ) , ( , ) , ( 2 2 0 2 0 0 0 0 0 R y x y x g extr y x f E Для решения таких задач предназначена команда |