利用python处理非线性规划问题

（例题和部分解答思路来自清风老师，代码由自己实现，仅供参考）

思路综述

局部优化算法
最小化 (Minimize)
描述：局部优化算法通常用于在给定初始点附近寻找最优解。SciPy中提供了多种方法，例如BFGS、L-BFGS-B、TNC等。
优点：速度较快，适合在已知局部最优初始点附近进行精细调整。
缺点：局部优化算法的结果高度依赖于初始点的选择。如果初始点选择不佳，算法可能会收敛到一个局部最优解，而不是全局最优解。

全局优化算法
a差分进化 (Differential Evolution)
描述：差分进化是一种基于群体的全局优化算法，适用于连续空间的非线性函数。它通过维护一个种群并利用变异、交叉和选择等操作来搜索最优解。
优点：适合处理复杂的、具有多个局部极小值的问题，能在较大的搜索空间内找到全局最优解。
缺点：计算成本较高，结果可能不稳定，运行多次可能会得到不同的结果。

综合思路
使用全局优化算法（例如，差分进化）来搜索整个搜索空间，找到一个较好的解或一个接近全局最优解的初始点。这一阶段的目的是探索更广泛的解空间，寻找可能的全局最优解。

以全局优化阶段得到的解作为初始点，使用局部优化算法（例如，BFGS或L-BFGS-B）在该点附近进行细致的优化。这一步是为了进一步改进解，充分利用局部优化算法在解的附近精细求解的能力。

基本例题

例题1

#使用全局优化算法结合局部最优
import numpy as np
from scipy.optimize import differential_evolution, minimize, NonlinearConstraint, Bounds

# 目标函数，需要最大化，因此取负
def objective(x):
    x1, x2 = x
    return -(x1**2 + x2**2 - x1 * x2 - 2 * x1 - 5 * x2)


# 线性约束
def constraint2(x):
    x1, x2 = x
    return 2 * x1 - 3 * x2 + 6

# 定义变量的边界（明确的范围而不是无界，取足够大的范围1e6）
bounds = Bounds([-1e6, -1e6], [1e6, 1e6])

# 不等式约束：
nonlinear_constraint1 = NonlinearConstraint(constraint1, 0, np.inf)
# 定义线性约束
linear_constraint1 = LinearConstraint([2, -3], -6, 1e6)  

# 进行差分进化全局优化
result_global = differential_evolution(objective, bounds=bounds, constraints=[nonlinear_constraint1, linear_constraint_1])

# 使用差分进化的结果作为初始猜测进行局部优化
result_local = minimize(objective, result_global.x, constraints=[nonlinear_constraint1, nonlinear_constraint2])

# 输出结果
print('Global optimization result:')
print(f'Optimal value: {-result_global.fun:.4f}')  # 目标函数取负值，因此应反转回去
print(f'Optimal solution: x1 = {result_global.x[0]:.4f}')
print(f'Optimal solution: x2 = {result_global.x[1]:.4f}')

print('\nLocal optimization result:')
print(f'Optimal value: {-result_local.fun:.4f}')  # 目标函数取负值，因此应反转回去
print(f'Optimal solution: x1 = {result_local.x[0]:.4f}')
print(f'Optimal solution: x2 = {result_local.x[1]:.4f}')

Global optimization result:
Optimal value: -1.0000
Optimal solution: x1 = 1.0000
Optimal solution: x2 = 0.0000

Local optimization result:
Optimal value: -1.0000
Optimal solution: x1 = 1.0000
Optimal solution: x2 = 0.0000

例题2

import numpy as np
from scipy.optimize import differential_evolution, minimize, NonlinearConstraint, Bounds,LinearConstraint

# 目标函数，需要最大化
def objective(x):
    x1, x2,x3 = x
    return -(x1*x2*x3)

A1=np.array([-1,2,2])
b1_=0
b1=np.inf

A2=np.array([1,2,2])
b2_=-np.inf
b2=72

A3=np.array([1,-1,0])
b3_=10
b3=10

# 定义变量的边界（明确的范围而不是无界）
bounds = Bounds([20, 10,-1e6], [30, 20,1e6])

# 不等式约束：约束函数应该大于等于0
linear_constraint1 = LinearConstraint(A1,b1_,b1)
linear_constraint2 = LinearConstraint(A2,b2_,b2)

