This article is participating in Python Theme Month. See the link for details

Do mathematical modeling in Python

In the past (a long time ago) I tried to learn mathematical modeling for a few days, AND I didn’t like the syntax of Matlab and Lingo, so I tried to do the examples from the book in Python, which resulted in this article.

In this article, I’ve written Python solutions to a few simple problems. I just want to show that mathematical modeling in Python is possible and elegant.

(As this is an old article written in the past, it was not carefully corrected before publishing. Please forgive me if there are any mistakes.)

Linear programming

Third party dependency library: Scipy.

Scipy.optimize. Linprog can solve linear programming problems:

linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='simplex', callback=None, options=None)
Copy the code

The standard form of the problem is:

        Minimize:     c^T * x
    
        Subject to:   A_ub * x <= b_ub
                      A_eq * x == b_eq
Copy the code

The following linear programming problems are solved:

\textrm{max } z=2x_1+3x_2-5x_3,\\ \textrm{s.t. } \left \{ \begin{array}{ll} x_1+x_2+x_3=7\\ 2x_1-5x_2+x_3\ge10\\ x_1+3x_2+x_3\le12\\ x_1,x_2,x_3\ge0 \end{array}\right.

Solution:

#! /usr/bin/python3

import numpy as np
from scipy import optimize


c = [-2, -3.5]
a = [[-2.5, -1], [1.3.1]]
b = [-10.12]
aeq = [[1.1.1]]
beq = [7]
bounds = [[0.None], [0.None], [0.None]]	# (0, None) means non-negative, this is a default value
result = optimize.linprog(c, a, b, aeq, beq, bounds)
print(result)
Copy the code

⚠️【 note 】 Here the problem is to find Max, the standard type is min, so when writing the c matrix, the values are written as negative in the problem; Similarly, the median value of a and b corresponding to the term >= should be negative. And then the final result is also going to be minus fun.

Output:

Fun: -14.571428571428571 message: 'Optimization terminated successfully.' nit: 2 Slack: Array ([3.85714286, 0.]) status: 0 SUCCESS: True x: array([6.42857143, 0.57142857, 0.])Copy the code

The optimal solution is x1=6.42857143,×2=0.57142857,×3=0x_1=6.42857143, X_2 =0.57142857, x_3=0x1=6.42857143,×2=0.57142857,×3=0, The corresponding optimal value z=14.571428571428571z=14.571428571428571 Z =14.571428571428571.

To get the values of fun and x, you can use result.fun and result.x directly.

Integer programming

The linear integer programming problem can be transformed into a linear programming problem.

For nonlinear integer programming, monte Carlo method can be used to obtain a satisfactory solution under a certain amount of computation.

Monte Carlo method (random sampling method)

Y =x2y=x^2y=x2, y=12−xy=12-xy=12−x and the XXX axis form a curved triangle in the first quadrant. A randomized experiment is designed to approximate the area of the image.

Solution: The idea of random experiment is as follows: generate 10^7 random points with uniform distribution in the rectangular region [0, 12] * [0, 9], calculate the frequency of random points falling on the curvilinear triangle, then the area of the curvilinear triangle is approximately the area of the above rectangle multiplied by the frequency.

import random


x = [random.random() * 12 for i in range(0.10000000)]
y = [random.random() * 9 for i in range(0.10000000)]

p = 0
for i in range(0.10000000) :if x[i] <= 3 and y[i] < x[i] ** 2:
            p += 1
    elif x[i] > 3 and y[i] < 12 - x[i]:
            p += 1
 
area_appr = 12 * 9 * p / 10 ** 7
print(area_appr)
Copy the code

It turned out to be around 49.5.

E.g. nonlinear integer programming is known as:

\textrm{max } z=x_1^2+x_2^2+3x_3^2+4x_4^2+2x_5^2-8x_1-2x_2-3x_3-x_4-2x_5\\ \textrm{s.t. } \left \{ \begin{array}{ll} 0 \le x_i \le 99, i=1,… ,5,\\ x_1+x_2+x_3+x_4+x_5 \le 400\\ x_1+2x_2+2x_3+x_4+6x_5 \le 800\\ 2x_1+x_2+6x_3 \le 200\\ x_3+x_4+5x_5 \le 200 \end{array}\right.

If you use enumeration, you have to calculate 100^5 = 10^10, which is too much work. Therefore, consider using Monte Carlo method to randomly calculate 10^6 points and get satisfactory points:

import time
import random

# Objective function
def f(x: list) - >int:
    return x[0] * *2 + x[1] * *2 + 3 * x[2] * *2 + 4 * x[3] * *2 + 2 * x[4] * *2 - 8 * x[0] - 2 * x[1] - 3 * x[2] - x[3] -2 * x[4]

# Constraint vector functions
def g(x: list) - >list:
    res = []
    res.append(sum(x) - 400)
    res.append(x[0] + 2 * x[1] + 2 * x[2] + x[3] + 6 * x[4] - 800)
    res.append(2 * x[0] + x[1] + 6 * x[2] - 200)
    res.append(x[2] + x[3] + 5 * x[4] - 200)
    return res

random.seed(time.time)

pb = 0
xb = []

for i in range(10 ** 6):
    x = [random.randint(0.99) for i in range(5)]		Generate a row of random integers on a five-column interval [0, 99]
    rf = f(x)
    rg = g(x)
    if all((a < 0 for a in rg)):		# if all elements of Rg are less than 0
        if pb < rf:
            xb = x
            pb = rf
            
print(xb, pb)
Copy the code

Assignment problem

Third party dependent libraries: Numpy, scipy.

Scipy.optimize. Linear_sum_assignment can solve the linear sum Assignment problem:

linear_sum_assignment(cost_matrix)
Copy the code

Note that I of element C[I, j] in the assignment matrix cost_matrix is worker and j is job.

Formally, let X be a boolean matrix where $X[i,j] = 1$ iff row i is assigned to column j. Then the optimal assignment has cost
$\min \sum_i \sum_j C_{i,j} X_{i,j}$
s.t. each row is assignment to at most one column, and each column to at most one row.

For the following assignment problem, the assignment matrix is known as

\left[ \begin{array}{cc} 3 & 8 & 2 & 10 & 3\\ 8 & 7 & 2 & 9 & 7\\ 6 & 4 & 2 & 7 & 5\\ 8 & 4 & 2 & 3 & 5\\ 9 & 10 & 6 & 9 & 10\\ \end{array} \right].

Solution:

>>> import numpy as np
>>> from scipy import optimize
>>> c = [[3.8.2.10.3], [8.7.2.9.7], [6.4.2.7.5], [8.4.2.3.5], [9.10.6.9.10]]
>>> c = np.array(c)
>>> optimize.linear_sum_assignment(c)
(array([0.1.2.3.4]), array([4.2.1.3.0]))	X15, x23, x32, x44, x51
>>> c[[0.1.2.3.4], [4.2.1.3.0]].sum(a)The result is substituted into cost_matrix to obtain the optimal value
21
Copy the code

Nonlinear programming

Nonlinear programming model

Third party dependent libraries: Numpy, scipy.

Scipy. Optimize integrates a series of programming functions that we’ve used many times before, including ways to solve nonlinear programming. For example, the minimize function can be used to solve many nonlinear programming problems with FMINcon solutions in Matlab.

minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)
Copy the code

Common parameters:

fun: : oThe minimum valueIs the objective function of,fun(x, *args) -> float, x isshape (n,)1 – D array
x0: Initial guess value,shape (n,)1 – D array
bounds: Sets a parameter range/constraint, tuple, as in((0, None), (0, None))
constraints:The constraint, put a series of dict tuples,({'type': TYPE, 'fun': FUN}, ...)
- TYPE:'eq'Means that the result of the function is equal to 0;'ineq'Means that the expression is greater than or equal to 0
- FUN: Constraint function

E.g. seek the following nonlinear programming:

\textrm{min } f(x)=x_1^2+x_2^2+x_3^2+8,\\ \textrm{s.t. } \left \{ \begin{array}{ll} x_1^2+x_2^2+x_3^2 \ge 0,\\ x_1+x_2^2+x_3^2 \le 20,\\ -x_1-x_2^2+2=0,\\ x_2+2x_3^2=3,\\ x_1,x_2,x_3 \ge 0. \end{array}\right.

import numpy as np
from scipy import optimize

f = lambda x: x[0] * *2 + x[1] * *2 + x[2] * *2 + 8

# Notice: eq ==; ineq >=
cons = ({'type': 'ineq'.'fun': lambda x: x[0] * *2 - x[1] + x[2] * *2},
        {'type': 'ineq'.'fun': lambda x: -x[0] - x[1] - x[2] * *3 + 20},
        {'type': 'eq'.'fun': lambda x: -x[0] - x[1] * *2 + 2},
        {'type': 'eq'.'fun': lambda x: x[1] + 2 * x[2] * *2 - 3})

res = optimize.minimize(f, (0.0.0), constraints=cons)

print(res)
Copy the code

Output:

Fun: 10.651091840572583 JAC: Array ([1.10433471, 2.40651834, 1.89564812]) message: 'Optimization terminated successfully.' nfev: 86 nit: 15 njev: 15 status: 0 success: True x: Array ([0.55216734, 1.20325918, 0.94782404])Copy the code

That is, Obtained when (x1, x2, x3) = (0.55216734, 1.20325918, 0.94782404) (x_1, x_2, x_3) = (0.55216734, 1.20325918, (x1, x2, x3) = 0.94782404) (0.55216734, 1.20325918, 0.94782404), The minimum y = 10.651091840572583 y = 10.651091840572583 y = 10.651091840572583.

Python solutions to unconstrained problems

Symbolic solution

Third party Dependence: Sympy.

E.g. seek the extremum of the multivariate function f(x, y)=x3−y3+3×2+3y2−9xf(x, y)=x ^ 3-y ^3 +3x ^2 +3y ^ 2-9xf (x, y)=x3−y3+3×2+3y2−9x

Thinking: to find the stagnation point, substitute into the Hessian matrix, positive definite value is minimum, negative definite maximum, may not be the extreme point.

Solution:

import sympy

f, x, y = sympy.symbols("f x y")

f = x**3 - y**3 + 3 * x**2 + 3 * y**2 - 9 * x
funs = sympy.Matrix([f])
args = sympy.Matrix([x, y])

df = funs.jacobian(args)        Student: First partial derivative
d2f = df.jacobian(args)         # Hessian matrix

stationaryPoints = sympy.solve(df)      # stagnation point

for i in stationaryPoints:
    a = d2f.subs(x, i[x]).subs(y, i[y])	# Hessian at stagnation
    b = a.eigenvals(multiple=True)      Find eigenvalues of Hessian matrix
    fv = f.subs(x, i[x]).subs(y, i[y])	# Function value at stagnation point
    
    if all((j > 0 for j in b)):
        print('point ({x}, {y}) is a minimum point corresponding to: f({x}, {y}) = {f}'.format(x=i[x], y=i[y], f=fv))
    elif all((j < 0 for j in b)):
        print('point ({x}, {y}) is the maximum point corresponding to: f({x}, {y}) = {f}'.format(x=i[x], y=i[y], f=fv))
    elif any((j < 0 for j in b)) and any((j > 0 for j in b)):
        print('Point ({x}, {y}) is not an extreme point'.format(x=i[x], y=i[y]))
    else:
        print('Unable to determine whether point ({x}, {y}) is an extreme point'.format(x=i[x], y=i[y]))
Copy the code

Output:

Point (-3, 0) is not an extreme point point (-3, 2) is a maximum point, the corresponding maximum is: f(-3, 2) = 31 point (1, 0) is a minimum point, the corresponding minimum is: f(1, 0) = -5 point (1, 2) is not an extreme pointCopy the code

Numerical solution

Third party dependent libraries: Numpy, scipy.

Using Python to find a numerical solution to an unconstrained extreme is similar to what we did in our nonlinear programming model. We still use the minimize function, but this time we do not even use constraints.

E.g. find the maximum value of the multivariate function f(x, y)=x3−y3+3×2+3y2−9xf(x, y)=x ^ 3-y ^3 +3x ^2 +3y ^ 2-9xf (x, y)=x3−y3+3×2+3y2−9x

import numpy as np from scipy import optimize f = lambda x: x[0]**3 - x[1]**3 + 3 * x[0]**2 + 3 * x[1]**2 - 9 * x[0] resMin = optimize.minimize(f, (0, 0)) # minimize resMax = optimize. Minimize (lambda x: -f(x), (0, 0)) print(" min: \n") print(resMin) print(" min: \n") print(resMax)Copy the code

Output:

Fun: -5.0HESS_INV: Array ([[8.34028325E-02, 3.27721596E-09], [3.27721596E-09, 1.00000000E]]) jAC: Array ([1.1920929e-07, 0.0000000e+00]) message: 'Optimization terminated successfully.' nfev: 20 nit: 4 njev: 5 status: 0 success: True x: array([1.00000000e+00, -5.40966234e-09]) fun: -30.99999999999847 hess_inv: Array ([[0.08280865, 0.00036445], [0.00036445, 0.16672048]]) JAC: ARRAY ([9.53674316E-07, -4.29153442E-06]) message: 'Optimization terminated successfully.' nfev: 172 nit: 11 njev: 43 status: 0 success: True x: Array ([2.99999994, 1.99999929])Copy the code

Find the zeros of functions and the solutions of equations

Third party Dependence: Sympy.

The sympay. Solve function can be used to obtain the symbolic solutions of the equations and each root can be converted to approximate values by using its EVALF method.

For polynomial f (x) = x3 – x2 + 2 x – 3 f (x) = x ^ 3 – x ^ 2 + 2 x – 3 f (x) = x3 – x2 + 2-3 x zero

>>> import sympy
>>> x = sympy.Symbol('x')
>>> s = sympy.solve(x**3 - x**2 + 2 * x -3)		Find the symbolic solution
>>> s
[1/3 + (-1/2 - sqrt(3)*I/2) * (65/54 + 5*sqrt(21) /18* * ()1/3) - 5/ (9* (-1/2 - sqrt(3)*I/2) * (65/54 + 5*sqrt(21) /18* * ()1/3)), 1/3 - 5/ (9* (-1/2 + sqrt(3)*I/2) * (65/54 + 5*sqrt(21) /18* * ()1/3)) + (-1/2 + sqrt(3)*I/2) * (65/54 + 5*sqrt(21) /18* * ()1/3), -5/ (9* (65/54 + 5*sqrt(21) /18* * ()1/3)) + 1/3 + (65/54 + 5*sqrt(21) /18* * ()1/3)]
>>> for i in s:
.    print(i.evalf())		# Approximate value
.
-0.137841101825493 - 1.52731225088663*I
-0.137841101825493 + 1.52731225088663*I
1.27568220365098

Copy the code

Find the solution of the following equations

\left\{\begin{array}{l} x^2+y-6=0\\ y^2+x-6=0\\ \end{array}\right.

>>> import sympy
>>> x, y = sympy.symbols('x y')
>>> s = sympy.solve((x**2+y-6, y**2+x-6))
>>> s
[{x: -3, y: -3}, {x: 2, y: 2}, {x: 6 - (1/2 - sqrt(21) /2) * *2, y: 1/2 - sqrt(21) /2}, {x: 6 - (1/2 + sqrt(21) /2) * *2, y: 1/2 + sqrt(21) /2}]
>>> for i in s:
.    print('({x}, {y})'.format(x=i[x].evalf(), y=i[y].evalf()))
.
(-3.00000000000000, -3.00000000000000)
(2.00000000000000.2.00000000000000)
(2.79128784747792, -1.79128784747792)
(-1.79128784747792.2.79128784747792)
Copy the code

Constrained extremum problem

See “Python Solutions to unconstrained problems” above.

To be continued

Add: numpy

Recommend a quality higher numpy Chinese document: www.numpy.org.cn/index.html.

Create a vector/matrix

With np. Array ()

>>> import numpy as np
>>> np.array([1.2.3])
array([1.2.3])
>>> np.array([[1.2.3], [4.5.6], [7.8.9]])
array([[1.2.3],
       [4.5.6],
       [7.8.9]])
Copy the code

Concatenation of matrices

Column merge/expansion:np.column_stack()
Line merge/expand:np.row_stack()

>>> import numpy as np
>>> A = np.array([[1.1], [2.2]])
>>> B = np.array([[3.3], [4.4]])
>>> A
array([[1.1],
       [2.2]])
>>> B
array([[3.3],
       [4.4]])
>>> np.column_stack((A, B))		Line #
array([[1.1.3.3],
       [2.2.4.4]])
>>> np.row_stack((A, B))			# column
array([[1.1],
       [2.2],
       [3.3],
       [4.4]])
Copy the code

Assignment of diagonal elements

Let’s say we have an NP.array X, and I want all of the diagonal values to be 0, and we can use nP.fill_diagonal (X, [0, 0, 0… ).

>>> D = np.zeros([3, 3])
>>> np.fill_diagonal(D, [1, 2, 3])
>>> D
array([[1., 0., 0.],
       [0., 2., 0.],
       [0., 0., 3.]])
Copy the code

To be continued

CDFMLR

mo4tech.com (Moment For Technology) is a global community with thousands techies from across the global hang out!Passionate technologists, be it gadget freaks, tech enthusiasts, coders, technopreneurs, or CIOs, you would find them all here.

Mathematical Modeling in Python | Python Theme Month

Do mathematical modeling in Python

Linear programming

Integer programming

Monte Carlo method (random sampling method)

Assignment problem

Nonlinear programming

Nonlinear programming model

Python solutions to unconstrained problems

Symbolic solution

Numerical solution

Find the zeros of functions and the solutions of equations

Constrained extremum problem

Add: numpy

Create a vector/matrix

Concatenation of matrices

Assignment of diagonal elements

Mathematical Modeling in Python | Python Theme Month

Do mathematical modeling in Python

Linear programming

Integer programming

Monte Carlo method (random sampling method)

Assignment problem

Nonlinear programming

Nonlinear programming model

Python solutions to unconstrained problems

Symbolic solution

Numerical solution

Find the zeros of functions and the solutions of equations

Constrained extremum problem

Add: numpy

Create a vector/matrix

Concatenation of matrices

Assignment of diagonal elements

Related Posts

Do you know how to use Lombok?

Maven builds the project to set the JDK version

In-depth understanding of thread communication