Матлаб. К. Ю. Петрова введение в matlab учебное пособие
 Скачать 2.57 Mb.
|
|
bandem тулбокса Optimization. Название программы связано с тем, что график функции Розенброка в окрестности точки минимума напоминает лежащий на боку банан. Найдем минимум с помощью команды fminunc, используя градиентный поиск. Поиск начинаем с точки х= – 1,9; у=2. >> f='100*(x(2)-x(1)^2)^2+(1-x(1))^2'; % описание функции >> GRAD='[100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2; 100*(2*x(2)-2*x(1)^2)]'; % описание градиента >>OPTIONS = optimset(OPTIONS,'gradobj','on'); % подключение опции градиентного поиска >> [x,fval] = fminunc({f,GRAD},[-1.9, 2],OPTIONS) x = 1.0000 1.0000 fval = 1.2371e-015 % результат Видим, что точка минимума успешно найдена. Пример. Другой известный тестовый пример, представляющий трудность для компьютерных алгоритмов оптимизации – функция Пауэла. Эта функция четырех аргументов: z =(a + 10b) 2 +5(c – d) 4 +(b – 2c) 4 + 10(a – d) 4 Она всегда положительна и достигает минимума при a = b = c = d = 0. Однако этот минимум очень плоский (в окрестности нуля функция растет крайне медленно), поэтому программы поиска экстремума могут остановиться (и выдать результат) довольно далеко от начала координат. Попытка аналитического отыскания экстремума функции Пауэла методом Ферма также наталкивается на трудности, поскольку приводит к необходимости решения системы из четырех кубических уравнений. Упражнение. Найдите минимум функции Пауэла средствами MATLAB. Для решения оптимизационных задач с ограничениями используются команды fminbnd и fmincon. Команда fminbnd (от bound – граница) предназначена для минимизации функции одной переменной на интервале 2 1 x x x . Например, x=fminbnd(inline('x^2+3*x+2'),0,1) даст ответ x = 0. Для минимизации функций от нескольких переменных с ограничениями более сложного вида применяют команду fmincon Список входных параметров функции довольно внушителен: x=fmincon(fun,x0, А,В,Аeq,Вeq,Lb,Ub,nonlcon,options,p1,p2,...). Первый параметр, как обычно, имя функции fun , второй – начальное приближение x0. Остальные параметры вплоть до опций options , 59 определяют ограничения. При необходимости можно указать дополнительные параметры p1, p2, ... . Если в задаче отсутствует какой-либо тип ограничений, соответствующие аргументы заменяются пустым массивом [ ] Рассмотрим подробнее типы ограничений. Матрицы A и B определяют ограничения вида B AX , матрицы Aeq и Beq определяют ограничения вида eq eq B X A . Векторы Ub и Lb определяют ограничения вида Ub X Lb . Нелинейные ограничения типа 0 ) ( X G описываются функцией nonlcon (от nonlinear constraints). Пример. Рассмотрим классическую задачу об экстремумах квадратичной формы f(X)=X T AX на сфере Х Т Х=1. Известно, что ее решения достигаются на собственных векторах симметричной матрицы А и равны ее собственным числам. Задав , 2 5 , 0 5 , 0 1 A получаем задачу об отыскании экстремума функции 2 2 2y xy x f при ограничении 1 2 2 y x Получим решение двумя способами, используя команды eig и fmincon. Спаособ 1.Найдем собственные векторы и собственные числа матрицы А, используя численный и символьный вариант ее задания: >> А=[1 1/2;1/2 2]; [v,d]=eig(А), [V,D]=eig(sym(А)), В результате получаем: v = -0.9239 0.3827 0.3827 0.9239 d = 0.7929 0 0 2.2071 V = [ 1, 1] [ 1+2^(1/2), 1-2^(1/2)] D = [ 3/2+1/2*2^(1/2), 0] [ 0, 3/2-1/2*2^(1/2)] Следовательно, экстремумы функции f определяются формулой 2 2 3 э f и равны 2,2 и 0,79. Максимум достигается в точке единичной окружности с координатами х=0,924, у=0,383; а минимум – в точке с координатами х=-0,924, у=0,383. Геометрически они соответствуют точкам касания единичной окружности с описанным и вписанным эллипсами (рис. 4.19). Рис. 4.19 x y x 2 + y 2 - 1 = 0 -2 -1 0 1 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 60 Способ 2. Для применения команды fmincon нужно описать функцию 2 2 2y xy x f и ограничение 0 1 2 2 y x g Для описания первой из них воспользуемся командой inline, для второй создадим m-файл nlcon.m: function [h,g] = nlcon(x) h = []; g = x(1)^2+x(2)^2-1; После этого набираем основную команду: >> f=inline('x(1)^2+x(1)*x(2)+2*x(2)^2'); x=fmincon(@(x) f(x),[1;1],[],[],[],[],[],[],[[],@(x)nlcon(x)]) x = -0.9239 0.3827 Мы получили координаты точки минимума. Чтобы получить координаты точки максимума, изменяем знак функции f: f1=inline('-(x(1)^2+x(1)*x(2)+2*x(2)^2)'); x=fmincon(@(x) f1(x),[1;1],[],[],[],[],[],[],[[],@(x)nlcon(x)]) x = 0.3827 0.9239 Видим, что оба способа решения привели к одинаковым результатам. Пример использования функции fmincon для подбора оптимальных параметров схемы моделирования, реализованной в SIMULINK, приводится в разд. 5.5.4. Все перечисленные команды могут возвращать один, два и более выходных аргументов: >> x=fminunc(@fff,[3 3]); >> [x,fval]=fminunc(@fff,[3 3]); >> [x,fval,flag]=fminunc(@fff,[3 3]); >> [x,fval,flag] ans = -0.0000 -0.0000 0.0000 1.0000 Первый выходной аргумент – точка минимума, второй – значение функции в этой точке, а третий – флаг, равный единице при успешном завершении алгоритма. Если флаг меньше нуля, то алгоритм не сошелся, а если он равен нулю, превышено максимальное количество вычисления функций. Отдельный класс оптимизационных задач образуют задачи линейного программирования, в которых и оптимизируемый критерий, и ограничения линейны. В них требуется найти экстремум критерия n n x c x c J 1 1 при наличии ограничений в виде неравенств , , 2 , 1 , 1 1 m i b x a x a i n in i Обычно эти условия записывают в матричной форме , b AX extr X c T Здесь b и c – векторы-столбцы, А – матрица размера m n. Ограничения типа равенств можно отдельно не рассматривать, так как они легко сводятся к неравенствам, например вместо 0 3 2 1 x x можно записать два неравенства 0 3 ; 0 3 2 1 2 1 x x x x Для решения задач линейного программирования разработаны различные численные методы, наиболее известным из которых является симплекс-метод. В MATLAB задачи линейного программирования решают с помощью функции linprog. В простейшем случае у нее три входных параметра: linprog(с, A. b). В этом случае решается задача минимизации выражения с T *x при условии A*x<=b (сравнение производится по всем строкам). Поясним ее применение на конкретном примере. 61 Пример. Задача о производстве стульев.Мебельная фабрика может выпускать стулья двух типов, стоимостью 6000 и 12000 рублей. Имеются следующие ресурсы: 440 погонных метров досок, 65 кв.м. обивочной ткани и 320 человеко-часов трудовых ресурсов. На изготовление одного стула требуются следующее количество ресурсов: Стул Расход досок Расход ткани Расход времени Первый 2 0.5 2 Второй 4 0.25 2.5 Ресурс 440 65 320 Требуется так спланировать производство стульев, чтобы общая цена продукции была максимальной. Решение. Перейдем к математической формулировке задачи. Обозначим через х количество стульев первого типа, через у – количество стульев второго типа. Тогда условия задачи сводятся к следующему: max 12 8 y x – оптимизируемый критерий. 440 4 2 y x – ограничение по расходу досок 65 25 0 5 0 y x – ограничение по расходу ткани 320 5 2 2 y x – ограничение по расходу времени. Матричная форма записи: , , 12 8 max, y x X c X c T 320 65 440 , 5 2 2 25 0 5 0 4 2 b A Приведем два способа решения этой задачи в MATLAB – графический и численный. Способ 1 (графический). Для графического решения построим на плоскости (x, y) три прямые, соответствующие ограничениям по трем ресурсам. По оси x будем откладывать количество стульев второго вида, по оси y количество стульев первого вида. Полученные прямые показаны на рис.4.20. Они, вместе с осями координат, задают область допустимых решений в виде неправильного пятиугольника. На том же рисунке показано семейство прямых 12 8 const y x Решение задачи дает крайняя правая прямая этого семейства, касающаяся многоугольника допустимых решений в точке с координатами (80, 60). Графики построены в МАТLAB с помощью следующей программы: x=0:0.2:300; y1=-2*x+220; y2=(-1/2)*x+130; y3=(-5/4)*x+160; plot(x,y1,x,y2,x,y3); grid; hold on for c=0:60:1460 y=-3/2*x+c/8; plot(x,y,'black');grid on; end 62 Рис. 4.20. Графическое решение. Способ 2 (численный). Для численного решения той же задачи применяем функцию linprog. >>X=linprog(-[8; 12],[2 4;0.5 0.25;2 2.5;],[440;65;320]) Optimization terminated successfully. X= 60.0000 80.0000 По условиям задачи требовалось найти максимум, поэтому, чтобы свести задачу к поиску минимума, первый параметр взят с коэффициентом -1. Мы получили то же решение, что и первым способом – максимальная прибыль будет получена, если выпустить 60 стульев первого типа и 80 стульев второго типа. Функции от матриц Нормы и числа обусловленности В разделе 1 были рассмотрены команды MATLAB для определения таких числовых характеристик матриц, как ранг, определитель, след, собственные числа. Из других команд этого ряда отметим svd, cond и norm, служащие для определения сингулярных чисел, чисел обусловленности и норм (векторных и матричных). По команде S=svd(A) вычисляются сингулярные числа матрицы А, т.е. положительные квадратные корни из собственных чисел матрицы А Т А. Они используются при определении ранга матрицы, при оценке ее обусловленности, при решении задач аппроксимации и редукции. Эта же команда с тремя выходными аргументами находит сингулярное разложение матрицы А. Заметим, 63 что для симметричных положительно определенных матриц сингулярные числа совпадают с собственными числами, в чем можно убедиться, сравнивая результаты команд svd и eig. Обе эти команды могут работать с символьными аргументами. Приведем простой пример: >>A =[2 4; 6 2]; >>S=svd(A) >>S= svd(sym(A)) А= 2 4 6 2 S=7.2361 2.7639 S=5+5^(1/2) 5-5^(1/2) Максимальное из сингулярных чисел матрицы определяет ее спектральную норму, а отношение максимального сингулярного числа к минимальному дает число обусловленности по этой норме. Так, в рассматриваемом примере норма матрицы равна , 2361 , 7 5 5 2 A а число обусловленности 618 , 2 2 5 3 5 5 5 5 2 (любители математики узнают здесь золотое сечение, увеличенное на единицу). В математике известен целый ряд векторных и матричных норм. Основные из них могут быть вычислены с помощью команды norm и ее модификаций. Три наиболее употребительные векторные нормы вычисляются с помощью следующих команд MATLAB: norm (х) и norm(х, 2) – вторая евклидова норма вектора 2 1 2 2 1 2 ) ( n x x X (именно ей мы пользуемся в обыденной жизни для измерения длины); norm(х, 1) - первая или модульная норма вектора n x x X 1 1 (мы пользуемся ей, определяя расстояние в городе с прямоугольной планировкой улиц); norm(x, inf) – бесконечная (чебышевская) норма вектора ) ( max i i x X (это габаритный размер параллелепипеда с диагональю Х). Эта же команда позволяет вычислять матричные нормы: norm(A) и norm(A, 2) – спектральная норма матрицы (равна ее наибольшему сингулярному числу) ), ( max 2 1 2 A A A T i i тот же результат можно получить как max(svd(A)); norm(A, 1) – первая (столбцовая) норма матрицы А (определяется столбцом с наибольшей суммой элементов) ) ( max 1 mj ij i a a A ; norm(A, inf) – бесконечная (строчная) норма матрицы А (определяется строкой с наибольшей суммой элементов) ) ( max 1 in i i a a A ; norm(A, 'fro') – фробениусова норма матрицы (равна евклидовой длине вектора, составленного из элементов матрицы) ). ( ) ( 2 1 2 1 2 2 12 2 11 A A tr a a a A T nm F Матричные нормы используются, в частности, при оценке степени обусловленности матриц. Как известно, числом обусловленности (condition number) называется величина ) ( cond 1 A A A Ee можно найти с помощью команды cond(A). Минимальное значение числа обусловленности μ=1 соответствует идеально обусловленной матрице. Чем больше μ,тем хуже обусловленность и тем большей погрешностью будет сопровождаться численное решение системы линейных алгебраических уравнений AX=b. Различным матричным нормам будут соответствовать различные числа обусловленности, требуемая модификация указывается вторым параметром команды cond. При решении задач аппроксимации, оптимального управления и многих других возникает необходимость вычисления норм функций одной переменной ), (t f y , заданных на интервале 0 T t В технических приложениях это чаще всего нормы входных и выходных сигналов некоторой системы. 64 Обычно рассматривают три классические нормы функций – модульную или первую, гильбертову или вторую (квадратичную) и чебышевскую (другие названия – равномерная, бесконечная): T dt t y y 0 1 ) ( , 2 / 1 0 2 2 ) ( T dt t y y , ) ( max 0 t y y T t В MATLAB нет специальных команд для вычисления этих норм, однако они легко могут быть найдены с помощью функции численного интегрирования trapz и команды выбора максимального элемента массива max: Вычисление модульной нормы функции y(t): N=trapz(t, abs(y)); Вычисление гильбертовой нормы функции y(t): N= sqrt(trapz(t, y.*y)); Вычисление чебышевской нормы функции y(t): N=max(abs(y)). При этом переменные t и y должны быть заданы массивами своих значений. Возможны и другие варианты, например, для вычисления интеграла можно использовать команду lsim с интегратором в качестве первого аргумента. Пример. Рассмотрим систему второго порядка с передаточной функцией Q(p)=1/(p+ ) 2 , представляющую собой последовательное соединение двух одинаковых апериодических звеньев. Весовая функция системы q(t) и ее энергия N 2 (квадрат гильбертовой нормы) на бесконечном интервале времениимеют вид 4 1 ; 0 3 2 2 dt q N te q at Приведем два способа получения нормы N в MATLAB, полагая α=1. а) Численный способ. Используем команды impulse и trapz: >>sys=tf(1,[1 2 1]); [q,t]=impulse(sys); % получение весовой функции системы и массива времени >>Nq=trapz(t, (q .* q) ); nq = sqrt(Nq) % Nq - интеграл на интервале [ 0, T ], где T = max(t); nq = 0.5000 % результат б) Символьный способ. Используем обратное преобразования Лапласа и символьное интегрирование: >>syms p; q=ilaplace(1/(p+1)^2) q =t*exp(-t) % весовая функция системы >>Mq=int(q^2,0, inf); mq=sqrt(Mq) % Mq – интеграл в пределах от нуля до бесконечности mq =1/2 % результат В обоих случаях результат совпадает с теоретическим. Упражнение. Дано апериодическое звено с передаточной функциейQ(p)=1/(p+α). Его весовая функция q(t) и квадрат ее гильбертовой нормы на интервале (0, Т) определяются формулами ). 1 ( 2 1 ; 2 0 2 2 T T t e dt q N e q Найдите средствами MATLAB норму N для значений , 2 , 1 T |