Нелинейное программирование средствами Python (Часть 1. Безусловное нелинейное программирование) -2


Введение


Закончив серию статей по линейному программированию средствами Python [1,2,3,4] логично перейти к рассмотрению задач нелинейного программирования. Следует заметить, что тема эта для нашего ресурса не новая. В публикации [5] рассмотрен алгоритм который заключается в формировании симплекса (simplex) и последующего его деформирования в направлении минимума, посредством трех операций. В публикации не использована функция optimize из модуля scipy. optimize популярной библиотеки для языка python, что не позволяет сделать вывод о сходимости алгоритма. В публикации [6] рассмотрен квази ньютоновский метод для не выпуклой функции, которая имеет непрерывные вторые производные. При этом автор утверждает, что неоспоримое преимущество алгоритма, конечно, состоит в том, что нет необходимости вычислять вторые производные.

Публикация [7] содержит описание квадратичной модели без линейных ограничений.

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

Постановка задачи нелинейного программирования


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

В обоих случаях нужно построить такую математическую модель, в которой на каждом этапе преобразования можно получить символьный результат в интерфейсах idle и jupyter-notebook и провести сравнение конечного результата с результатом, полученным с использованием модуля scipy optimize minimize.

Стандартная форма задачи нелинейного программирования


В стандартной форме задачу НП можно записать в следующем виде:

(1)
при ограничениях:

Если ввести обозначение:
то задача НП (1) примет более простой вид: (2)
при ограничениях:

Множество – будем называть множеством допустимых решений задачи НП. Задачу (2) можно переписать следующим образом:

Допустимое решение называется оптимальным решением задачи НП, если

Заметим, что последняя запись означает, что для любого , т.е. наибольшее значение функции f на множестве M. Число будем называть значением задачи НП.

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

Выпуклые и вогнутые функции


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

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

Пусть выпуклое множество вещественная функция, заданная на множестве M.

Функция f называется выпуклой, если для любой пары точек и выполняется неравенство .
Функция f называется вогнутой, если для любой пары точек и выполняется неравенство .

Теорема о выпуклости множества допустимых решений


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

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

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

Безусловный метод нелинейного программирования


Если в задаче нелинейного программирования (1) нет ограничений на переменные, т. е. то такая задача называется задачей безусловной оптимизации и имеет вид (3)
или в векторной форме

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

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


Если – оптимальное решение задачи (3), то все частные производные функции f в точке равны нулю, т. е. (4)
Необходимое условие оптимальности вида (4) часто называют условием стационарной точки.

Замечание о достаточности условий оптимальности


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

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

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

Листинг программы
from sympy import *
x,y,z=symbols(' x y z' )
f=x**2+2*y**2+3*z**2-2*x*y-2*x*z
print('Анализируемая функция f  для переменных x,y,z -\n f= ', f )
print('Необходимые условия экстремума')
fx=f.diff(x)
print('df/dx =',fx,'=0')
fy=f.diff(y)
print('df/dy =',fy,'=0')
fz=f.diff(z)
print('df/dz =',fz,'=0')
sols=solve([fx,fy,fz],x,y,z)
print('Стационарные точки M(x,y,z) -',sols)
print('Вторые производные в точке М')
fxx=f.diff(x,x).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fxx=',fxx)
fxy=f.diff(x,y).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fxy=',fxy)
fxz=f.diff(x,z).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fxz=',fxz)
fyx=f.diff(y,x).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fyx=',fyx)
fyy=f.diff(y,y).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fyy=',fyy)
fyz=f.diff(y,z).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fyz=',fyz)
fzx=f.diff(z,x).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fzx=',fzx)
fzy=f.diff(z,y).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fzy=',fzy)
fzz=f.diff(z,z).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fzz=',fzz)
print('Расчёт определителей матрицы Гессе \n («разрастаются» из левого верхнего угла)')
d1=fxx
print('Определитель М1- d1=',d1)
M2=Matrix([[fxx,fxy],[fyx,fyy]])
print('Матрица М2-',M2)
d2=M2.det()
print('Определитель М2- d2=',d2)
M3=Matrix([[fxx,fxy,fxz],[fyx,fyy,fyz],[fzx,fzy,fzz]])
print('Матрица М3-',M3)
d3=M3.det()
print('Определитель М3- d3=',d3)
print ('Достаточные условия экстремума')
if  d1>0 and d2>0 and d3>0:
         print('При d1=%s,>0, d2=%s>0, d3=%s>0, минимум f в точке М(%s,%s,%s)'%(str(d1),str(d2),str(d3),str(sols[x]),str(sols[y]),str(sols[z])))
elif  d1<0 and d2>0 and d3<0:
         print('При d1=%s,<0, d2=%s>0, d3=%s<0,максимум f в точке М(%s,%s,%s)'%(str(d1),str(d2),str(d3),str(sols[x]),str(sols[y]),str(sols[z])))
elif  d3!=0:
         print('Седло в точке М(%s,%s,%s)'%(str(sols[x]),str(sols[y]),str(sols[z])))
else:
          print('Нет экстремума в точке М(%s,%s,%s)'%(str(sols[x]),str(sols[y]),str(sols[z])))