#等式约束：
linear_constraint3=LinearConstraint(A3,b3_,b3)


# 进行差分进化全局优化
result_global = differential_evolution(objective, bounds=bounds, constraints=[linear_constraint1,linear_constraint2,linear_constraint3])

# 使用差分进化的结果作为初始猜测进行局部优化
result_local = minimize(objective, result_global.x, constraints=[linear_constraint1,linear_constraint2,linear_constraint3])

# 输出结果
print('Global optimization result:')
print(f'Optimal value: {result_global.fun:.4f}')  
print(f'Optimal solution: x1 = {result_global.x[0]:.4f}')
print(f'Optimal solution: x2 = {result_global.x[1]:.4f}')
print(f'Optimal solution: x3 = {result_global.x[2]:.4f}')

print('\nLocal optimization result:')
print(f'Optimal value: {result_local.fun:.4f}')  
print(f'Optimal solution: x1 = {result_local.x[0]:.4f}')
print(f'Optimal solution: x2 = {result_local.x[1]:.4f}')
print(f'Optimal solution: x3 = {result_global.x[2]:.4f}')

Optimal value: -3445.6051
Optimal solution: x1 = 22.5850
Optimal solution: x2 = 12.5850
Optimal solution: x3 = 12.1225

Local optimization result:
Optimal value: -3445.6051
Optimal solution: x1 = 22.5850
Optimal solution: x2 = 12.5850
Optimal solution: x3 = 12.1225

典型例题

例一：工地选址问题(2)

建模

代码

import numpy as np
from math import sqrt
from scipy.optimize import differential_evolution, minimize, LinearConstraint, Bounds

# 工地坐标（p_j, q_j） j=1,2,...,6
loc_wok = [[1.25, 1.25], [8.75, 0.75], [0.5, 4.75], [5.75, 5], [3, 6.5], [7.25, 7.25]]

min_x=np.array(loc_wok).min(axis=0)[0]
min_y=np.array(loc_wok).min(axis=0)[1]
max_x=np.array(loc_wok).max(axis=0)[0]
max_y=np.array(loc_wok).max(axis=0)[0]

# 工地日需求量
need_wok = [3, 5, 4, 7, 6, 11]

# 最大供应量
max_sup = 20

# 目标函数，需要最小化
def objective(vars):
    result = 0
    a = [vars[12], vars[14]]
    b = [vars[13], vars[15]]
    for i in range(2):
        for j in range(6):
            result += (vars[i*6+j] * sqrt((a[i] - loc_wok[j][0])**2 + (b[i] - loc_wok[j][1])**2))
    return result

# 不等式约束
index1 = np.zeros((2, 16))
index1[0, :6] = 1
index1[1, 6:12] = 1

unequal_constraint = [LinearConstraint(index1[i], 0, max_sup) for i in range(2)]

# 等式约束
index2 = np.zeros((6, 16))
for i in range(6):
    index2[i, i] = index2[i, 6 + i] = 1

equal_constraint = [LinearConstraint(index2[i], need_wok[i], need_wok[i]) for i in range(6)]

# 定义变量的边界（减少搜索范围）
lowBound = [0] * 12 + [min_x,min_y]*2
highBound = [max(need_wok)] * 12 + [max_x,max_y]*2
bounds = Bounds(lowBound, highBound)

# 多阶段优化：先进行快速全局优化，然后进行局部优化

# 快速全局优化替代差分进化
result_global = differential_evolution(objective, bounds=bounds, constraints=equal_constraint + unequal_constraint, maxiter=500, popsize=10)

# 使用差分进化的结果作为初始猜测进行局部优化
result_local = minimize(objective, result_global.x, bounds=bounds, constraints=equal_constraint + unequal_constraint)

# 输出结果
print('Global optimization result:')
print(f'Optimal value: {result_global.fun:.2f}')
for i in range(16):
    print(f'Optimal solution: x{i+1} = {result_global.x[i]:.2f}')

print('\nLocal optimization result:')
print(f'Optimal value: {result_local.fun:.4f}')
for i in range(16):
    print(f'Optimal solution: x{i+1} = {result_local.x[i]:.2f}') 

结果

