Ход работы

Задание на лабораторную работу

1. Изучить материал теоретической части.

2. Синтезировать систему нечеткого логического вывода для функции, предложенной в практической части.

3. Построить трехмерное изображение функции, использую программный код, представленный в первом пункте практической части.

4. Сравнить полученную поверхность с поверхностью «входы – выход», которая соответствует синтезированной нечеткой системе. Сделать вывод.

Отчет должен содержать цель работы, задание и практическую часть, включающую:

1. скриншоты типа рис. 2, 4, 6, 7 с пояснениями;

2. качественную оценку сравнения поверхностей, полученных разными способами;

3. выводы, наличие которых определено заданием.

Практическая работа №2
«Изучение пакета прикладных программ NeuralNetworkToolbox в среде системы MATLAB»

Цель: изучить принципы построения нейронных сетей в пакете прикладных программ NeuralNetworkToolbox в среде системы MATLAB.

Постановка задачи: Создать и обучить нейронную сеть выполнению операции у = х₁² + х₂, если заданы последовательности входа

Р1 = [1 0.5 0 1; -2 0 0.5 1] и цели Т1 = [-1 0.25 0.5 2].

Ход работы

Сформируем последовательности входов и целей в рабочей области GUI-интерфейса NNTool, используя окно CreateNewData (рис. 2.1).

Рис. 2.1. Ввод входных данных.

Выберем нейронную сеть типа feed-forwardbackprop с прямой передачей сигнала и с обратным распространением ошибки. Схема этой сети показана на рис. 2.2.

Рис. 2.2. Структура нейронной сети типа

Выполним инициализацию сети, для чего выберем закладку Initialize, откроется диалоговая панель, показанная на рис. 2.3. Диапазоны значений исходных данных выберем по входам из ниспадающего меню Getfrominput. Для ввода установленных диапазонов и инициализации весов надо воспользоваться кнопками SetRanges (Установить диапазоны) и InitializeWeights (Инициализировать веса). Если требуется вернуться к прежним диапазонам, то следует выбрать кнопки RevertRanges (Вернуть диапазоны) и RevertWeights (Вернуть веса).

Рис. 2.3. Инициализация сети.

Затем выполняется обучение сети, для чего выбирается закладка Train и открывается диалоговая панель, показанная на рис. 2.4.

Рис. 2.4. Обучение сети. Информация об обучающих последовательностях.

На рис. 2.5. представлены параметры обучения выбранной сети.

Рис. 2.5. Обучение сети. Параметры обучения.

Применяя эти закладки, можно установить имена последовательностей входа и цели, а также параметров процедуры обучения. Теперь можно приступить к обучению сети (кнопка TrainNetwork).

Рис. 6. Результат обучения сети (график).

Качество обучения сети с прямой передачей сигнала на выбранной обучающей последовательности поясняется на рис. 2.6. Соответствующие веса и смещения можно увидеть, если выбрать закладку Weights (рис. 7).

Рис. 7. Веса и смещения сети.

Для удобства работы можно экспортировать созданную нейронную сеть в рабочую область системы MATLAB и получить информацию о весах и смещениях непосредственно в рабочем окне системы.

Результаты обучения можно просмотреть в окне Network/DataManager, выбрав кнопку Manager (рис. 2.7).

Проведем тестирование нейронной сети. Для этого создадим дополнительный набор данных Dtest=[0.25 0.5 0.5 0.75 1; -2 -1 0 0.4 1] и набор желаемого результата Ttest=[-1.9375 -0.75 0.25 0.9625 2]. На рисунке 9 показаны выходные значения и абсолютная погрешность вычисления.

В результате настройки будут установлены следующие значения весов и смещений, которые можно увидеть, выбрав закладку Weights. Для данной сети вектора весов равны IW{1, 1} = [-2.0959 -2.149], LW=[-14.172], а смещения b{1} = [0.55339], b{2} = [-6.2327]. Таким образом, линия переключения, разделяющая плоскость на 2 области, описывается следующим образом:

L: -14.172(-2.0959p1 – 2.149р2 +0, 55339)-6.2327 = 0.

Рис. 2.8. Запуск алгоритма на тестирующих данных.

Рис. 2.9. Результаты тестирования.

Варианты заданий

Практическая работа №3
«Оптимизация многоэкстремальных функций нескольких переменных»

Цель: получить практический опыт использования, настройки, оптимизации генетических алгоритмов для нахождения глобального минимума/максимума многоэкстремальных функций. Сравнить ГА с классическими алгоритмами поиска экстремумов.

Необходимые навыки работы:

· работа с файлами (открыть, сохранить, создать m - файл);

· базовые навыки работы в консоли: работа с массивами, матрицами, построение двумерных, трехмерных графиков функции;

· навыки работы с графическим интерфейсом генетических алгоритмов (gatool);

· умение писать и вызывать собственные функции для исполнения виде m-файлов.

Задание к лабораторной работе:

1. Построить двумерный, трехмерный график заданной функции.

2. Определить интервал поиска глобального экстремума таким образом, чтобы в него попадало несколько локальных максимумов.

3. Найти глобальный минимум/максимум функции 1-ой, 2-х, 20-ти, 40-ка переменныхс помощью инструментария генетический алгоритмов.

4. Определить критический размер популяции.

Порядок выполнения работы:

1. Изучить порядок работы с инструментом gatool.

2. Рассмотреть предложенный пример по поиску глобального минимума функции.

3. Выбрать свою функцию и реализовать соответствующую ей функцию пригодности в m-файле в соответствии с примером для дальнейшей оптимизации.

4. Построить одно- и двумерный график заданной функции, указать глобальный экстремум и выбрать интервал поиска.

5. Найти экстремум функции 2-х переменных классическим методом

6. Запустить gatoolс значениями по умолчанию и найти экстремум функции 1-ой, 2-х переменных.

7. Варьируя настройки размера популяции, критериев остановки, настройки функции масштабирования, пригодности, отбора, мутации, скрещивания определить оптимальные параметры настройки ГА для функции 2-х переменных.

8. Найти критический размер популяции, при котором ГА перестает находить глобальный минимум.

9. Найти экстремум функции 20-ти, 40-ка переменных.

10. Найти минимум, используя гибридную функцию.

