Министерство образования Республики Беларусь
Учреждение образования «Белорусский государственный университет информатики и радиоэлектроники»
Факультет информационных технологий и управления
Кафедра систем управления
ПОЯСНИТЕЛЬНАЯ ЗАПИСКА
к курсовому проекту по курсу:
«Математические основы теории систем»
Выполнил: Руководитель:
студент гр. 902402 Павлова А.В.
Потягевич Д.М.
Минск 2012
Содержание
1 МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ ЛИНЕЙНЫХ СИСТЕМ 3
1.1 Дифференциальное уравнение системы. Характеристическое уравнение и его корни. 3
1.2 Разложение передаточной функции на сумму простых слагаемых. Вычисление импульсной переходной характеристики ω(t) с помощью обратного преобразования Лапласа и переходной характеристики h(t) 4
1.3 Построение ЛАЧХ и ЛФЧХ 6
1.4 Уравнение состояния в нормальной форме,схема моделирования 8
1.5 Уравнение состояния в канонической форме, схема моделирования 10
1.6 Решение уравнения состояния в нормальной и канонической формах 13
1.7 Проверка: одинаково ли значение коэффициента усиления по передаточной функции, переходной характеристике, моделям в пространстве состояний, аналитической записи импульсной переходной характеристики 16
2 ЛИНЕЙНОЕ ПРОГРАММИРОВАНИЕ 18
2.1 Математическая модель задачи. Нахождение оптимального плана х* и экстремального значения функции 18
2.2 Построение и решение задачи, двойственной к исходной. Сравнение решения прямой и двойственной задач 20
2.3 Получение целочисленного решения путем введения дополнительных ограничений по методу Гомори 24
3 НЕЛИНЕЙНОЕ ПРОГРАММИРОВАНИЕ 26
3.1 Построение ОДЗП, выбор начальной точки поиска 26
3.2 Нахождение экстремального значения функции F(x) без учета ограничений на переменные 27
3.2.1 Метод наискорейшего спуска 27
3.2.2 Метод Ньютона-Рафсона 30
3.3 Нахождение экстремального значения функции F(x) с учетом системы ограничений задачи 32
3.3.1 Метод допустимых направлений Зойтендейка 32
3.3.2 Метод линейных комбинаций 35
3.3.3 Условия теоремы Куна-Таккера 40
4 ТЕКСТЫ ПРОГРАММ В СРЕДЕ MATLAB 44
4.1 Математическое описание линейных систем 44
4.2 Линейное программирование 48
4.3 Нелинейное программирование 49
1 Математическое описание линейных систем
1.1 Дифференциальное уравнение системы.Характеристическое уравнение и его корни
Передаточная функция системы имеет вид:
(1.1)
Передаточная функция системы W(s) –это отношение изображения по Лапласу сигнала на выходе к изображению по Лапласу сигнала на входе при нулевых начальных условиях.
Чтобы перейти от передаточной функции к дифференциальному уравнению системы, нужно перейти из области изображений по Лапласу во временную область. Из (1.1) следует:
Для перехода во временную область воспользуемся формальными правилами:
Тогда дифференциальное уравнение, соответствующее передаточной функции, имеет вид:
(1.2)
где - выходной сигнал, - входной сигнал.
Характеристическое уравнение системы определяется знаменателем передаточной функции:
Один из корней уравнения можно найти подбором, это будет λ1=-3, а затем понизить порядок уравнения и решить его:
Итак,
,
тогда корни характеристического уравнения имеют вид:
Передаточная функция системы в форме нулей и полюсов (zpkформе) имеет вид:
1.2 Разложение передаточной функции на сумму простых слагаемых. Вычисление импульсной переходной характеристики ω(t) спомощью обратного преобразования Лапласа и переходной характеристики h(t)
Получим разложение передаточной функции на сумму простых слагаемых:
Найдём a, b, c.
.
Следовательно,
Получим систему уравнений:
Решая эту систему, получим корни:
Передаточная функция:
Импульсная переходная характеристика w(t) – это процесс изменения сигнала на выходе при подаче на вход δ-функции. Ее можно найти в результате обратного преобразования Лапласа, примененного к каждому слагаемому передаточной функции.
В соответствии с таблицами соответствия , тогда:
Переходная характеристика h(t) – это процесс изменения сигнала на выходе при подаче на вход единичного ступенчатого воздействия. Ее можно вычислить следующим образом:
,
таким образом аналитическая форма переходной характеристики имеет вид:
.
Переходную характеристику можно также вычислить следующим образом: так как на вход подается единичное ступенчатое воздействие, а преобразование по Лапласу 1(t) это , то , получим такой же результат.
Выполним проверку:
-верно;
- верно.
1.3 Построение ЛАЧХ и ЛФЧХ
При определении частотных характеристик подразумевается, что на входе и выходе системы сигналы являются гармоническими.
Амплитудно-частотная характеристика (АЧХ) показывает, как изменяется отношение выходного сигнала к входному в зависимости от частоты.
Фазочастотная характеристика (ФЧХ) показывает изменение сдвига фаз между входным и выходным сигналами в зависимости от частоты.
ЛАЧХ строится в двойных логарифмических шкалах. По одной логарифмической оси откладывается круговая частота ω, по другой значение L(ω)=20lgK, выраженное в децибелах. Асимптотическая ЛАЧХ состоит из отрезков прямых линий с наклонами кратными 20дБ/дек.
Преобразуем передаточную функцию к следующему виду:
Теперь она представляет собой произведение трех апериодических и одного форсирующего звена с постоянными времени ; Коэффициент усиления К=8. Сопрягающие частоты звеньев равны
; ; ; .
Далее необходимо правильно разметить оси и отметить на оси ω сопрягающие частоты. ЛАЧХ приведена на рис. 1.1, а.
Рисунок 1.1 –Асимптотические ЛАЧХ (а) и ЛФЧХ (б)
Первая линия в области низких частот проводится через точку ω=1;L(ω)=20lgK. Через эту точку проводится первый наклон, и он равен , где - количество интегрирующих звеньев. Так как есть дифференцирующее звено (=-1), то первый наклон будет +1. Он идет с наклоном +1доL(ω)=20lgK=20lg8=18 частоты . Эта частота относится к апериодическому звену Следовательно, наклон изменится на -1. ЛАЧХ параллельна оси частот. Этот наклон будет идти до сопрягающей частоты . Так как это частота относится к апериодическому звену, то наклон изменится на -1. После частоты наклон изменится на -1и станет равным -2. Частота, при которой частотная характеристика пересечет ось частот, называется частотой среза,
Фазочастотная характеристика (рис. 1.1, б) построена в соответствии с выражением
Значения каждого из слагаемых определяются приближенно для значений ω → 0, ω → ∞, . В этих точках
1.4 Уравнение состояния в нормальной форме,схема моделирования
Кроме входных и выходных переменных при описании систем выделяют переменные х, связанные с внутренней структурой устройства,- переменные состояния. Тогда систему можно описать с помощью уравнений состояния.
Нормальная форма уравнения состояния имеет вид:
(1.3)
гдеА – матрица Фробениуса.
Матрица Фробениуса – это квадратная матрица, размер которой определяется порядком дифференциального уравнения. Элементы, стоящие над главной диагональю – единицы; элементы нижней строки - это коэффициенты левой части дифференциальногоуравнения, взятые с противоположным знаком, все остальные элементы – нулевые.
Согласно (1.2) дифференциальное уравнение системы имеет вид:
где и - коэффициенты уравнения.
Для систем с одним входом и одним выходом D – одноэлементная матрица, В – вектор-столбец, состоящий из 3 элементов, которые определяются следующим образом:
;
Матрица С – вектор-строка, состоящая из 3 элементов, первый элемент единица, остальные нули:
.
Подставим рассчитанные матрицы в систему (1.3), получим:
От матричной форме перейдем к скалярной:
Схема модели приведена на рис.1.2:
Рисунок 1.2 – Схема модели, соответствующая уравнениям состояния в нормальной форме
1.5 Уравнение состояния в канонической форме,
схема моделирования
Запишем уравнения состояния в канонической форме. Чтобы перейти к канонической форме, введем новую переменную q, которая связана с переменной состояния x следующим образом: .
М – модальная матрица, которая имеет вид:
где - характеристические числа матрицы Фробениуса А. Модальная матрица имеет такой вид, так как матрица А имеет форму Фробениуса и все корни характеристического уравнения различны.
При подстановке q вместо х в нормальную форму уравнений состояния (1.3) получим уравнения состояния системы в канонической форме:
(1.4)
Здесь - диагональная матрица:
где - матрица, обратная модальной, определяемая выражением:
.
Здесь -матрица, присоединенная к М, т.е. транспонированная матрица алгебраических дополнений.
Подставим найденные значения в (1.4), получим:
Схема модели, соответствующая полученной системе, приведена на рис.1.3. Для нее характерно параллельное соединение интеграторов, выходы которых определяются переменными состояния .
Рисунок 1.3 – Схема модели, соответствующая уравнениям состояния в канонической форме
Получим уравнения состояния в нормальной форме, с помощью Matlaba.
ss(W) – команда, выводящая матрицы A, B, C, D
Получим:
Рисунок 1.4 – Схема модели, соответствующая уравнениям состояния в нормальной форме
Рисунок 1.5 – Переходная характеристика
1.6 Решение уравнения состояния в нормальной и канонической формах
Найдем решение y(t) для системы уравнений в нормальной форме, если начальные условия имеют вид y(0)=1.8; . Сигнал u(t)=18∙1(t).переходя к начальным условиям для х, в соответствии с принятыми ранее обозначениями, получим:
Решение уравнения состояния складывается из двух составляющихx(t)=x1(t)+x2(t) – свободной и вынужденной.
Свободная составляющая – это общее решение дифференциального уравнения системы с нулевой правой частью. Оно не зависит от внешнего воздействия и характеризует естественное поведения системы.
Вынужденная составляющая – это частное решение дифференциального уравнения с ненулевой правой частью. Оно зависит от сигнала u(t) и характеризует поведения системы под его воздействием.
Решение уравнения состояния имеет вид:
где еAt – фундаментальная матрица или матрица перехода.
Она вычисляется по следующей формуле:
еAt= γ0E+ γ1A+ γ2A2,
гдеγ0, γ1, γ2– неизвестные коэффициенты.
Вычислить их можно, решая матричное уравнение:
Для рассматриваемого примера:
Перемножая матрицы, получаем систему уравнений следующего вида:
γ0 -3γ1+9γ2= е-3t,
γ0 -5γ1+25γ2= е-5t,
γ0 -10γ1+100γ2= е-10t.
Решение данной системы уравнений имеет вид:
γ0=3.571e-3t – 3e-5t+0.429е-10t,
γ1=1.071e-3t -1,3e-5t+0,229-10t,
γ2=0,071e-3t-0,1e-5t+0,029e-10t.
=>
Итак,
Так как y=x1, то свободная составляющая выходного сигнала будет равна 6.429e-3t– 5.4e-8t+0.771e-10t. Определим вынужденную составляющую при входном сигнале u(t)=18∙1(t). Сигнал на выходе при поступлении на вход 1(t) уже вычислен – это переходная характеристика системы (1.1). Чтобы получить вынужденную составляющую сигнала в нашем случае, умножим переходную характеристику на 18.
Таким образом, сигнал на выходе системы будет следующим:
y(t)=6.429e-3t– 5.4e-5t+0.771e-10t+18( )=
=1549.287e-3t-2165.4e-5t+617.919e-10t.
Выполним проверку:
y(0)=1549.287-2165.4+617.919=1.8 - верно;
y(∞)=0 - верно.
Найдем решение уравнений состояния, представленных в канонической форме (1.4). Каждое из дифференциальных уравнений первого порядка зависит только от одной переменной, и его решение в общем виду имеет вид:
Определим начальные условия q(0) для вектора q(t).
Так как , то
Найдем выражения для , и :
=
В результате получим:
Выполним проверку:
y(0)=1549.287-2165.4+617.914=1.8 - верно;
y(∞)=0 - верно.
Решения нормальных и канонических уравнений состояния совпадают.
1.7 Проверка: одинаково ли значение коэффициента усиления по передаточной функции, переходной характеристике,моделям в пространстве состояний, аналитической записи импульсной переходной характеристики
Проверим значение коэффициента усиления по:
передаточной функции:
переходной характеристике:
- верно.
моделям в пространстве состояний (в установившемся режиме на входах интеграторов нули и u=1):
- каноническая форма:
- нормальная форма:
К=0
аналитической записи импульсной переходной характеристики:
;проверяем:
.
Мы видим, что значение коэффициента усиления одинаково.
2 Линейное программирование
2.1 Математическая модель задачи. Нахождение оптимального плана х* и экстремального значения функции
Найти максимальное значение функции F(x)=-1x1-5x2-3x3 при следующих ограничениях:
{█(〖-5x〗_1 〖+4x〗_2 〖-1x〗_3=-3@〖-5x〗_1+〖0x〗_2+〖5x〗_3≥-9@〖0x〗_1-〖5x〗_2-〖4x〗_3≤-42@〖-4x〗_1 〖+5x〗_2+〖2x〗_3≤9@x_1,2,3≥0)┤
Домножим второе ограничение на (-1) и введем в ограничения дополнительные переменные x4, x5, x6, и искусственную переменную R следующим образом:
{█(〖-5x〗_1 〖+4x〗_2 〖-1x〗_3+R=-3@〖5x〗_1+〖0x〗_2-〖5x〗_3+x_4=9@〖0x〗_1-〖5x〗_2-〖4x〗_3+x_5=-42@〖-4x〗_1 〖+5x〗_2+〖2x〗_3+x_6=9)┤
Пусть R, x4, x5, x6 – базисные переменные, а x1, x2, x3 – небазисные. Функция цели
F(x)=F(x)+M∙∑▒R=〖-1x_1-5x〗_2-3x_3+M∙(-3+5x_1-〖4x〗_2+1x_3).
В первой симплекс-таблице (табл. 2.1) коэффициенты при небазисных переменных в F- и М-строках меняют знаки на противоположные, т.к. осуществляется максимизация функции. Свободный член в М-строке берется с противоположным знаком.
Таблица 2.1
БП Своб. члены НП
x1 x2 x3
R -3 -5 4 -1
x4 9 5 0 -5
x5 -42 0 -5 -4
x6 9 -4 5 2
F 0 1 5 3
M 3 5 -4 1
Решение, соответствующее таблице 2.1, не является допустимым, т.к. есть отрицательный свободный член.
Выберем ведущий столбец и строку: наибольший по модулю отрицательный свободный член находится в x5-строке, в этой строке наибольший по модулю отрицательный элемент соответствует столбцу x2, следовательно, столбец x2 – ведущий. На пересечении ведущей строки и ведущего столбца будет ведущий элемент. После пересчета получим таблицу 2.2.
Таблица 2.2
БП Своб. члены НП
x1 X5 x3
R -183/5 -5 4/5 -21/5
X4 9 5 0 -5
X2 42/5 0 -1/5 4/5
x6 -33 -4 1 -2
F -42 1 1 -1
M 183/5 5 -4/5 21/5
Сначала в новой таблице заполняются строка и столбец, которые в предыдущей таблице были ведущими. Элемент на месте ведущего равен обратой величине элемента предыдущей таблицы. Элементы строки делятся не ведущий элемент, а элементы столбца так же делятся на ведущий элемент, но берутся с противоположным знаком.
Пересчет остальных элементов производится по правилу прямоугольника: прямоугольник строится по старой таблице таким образом, что одну из его диагоналей образует пересчитываемый (a_ji) и ведущий (a_rl) элементы. Вторая диагональ определяется однозначно. Для нахождения нового элемента (a_ji^') из элемента a_ji вычитается произведение элементов противоположной диагонали:
a_jl^'=a_jl-(a_jl∙a_ri)/a_rl (2.1)
Далее по аналогии находим ведущий элемент и составляем симплекс-таблицы до тех пор, пока не получим допустимое и оптимальное решение (столбец свободных членов и F-строка не должны содержать отрицательных элементов):
Таблица 2.3 Таблица 2.4
БП Своб. члены НП
X5 X4
X1 24/5 -2/23 21/230
X3 3 -2/23 -5/46
X2 6 -3/23 2/23
X6 -39/5 11/23 17/115
F -219/5 1 -1/5
M 0 0 0
БП Своб. члены НП
X5 X3
X1 183/25 -4/25 21/25
X4 -138/5 4/5 -46/5
X2 42/5 -1/5 4/5
X6 -93/25 9/25 34/25
F -1233/25 29/25 -46/25
M 0 0 0
Для нахождения ведущей строки найдем максимальный по модулю отрицательный свободный член (-39/5). Ведущая строка - X6. В строке X6 так же найдем максимальный по модулю отрицательный свободный член (17/115). Столбец X4 - ведущий. Так как в строке с отрицательным свободным членом нет отрицательных элементов, то система ограничений не совместна и задача не имеет решения
2.2 Построение и решение задачи, двойственной к исходной. Сравнение решения прямой и двойственной задач
Рассмотрим соотношение прямой и двойственной задач:
(2.2)
Число переменных двойственной задачи совпадает с числом ограничений прямой задачи.
Исходная задача:
F(x)=-1x1-5x2-3x3 (max)при следующих ограничениях:
{█(〖-5x〗_1 〖+4x〗_2 〖-1x〗_3=-3@〖-5x〗_1+〖0x〗_2+〖5x〗_3≥-9@〖0x〗_1-〖5x〗_2-〖4x〗_3≤-42@〖-4x〗_1 〖+5x〗_2+〖2x〗_3≤9@x_1,2,3≥0)┤
Так как требуется найти максимум целевой функции, то неравенства в системе ограничений должны быть вида <. Второе неравенство ограничений прямой задачи умножим на (-1):
{█(〖-5x〗_1 〖+4x〗_2 〖-1x〗_3=-3@〖5x〗_1+〖0x〗_2-〖5x〗_3≤9@〖0x〗_1-〖5x〗_2-〖4x〗_3≤-42@〖-4x〗_1 〖+5x〗_2+〖2x〗_3≤9)┤
тогда
A=[■(-5&4&-1@5&0&-5@0&-5&-4@-4&5&2)] , B=[█(-3@9@-42@-9)], C^T=[■(-1&-5&-3)]
Двойственная задача будет иметь 4 переменные, так как прямая содержит 4 ограничения. В соответствии с (2.2) запишем двойственную задачу в виде:
F(y)=B^T Y=〖-3y〗_1+〖9y〗_2-42y_3+〖9y〗_4 (min),
Тогда условие A^T Y≥C пример следующий вид:
[■(-5&5&0&-4@4&0&-5&5@-1&-5&-4&2)]∙[█(y_1@y_2@y_3@y_4 )]≥[█(-1@-5@-3)],
следовательно:
{█(〖-5y〗_1+〖5y〗_2+〖0y〗_3-4y_4≥-1@〖4y〗_1-0y_2 〖-5y〗_3+5y_4≥-5@〖-1y〗_1-5y_2 〖-4y〗_3+〖2y〗_4≥-3@y_2,3,4≥0)┤.
Так как в прямой задаче есть ограничение равенство, то на у1, соответствующая переменная двойственной задачи, не накладывается условие неотрицательности.
Введем дополнительные переменные:
{█(〖5y〗_1-〖5y〗_2+〖0y〗_3+4y_4 〖+y〗_5=1@〖-4y〗_1-0y_2 〖+5y〗_3-5y_4 〖+y〗_6=5@〖1y〗_1+5y_2 〖+4y〗_3-〖2y〗_4 〖+y〗_7=3)┤
Составим симплекс-таблицу:
Таблица 2.5
БП Своб. члены НП
y1 y2 y3 y4
y5 1 5 -5 0 4
y6 5 -4 0 5 -5
y7 3 1 5 4 -2
F 0 -3 9 -42 9
Таблица 2.6
БП Своб. члены НП
Y1 Y2 Y3 Y7
Y5 -5 7 5 8 2
Y6 25/2 -13/2 -25/2 -5 -5/2
Y4 3/2 -1/2 -5/2 -2 -½
F -27/2 3/2 63/2 -24 9/2
Для нахождения ведущей строки найдем максимальный по модулю отрицательный свободный член (-5). Ведущая строка - X5. В строке X5 так же найдем максимальный по модулю отрицательный свободный член (2). Столбец X7 - ведущий.Так как в строке с отрицательным свободным членом нет отрицательных элементов, то система ограничений не совместна и задача не имеет решения Запишем соответствие между переменными прямой и двойственной задач:
Исходные переменные Дополнительные переменные
прямой задачи прямой задачи
x1 x2 x3 R x4 x5 x6
(2.3)
y5 y6 y7 y1 y2 y3 y4
Дополнительные переменные Исходные переменные
двойственной задачи двойственной задачи
Соответствие не означает равенство. Оптимальное решение прямой задачи определяется коэффициентами F-строки. Переменные прямой задачи приравниваются к коэффициентамприсоответствующим им небазисных переменных в F-строке оптимальной симплекс-таблицы двойственной задачи.
2.3 Получение целочисленного решения путем введения дополнительных ограничений по методу Гомори
При решении задачи не было получено решения, т.к. система ограничений не совместна.
3 Нелинейное программирование
3.1 Построение ОДЗП, выбор начальной точки поиска
Целевая функция имеет вид:
Построим ОДЗП:
Пусть x1=3, тогда x2=5
Пусть x1=0, тогда x2=-4
Пусть x1=3, тогда x2=5
Пусть x1=0, тогда x2=6
Рисунок 3.1 - ОДЗП
Внутри области допустимых значений выбираем точку x0, которая в дальнейшем будет являться начальной в процессе поиска экстремума:
x0=(1;3).
3.2 Нахождение экстремального значения функции F(x) без учета ограничений на переменные
3.2.1 Метод наискорейшего спуска
В методе наискорейшего спуска (подъема) очередная точка при поиске максимума функции вычисляется по формуле:
где направление движения задается вектором антиградиента функции F(x), вычисленном в точке , а величина шага перемещения определяется числовым параметром .
Найдем градиент :
.
На первом шаге движение осуществляется из точки вдоль вектора
- в новую точку :
Величина шага на любом шаге выбирается из условия обеспечения экстремума функции в рассматриваемом направлении. Подставляя координаты точки в функцию , получим:
Из условия:
;
найдем :
;
В результате после первого шага координаты очередной точки получаются равными:
Вычисляем
.
На втором шаге движение осуществляется в направлении вектора
- :
Подставив полученные выражения для x12 и x22 в функцию цели и преобразовав получаем
;
Из условия:
;
найдем :
;
Тогда:
В результате после второго шага координаты очередной точки получаются равными:
Вычисляем :
На третьем шаге движение осуществляется в направлении вектора
- :
;
Из условия:
;
найдем :
;
В результате после третьего шага координаты очередной точки получаются равными:
После проведения четвертой итерации получаем точку
На четвертой итерации закончим вычисления, значение функции цели:
Рисунок 3.2 – Графическая интерпретация метода наискорейшего спуска
Как видим из графической интерпретации, происходит приближение к экстремальному значению функции без учета ограничений.
3.2.2 Метод Ньютона-Рафсона
Условие задания:
=[1;3]
Данный метод дает решение задачи за 1 шаг. Очередная точка поиска вычисляется в соответствии с выражением:
где – матрица Гессе функции ; – обратная по отношению к матрица.
Градиент F(x):
;
.
где det H – определитель матрицы H; AdjH – присоединенная к H матрица (транспонированная матрица алгебраических дополнений).
Найдем определитель матрицы Гессе:
Найдем транспонированную матрицу алгебраических дополнений AdjH:
Теперь найдем матрицу обратную по отношению к - матрицу :
тогда:
Следовательно, в точке функция F(x) достигает максимального значения:
3.3 Нахождение экстремального значения функции F(x) с учетом системы ограничений задачи
3.3.1 Метод допустимых направлений Зойтендейка
Условие задания:
.
Тогда координаты очередной точки:
Здесь решение совпадает с первой итерацией метода наискорейшего спуска (Подъёма), тогда:
Определяем интервал допустимых значений для 0, при котором точка x1 будет принадлежать ОДЗП. Для этого подставим координаты точки x1 в ограничения задачи:
=>
Тогда:
Находим величину ,которая обеспечит экстремум функции F(x). Воспользуемся уже найденным = , но т.к. оно не входит в наш интервал, то = При этом очередная точка поисковой траектории оказывается на границе области. Координаты точки и значение градиента функции в этой точке определяются выражениями
Движение в направлении антиградиента выводит за пределы ОДЗП, поэтому очередную точку поиска вычисляем исходя из выражения:
где - новое направление, которое составляет минимальный острый угол с вектором градиента и направлено либо внутрь, либо по границе ОДЗП. При этом очередная точка должна принадлежать ОДЗП, а функция цели при переходе к очередной точке должна уменьшаться максимальным образом.
Направление находим, как решение задачи:
Найдем направление очередного шага: т.к. x1 лежит на , то условие (где - вектор коэффициентов при переменных в первом ограничении, на котором находится точка x1) запишется:
При движении из точки x1 в точку x2 следует двигаться по граничной прямой в направлении .
Координаты точки x2 определяются выражением:
или
Находим интервал изменения , при котором точка принадлежит ОДЗП, причем ограничение отбросим:
Получим интервал:
Найдем такое 1 , которое обеспечит максимум F(x) в направлении . Для этого координаты точки x2 подставляются в функцию F(x),тогда:
Значение 1 не принадлежит ранее найденному интервалу, поэтому для расчета координат точки принимается = :
Вычисляются составляющие вектора градиента в точке x2:
Вектор градиента не перпендикулярен вектору S1, но при выборе дальнейшего направления из условия:
И при движении вдоль оси ординат не происходит улучшения экстремального значения функции, следовательно, мы находимся в точке экстремума функции с учетом ограничений. А значение функции цели в этой точке равно:
Рисунок 3.3 - Графическая интерпретация метода допустимых направлений Зойтендейка
3.3.2 Метод линейных комбинаций
Условие задания:
Вычислим градиент функции F(x):
;
.
На следующем этапе вычислим значение градиента в точке x0:
Суть метода линейных комбинаций заключается в линеаризации функции F(x) и замене ее линейной функцией в соответствии с выражением:
Решаем задачу линейного программирования при следующих ограничениях:
Процедура решения задачи иллюстрируется следующей симплекс таблицей:
Таблица 3.1
Таблица 3.2
Получено оптимальное допустимое решение, которое имеет вид:
Произведем корректировку найденного решения в соответствии с выражением:
Найдем значение , которое доставляет экстремальное значение функции F(x1):
Определяем интервал допустимых значений для 0, при котором точка x1 будет принадлежать ОДЗП. Для этого подставим координаты точки x1 в ограничения задачи:
=>
Тогда:
Величина , не входит в наш интервал, тогда =1. Координаты точки и значение градиента функции в этой точке определяются выражениями
Линеаризуем функцию F(x) относительно точки x1 и заменим ее линейной функцией w(x1):
Решаем задачу линейного программирования при следующих ограничениях:
Процедура решения задачи иллюстрируется следующей симплекс таблицей:
Таблица 3.3
Таблица 3.4
Решение данной таблицы допустимо и оптимально, следовательно решение имеет вид:
Так как полученная точка является точкой предыдущего шага, следовательно, нет продвижения к точке экстремума, мы находимся в точке экстремума задачи с учетом ограничений.
Значение функции цели в этой точке равно:
Рисунок 3.4 - Графическая интерпретация метода линейных комбинаций.
3.3.3 Условия теоремы Куна-Таккера
Условие задания:
Составляем функцию Лагранжа:
здесь – левые части ограничений, приведенных к нулевой правой части; – неопределенные множители Лагранжа.
Точка экстремума является седловой точкой с минимумом по x и максимумом по , поэтому ограничения приведены к виду :
Условия теоремы Куна – Таккера записываем следующим образом:
Частные производные функции Лагранжа определяются выражениями:
Для того, чтобы вышеуказанные выражения имели вид равенств, введем в них дополнительные переменные:
Решение этой системы из четырех алгебраических уравнений, содержащих восемь неизвестных, можно найти с помощью симплекс-процедуры. На первом шаге в базис включаются все введенные дополнительные переменные. Строка для функции цели отсутствует.
Проведем симплекс – преобразования и найдем допустимое базисное решение. Если решение удовлетворяет условию оптимальности:
то оно является оптимальным.
Таблица 3.5
Таблица 3.6
Таблица 3.7
Полученное решение удовлетворяет следующим условиям:
Следовательно, координаты точки экстремума x*:
;
А функция цели в этой точке имеет значение:
4 Тексты программ в среде MATLAB
4.1 Математическое описание линейных систем
>> w=tf([1200 0],[1 18 95 150])
Transfer function:
1200 s
-------------------------
s^3 + 18 s^2 + 95 s + 150
>> pole(w)
ans =
-10.0000
-5.0000
-3.0000
>> zpk(w)
Zero/pole/gain:
1200 s
------------------
(s+10) (s+5) (s+3)
>> ss(w)
a =
x1 x2 x3
x1 -18 -5.938 -2.344
x2 16 0 0
x3 0 4 0
b =
u1
x1 8
x2 0
x3 0
c =
x1 x2 x3
y1 0 9.375 0
d =
u1
y1 0
Continuous-time model.
>> ch=[1200 0]
ch =
1200 0
>> zn=[1 18 95 150]
zn =
1 18 95 150
>> [x]=residue(ch,zn)
x =
-342.8571
600.0000
-257.1429
>> [c]=residue(ch,[zn,0])
c =
34.2857
-120.0000
85.7143
0
>>step(w)
Рисунок 4.1 – Переходная характеристика h(t)
>>impulse(w)
Рисунок 4.2 – Импульсная переходная характеристика w(t)
>>nyquist(w)
Рисунок 4.3 - АФЧХ системы
>>margin(w)
Рисунок 4.4 - ЛАЧХ и ЛФЧХ системы
>> M=[1 1 1; -3 -5 -10; 9 25 100]
M =
1 1 1
-3 -5 -10
9 25 100
>> inv(M)
ans =
3.5714 1.0714 0.0714
-3.0000 -1.3000 -0.1000
0.4286 0.2286 0.0286
>> B=[0;1200;-21600]
B =
0
1200
-21600
>> M^-1*B
ans =
-257.1429
600.0000
-342.8571
4.2 Линейное программирование
>> F=[1 5 3];
A=[5 0 -5;0 -5 -4; -4 5 2; -1 0 0; 0 -1 0 ; 0 0 -1 ];
B=[9;-42;9;0;0;0];
Aeq=[-5 4 -1];
Beq=[-3];
x=linprog(F,A,B,Aeq,Beq)
Exiting: One or more of the residuals, duality gap, or total relative error
has stalled:
the primal appears to be infeasible and the dual unbounded since
the dual objective > 1e+10
and the primal objective > -1e+6.
Решение двойственной задачи:
>> F=[-3 9 -42 9];
>> A=[5 -5 0 4; -4 0 5 -5; 1 5 4 -2; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1];
B=[1;5;3;0;0;0];
>> x=linprog(F,A,B)
Exiting: One or more of the residuals, duality gap, or total relative error
has grown 100000 times greater than its minimum value so far:
the dual appears to be infeasible (and the primal unbounded).
(The primal residual < TolFun=1.00e-008.)
4.3 Нелинейное программирование
Изображение фигуры
[x1,x2]=meshgrid([-5:0.1:5]);
F=6*x1.^2+6*x2.^2-3*x1.*x2+8*x2;
figure(1)
plot3(x1,x2,F)
grid on
figure(2)
contour(x1,x2,F)
grid on
figure(3)
meshc (x1,x2,F)
grid on
Рисунок 4.1 – Изображения функции цели
Построение линий уровня:
[x1,x2]=meshgrid([-2:0.1:9]);
F=6*x1.^2+6*x2.^2-3*x1.*x2+8*x2;
figure;
cl=[-1 15.5 130 300 500 ];
[c,h]=contour3(x1,x2,F,cl,'r');
clabel(c,h);
view(0,90);
Рисунок 4.2 - Линии уровня
Решение задачи нелинейного программирования без ограничений в Matlab:
F=inline('6*x(1)^2+6*x(2)^2-3*x(1)*x(2)+8*x(2)');
x0=[1 3];
options=optimset('LargeScale','off');
[x,fval,exitflag,output]=fminunc(F,x0,options);
x
fval
Local minimum found.
Optimization completed because the size of the gradient is less than
the default value of the function tolerance.
<stopping criteria details>
x =
-0.1778 -0.7111
fval =
-2.8444
Решение задачи нелинейного программирования с ограничениями в Matlab:
F=6*x1.^2+6*x2.^2-3*x1.*x2+8*x2;
A=[7 -4;1 4];
B=[-4 36];
x0=[1 3];
H=[12 -3;-3 15];
f=[0 8];
vlb=zeros(2,1);
vub=[];
x=quadprog(H,f,A,B,[],[],vlb,vub);
z=f*x+0.5*x'*H*x;
x
z
Warning: Large-scale algorithm does not currently solve this
problem formulation,
using medium-scale algorithm instead.
> In quadprog at 293
Optimization terminated.
x =
0
1.0000
z =
14.0000