Global optimization result:
Optimal value: 89.57
Optimal solution: x1 = 2.98
Optimal solution: x2 = 0.03
Optimal solution: x3 = 3.99
Optimal solution: x4 = 0.10
Optimal solution: x5 = 5.99
Optimal solution: x6 = 2.93
Optimal solution: x7 = 0.02
Optimal solution: x8 = 4.97
Optimal solution: x9 = 0.01
Optimal solution: x10 = 6.90
Optimal solution: x11 = 0.01
Optimal solution: x12 = 8.07
Optimal solution: x13 = 3.00
Optimal solution: x14 = 6.50
Optimal solution: x15 = 6.14
Optimal solution: x16 = 5.17

Local optimization result:
Optimal value: 89.5715
Optimal solution: x1 = 2.98
Optimal solution: x2 = 0.03
Optimal solution: x3 = 3.99
Optimal solution: x4 = 0.10
Optimal solution: x5 = 5.99
Optimal solution: x6 = 2.93
Optimal solution: x7 = 0.02
Optimal solution: x8 = 4.97
Optimal solution: x9 = 0.01
Optimal solution: x10 = 6.90
Optimal solution: x11 = 0.01
Optimal solution: x12 = 8.07
Optimal solution: x13 = 3.00
Optimal solution: x14 = 6.50
Optimal solution: x15 = 6.14
Optimal solution: x16 = 5.17

例二：飞行管理问题

建模

代码

import numpy as np
from math import sqrt, cos, sin, radians
from scipy.optimize import differential_evolution, minimize, NonlinearConstraint, Bounds

# 工地坐标与初始角度
loc = np.array([[150, 140], [85, 85], [150, 155], [145, 50], [130, 150], [0, 0]])
angle = np.array([243, 236, 220.5, 159, 230, 52])
v = 800  # 速度

def objective(vars):
    return np.sum(vars**2)

# 非线性约束函数生成器
def create_func(i, j):
    def func(x):
        A1 = cos(radians(angle[j] + x[j])) - cos(radians(angle[i] + x[i]))
        B1 = sin(radians(angle[j] + x[j])) - sin(radians(angle[i] + x[i]))
        t0 = -((loc[j][1] - loc[i][1]) * B1 + (loc[j][0] - loc[i][0]) * A1) / (v * (A1**2 + B1**2))
        
        t = max(0, t0)
        
        x_j = loc[j][0] + v * t * cos(radians(angle[j] + x[j]))
        x_i = loc[i][0] + v * t * cos(radians(angle[i] + x[i]))
        y_j = loc[j][1] + v * t * sin(radians(angle[j] + x[j]))
        y_i = loc[i][1] + v * t * sin(radians(angle[i] + x[i]))
        
        return sqrt((x_j - x_i)**2 + (y_j - y_i)**2)
    
    return func

# 构建非线性约束
Constraints = []
for i in range(6):
    for j in range(6):
        if i != j:
            Constraints.append(NonlinearConstraint(create_func(i, j), 8, 1e6))

# 定义变量的边界
lowBound = [-30] * 6
highBound = [30] * 6
bounds = Bounds(lowBound, highBound)

# 全局优化
result_global = differential_evolution(objective, bounds=bounds, constraints=Constraints)

# 局部优化
result_local = minimize(objective, result_global.x, bounds=bounds, constraints=Constraints)

# 输出结果
print('Global optimization result:')
print(f'Optimal value: {result_global.fun:.2f}')
for i in range(6):
    print(f'Optimal solution: x{i+1} = {result_global.x[i]:.2f}')

print('\nLocal optimization result:')
print(f'Optimal value: {result_local.fun:.2f}')
for i in range(6):
    print(f'Optimal solution: x{i+1} = {result_local.x[i]:.2f}')

结果

Global optimization result:
Optimal value: 6.95
Optimal solution: x1 = 0.00
Optimal solution: x2 = -0.00
Optimal solution: x3 = 2.06
Optimal solution: x4 = -0.50
Optimal solution: x5 = -0.00
Optimal solution: x6 = 1.57

Local optimization result:
Optimal value: 6.95
Optimal solution: x1 = 0.00
Optimal solution: x2 = -0.00
Optimal solution: x3 = 2.06
Optimal solution: x4 = -0.50
Optimal solution: x5 = -0.00
Optimal solution: x6 = 1.57