11. Экспортировать настройки gatoolв рабочее пространство Matlaba, и найти экстремум из командной строки.

12. Определить точность найденных экстремумов, результаты занести в таблицу.

13. Задействовать субпопуляцию, изменяя настройки миграции определить её влияние на конечный результат.

14. Для одного из случаев поиска минимума построить все возможные графики, отражающие ход процесса поиска, и объяснить их значение.

15. Изменить критерий остановки ГА на меньшее количество поколений = 10 и построить генеалогическое дерево.

16. В каждом случае поиска минимума провести его не менее 3-х раз, замерить среднее время выполнения, занести в таблицу усредненные результаты.

17. Повторить п. №7, варьируя настройки ГА следующим образом:

a. увеличить размер популяции в двое, уменьшить – вдвое;

b. изменить критерии остановки, увеличив, уменьшив максимальное число поколений;

c. опробовать все возможные значения для функции отбора и отметить их влияние на результат;

d. изменить тип скрещивания на одноточечный, двухточечный, разрозненный, центральный, эвристический;

e. использовать следующие настройки отбора: стохастические, остаточные, метод рулетки, турнир, унифицированные;

f. установки мутации: Гауссовская, универсальная;

18. Результаты занести в таблицу:

Варианты заданий:

1. Функция Гаусса четвертой степени (Ди Янг 4) (Gaussianquartic (DeJong4)):

, где rand(1) - функция, возвращающая случайную величину в интервале(0; 1).

Минимум равен 0 в точке, где x_i=0.

