Матлаб. К. Ю. Петрова введение в matlab учебное пособие
 Скачать 2.57 Mb.
|
|
fmincon. Для ее применения нужно создать две вспомогательные функции, содержащие описание минимизируемого критерия и ограничений. Оформим эти функции в виде m-файлов с именами energy и nlcon : function energy function nlcon function E = energy (x) %open('sys'); x0=x(1); y0=x(2); set_param('sys_old/x','initial',num2str(x0,15)); set_param('sys_old/y','initial',num2str(y0,15)); [tout,xx,yout]=sim('sys_old',[0,40],[],[]); E=trapz(tout,yout.^2); function [h,g] = nlcon(x) h = []; R=1; %R=2.29; g = x(1)^2+x(2)^2-R^2; Запускаем процесс оптимизации, набирая в рабочем окне команду fmincon : >> [X,E]=fmincon('energy',[1;1],[],[],[],[],[],[],'nlcon') По окончании итерационного процесса получаем сообщение об его успешном завершении и значения выходных параметров (вектор начальных условий и минимальную энергию для R=1). Чтобы решить задачу о максимальной энергии, нужно изменить знак перед командой trapz в последней строке функции energy и повторить вычисления. В результате получаем: R=1 norm(X) E=trapz(tout,yout.^2) R=1 norm(X) E=-trapz(tout,yout.^2) X = 0.7198 -0.6942 1 E =1.2637 X = 0.6844 0.7291 1 E = -2.4118 Найденные значения близки к полученным первым способом, но более точные. 106 Выполняя аналогичные вычисления для R=2.29, также получаем уточненные значения оптимальных начальных условий: R=2.29 norm(X) E=-trapz(tout,yout.^2) R=2.29 norm(X) E=trapz(tout,yout.^2) X = 1.7201 1.5118 2.2900 E = -51.8823 X = 1.2100 -1.9442 2.2900 E = 6.5955 Обратим внимание на два нюанса в оформлении функции energy. Если во второй строке убрать знак комментария (%), то вызов модели sys.mdl и перерисовка схемы моделирования будет происходить на каждом шаге итерационного процесса, что сильно увеличит время решения. Проще открыть эту модель вручную до начала решения, тогда команда open('sys') не потребуется. Второй нюанс касается синтаксиса команды num2str(x0,15) . Ее второй аргумент задает количество разрядов числа x0 , преобразуемых в строку. Если его опустить, то по умолчанию оно может оказаться равным 4, что полностью исказит результаты оптимизации (авторы столкнулись с этим явлением и долго не могли понять, в чем дело). С подобными тонкостями нередко приходится встречаться при моделировании в MATLAB и других пакетах, они требуют от пользователя хорошего понимания сущности решаемой задачи и применяемого математического аппарата с одной стороны, и знания языка программирования и понимания процесса компьютерной реализации используемого алгоритма, с другой стороны. Задачи и упражнения 1. Нелинейные осцилляторы. Выполнить моделирование в MATLAB приводимых ниже нелинейных дифференциальных уравнений и определить области их устойчивости в зависимости от начальных условий и значений параметров: а) Уравнение математического маятника 0 sin x a x ; б) Уравнение Ван-дер-Поля 0 ) 1 ( 2 x x x x ; в) Уравнение Дюффинга 0 3 2 x x x Указание: В случае уравнения Дюффинга при 0 – решения колебательные. При 0 и при 0 x происходит потеря устойчивости. 2. Цифровая модель. Построить теоретически и в пакете MATLAB цифровую модель непрерывной системы второго порядка ; 4 2 3 u y y y начальные условия – нулевые, шаг h = 0,1. Получить и сравнить графики переходных функций исходной системы и цифровой модели. Теоретическое решение. Корни характеристического уравнения непрерывной системы 1 , 0 ; 2 ; 1 2 1 h p p Корни характеристического уравнения цифровой модели находим по формуле : h p i i e z 2 , 0 2 1 , 0 1 ; e z e z Характеристическое уравнение цифровой модели 0 741 , 0 724 , 1 ) ( ) )( ( 2 3 , 0 2 , 0 1 , 0 2 2 1 z z e z e e z z z z z Ему соответствует разностное уравнение 741 , 0 724 , 1 2 2 1 n n n n x y y y Для определения коэффициента приравниваем установившиеся реакции непрерывной системы и цифровой модели на единичный входной сигнал х = 1. Для непрерывной системы ; 2 уст y для цифровой модели 017 , 0 741 , 0 724 , 1 1 уст n y , откуда 034 , 0 107 Итак, цифровая модель имеет вид 2 2 1 034 , 0 741 , 0 724 , 1 n n n n x y y y Переходные функции исходной системы и цифровой модели: 2 2 4 ) ( 2 t t e e t y , 2 ) 819 , 0 ( 2 ) 905 , 0 ( 4 n n n y Решение в MATLAB. Приведем варианты решения с помощью команд c2d, impinvar и bilinear: >> sys= tf(4,[1 3 2]) >> sysd = c2d(sys,0.1,'imp') >> [bz,az] = impinvar(4,[1 3 2],10) >> [numd,dend] = bilinear(4,[1 3 2],10) 4 ------------- s^2 + 3 s + 2 0.3444 z ---------------------- z^2 - 1.724 z + 0.7408 bz = 0 0.0344 az =1.0000 -1.7236 0.7408 numd = 0.0087 0.0173 0.0087 dend = 1.0000 -1.7229 0.7403 Видно, что с теоретическим решением совпадает результат команды impinvar. 3. Газопровод. Модель газопровода, условно разбитого на шесть участков, описывается системой дифференциальных уравнений u X X 0 0 0 0 0 1 1 1 0 0 0 0 1 2 1 0 0 0 0 1 2 1 0 0 0 0 1 2 1 0 0 0 0 0 2 1 0 0 0 0 1 2 , 1 0 0 0 0 0 0 0 0 0 0 1 X Y Здесь u – давление газа на входе, x i – давление на отдельных участках газопровода. Вектор Y описывает места установки измерительных датчиков. Требуется: а) выполнить моделирование системы в SIMULINK при скачкообразном и синусоидальном изменении входного сигнала; б) найти модальное, сопровождающее и сбалансированное представления это системы и сравнить графики их весовых функций. 4. Ограниченная задача двух тел. Движение спутника вокруг Земли в плоскости орбиты характеризуется дифференциальными уравнениями: 3 2 2 ) ( y x kx x , 3 2 2 ) ( y x ky y Гравитационная постоянная k=3,9870 5 10 2 3 c km Требуется составить схему и выполнить моделирование в SIMULINK при начальных условиях 388 , 6478 ; 0 0 0 y x км, 89 , 7 0 x км/c, 0 0 y . Получить графики координат спутника, его траектории и высоты относительно поверхности Земли , 3 2 2 2 R z y x h где радиус Земли 3 R =6378,388 км. 5. Управление луноходом. Луноход-1 (1970-71г.), имел вес 700кг, 2 скорости 1км/ч и 2км/ч, прошел 11 км за 10 месяцев. Управление им осуществлял экипаж, находящийся на Земле, задержка сигнала в контуре управления составляла порядка трех секунд. Упрощенная математическая модель системы управления характеризуется уравнением: зад y a T t ay t y ) ( ) ( 108 Этому уравнению при a=1 соответствует структурная схема на интеграторе и элементе задержки ЭЗ, показанная на рис.5.16. Требуется выполнить ее моделирование в SIMULINK при 1 зад y , получить графики y(t) и ( ) y f y , определить период колебаний. Найти значение а, при котором решение носит периодический характер. Построить область устойчивости на плоскости (а, Т). y зад u y – 1 Э.З. 1 Рис. 5.16. 109 16.06.2005 Файл mtl_M5_checked ПРОГРАММИРОВАНИЕ В MATLAB Типы данных MATLAB поддерживает большое количество типов данных, которые, по сути, являются классами. Для того чтобы иметь общее представление о возможностях пакета, на рис. 5.1 приведено полное дерево классов пакета. Однако более чем в 90% приложений пользователь имеет дело лишь с классами double, cell и пользовательскими классами. Все классы MATLAB являются наследниками массива (array). Таким образом, MATLAB является по своей структуре типичным векторным процессором. NUMERIC численный char символьный logical логический Java class structure структура function handle дескриптор функции ARRAY массив cell массив ячеек int8, uint8, int16, uint16, int32, int32, int64, uint64 single одинарной точности double двойной точности user class пользовательские классы Рис. 6.1 Класс double array является, в некотором смысле, «классом по умолчанию». Даже если ввести в командном окне MATLAB скалярную переменную, то она будет восприниматься как массив размером 1x1: >> c=3; size (c) ans =1 1 110 Еще один важный тип – «массив ячеек» (cell array). Элементами подобного массива может быть что угодно: числа, строки, массивы, другие массивы ячеек, экземпляры класса и т.п. Задается он при помощи фигурных скобок { }: >> c={1,2 3; 'scrambled eggs',[12, 13, 14], -0.3} c = [ 1] [ 2] [ 3] 'scrambled eggs' [1x3 double] [-0.3000] Употребление точки с запятой и запятой такое же, как и при задании массивов. Чтение и запись отдельных элементов производится так же, как и в случае с массивами типа double, однако вместо круглых используются фигурные скобки: c{1,2} c{2,1} c{2,2}=18 ans = 2 ans = scrambled eggs c = [ 1] [ 2] [ 3] 'scrambled eggs' [18] [-0.3000] Команды size и length в случае с массивами ячеек ведут себя так же, как и с обычными: >> size(c) ans = 2 3 >> length(c) ans = 3 >> c1={1,2 ,3} c1 = [1] [2] [3] >> size(c1) ans = 1 3 >> length(c1) ans = 3 >> c2={3,; 4 ; 7} c2 = [3] [4] [7] >> size(c2) ans = 3 1 >> length(c2) ans = 3 Операция «вырезка» из массива ячеек ведет себя довольно неожиданно, поэтому на начальном этапе знакомства с пакетом лучше ее избегать. Для взаимного преобразования между массивами double и cell array, помимо присваивания в цикле, существуют стандартные средства – команды num2cell, mat2cell и cell2mat. Команда num2cell, как явствует из названия, предназначена для преобразования численных массивов в массивы ячеек. Цифру 2 (two) в обозначении операции следует читать как предлог to (к). >> a=[1 2; 3 4] a = 1 2 3 4 >> num2cell(a) ans = [1] [2] [3] [4] >>n1=num2cell(a,1) n1 = [2x1 double] [2x1 double] >> n1{1} ans = 1 3 >>n2=num2cell(a,2) n2 = [1x2 double] [1x2 double] Простейший вариант вызова команды num2cell – с одним аргументом. Результат в этом случае очевиден. При вызове с двумя аргументами первый содержит исходный массив, а второй – номер измерения, по которому матрицу предстоит «разрезать» на полосы. Так num2cell(a,1) «режет» матрицу по вертикали, а num2cell(a,2) – по горизонтали. Действие команды cell2mat обратно действию num2cell, причем вызов осуществляется только с одним входным аргументом: >> cell2mat(num2cell(a)) ans = 1 2 3 4 >> cell2mat(num2cell(a))-a ans = 0 0 0 0 Команда mat2cell наиболее интересна из всех перечисленных выше. На входе она принимает три аргумента – исходный массив и два массива, содержащих размеры вертикальных и горизонтальных блоков. Линии, по которым следует разбить массив, показаны на рисунке. Видно, 111 что размеры вертикальных блоков – 1 и 2, а горизонтальных – 2, 3 и 1. Эти значения – [1 2] и [2 3 1] как раз и являются вторым и третьим аргументами команды mat2cell. Листинг программы и значения массивов в ячейках после применения команды mat2cell показаны ниже. >> b=[ 1 2 3 4 5 6; 7 8 9 10 11 12; 13,14,15, 16, 17, 18 ] b = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 >> b1=mat2cell(b,[1 2],[2 3 1]) b1 = [1x2 double] [1x3 double] [ 6] [2x2 double] [2x3 double] [2x1 double] >> b1{1,1} ans = 1 2 >> b1{1,2} ans = 3 4 5 >> b1{1,3} ans = 6 >> b1{2,1} ans = 7 8 13 14 >> size(b) ans = 3 6 >> b1{2,2} ans = 9 10 11 15 16 17 >> b1{2,3} ans = 12 18 Использование структур и пользовательских классов Скорее всего, Вам не придется писать собственные классы средствами MATLAB . Тем не менее, многие стандартные тулбоксы реализованы с их применением. Поэтому полезно владеть как минимум основами использования классов в MATLAB . Первый вопрос, который возникает – как определить, к какому типу или классу принадлежит объект. Для этого существуют функции class и isa (is a ….). Первая из них возвращает имя класса, которому принадлежит объект, а вторая возвращает 1 или 0, в зависимости от того, принадлежит ли объект указанному в качестве второго параметра классу. Покажем их работу на примере уже рассмотренных типов данных double и cell array: >> a=[1 2 3] a = 1 2 3 >> b={4,5,6} b = [4] [5] [6] >> class(a) ans = double >> class(b) ans = cell >> isa(a,'double') ans = 1 >> isa(a,'cell') ans = 0 >> isa(b,'double') ans = 0 >> isa(b,'cell') ans = 1 >> class(class(a)) ans =char Напомним, что все пользовательские классы являются наследниками структуры (struct). Структура в пакете MATLAB создается командой struct, которая в качестве входных параметров принимает набор пар имени и значения поля: struct (имя1, значение1, имя2, значение2, …..). Обращение к полям структуры осуществляется при помощи оператора «точка». Значением поля структуры может быть объект любого класса, в том числе массив чисел или массив ячеек: >> d = struct('strings',{{'hello','yes'}},'lengths',[5 3]) d = strings: {'hello' 'yes'} lengths: [5 3] >> d.strings ans = 'hello' 'yes' >> d.strings{1}='no' d = strings: {'no' 'yes'} lengths: [5 3] Обратите внимание, что при объявлении структуры, у которой одно из полей является массивом ячеек, используются одни «лишние» скобки. Если их не поставить, то будет создан целый массив структур. При этом функция struct демонстрирует следующее поведение. Если лишь одно поле f1 в качестве значения имеет массив ячеек, то создается массив структур, в котором поле f1 каждого из элементов принимает поочередно значения из массива ячеек, а остальные поля 112 фиксированы. Если два или более полей заданы массивами ячеек, то эти массивы должны быть одинаковой длины и команда struct создает массив структур, в котором элементы поочередно принимают соответствующие значения. Если же количество ячеек разное, то выдается сообщение об ошибке: >> d = struct('strings',{'hello','yes'},'lengths',[5 3]) d = 1x2 struct array with fields: strings lengths >> d.strings ans = hello ans = yes >> d.lengths ans = 5 3 ans = 5 3 >> size(d) ans = 1 2 >> d(1) ans = strings: 'hello' lengths: [5 3] >> u=struct('f1',{1,2},'f2',{3,4}) u = 1x2 struct array with fields: f1 f2 >> u(1) ans = f1: 1 f2: 3 >> u(2) ans = f1: 2 f2: 4 >> d(2) ans = strings: 'yes' lengths: [5 3] >> u=struct('f1',{1,2},'f2',{3,4,5}) ??? Error using ==> struct Array dimensions of input 4 must match those of input 2 or be scalar. >> u=struct('f1',{1,2},'f2',{3,4},'f3',{5,6}) u = 1x2 struct array with fields: f1 f2 f3 >> u(1) ans = f1: 1 f2: 3 f3: 5 >> u(2) ans = f1: 2 f2: 4 f3: 6 Для работы со структурами имеются следующие функции: isfield, getfield, setfield, rmfield, fieldnames. Команда isfield (s,fldname) возвращает 1 если у структуры s есть поле fldname и 0 в противном случае Команда getfield (s,fldname) возвращает значение поля fldname структуры s, т.е. является синонимом для s.fldname. Функция setfield (s,fldname,fldvalue) возвращает копию объекта s, у которой полю fldname присвоено значение fldvalue. Сам объект s при этом не изменяется, т.е. это неравносильно присваиванию s.fldname=fldvalue. Функция rmfield(s,fldname) предназначена для удаления поля fldname у объекта s. Для того, чтобы добавить новое поле достаточно просто присвоить ему значение: s.new_field_name=some_value. Функция fieldnames(s) возвращает массив ячеек, содержащий имена полей структуры s. Ниже приведены примеры использования этих команд. >>u=struct('f1',1,'f2',0.1) u = f1: 1 f2: 0.1000 >> isfield(u,'f1') ans = 1 >> isfield(u,'f3') ans = 0 >> getfield(u,'f1') ans = 1 >>setfield(u,'f1',-1) ans = f1: -1 f2: 0.1000 >> u u = f1: 1 f2: 0.1000 >>u1=setfield(u,'f1',-1) u1 = f1: -1 f2: 0.1000 >> fieldnames(u) ans = 'f1' 'f2' >>rmfield(u1,'f2') ans = f1: -1 >> u1 u1 = f1: -1 f2: 0.1000 >>u2=rmfield(u1,'f2') u2 = f1: -1 >> fieldnames(u2) ans = 'f1' >> u2.f1=19 u2 = f1: 19 Поскольку пользовательские классы являются наследниками |