r=f.subs({x:sols[x],y:sols[y],z:sols[z]})          
print('Значение %s функции в точке М(%s,%s,%s)'%(str(r),str(sols[x]),str(sols[y]),str(sols[z])))         


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

Результат работы программы (в интерфейсах IDLE и jupyter-notebook результаты совпадают)


Анализируемая функция f для переменных x,y,z — f= x**2 — 2*x*y — 2*x*z + 2*y**2 + 3*z**2
Необходимые условия экстремума

df/dx = 2*x — 2*y — 2*z =0
df/dy = -2*x + 4*y =0
df/dz = -2*x + 6*z =0

Стационарные точки M(x,y,z) — {z: 0, y: 0, x: 0}
Вторые производные в точке М

fxx= 2
fxy= -2
fxz= -2
fyx= -2
fyy= 4
fyz= 0
fzx= -2
fzy= 0
fzz= 6

Расчёт определителей матрицы Гессе («разрастаются» из левого верхнего угла)
Определитель М1- d1= 2
Матрица М2- Matrix([[2, -2], [-2, 4]])
Определитель М2- d2= 4
Матрица М3- Matrix([[2, -2, -2], [-2, 4, 0], [-2, 0, 6]])
Определитель М3- d3= 8
Достаточные условия экстремума
При d1=2,>0, d2=4>0, d3=8>0, минимум f в точке М(0,0,0)
Значение 0 функции в точке М(0,0,0)

Программа позволяет контролировать все этапы определения экстремума и кроме максимума и минимума определяет «седло» и отсутствие экстремума по следующим условиям:

if  d1>0 and d2>0 and d3>0:
         print('При d1=%s,>0, d2=%s>0, d3=%s>0, минимум f в точке М(%s,%s,%s)'%(str(d1),str(d2),str(d3),str(sols[x]),str(sols[y]),str(sols[z])))
elif  d1<0 and d2>0 and d3<0:
         print('При d1=%s,<0, d2=%s>0, d3=%s<0,максимум f в точке М(%s,%s,%s)'%(str(d1),str(d2),str(d3),str(sols[x]),str(sols[y]),str(sols[z])))
elif  d3!=0:
         print('Седло в точке М(%s,%s,%s)'%(str(sols[x]),str(sols[y]),str(sols[z])))
else:
          print('Нет экстремума в точке М(%s,%s,%s)'%(str(sols[x]),str(sols[y]),str(sols[z])))

Сравним эти результаты с результатами, полученными при помощи scipy optimize. Сначала применим метод простого перебора.

from scipy import optimize:
def  f(x):
         return x[0]**2 +2* x[1]**2+3*x[2]**2-2*x[0]*x[1]-2*x[0]*x[2]
result=optimize.brute(f,((-5,5),(-5,5),(-5,5)))
print(result)
Получим
[ 3.11324802e-05   1.21277436e-05   1.62810432e-05]

Метод даёт неплохие результаты, но иногда интересно узнать о количестве итераций. Используем генетический алгоритм дифференциальная эволюция:

from scipy import optimize
def  f(x):
         return x[0]**2 +2* x[1]**2+3*x[2]**2-2*x[0]*x[1]-2*x[0]*x[2]
print(optimize.differential_evolution(f,((-5,5),(-5,5),(-5,5))))
Получим
     fun: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 8914
     nit: 197
 success: True
       x: array([ 0.,  0.,  0.])

Результаты отличные, но большое количество итераций, поэтому ещё рассмотрим метод optimize.minimize не требующий дополнительных параметров:

 from scipy import optimize
def  f(x):
         return x[0]**2 +2* x[1]**2+3*x[2]**2-2*x[0]*x[1]-2*x[0]*x[2]
print(optimize.minimize(f,[2,2,2]))
Получим
     fun: 1.732026196804416e-14
 hess_inv: array([[ 2.9998198 ,  1.49965161,  0.99983175],
       [ 1.49965161,  0.99932612,  0.49967456],
       [ 0.99983175,  0.49967456,  0.49984284]])
      jac: array([  1.55507804e-07,  -2.38747034e-07,  -1.50008467e-07])
  message: 'Optimization terminated successfully.'
     nfev: 40
      nit: 6
     njev: 8
   status: 0
  success: True
        x: array([ -1.75716056e-07,  -1.54995367e-07,  -9.10240104e-08])

По числу итераций это лучший результат.

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

Вывод


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

Проведено сравнение результатов работы разработанной программы с результатами, полученными с помощью популярной библиотеки scipy optimize.

Ссылки


  1. Расширение аналитических возможностей метода линейного программирования средствами Python.
  2. Решение прямой и двойственной задачи линейного программирования средствами Python.
  3. Решение закрытой транспортной задачи с дополнительными условиями средствами Python.
  4. Решение задач линейного программирования с использованием Python.
  5. Метод оптимизации Нелдера — Мида. Пример реализации на Python.
  6. Метод BFGS или один из самых эффективных методов оптимизации. Пример реализации на Python.
  7. Метод оптимизации Trust-Region DOGLEG. Пример реализации на Python.




К сожалению, не доступен сервер mySQL