Symbolic solution of nonlinear programming problems

Introduction

With the advent of the SymPy library for solving mathematical problems, additional features have appeared that allow you to display the results in symbolic form.

A detailed description of the use of symbolic calculations is given in the publication [1] under the title “Introduction to Scientific Python” in the section “Character Computing”.

Expanding the scope of the use of symbolic computing to solve individual nonlinear programming problems, I hope will contribute to the popularization of Python, including as an alternative to expensive mathematical packages.

Formulation of the problem

To give examples of symbolic calculations for the unconditional extremum of a differentiable nonlinear function of the target with the determination of sufficient conditions for the existence of an extremum in the Hessian matrix. Consider also the problem of conditional nonlinear programming with linear constraints using Lagrange multipliers.

In order to determine the terminology, I will give the following definition [2]. The nonlinear programming task (NP problem) is the problem of finding the maximum (minimum) of the nonlinear function of many variables when there are (there are not) restrictions on the variables such as equalities or inequalities.

Symbolic calculation of the unconditional extremum of a differentiable function of three variables

Despite the complexity of the tasks to be solved with a symbolic solution, everything becomes simple and intuitive. Consider the listing of the first example.

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])))   

Result

The analyzed function f for the variables 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

Necessary conditions for the extremum

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

Stationary points M (x, y, z) - {x: 4, y: 1, z: 2}

Second derivatives at the point M

fxx = -2
fxy = 1
fxz = -2
fyy = -10
fyz = 2
fzz = -6

Calculation of the determinants of the Hessian matrix (“grow” from the upper left corner)

Determinant M1-d1 = -2
Matrix M2-Matrix ([[- 2 , 1], [1, -10]])
Determinant M2-d2 = 19
Matrix M3-Matrix ([[[- 2, 1, -2], [1, -10, 2], [-2, 2, - 6]])
Determinant M3 - d3 = -74

Sufficient conditions for the extremum
For d1 = -2, <0, d2 = 19> 0, d3 = -74 <0, maximum f at point M (4,1,2)
Value of 51 functions at point M (4,1,2)

There are a lot of explanations in the program, but this is a positive moment both in training and in research. The number of variables in the goal function can be either increased or decreased as suggested in the publication [3].

The program determines the extreme minimum or maximum, as well as the "saddle". If the target function is not a differentiable program, it displays a message and stops working. For example, changing the function of the target we get: The

analyzed function f for the variables x, y, z - f = x ** 2 - 2 * x * y - 2 * x * z + 2 * y ** 2 + 3 * z ** 2

Necessary conditions for the extremum

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

Stationary points M (x, y, z) - {z: 0, y: 0, x: 0} The

second derivatives at the point M

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

Calculation of the determinants of the Hessian matrix (“grow” from the upper left corner)

Determinant M1-d1 = 2
Matrix M2-Matrix ([[[2, -2], [-2, 4]]))
Determinant M2-d2 = 4
Matrix M3-Matrix ([[2, -2, -2], [-2, 4, 0], [- 2, 0, 6]]) The
determinant M3 is d3 = 8

Sufficient conditions for the extremum
For d1 = 2,> 0, d2 = 4> 0, d3 = 8> 0, the minimum f at the point M (0,0,0)

Value 0 functions at point M (0,0,0)

Symbolic determination of the conditional extremum of a differentiable function by the method of Lagrange multipliers

Let us demonstrate the Lagrange method by the example of a problem taken from a publication [4] by which we can compare the results:

There are two ways of producing a certain product. The production costs for each method depend on the produced x and y as follows: g (x) = 9 * x + x ** 2, g (y) = 6 * y + y ** 2. For a month it is necessary to produce 3 × 50 units of products, distributing it between the two methods so as to minimize total costs.

Consider a symbolic solution to the problem, which explains the method.

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]))

Result

The analyzed function f for the variables x, y:

f = x ** 2 + 9 * x + y ** 2 + 6 * y

Restrictions: x + y - 150 = 0

Auxiliary Lagrange function:

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

Stationary point M (x, y):
74.25, 75.75

Instead of conclusions

The result coincides with that given in [4], but the main thing, in my opinion, is different. Using symbolic computing as an example, it’s easier to understand any method and use it. But for mathematical packages, this idea is not new in both Mathcad and Matlab, symbolic calculations have long been widely used. But this is in Python!

Thanks to everyone who read the article.

References

1. Introduction to scientific Python.
2. Theme 6. Nonlinear programming.
3. Extremums of functions of two and three variables.
4. The method of Lagrange multipliers. An example of a solution.