2. ФункцияРастригина (Rastrigin's function):

.

Минимум равен 0 в точке, где x_i=0.

Локальный минимум в точке, где одна координата равна 1, а остальные равны 0.

3. Функция Швефела (Schwefel's (Sine Root) Function):

, где .

Минимум равный 0 в точке, где x_i=420.9687.

Локальный минимум в точке, где одна координата равна -302.5232, а остальные равны 420.9687. Т.к. локальный минимум находится далеко от глобального, то есть опасность, что алгоритм " собьется с пути". Угол под знаком синуса получается в радианах.

4. Функция Гривангка (Griewangk's Function):

, где .

Минимум равный 0 в точке, где x_i=0.0.

При n=10 есть ещё 4 локальных минимума равных 0, 0074 приблизительно в точке (+-Pi, +-Pi*sqrt(2), 0,..., 0).

5. ФункцияЭкклея (Ackley's Function):

, где .

Минимум равный 0 в точке, где x_i=0.0.

Оптимизация многоэкстремальных функций нескольких переменных с помощью ГАУСА

В типичной задаче оптимизации существует набор переменных, влияющих на процесс, и формула или алгоритм, который использует эти переменные, чтобы построить полную модель этого процесса. При этом задача заключается в том, чтобы найти такие значения переменных, которые некоторым образом оптимизируют модель. Если моделью является формула, то обычно ищут максимум или минимум функции, которую данная формула представляет. Существуют много математических методов, которые решают (и очень быстро) задачи оптимизации в том случае, если это задачи с " хорошим поведением". Однако традиционные методы часто терпят крах, если задача " ведет себя" недостаточно хорошо.

Реальные прикладные задачи оптимизации очень сложны. Современные методы оптимизации далеко не всегда справляются с решением реальных задач без помощи человека. Нет, пока такой теории, которая учла бы любые особенности функций, описывающих постановку задачи. Рассмотрим на простом примере основные моменты оптимизации функции двух переменных при помощи генетических алгоритмов с использованием MatlabGeneticAlgorithmTool. Взятая нами целевая функция является кусочно-непрывной функцией, то есть имеются участки гладких функций с встроенными между ними разрывами одним примечательным участком, где функция является всюду недифференцируемой. Дополнительно, что характерно для большинства задач оптимизации, следует отметить, что взятая нами функция является гладкой в области точки минимума. Наша задача - найти минимум следующей функции: .

Первым шагом будет составление целевой функции ввидеm-файла. В данном случае она будет совпадать с уравнением минимизируемой функции.

Требуемый для расчета данной функции М-файл должен формировать вектор-строку размерности i в соответствии с имеющимися переменными x_i, а также возвращать скаляр со значением функции в точке х. для формирования m-файла следует проделать следующие шаги:

1. Выбрать команду New в меню MATLAB File.

2. Выбрать команду m-File. В этом случае откроется редактор нового М-файла.

3. В М-файле ввести следующий текст:

function scores = rastriginsfcn(pop)

scores = 10.0 * size(pop, 2) + sum(pop.^2 - 10.0 * cos(2 * pi.* pop), 2);

где:

.* - умножение 2-х массивов. Т.е. для массивов А, В запись вида А.* В будет означать для каждого элемента C(i, j) = A(i, j) * B(i, j). Массивы должны иметь одинаковую размерность, за исключением, если один из массивов имеет размерность 1х1.

.^ - соответственно - возведение массива в степень.

pop- параметры передаваемый функции. В нашем случае матрица текущей популяции для каждой строки который необходимо вернуть значение пригодности, имеющая количество столбцов равного количеству генов в хромосоме (количеству переменных искомой функции), а количество строк - равному размеру популяции.

Или другой вариант записи этой же функции:

function scores = rastriginsfcn2(pop)

for j = 1: size(pop, 2)

f(j) = 0.0;

x = pop(j,:);

scores(j, 1) = 10.0 * size(x, 2) + sum(x.^2 - 10.0 * cos(2 * pi.* x));

end

4. Сохранить М-файл в текущей рабочей директории MATLAB

Далее необходимо построить трехмерный график данной функции. Это делается с помощью команды[-5 5; -5 5])., в главном окне MatLab'а. В результате получим график изображенный на 3.1.

Рисунок 3.1. - Функция Растригина

Для работы с инструментарием генетического алгоритма используются следующие функции:

ga - Поиск минимума функции при помощи генетического алгоритма

gaoptimget - Получить значения структуры опций генетического алгоритма

gaoptimset - Создать структуры опций генетического алгоритма

gatool – Открыть графический инструментарий генетического алгоритма

Воспользуемся последней, в результате откроется окно (.2 - Графический интерфейс GeneticAlgorithmTool), в котором зададим необходимые нам параметры.

Рисунок 3.2 - Графический интерфейс GeneticAlgorithmTool

· Функция пригодности (Fitnessfunction) = @rastriginsfcn - та целевая функция, которую необходимо оптимизировать.

· Количество переменных (Numberofvariables) = 2 оптимизируемой функции.

· Plotinterval = 1 - устанавливает количество поколений, между обновлениями графиков.

· Bestfitness = 1 - отображает в осях значение - поколение - лучшее и среднее значение на каждом шаге итерации.

Нажмем кнопку старт, в результате чего будет отображен ход процесса поиска и найденное конечное значение минимума: (0, 01224; 0, 00847).

Замер времени поиска можно выполнить следующим образом:

· Кликнуть мышкой на команду «ExporttoWorkspace».

· Вдиалоговомбоксе «Export to Workspace dialog box» ввестиимядляструктурызадачи, например, «my_problem» вполе «Export problems and options to a MATLAB structure named».

· В консоле MATLAB вызвать функцию ga с аргументом my_problem:

tic; [x fval] = ga(my_problem); toc

· В результате получим в главном окне время выполнения задачи:

Elapsedtimeis 0.163189 seconds.

Вопросы и задания к защите:

1. Какова схема работы ГА?

2. В чем преимущество ГА по сравнению с другими алгоритмами поиска.

3. На каком этапе генетического алгоритма в нескольких случайно выбранных особях нового поколения изменяются некоторые гены?

4. По какому принципу выбирается индивидуум для скрещивания?

5. Какие этапы включает в себя генетический алгоритм?

6. На сколько точные решения получаются в результате работы генетического алгоритма и с чем это связано?

7. На основе какой функции строится функция приспособленности?

8. Приспособленность особи на фенотипическом уровне определяет …

9. Приспособленность особи на генотипическом уровне определяет …

10. Каковы сильные и слабые стороны ГА?

11. Что обозначают красные, синие, зеленые линии на генеалогическом дереве?

12. В какой зависимости находится время выполнения от количества особей в поколении? От количества переменных функции?

13. Почему среднее расстояние между особями одного поколения сокращается?

14. Почему совместно с ГА применяют другие методы оптимизации (такие как градиентный спуск)?

Практическая работа №4
«Оптимизация ПИД регулятора САР с помощь генетических алгоритмов»

Цель: изучить примеры практического применения генетических алгоритмов. Сравнить полученные результаты с классическими методами настройки регуляторов.

Необходимые навыки работы:

· работа с файлами (открыть, сохранить, создать m - файл);

· базовые навыки работы в консоли: работа с массивами, матрицами, построение двумерных, трехмерных графиков функции;

· навыки работы с графическим интерфейсом генетических алгоритмов (gatool);

· умение писать и вызывать собственные функции для исполнения виде m-файлов;

· основы отображения информации в графическом виде (функция plot).

Задание к лабораторной работе:

1. Для полученной модели системы выполнить линеаризацию.

2. Настроить ПИД-регулятор стандартными методами и с помощью генетических алгоритмов.

3. Сделать выводы о качестве настройки.

4. Повторить пункты 1-3 для не линеаризованной системы.

Порядок выполнения работы:

1. Получить вариант автоматической системы для моделирования.

2. Линеаризовать систему.

3. Определить устойчивость.

4. Рассчитать настройки ПИД-регулятора классическим методом

5. Составить схему линейной системыс ПИД-регулятором в simulink.

6. Провести моделирование и зафиксировать показатели качества ПХ и её график.

7. Добавитьнасхемублок “simulink Response Optimization/ Signal Constraint “.

8. Оптимизировать систему, используя градиентный спуск, симплекс метод.

9. Зафиксировать полученные значения и определить показатели качества ПХ.

10. Оптимизировать систему с использованием “GeneticAlgorithmtoolbox”.

11. Зафиксировать полученные значения и определить показатели качества ПХ.

12. Собрать нелинейную систему и повторить для нее пункты 8-11.

13. Свести данные в таблицу.

14. Сравнить результаты оптимизации, полученные разными методами.

15. Сделать выводы о проделанной работе.

Отчет должен содержать:

1. Схема линейной системы.

2. График ПХ для разомкнутой системы.

3. ПХ дляполученная после каждого метода оптимизации

4. Схема нелинейной системы, и её графики аналогично линейной.

5. Сводная таблица результатов.

6. Выводы.

Расчет коэффициентов классических регуляторов методом генетических алгоритмов.

Получившие распространение имитационные методы моделирования в задачах управления предъявляют качественно новые требования к решению задач параметрической оптимизации. На замену аналитическим кос­венным приемам вычисления оптимальных настроечных параметров регуляторов все активнее приходят численные алгоритмы оптимизации.

Опыт исследования алгоритмов управления показал, что для простых одноконтурных систем автоматического регулирования (АСР) с линейными регуляторами задачи оптимизации, как правило, являются одноэкстремальными. Однако для сложных многоконтурных систем управле­ния и систем управления с нейроконтроллерами характерно наряду с глобальным наличие большого числа локальных экстремумов. Кроме того, локальные экстремумы появляются и при введении ограничений на пространство поиска.

Для решения одноэкстремальных задач оптимизации существует достаточное число градиентных и численных алгоритмов. Одним из таких алгоритмов является метод деформируемого многогранника Нелдера-Мида. При оптимизации одноконтурной АСР с ПИ-регулятором такой алгоритм устойчиво находит оптимальные значениянастроечных параметров и для функции цели вида

где - интеграл по модулю регулируемой величины на интервале времени переходного процесса ; - соответственно, заданная и текущая степень затухания переходного процесса регулирования; - масштабный коэффициент, учитывающий вес штрафной функции.

Применениеметода деформируемого многогранника для оптимизации двухконтурной АСР с дифференциатором и АСР с нейроконтроллерами различ­ной структуры приводит к неоднозначности решения. В каждом случае результаты зависят от выбранных начальных координат. Из этого следует вывод о многоэкстремальности подобных задач, и для их решения требу­ются методы глобальной оптимизации, одними из которых и является генетические алгоритмы.

В этом разделе рассматривается метод настройки параметров ПИД-регулятора для систем автоматического управления с применением генетического алгоритма, сравнение с аналогичными регуляторами, полученными стандартными методами проектирования. Расчёт коэффициентов этих регуляторов осуществляется по классическим формулам исходя из требуемых критериев качества переходных процессов.

Передаточная функция ПИД-регулятора представлена формулой:

.

Видно, что необходимо определить три коэффициента регулятора: K_P, K_I, K_D, которые можно рассчитать воспользовавшись любым классическим подходом синтеза ПИД регулятора (с использованием диаграмм Боде или методом корневого годографа), но при этом необходимо чтобы система была устойчивой и имела запас по фазе и амплитуде, следовательно, необходимо построение корневого годографа, ЛАЧХ, ФЧХ, что соответственно занимает некоторое время, и при изменении параметров системы необходимо производить построение ЛАЧХ и ФЧХ заново. Таких же результатов можно достичь при помощи ГА, для которого не нужно определять запасы устойчивости и строить диаграммы. Для нахождения коэффициентов методом ГА можно использовать как устойчивый, так и неустойчивый объект управления. Коэффициенты регулятора подбираются самим ГА таким образом, что система станет устойчивой. Если система изначально устойчива, то проблем при использовании ГА не возникает, он всегда подберёт коэффициенты такими, чтобы переходный процесс (ПП) удовлетворял требуемым показателям качества регулирования.

При синтезе коэффициентов регулятора для устойчивого объекта метод ГА определит коэффициенты, при которых система приобретет требуемую динамику.Рассмотрим пример получения коэффициентов ПИД-регулятора с использованием метода генетических алгоритмов в среде MATLAB для устойчивого объекта системы заданной передаточной функцией . Для этого соберем её в simulink (см. рисунок 4.1). В представленной модели оптимизация будет проводиться по минимизации среднеквадратичной ошибки, но имеется возможность задать любой другой критерий.Для этого требуется лишь немного откорректировать функцию пригодности.

Рисунок 4.1 - Модель оптимизируемой системы

Итак, сначала снимем ПХ разомкнутой системы без регулятора. Получим характеристику, представленную на Рисунок4.2.

Следующим шагом будет добавление ПИД -регулятора. И создание функции пригодности, которая для каждой особи будет запускать объект на моделирование, и получат конечное значение интегратора на выходе " Out1". Вот её код:

functionf = PID(x)

%Устанавливаем значение параметров ПИД регулятора исходя из

%полученныххромосом

set_param ('GA_OPT/PID Controller', 'p', num2str(x(1)));

set_param ('GA_OPT/PID Controller', 'i', num2str(x(2))); %

set_param ('GA_OPT/PID Controller', 'd', num2str(x(3))); %

%Запускаем моделирование в simulink

[t, z, y] = sim('GA_OPT', [0 30]);

%Получаем последнее значение интегратора суммирующего квадрат %ошибки и присваиваем его значение функции пригодности

f = y(length(y));

end

Заодно напишем функцию, отображающую ПХ для лучшей хромосомы в популяции:

function state = PIDPlot(options, state, flag, locations)

% получимномерлучшейособи в популяциисохранимеё в х

[unused, i] = min(state.Score);

x = state.Population(i,:);

% настроимсоответствующе ПИД-регулятор

set_param ('GA_OPT/PID Controller', 'p', num2str(x(1)));

set_param ('GA_OPT/PID Controller', 'i', num2str(x(2)));

set_param ('GA_OPT/PID Controller', 'd', num2str(x(3)));

% и запустимнамоделированиесимулинк

[t, z, y] = sim('GA_OPT', [0 30]);

% пополученнымпараметрампостроимграфик

grid on;

hold all;

plot(t, z(:, 4));

holdoff;

Настроив и запустив gatool получим оптимальные параметры ПИД-регулятора:

K_P = 8, 8355; K_I = 1, 487; K_D = 9, 6923.

В результате чего суммарная квадратичная ошибка за 50 секунд моделирования составила 0, 7612.

А также график отражающий оптимизацию ПХ с течением поколений (рисунок 4.3).

Рисунок 4.2 - ПХ системы без регулятора

Рисунок 4.3 - Изменение ПХ в процессе настройки ПИД - регулятора

По графикам видно, что время переходного процесса сократилось примерно с 13 до 5 с, а выброс уменьшился с 50% до 18%.

Теперь воспользуемся методом градиентного спуска для настройки параметров регулятора. Получим: K_P = 12, 3459; K_I = 2, 4111; K_D = 14, 5429 и суммарную квадратичную ошибку за 50 секунд моделирования 0.771.

Таким образом, генетический алгоритм справился с поставленной задачей не хуже градиентного спуска и симплекс метода.

Анализируя полученные результаты моделирования можно сделать вывод, что данный подход получения коэффициентов классических регуляторов даёт хорошие результаты как для устойчивых, так и не устойчивых объектов.

Главное достоинство метода ГА в том, что его применение не требует от разработчика приводить описание объекта управления или контроллера к определенному виду, например, к описанию в матрицах состояния, что необходимо для нахождения коэффициентов регулятора состояния. Единственное требование, которое необходимо соблюдать - это вычисление критерия качества при любых настройках регулятора.

Очевидно, что в случае применения ГА модель объекта может быть существенно нелинейной, может иметь шумы и ограничения, может иметь любой порядок и параметрические возмущения. Все это делает невозможным рассчет классических регуляторов, без внесения существенных упрощений и допущений. Что на практике совершенно не гарантирует применимость полученных результатов.

Преимущества метода особенно заметны в условиях ограниченного времени на проектирование системы управления.

Практическая работа №5
«Решение задач сводящихся к задаче коммивояжера. Разработка собственных функций создания, мутации и скрещивания».

Цель: изучить примеры практического применения генетических алгоритмов. Получить навыки модернизации стандартных функций мутации, скрещивания и создания популяции для конкретных приложений. Научится графически интерпретировать конечный результат и использовать нелинейные ограничения.

Необходимые навыки работы:

· работа с файлами (открыть, сохранить, создать m - файл);

· базовые навыки работы в консоли: работа с массивами, матрицами, построение двумерных, трехмерных графиков функции;

· навыки работы с графическим интерфейсом генетических алгоритмов (gatool);

· умение писать и вызывать собственные функции для исполнения виде m-файлов;

· основы отображения информации в графическом виде (функция plot).

Задание к лабораторной работе:

1. Сформулировать задачу в терминах генетического алгоритма.

2. Реализовать собственные функции мутации, скрещивания и задания популяции.

3. Решить задачу, изменяя настройки генетического алгоритма для достижения наилучшего результата.

4. Построить графики.

Порядок выполнения работы:

1. Рассмотреть предоставленный пример.

2. Определить особенности задачи не позволяющие воспользоваться стандартными функциями инструментария генетического алгоритма.

3. Составить список разрабатываемых функций.

4. Составить алгоритм их работы, используя отличные от примера варианты кроссовера и мутации.

5. Реализовать полученные алгоритмы в виде отдельныхm-файлов.

6. Изменить настройки инструментария в соответствии с написанными функциями.

7. Решить поставленную задачу.

8. Определить корректность найденного решения и в зависимости от него внести коррективы в соответствующие алгоритмы.

9. Написать функцию для наглядного графического отображения лучшей особи в текущем поколении.

10. Сделать выводы о проделанной работе.

Отчет должен содержать:

1. Алгоритм работыизмененных функций

2. Код написанных функций, с комментариями.

3. Графическое изображение лучшей особи в начальной популяции.

4. Графики, отражающие ход процесса поиска, а также данные о его результате причине остановке.

5. Графическое изображение найденного решения.

6. Выводы.

Варианты заданий:

1. Решить задачу коммивояжера для 10 и 100 городов, расстояния между городами сгенерировать случайным образом.

2. Решить задачу о дыропробивном процессе, если требуется пробить 50 дырок, количество разных отверстий – 4, время смены инструмента на соседний с ним 2 с. Расположение отверстий сгенерировать случайным образом.

3. Решить задачу о ранце. Требуется загрузить в рюкзак: яблоки (полезность 8, вес 200г.), апельсины(10; 200), мясо – (13; 600), тушенка – (12; 500), спички (20; 50), карта (15; 50), квас (21; 1500), палатка (30; 3333), яйца (13; 434), пылесос (1; 2634), лопата (28; 1000), компас (34; 100), ноутбук (5; 3444). Максимальная снаряженная масса рюкзака 5000 г.

4. Решить Диофантово уравнение a+2b+3c+4d=30, где a, b, c и d - некоторые положительные целые числа.

Решение классической задачикоммивояжера и сводящихся к ней.

Задача коммивояжера. Задача коммивояжера (в дальнейшем сокращённо - ЗК) является одной из знаменитых задач теории комбинаторики. Она была поставлена в 1934 году, и об неё, как об Великую теорему Ферма обламывали зубы лучшие математики. В своей области (оптимизации дискретных задач) ЗК служит своеобразным полигоном, на котором испытываются всё новые методы.

Постановка задачи следующая.

Коммивояжер (бродячий торговец) должен выйти из первого города, посетить по разу в неизвестном порядке города 2, 1, 3..n и вернуться в первый город. Расстояния между городами известны. В каком порядке следует обходить города, чтобы замкнутый путь (тур) коммивояжера был кратчайшим?

Чтобы привести задачу к научному виду, введём некоторые термины. Итак, города перенумерованы числами jÎ Т=(1, 2, 3..n). Тур коммивояжера может быть описан циклической перестановкой t=(j₁, j₂,.., j_n, j₁), причём все j₁..j_n – разные номера; повторяющийся в начале и в конце j₁, показывает, что перестановка зациклена. Расстояния между парами вершинС_ijобразуют матрицу С. Задача состоит в том, чтобы найти такой тур t, чтобы минимизировать функционал:

Задача о производстве красок. Имеется производственная линия для производства n красок разного цвета; обозначим эти краски номерами 1, 2… n. Всю производственную линию будем считать одним процессором.. Будем считать также, что единовременно процессор производит только одну краску, поэтому краски нужно производить в некотором порядке Поскольку производство циклическое, то краски надо производить в циклическом порядке p=(j₁, j₂,.., j_n, j₁). После окончания производства краски i и перед началом производства краски jнадо отмыть оборудование от краски i. Для этого требуется время C[i, j]. Очевидно, что C[i, j] зависит как от i, так и от j, и что, вообще говоря, C[i, j]≠ C[j, i]. При некотором выбранном порядке придется на цикл производства красок потратить время

Где t_k- чистое время производства k-ой краски (не считая переналадок). Однако вторая сумма в правой части постоянна, поэтому полное время на цикл производства минимизируется вместе с общим временем на переналадку.

Таким образом, ЗК и задача о минимизации времени переналадки – это просто одна задача, только варианты ее описаны разными словами.

Задача о дыропробивном прессе. Дыропробивной пресс производит большое число одинаковых панелей – металлических листов, в которых последовательно по одному пробиваются отверстия разной формы и величины. Схематически пресс можно представить в виде стола, двигающегося независимо по координатам x, y, и вращающегося над столом диска, по периметру которого расположены дыропробивные инструменты разной формы и величины. Каждый инструмент присутствует в одном экземпляре. Диск может вращаться одинаково в двух направлениях (координата вращения z). Имеется собственно пресс, который надавливает на подвешенный под него инструмент тогда, когда под инструмент подведена нужная точка листа.

Операция пробивки j-того отверстия характеризуется четверкой чисел (xj, yj, zj, tj),, гдеxj, y_j- координаты нужного положения стола, z_j - координата нужного положения диска и t_j- время пробивки j-того отверстия.

Производство панелей носит циклический характер: в начале и конце обработки каждого листа стол должен находиться в положениях (x₀, y₀) диск в положении z₀причем в этом положении отверстие не пробивается. Это начальное состояние системы можно считать пробивкой фиктивного нулевого отверстия. С параметрами (x₀, y₀, z₀, 0).

Чтобы пробить j-тое отверстие непосредственно после i-того необходимо произвести следующие действия:

1. Переместить стол по оси x из положения x_i в положение x_j, затрачивая при этом время t^(x)(|x_i-x_j|)=t_{i, j}^(x)

2. Проделать то же самое по оси y, затратив время t_{i, j}^(y)

3. Повернуть головку по кратчайшей из двух дуг из положения z_i в положение z_j, затратив время t_{i, j}^(z).

4. Пробить j-тое отверстие, затратив время t_j.

Конкретный вид функций t^(x), t^(y), t^(z) зависит от механических свойств прессаи достаточно громоздок. Явно выписывать эти функции нет необходимости

Действия 1-3 (переналадка сi-того отверстия j-тое) происходит одновременно, и пробивка происходит немедленно после завершения самого длительного из этих действий. Поэтому

С[i, j] = max(t^(x), t^(y), t^(z)).

Теперь, как и в предыдущем случае, задача составления оптимальной программы для дыропробивного пресса сводится к ЗК (здесь - симметричной).

Задача о ранце. Имеется «ранец» (машина, самолет), с известнойгрузоподъемностью. Заданнабор наименований предметов, которые желательно взять с собой. Указаны вес и «полезность» (объем) каждого предмета. Нужно выбрать набор предметов, загружаемых в ранец, так чтобы суммарная полезность набора оказалась наибольшей. Задача имеет множество применений помимо собственно туристского. Например, нужно загрузить боевой самолет ракетами и бомбами так, чтобы нанести наибольший урон противнику. Или загрузить транспортное средство так, чтобы доставить груз наибольшей прибыльности.

Рассмотрим классическую задачу коммивояжера для 40 городов на плоскости. Пусть каждый отдельный город будет представлен целым положительным числом, а хромосома будет состоять из 40 не повторяющихся номеров, символизирующих порядок обхода. Вынесем координаты каждого города в отдельную матрицу locations размерностью общее количество городов на количество координат города, т.е. 40х2. Тогда номер города будет соответствовать ее строке. Готовый инструментарий генетического алгоритма не позволяет решить такую задачу без изменений по следующим причинам:

· Стандартная функция инициализации популяции не отслеживает уникальность каждого гена в хромосоме.

· При использовании встроенной функции кроссовера может возникнуть ситуация, при которой часть городов не попадет в потомка, а часть удвоится.

· Функция мутации должна лишь переставлять порядок обхода городов.

В связи с этим необходимо написать собственные функции мутации, кроссовера и создания популяции. Опишем алгоритм их работы.

Мутация – будет осуществлять простой обмен двух генов в хромосоме местами.

Кроссовер – выделяет часть хромосомы и и получает новую путем отражения данной части.

Создание поколения будет проходить случайным образом с использованием генератора случайных чисел и последующего их округления до целых.

Не смотря на то, что функции мутации и кроссовера получаются достаточно примитивными, они вполне подойдут для успешно решения задачи.

Функция пригодности должна рассчитывать общее сумму расстояний между городами, так как набор городов постоянен, то имеет смысл заранее рассчитать расстояния между ними, которое затем просто складывать. Из курса геометрии известно, что расстояние между 2-мя точками рассчитывается по след формуле: . Сохраним их в матрице distances. И реализуем функцию пригодности (см. приложение).

Далее остается добавить возможность графической интерпретации найденного решения и установить опции в соответствии с новыми функциями:

Fitness function: @(x) fitness(x, distances);

Creation function: @create

Crossover function: @crossover

Mutation function: @mutate

Plot function: @my_plot

Generations: 500

Numberofvariables: 40

После запуска получим следующие графики:

Рисунок 5.1 - Первоначальные данные (слева), кратчайший путь (справа)

Рисунок 5.2 - График средних и лучших значений функции приспособленности в каждом поколении

Т.е. в результате работы генетического алгоритма расстояние сократилось с 12000 до 3627.

Вопросы и задания к защите:

1. По какому алгоритму работает функция мутации? Кроссовера?

2. Как создается начальная популяция?

3. Каково возможное количество решений данной задачи?

4. Сколько времени потребовалось бы для решения аналогичной задачи методом перебора?

5. Имеет ли смысл использовать в данных задачах гибридные методы?

6. Какие еще методы существуют для решения задач подобных коммивояжеру? В чем их преимущества и недостатки?

7. Какие методы отбора в генетическом алгоритме вы знаете?

8. При каком методе отбора обязательно будут выживать лучший или лучшие члены популяции?

9. В чем различие между одноточечным и двухточечным кроссовером?

10. Назовите параметры скрещивания, которые настраиваются перед началом работы генетического алгоритма.

11. В чем основное различие между целевой функции и функцией приспособленности?

12. Какова причина того, что решения всегда разные?

13. Какие типы задач, решаются с помощью генетических алгоритмов?

Практическая работа №6 Оптимизации отклика нелинейных систем MATLAB Simulink

Рассмотрим применение с пакетом расширения Simulink некоторых наиболее важных пакетов расширения, входящих в состав инструментального ящика BlockSet.

Для динамической оптимизации систем управления в ходе их имитационного моделирования служит пакет прикладных программ построения нелинейных систем управления Nonlinear Control Design (NCD), входящий в состав инструментальных средств BlockSets Simulink 1. Он обеспечивает:

· легкую настройку переменных;

· указание неопределенных параметров систем;

· интерактивную оптимизацию;

· моделирование методом Монте-Карло;

· поддержку проектирования как одномерных, так и многомерных систем управления;

· моделирование подавления помех;

· моделирование процессов слежения;

· •моделирование объектов с запаздыванием и решение других задач управления.

Главное назначение пакета – оптимизация отклика нелинейных систем при наличии заданных ограничений на переходные процессы в таких системах.

В версиях MATLAB 6.0-7.0 используются, соответственно, версии данного пакета 1.0/6.0/7.0. Однако обновления этого пакета шли в основном в части совершенствования алгоритмов оптимизации. Поэтому излагаемый ниже материал в равной мере относится ко всем версиям этого пакета. В версии MATLAB 7 + Simulink 6 и в последующих версиях пакет Nonlinear Control Design переименован в Simulink Response Optimization. Особенности этого пакета расширения описаны в конце данного урока.

Состав блоков пакетов

Набор блоков пакетов оптимизации моделируемых нелинейных систем содержит всего три блока: CRMS (Continue RMS), DRMS (Discrete RMS) и NCD Output (в Simulink Response Optimization блок Signal Constraint). Блок CRMS реализует математическую зависимость

где u (t) – входной сигнал блока, y (t) – его выходной сигнал. Для стационарных случайных процессов с нулевым математическим ожиданием выходной сигнал блока при t → ∞ является среднеквадратическим отклонением.

Блок DRMS, по сути, реализует ту же зависимость, что и блок CRMS, но для сигналов, определенных в дискретные моменты времени:

.

Рассматриваемые блоки могут применяться, в частности, в системах моделирования, где качество функционирования целесообразно оценивать интегральным квадратичным критерием или стандартным отклонением ошибки.

Блок NCD Output является основным в рассматриваемом наборе блоков. Он имеет свое рабочее окно и меню и позволяет в интерактивном режиме выполнить следующие операции:

• задать требуемые ограничения во временной области на любой сигнал оптимизируемой системы;

• указать параметры, подлежащие оптимизации;

• указать неопределенные параметры;

• провести параметрическую оптимизацию системы с учетом заданных ограничений.

10.1.3. Демонстрация работы блоков пакета NCD

Для запуска демонстрационного примера в режиме командной строки MATLAB введем команду: > > rmsdemo

В результате выполнения этой команды появится рабочее окно Simulink с моделью, содержащей источник синусоидального сигнала с единичной амплитудой, к выходу которого подключены блоки CRMS и DRMS (рис. 10.1). С помощью блока Mux их выходы подключаются к входу виртуального осциллографа.

Активизируем блок осциллографа Scope (двойным щелчком мыши) и запустим процесс моделирования. Его результат отображен на осциллограммах рис. 10.1, где сплошная линия - выход блока CRMS, а ступенчатая - выход блока DRMS. Они также отличаются цветом линий. С течением времени сигналы на выходах обоих блоков стремятся к одному и тому же установившемуся значению – действующему значению синусоиды с единичной амплитудой, равному 1/V2 ^0.707. Для блока CMRS это приближение получается плавным, а в случае блока DRMS - дискретным.

10.2. Оптимизация нелинейных систем с помощью пакета NCD

Настройка параметров PID-регулятора

Пример использования блоков NCD Blockset содержится в файле ncddemo1. Запустив этот файл, в окне Simulink мы увидим модель системы (рис. 4.40).

Рисунок. 1.40. Моделирование PID-регулятора

Основными элементами замкнутой системы PID-регулятора являются:

· объект регулирования (блок Plant & Actuator);

· PID-регулятор (Controller) на основе дифференцирующего устройства и интегратора;

· цепь обратной связи и узел сравнения.

Кроме того, в модель входит блок задающего воздействия (в виде единичного скачка) Step и блок NCD Output (NCD OutPort 1). Раскроем блок Plant & Actuator двойным щелчком мыши. В состав объекта регулирования (см. его окно на рис. 10.8 справа) входит нелинейность с уровнями ограничения –2, 2 (блок Limit), блок динамического ограничения коэффициента усиления (Rate), осуществляющий ограничение величины

диапазоном [-0.8, 0.8], а также линейное динамическое звено предаточной функции.

где коэффициент a ₂ может принимать значения в диапазоне [40, 50] с номинальным значением a ₂ = 43, а коэффициент a ₁ – в диапазоне от [0, 5, 3] с номинальным значением a ₁ = 1, 5.

Постановка задачи оптимизации в данном случае такова: при заданной структуре объекта управления и известных неопределенностях его параметров найти значения коэффициентов Kp, Ki и Kd регулятора, при которых в представленной замкнутой структуре переходный процесс будет иметь параметры, заданные по умолчанию.

Отметим, что моделирование в данном случае будет проводиться для номинальных значений коэффициентов a ₁ и a ₂ и при начальных значениях параметров оптимизации Kp = 0, 63, Ki = 0, 0504, Kd = 1, 96810. Такие значения выбраны в соответствии с методикой настройки ПИД-регуляторов Зиглера–Николса (Ziegler-Nichols method), согласно которой:

коэффициенты Ki и Kd устанавливаются равными нулю, а коэффициент Kp увеличивается до тех пор, пока система не потеряет устойчивость;

предельное значение Kp обозначается как Ku, а период автоколебаний – как Pu;

задаются следующие значения коэффициентов регулятора: Kp = 3·Ku/5, Ki = 6·Ku/(5ЧPu), Kd = 3·Ku·Pu/40.

Запустим процесс оптимизации, нажав кнопку Start. Полученный результат отражен на рис. 10.10. Найденные оптимальные значения (в командном режиме работы MATLAB) при этом равны:

> > Kp

Kp = 1.3366

> > Ki

Ki = 0.1548

> > Kd

Kd = 10.3318

и существенно отличаются от начальных.

Проведем теперь оптимизацию, учитывая интервальную неопределенность коэффициентов a ₁ и a ₂, то есть установив сначала в окне задания неопределенных переменных флажки Constrain lower simulation и Constrain upper simulation, а затем повторно запустив процесс моделирования. Результат для этого случая показан на рис. 1.41.

Рисунок. 1.41. Иллюстрация процесса оптимизации при неопределенных параметрах a1 и a2

В командном режиме вычислений находим:

> > Kp

Kp = 1.6740

> > Ki

Ki = 0.1273

> > Kd

Kd = 10.2597.

Оптимизация регулятора привода

По выше приведенным методом была проведена оптимизация модели привода. Для сравнения переходных характеристик была составлена Simulink модель, в которую входят две модели привода (рисунок 4.42). Они различаются настройками регуляторов, первый регулятор настроен в ручную, другой настраивается при помощи блока из пакета NCD.

Рисунок.1.43. Simulink модель привода

Была произведена атематическая настройка регулятора скорости, путем задания границ приходного процесса (рисунок 1.44)

Рисунок 1.44. Процесс подпорки параметров регулятора.

После завершения моделирования выводится числовые параметры регуляторов (рис 1.45 а). После проведение моделирования можно констатировать, что переходный процесс с автоматически настроенным регулятором имеет меньшую калебательность (рис 1.45 б).

Рисунок 1.45. Итоги моделирования

Инициализация начальных данных:

%%

% Путешественникпосещаетгорода в России. Файл

% russia.matсодержиткоординатыграницстраны в переменных x и y,

% и геометрическиупрощенныекоординаты в переменныххх и уу.

load('russia.mat', 'x', 'y', 'xx', 'yy');

plot(x, y, 'Color', 'red'); hold on;

%%

% Сгенерируемслучайнымобразомрасположениегородоввнутриполигона

% имеющегоформуграницРоссии. Дляэтоговоспользуемсяфункцией INPOLYGON,

% котораяпозволяетубедитсялежитлиданнаяточкавнутриполигона.

cities = 40;

locations = zeros(cities, 2);

n = 1;

while (n < = cities)

xp = rand*1200;

yp = rand*900;

ifinpolygon(xp, yp, xx, yy)

locations(n, 1) = xp;

locations(n, 2) = yp;

n = n+1;

end

end

plot(locations(:, 1), locations(:, 2), 'bo');

%%

% Синиекружкипредставляютсобойгородакоторыепутешественникдолженпосетить.

% Имеяматрицу locations, представляющуюсобойкоординатывсехгородов,

% заранеерассчитаемматрицу distance представляющуюрасстояниямеждугородами.

distances = zeros(cities);

for count1=1: cities,

for count2=1: count1,

x1 = locations(count1, 1);

y1 = locations(count1, 2);

x2 = locations(count2, 1);

y2 = locations(count2, 2);

distances(count1, count2)=sqrt((x1-x2)^2+(y1-y2)^2);

distances(count2, count1)=distances(count1, count2);

end;

end;

%%

% Инструментарий GA вызываетфункциюпригодноститолько с однимпараметром 'x',

% нонашафункцияпригодностиимеет 2 аргумента: x, distances.

% Дляпередачиейвторогопараметра - матрицы 'distances'

% воспользуемсяанонимнойфункцией.Таккакдистанциимыужеопределили

% и онинаходятся в матрице distances в рабочемпространствеMatLab'а,

% тоанонимнаяфункциябудетполучатьаргумент 'x', и простовызывать

% нашуфункциюпригодностистворомпараметром 'distances'.

FitnessFcn = @(x) fitness(x, distances);

%%

% Аналогичнымобразомвоспользуемсяанонимнойфункцией и дляпостроения

% гарфика, отражающегопорядокобходагородов, полученного с помощью

% лучшейтекущейособи.Здесьнам в дополнениепотребуетсяматрица locations.

my_plot = @(options, state, flag) tr_plot(options,...

state, flag, locations);

%% Настройкигенетическогоалгоритма

% Сначалаопределимтиппопуляции и интервалначальныхзначенийгенов

options = gaoptimset('PopulationType', 'custom', 'PopInitRange',...

[1; cities]);

%%

% Сделаемнастройки в соответствии с написанныминамифункциямимутации,

% скрещивания, создания и отображения, а такжеизменим

% некоторыекритерииустановки

options = gaoptimset(options, 'CreationFcn', @create,...

'CrossoverFcn', @crossover,...

'MutationFcn', @mutate,...

'PlotFcn', my_plot,...

'Generations', 500, 'PopulationSize', 60,...

'StallGenLimit', 200, 'Vectorized', 'on');

%%

% Теперьможнозапуститьгенетическийалгоритмнарешениезадачи

numberOfVariables = cities;

[x, fval, reason, output] = ga(FitnessFcn, numberOfVariables, options)

Функциясоздания:

function pop = create(NVARS, FitnessFcn, options)

% NVARS – количествопеременныхфункцийпригодности

totalPopulationSize = sum(options.PopulationSize); %определяемразмерпопуляции

n = NVARS; %определяем количество переменных функции пригодности

pop = zeros(totalPopulationSize, n); %создаемпустойвектор

for i = 1: totalPopulationSize

pop(i,:) = randperm(n); % изаполняемегослучайнымичислами

end

Функцияскрещивания:

functionxoverKids= crossover(parents, options, NVARS,...

FitnessFcn, thisScore, thisPopulation)

nKids = length(parents)/2; %детейбудетв 2 разаменьшеродителей

xoverKids = zeros(nKids, NVARS); %новыйвектор

index = 1;

for i=1: nKids

parent = thisPopulation(parents(index),:);

index = index + 2; % каждый второй родитель

% кроссовер будет заключатся в перевороте секции одного родителя

p1 = ceil((length(parent) -1) * rand); %выбираемконечнуюточку

p2 = p1 + ceil((length(parent) - p1- 1) * rand); %начальную

child = parent;

child(p1: p2) = fliplr(child(p1: p2)); %ипереворачиваемполученнуюсекцию

xoverKids(i,:) = child;

end

Функциямутации:

functionmutationChildren = mutate(parents, options, NVARS,...

FitnessFcn, state, thisScore, thisPopulation, mutationRate)

mutationChildren = zeros(length(parents), NVARS);

for i=1: length(parents)

parent = thisPopulation(parents(i),:);

p = ceil(length(parent) * rand(1, 2)); %получаем случайно заполненную матрицу 1х2 приводим в масштаб кол-ва переменных и округляем до целых

child = parent;

child(p(1)) = parent(p(2)); %производимобменчастями

child(p(2)) = parent(p(1));

mutationChildren(i,:) = child;

end

Функцияпригодности:

functionscores = fitness(x, distances) %функция пригодности вызывается с уже заполненной матрицей дистанций для всех городов.

scores = zeros(length(x), 1);

for j = 1: size(x, 1)

p = x(j,:);

f = distances(p(end), p(1));

for i = 2: length(p)

f = f + distances(p(i-1), p(i));

end

scores(j) = f;

end

Список использованной литературы

1. Батыршин И.З., Недосекин А.О., Стецко А.А., Тарасов В.Б., Язенин А.В., Ярушкина Н.Г. Нечеткие гибридные системы. Теория и практика / Под редакцией Н.Г. Ярушкиной. – М.: Физматлит, 2007. – 208с.
2. Заде Л.А. Понятие лингвистической переменной и его применение к принятию приближенных решений. – М.: Мир, 1976.
3. Леоненков А.В. Нечеткое моделирование в среде MATLAB и fuzzyTECH. – СПб.: БХВ-Петербург, 2005. – 736с.: ил.
4. Ярушкина Н.Г. Основы теории нечетких и гибридных систем. – М.: Финансы и статистика, 2004. – 320с.: ил.