Введение
С появлением библиотеки SymPy для решения математических задач появились дополнительные возможности, позволяющие отображать результаты в символьной форме.
Подробное описание использования символьных вычислений приведено в публикации [1] под названием «Введение в научный Python» в разделе «Символьные вычисления».
Расширение области применения символьных вычислений на решение отдельных задач нелинейного программирования надеюсь будет способствовать популяризации Python в том числе и как альтернатива дорогостоящих математических пакетов.
Постановка задачи
Привести примеры символьных вычислений для безусловного экстремума дифференцируемой нелинейной функции цели с определением достаточных условий существования экстремума по матрице Гессе. Рассмотреть так же задачу условного нелинейного программирования с линейными ограничениями при помощи множетелей Лагранжа.
Для того, чтобы определиться с терминологией приведу следующее определение [2]. Задачей нелинейного программирования (задачей НП) называется задача нахождения максимума (минимума) нелинейной функции многих переменных, когда на переменные имеются (не имеются) ограничения типа равенств или неравенств.
Символьное вычисление безусловного экстремума дифференцируемой функции трёх переменных
Несмотря на сложность решаемых задач при символьном решении всё становится простым и наглядным. Рассмотрим листинг первого примера.
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
f=-x**2-5*y**2-3*z**2+x*y-2*x*z+ 2*y*z+11 *x+2*y+18*z+10
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')
try:
sols=solve([fx,fy,fz],x,y,z)
except:
print(' Функция не дифференцируема')
raise SystemExit(1)
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)
fyy=f.diff(y,y).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fyy=',fyy)
fzy=f.diff(z,y).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fyz=',fzy)
fzz=f.diff(z,z).subs({x:sols[x],y:sols[y],z:sols[z]})
print('fzz=',fzz)
fyx=fxy;fzx=fxz;fyz=fzy# равенства из условия симметричности матрицы Гессе
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])))
Результат
Анализируемая функция f для переменных x,y,z — f= -x**2 + x*y — 2*x*z + 11*x — 5*y**2 + 2*y*z + 2*y — 3*z**2 + 18*z + 10
Необходимые условия экстремума
df/dx = -2*x + y — 2*z + 11 =0
df/dy = x — 10*y + 2*z + 2 =0
df/dz = -2*x + 2*y — 6*z + 18 =0
Стационарные точки M (x, y, z) — {x: 4, y: 1, z: 2}
Вторые производные в точке М
fxx= -2
fxy= 1
fxz= -2
fyy= -10
fyz= 2
fzz= -6
Расчёт определителей матрицы Гессе
(«разрастаются» из левого верхнего угла)
Определитель М1- d1= -2
Матрица М2- Matrix([[-2, 1], [1, -10]])
Определитель М2- d2= 19
Матрица М3- Matrix([[-2, 1, -2], [1, -10, 2], [-2, 2, -6]])
Определитель М3- d3= -74
Достаточные условия экстремума
При d1=-2,<0, d2=19>0, d3=-74<0, максимум f в точке М(4,1,2)
Значение 51 функции в точке М (4,1,2)
В программе много пояснений, но это и является положительным моментом как при обучении так и при исследовании. Количество переменных в функции цели можно как увеличить, так и уменьшить как это предлагается в публикации [3].
Программа определяет экстремум минимум или максимум, а также «седло». Если функция цели не дифференцируемая программа выводит сообщение и прекращает работу. Например, изменив функцию цели получим:
Анализируемая функция 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
fyy= 4
fyz= 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)
Символьное определение условного экстремума дифференцируемой функции по методу множителей Лагранжа
Продемонстрируем метод Лагранжа на примере задачи взятой из публикации [4] по которой можно сравнить результаты:
Имеется два способа производства некоторого продукта. Издержки производства при каждом способе зависят от произведенных x и y следующим образом: g(x) = 9*x + x**2, g(y)=6*y + y**2. За месяц необходимо произвести 3×50 единиц продукции, распределив ее между двумя способами так, чтобы минимизировать общие издержки.
Рассмотрим символьное решение задачи, которое и объясняет метод.
from sympy import *
x,y,w=symbols(' x y w' )
g=9*x+x**2+6*y+y**2
print('Анализируемая функция f для переменных x,y :n f= ', g)
q=x+y-150
print('Ограничения: ', q,'=0')
f=9*x+x**2+6*y+y**2+w*(x+y-150)
print('Вспомогательная функция Лагранжа :n ',f)
fx=f.diff(x)
print('df/dx =',fx,'=0')
fy=f.diff(y)
print('df/dy =',fy,'=0')
fw=f.diff(w)
print('df/dw =',fw,'=0')
sols=solve([fx,fy,fw],x,y,w)
print('Стационарная точка M(x,y) :n',float(sols[x]),',',float(sols[y]))
Результат
Анализируемая функция f для переменных x,y:
f= x**2 + 9*x + y**2 + 6*y
Ограничения: x + y — 150 =0
Вспомогательная функция Лагранжа:
w*(x + y — 150) + x**2 + 9*x + y**2 + 6*y
df/dx = w + 2*x + 9 =0
df/dy = w + 2*y + 6 =0
df/dw = x + y — 150 =0
Стационарная точка M(x,y):
74.25, 75.75
Вместо выводов
Результат совпадает с приведенным в [4], но главное, по моему мнению, другое. На примере символьных вычислений легче понять любой метод и им воспользоваться. Но для математических пакетов мысль эта не новая и в Mathcad и в Matlab символьные вычисления давно и широко используются. Но это же на Python!
Спасибо всём, кто прочитал статью.
Ссылки
- Введение в научный Python.
- Тема 6. Нелинейное программирование.
- Экстремумы функций двух и трёх переменных.
- Метод множителей Лагранжа. Пример решения.
Автор: Юрий Тараненко
