数学优化（Mathematical Optimization），也称为数值优化或数理优化，是关于选择最优元素（通常是最大化或最小化某一特定函数的值）的一门数学学科。这些元素来自某个定义域，该定义域受到一系列约束条件的限制。数学优化在工程、经济学、运营管理、物流和医学等多个领域有广泛应用。
目标函数（Objective Function）：
决策变量（Decision Variables）：
约束条件（Constraints）：
线性规划（Linear Programming, LP）：
非线性规划（Non-linear Programming, NLP）：
混合整数规划（Mixed-Integer Programming, MIP）：
动态规划（Dynamic Programming, DP）：
设想一个简单的问题，一个公司生产两种产品 (A) 和 (B)，利润分别是 (20) 和 (30) 元，每种产品需要的生产时间分别是 (1) 小时和 (2) 小时，而生产时间的总约束是 (8) 小时。公司希望最大化利润。
from scipy.optimize import linprog
# 定义目标函数系数（利润）
c = [-20, -30] # 因为 linprog 函数默认是最小化问题，所以这里用负数表示最大化
# 定义不等式约束，A_ub * x <= b_ub
A = [[1, 2]] # 每种产品的时间约束
b = [8] # 总时间限制
# 定义变量的非负约束
x_bounds = (0, None)
y_bounds = (0, None)
# 求解线性规划
res = linprog(c, A_ub=A, b_ub=b, bounds=[x_bounds, y_bounds], method='highs')
print('Optimal value:', -res.fun) # 恢复成正数表示最大化
print('X:', res.x)Optimal value: 160.0
X: [8. 0.]在这个例子中，数学优化帮助我们确定了在生产时间有限的情况下，应该生产多少种产品来最大化总利润。这就是数学优化在实际应用中的一个经典示例。
import numpy as np
import matplotlib.pyplot as plt
# 提取结果
optimal_value = -res.fun # 最优目标函数值（转正数）
optimal_x = res.x # 最优变量值
# 设置图表数据
x = np.linspace(0, 10, 10)
y = (8 - x) / 2 # 从约束方程 y = (8 - x) / 2 得到
plt.figure(figsize=(10, 6))
# 绘制约束条件的可行域
plt.plot(x, y, label=r'$x + 2y \leq 8$', color='r', linestyle='-')
plt.fill_between(x, 0, y, where=(y >= 0), color='gray', alpha=0.3)
# 绘制目标函数等高线（等价值线）
contour_lines = np.arange(100, 240, 20)
for value in contour_lines:
plt.plot(x, (value / 30) - (20 / 30) * x, linestyle='--', label=f'$20x + 30y = {value}$')
# 绘制最优点
plt.plot(optimal_x[0], optimal_x[1], 'bo') # 最优点
plt.text(optimal_x[0], optimal_x[1], f' Optimal Point\n({optimal_x[0]:.2f}, {optimal_x[1]:.2f})', fontsize=12, color='blue')
# 添加图表元素
plt.xlabel('Product A')
plt.ylabel('Product B')
plt.title('Linear Programming Optimization')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.legend(loc='upper right')
plt.grid(True)
plt.show()
之前的例子没有考虑到 Product A 和 Product B 的取值需要是非负整数。如果我们要在这个条件下进行优化，就需要进行整数规划。在这种情况下，我们需要使用混合整数线性规划（Mixed-Integer Linear Programming, MILP）来求解。
使用 Gurobipy 库可以解决这个问题。
安装 Gurobipy
conda install -c gurobi gurobi
import numpy as np
import matplotlib.pyplot as plt
from gurobipy import Model, GRB
# 创建模型
model = Model()
# 添加变量
x = model.addVar(vtype=GRB.INTEGER, name="x") # Product A
y = model.addVar(vtype=GRB.INTEGER, name="y") # Product B
# 设置目标函数：最大化利润
model.setObjective(20 * x + 30 * y, GRB.MAXIMIZE)
# 添加约束：时间限制
model.addConstr(x + 2 * y <= 8, 'time_limit')
# 优化模型
model.optimize()
# 提取结果
optimal_value = model.objVal
optimal_x = [v.x for v in model.getVars()]
# 打印结果
for v in model.getVars():
print(f'{v.VarName}: {v.x}')
print(f'Optimal value: {optimal_value}')
# 画图
x_vals = np.arange(0, 9, 1) # Product A 的整数取值范围
y_vals = np.arange(0, 5, 1) # Product B 的整数取值范围
X, Y = np.meshgrid(x_vals, y_vals)
Z = 20 * X + 30 * Y # 目标函数
plt.figure(figsize=(10, 6))
contour = plt.contourf(X, Y, Z, levels=10, cmap='Blues')
plt.colorbar(contour)
plt.plot(optimal_x[0], optimal_x[1], 'bo') # 最优点
plt.text(optimal_x[0], optimal_x[1], f' Optimal Point\n({int(optimal_x[0])}, {int(optimal_x[1])})', fontsize=12, color='blue')
# 绘制约束线和可行域
plt.plot(x_vals, (8 - x_vals) / 2, label=r'$x + 2y \leq 8$', color='red', linestyle='-')
plt.fill_between(x_vals, 0, (8 - x_vals) / 2, where=(8 - x_vals) / 2 >= 0, color='gray', alpha=0.3)
# 设置Y轴整数刻度
plt.yticks(np.arange(0, max(y_vals)+1, 1)) # 设置Y轴仅显示整数刻度
# 添加图表元素
plt.xlabel('Product A')
plt.ylabel('Product B')
plt.title('Integer Linear Programming Optimization')
plt.legend(loc='upper right')
plt.grid(True)
plt.show()Set parameter Username
Academic license - for non-commercial use only - expires 2025-04-02
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (mac64[arm] - Darwin 24.1.0 24B83)
CPU model: Apple M3 Pro
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 1 rows, 2 columns and 2 nonzeros
Model fingerprint: 0xd214e7d0
Variable types: 0 continuous, 2 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+00]
Objective range [2e+01, 3e+01]
Bounds range [0e+00, 0e+00]
RHS range [8e+00, 8e+00]
Found heuristic solution: objective 160.0000000
Presolve removed 1 rows and 2 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 12 available processors)
Solution count 1: 160
Optimal solution found (tolerance 1.00e-04)
Best objective 1.600000000000e+02, best bound 1.600000000000e+02, gap 0.0000%
x: 8.0
y: -0.0
Optimal value: 160.0
约束条件的可行域
目标函数的等高线
最优点
Optimal Value（最优值）：
optimal_value 为 160，即表示在给定约束条件下所能达到的最大目标函数值（利润）。Optimal X（最优解）：
optimal_x 为 [8., 0.]，即为达到最优目标所需要的产品 A 和产品 B 的生产数量。通过上述图示和解释，我们可以清晰、准确地理解在整数规划条件下优化问题的求解过程及最优值和最优解的实际含义。
上面介绍了数学优化的好处，现在让我们看看它的实际应用，以及优化与机器学习如何协同工作
线性回归模型中的参数求解过程——其实就在最小化损失函数。
# 线性回归示例
import numpy as np
from sklearn.linear_model import LinearRegression
# 生成数据
X = np.array([[1, 1], [1, 2], [2, 2], [2, 3]])
y = np.dot(X, np.array([1, 2])) + 3
# 创建线性回归模型并训练
model = LinearRegression().fit(X, y)
# 优化结果
print(model.coef_, model.intercept_)[1. 2.] 3.0000000000000018混合整数规划（Mixed-Integer Programming, MIP）是一种数学优化技术，包含二类变量：连续变量（可以取任意实数值）和离散变量（只能取整数值）。MIP在许多实际应用中都有广泛的应用，尤其在那些需要对决策变量进行整数限定的问题领域。以下是一些常见的应用领域和具体实例：
应用领域
生产与供应链管理：
金融与投资：
能源管理：
设施选址：
健康管理：
具体实例：生产计划问题
假设一个工厂生产两种产品：产品 A 和产品 B。工厂有两个机器，分别用来生产这两种产品。每个产品的生产时间和利润如下表所示：
| 产品 | 机器1时间（小时） | 机器2时间（小时） | 利润（美元） |
|---|---|---|---|
| A | 2 | 1 | 30 |
| B | 1 | 2 | 20 |
机器 1 每周的总可用时间为 40 小时，机器 2 每周的总可用时间为 60 小时。目标是确定产品A和B的生产数量，使得总利润最大化。在这里，生产数量是整数变量。
混合整数规划模型如下：
from gurobipy import Model, GRB
# 创建模型
model = Model()
# 添加变量
x = model.addVar(vtype=GRB.INTEGER, name="x") # 产品A
y = model.addVar(vtype=GRB.INTEGER, name="y") # 产品B
# 设置目标函数：最大化利润
model.setObjective(30 * x + 20 * y, GRB.MAXIMIZE)
# 添加约束
model.addConstr(2 * x + y <= 40) # 机器1时间约束
model.addConstr(x + 2 * y <= 60) # 机器2时间约束
# 优化模型
model.optimize()
# 提取结果
for v in model.getVars():
print(f'{v.VarName}: {v.x}')
print(f'Optimal value: {model.objVal}')Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (mac64[arm] - Darwin 24.1.0 24B83)
CPU model: Apple M3 Pro
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Model fingerprint: 0x8826119c
Variable types: 0 continuous, 2 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+00]
Objective range [2e+01, 3e+01]
Bounds range [0e+00, 0e+00]
RHS range [4e+01, 6e+01]
Found heuristic solution: objective 600.0000000
Presolve time: 0.00s
Presolved: 2 rows, 2 columns, 4 nonzeros
Variable types: 0 continuous, 2 integer (0 binary)
Root relaxation: objective 7.333333e+02, 2 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 733.33333 0 2 600.00000 733.33333 22.2% - 0s
H 0 0 730.0000000 733.33333 0.46% - 0s
0 0 733.33333 0 2 730.00000 733.33333 0.46% - 0s
Explored 1 nodes (2 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 12 (of 12 available processors)
Solution count 2: 730 600
Optimal solution found (tolerance 1.00e-04)
Best objective 7.300000000000e+02, best bound 7.300000000000e+02, gap 0.0000%
x: 7.0
y: 26.0
Optimal value: 730.0在这个实例中，MIP 帮助我们找到最大化利润的最佳生产数量。
假设一个城市的交通网络如下图所示，每条边（道路）上标出的是该路段的道路容量（即最多允许通过的车辆数）。我们的目标是从源节点 S 到目标节点 T 的流量最大化。
graph TD
S -->|4| A
S -->|5| B
A -->|3| C
A -->|2| D
B -->|3| D
B -->|4| E
C --->|3| T
D --->|2| T
E -->|2| F
F -->|3| T
路径上的数值表示该路径的最大容量，例如 S → A = 4 指 S 点到 A 点最大流量为 4.
下面是使用 Gurobi 进行交通流量优化建模和求解的代码。
from gurobipy import Model, GRB
# 创建模型
model = Model("MaxFlow")
# 边及其容量
edges = {'S': {'A': 4, 'B': 5},
'A': {'C': 3, 'D': 2},
'B': {'D': 3, 'E': 4},
'C': {'T': 3},
'D': {'T': 2},
'E': {'F': 2},
'F': {'T': 3}}
# 添加变量
flows = {}
for u in edges:
for v in edges[u]:
flows[u, v] = model.addVar(name=f"flow_{u}_{v}", ub=edges[u][v])
# 设置目标函数：最大化从 S 到 T 的流量
model.setObjective(sum(flows['S', v] for v in edges['S']), GRB.MAXIMIZE)
# 添加节点平衡约束
# 所有中间节点的流量平衡：输入等于输出
nodes = set(edges.keys()).union({v for u in edges for v in edges[u]})
for node in nodes:
if node not in ['S', 'T']:
inflow = sum(flows[u, node] for u in edges if (u, node) in flows)
outflow = sum(flows[node, v] for v in edges[node]) if node in edges else 0
model.addConstr(inflow == outflow, name=f"node_{node}_balance")
# 优化模型
model.optimize()
# 打印结果
if model.status == GRB.OPTIMAL:
print('\nOptimal flow value: ', model.objVal)
for (u, v), flow_var in flows.items():
print(f"{u} -> {v}: {flow_var.x}")Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (mac64[arm] - Darwin 24.1.0 24B83)
CPU model: Apple M3 Pro
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 6 rows, 10 columns and 15 nonzeros
Model fingerprint: 0xcaf27f06
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [1e+00, 1e+00]
Bounds range [2e+00, 5e+00]
RHS range [0e+00, 0e+00]
Presolve removed 6 rows and 10 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration Objective Primal Inf. Dual Inf. Time
0 7.0000000e+00 0.000000e+00 0.000000e+00 0s
Solved in 0 iterations and 0.00 seconds (0.00 work units)
Optimal objective 7.000000000e+00
Optimal flow value: 7.0
S -> A: 3.0
S -> B: 4.0
A -> C: 3.0
A -> D: 0.0
B -> D: 2.0
B -> E: 2.0
C -> T: 3.0
D -> T: 2.0
E -> F: 2.0
F -> T: 2.0ub）设为边的容量。例如 flow_S_A 表示从 S 到 A 的流量，它的上限为 4。S 的总出流量，即从 S 到所有邻接节点（即 S → A 和 S → B）的流量之和。A 的约束是 flow_S_A == flow_A_C + flow_A_D。S 和 T 将没有这一平衡约束，因为 S 是源点，而 T 是目的点。model.optimize() 来求解最大流量问题。各条路径上的流量如下图所示，最大的流量即为 7。
graph TD
S -->|3| A
S -->|4| B
A -->|3| C
A .->|0| D
B -->|2| D
B -->|2| E
C --->|3| T
D --->|2| T
E -->|2| F
F -->|2| T
这个实例展示了如何使用混合整数规划方法处理最大流问题，确定从源节点到目标节点的最大流量，同时满足所有边的容量约束。
数学优化在改进机器学习模型的各个方面起着至关重要的作用，包括超参数调优、特征选择和构建更稳健的预测模型。以下是这些方面的详细解释以及如何使用数学优化来实现它们的改进。
X_train 和 y_train 被假定是已经分割好的训练数据集。为了完整地说明这些概念，我们将使用一个流行的内置数据集，并展示如何加载、预处理和分割数据集。在这次演示中，我们将使用 scikit-learn 中的 iris 数据集。
我们将首先展示如何从 scikit-learn 加载数据集，并进行数据预处理和分割。例如，对于 iris 数据集，以下是代码示例：
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import pandas as pd
# 加载 iris 数据集
iris = load_iris()
X = iris.data # 特征矩阵
y = iris.target # 目标变量
# 数据标准化
scaler = StandardScaler()
X = scaler.fit_transform(X)
# 分割数据集为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
# 将数据转换为 pandas DataFrame 以便后续使用
X_train = pd.DataFrame(X_train, columns=iris.feature_names)
X_test = pd.DataFrame(X_test, columns=iris.feature_names)
y_train = pd.Series(y_train)
y_test = pd.Series(y_test)
print("训练集特征维度:", X_train.shape)
print("训练集标签维度:", y_train.shape)
print("测试集特征维度:", X_test.shape)
print("测试集标签维度:", y_test.shape)训练集特征维度: (105, 4)
训练集标签维度: (105,)
测试集特征维度: (45, 4)
测试集标签维度: (45,)超参数调优是优化机器学习模型性能的关键部分，因为超参数对模型的训练过程和最终性能有重大影响。
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestClassifier
# 定义参数网格
param_grid = {
'n_estimators': [100, 200, 300],
'max_depth': [10, 20, 30]
}
# 创建模型
model = RandomForestClassifier()
# 创建网格搜索实例
grid_search = GridSearchCV(estimator=model, param_grid=param_grid, cv=5)
grid_search.fit(X_train, y_train)
# 最优参数
print(f'Best Parameters: {grid_search.best_params_}')Best Parameters: {'max_depth': 10, 'n_estimators': 100}from sklearn.model_selection import RandomizedSearchCV
from sklearn.ensemble import RandomForestClassifier
from scipy.stats import randint
# 定义参数分布
param_dist = {
'n_estimators': randint(100, 500),
'max_depth': randint(10, 50)
}
# 创建模型
model = RandomForestClassifier()
# 创建随机搜索实例
random_search = RandomizedSearchCV(estimator=model, param_distributions=param_dist, n_iter=50, cv=5)
random_search.fit(X_train, y_train)
# 最优参数
print(f'Best Parameters: {random_search.best_params_}')Best Parameters: {'max_depth': 41, 'n_estimators': 208}特征选择有助于减少模型的复杂性，提高模型的泛化能力和性能。
from sklearn.linear_model import Lasso
# 创建Lasso模型
model = Lasso(alpha=0.1)
model.fit(X_train, y_train)
# 选择非零特征
selected_features = X_train.columns[model.coef_ != 0]
print(f'Selected Features: {selected_features}')Selected Features: Index(['petal length (cm)', 'petal width (cm)'], dtype='object')from sklearn.feature_selection import RFE
from sklearn.ensemble import RandomForestClassifier
# 创建模型
model = RandomForestClassifier()
# 创建RFE实例
selector = RFE(estimator=model, n_features_to_select=5, step=1)
selector = selector.fit(X_train, y_train)
# 选择特征
selected_features = X_train.columns[selector.support_]
print(f'Selected Features: {selected_features}')Selected Features: Index(['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)',
'petal width (cm)'],
dtype='object')数学优化有助于构建更稳健的预测模型，增强其泛化能力，处理数据噪声和异常值。
# 使用随机森林模型
model = RandomForestClassifier()
model.fit(X_train, y_train)
# 使用梯度生成对抗样本（这里模拟一个简单的对抗样本生成过程，但实际效果存在局限）
def adversarial_training(model, X_train, y_train, epsilon=0.1):
# 模拟对抗样本生成（注意：深度模型通常用梯度生成对抗样本，这里简化处理）
perturbation = epsilon * np.sign(np.random.randn(*X_train.shape))
X_adv = X_train + perturbation
X_combined = np.vstack([X_train, X_adv])
y_combined = np.hstack([y_train, y_train])
# 合并对抗样本和原样本进行训练
model.fit(X_combined, y_combined)
return model
# 应用对抗训练
model = adversarial_training(model, X_train.values, y_train.values)
# 验证模型
acc = model.score(X_test, y_test)
print("模型在测试集上的准确率: ", acc)模型在测试集上的准确率: 1.0from sklearn.model_selection import cross_val_score
# 创建模型
model = RandomForestClassifier()
# 交叉验证评估模型
scores = cross_val_score(model, X, y, cv=5)
print(f'Cross-Validation Scores: {scores}')
print(f'Mean Score: {scores.mean()}')Cross-Validation Scores: [0.96666667 0.96666667 0.93333333 0.96666667 1. ]
Mean Score: 0.9666666666666668数学优化在机器学习模型的改进中扮演着重要角色，通过优化超参数、选择最佳特征和构建更稳健的模型，可以显著提升模型的性能和泛化能力。这些方法的具体实现通过示例代码得到了展示，帮助理解和应用数学优化技术来改进机器学习模型。
我们可以使用不同的优化方法解决复杂的问题。比如在生产计划中，我们希望尽可能减少生产成本。我们可以使用线性规划来解决这个问题。
from scipy.optimize import linprog
# 定义目标函数系数
c = [29, 45]
# 定义不等式约束，A_ub * x <= b_ub
A = [[-5, 20], [10, -50]]
b = [-15, -30]
# 求解线性规划
res = linprog(c, A_ub=A, b_ub=b)
print('Optimal value:', res.fun, '\nX:', res.x)Optimal value: 1053.0
X: [27. 6.]哦，所以你认为你可以优化：不同的方法如何帮助你做出更好的决策
非线性规划（Nonlinear Programming, NLP）是数学优化中的一类重要问题，其中目标函数或约束条件至少有一个是非线性的。解决非线性规划问题的经典方法有很多，以下是一些常用的方法及其简要介绍：
梯度下降法是一种迭代优化算法，主要用于寻找连续可微函数的局部最小值。
基本步骤
优缺点
import numpy as np
def gradient_descent(f, grad_f, x0, lr=0.01, max_iter=1000, tol=1e-6):
x = x0
for i in range(max_iter):
grad = grad_f(x)
x_new = x - lr * grad
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
return x
# 示例
f = lambda x: x[0]**2 + x[1]**2
grad_f = lambda x: np.array([2*x[0], 2*x[1]])
x0 = np.array([2, 3])
opt_x = gradient_descent(f, grad_f, x0)
print(f'Optimal x: {opt_x}')Optimal x: [2.75597916e-05 4.13396874e-05]牛顿法是一种基于二阶导数（Hessian矩阵）的优化方法，用于寻找可微可导函数的极小值。
基本步骤
优缺点
import numpy as np
def newton_method(f, grad_f, hessian_f, x0, max_iter=100, tol=1e-6):
x = x0
for i in range(max_iter):
grad = grad_f(x)
hessian = hessian_f(x)
x_new = x - np.linalg.inv(hessian).dot(grad)
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
return x
# 示例
f = lambda x: x[0]**2 + x[1]**2
grad_f = lambda x: np.array([2*x[0], 2*x[1]])
hessian_f = lambda x: np.array([[2, 0], [0, 2]])
x0 = np.array([2, 3])
opt_x = newton_method(f, grad_f, hessian_f, x0)
print(f'Optimal x: {opt_x}')Optimal x: [0. 0.]内点法是一种用于处理具有约束条件的非线性规划问题的优化方法。内点法通过引入障碍函数来避免逐出可行域。
基本步骤
优缺点
信赖域法是一类迭代优化算法，通过在每一步迭代中求解一个包含约束的小规模子问题来更新变量。
基本步骤
优缺点
import numpy as np
from scipy.optimize import minimize
# 目标函数
def rosenbrock(x):
return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
x0 = np.array([1.2, 1.2])
result = minimize(rosenbrock, x0, method='trust-constr')
print(f'Optimal x: {result.x}, Optimal value: {result.fun}')Optimal x: [0.99999552 0.99999104], Optimal value: 2.005102033632681e-11以上是一些经典的非线性规划优化方法，每种方法都有其优缺点和适用范围。选择合适的优化方法需要根据具体问题的性质、规模和目标函数的特点来决定。通过这些优化技术，可以高效解决实际中的复杂非线性规划问题。
高级优化算法如遗传算法、模拟退火和粒子群优化等在处理复杂的非线性、多峰、多约束优化问题时非常有用。它们通常基于生物进化、物理退火或群体智能等自然现象的启发。
遗传算法是一种基于自然选择和遗传机制的优化算法。它通过模拟生物进化过程中的选择、交叉和突变操作来发现最优解。
基本步骤
示例代码
import numpy as np
from deap import base, creator, tools, algorithms
def evaluate(individual):
return sum(individual),
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)
toolbox = base.Toolbox()
toolbox.register("attr_float", np.random.uniform, -10, 10)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=10)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", evaluate)
population = toolbox.population(n=100)
NGEN = 50
CXPB, MUTPB = 0.5, 0.2
for gen in range(NGEN):
offspring = toolbox.select(population, len(population))
offspring = list(map(toolbox.clone, offspring))
for child1, child2 in zip(offspring[::2], offspring[1::2]):
if np.random.rand() < CXPB:
toolbox.mate(child1, child2)
del child1.fitness.values
del child2.fitness.values
for mutant in offspring:
if np.random.rand() < MUTPB:
toolbox.mutate(mutant)
del mutant.fitness.values
invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
fitnesses = map(toolbox.evaluate, invalid_ind)
for ind, fit in zip(invalid_ind, fitnesses):
ind.fitness.values = fit
population[:] = offspring
best_ind = tools.selBest(population, 1)[0]
print(f'Best individual is {best_ind}, {best_ind.fitness.values}')Best individual is [17.74946320369606, 12.628846568393326, 12.814644772889038, 11.284328141367586, 13.757849401094125, 11.719117378960084, 14.833331027352905, 13.90895408387577, 13.031288687132951, 9.973760189588456], (131.7015834543503,)模拟退火是一种基于物理退火过程的优化算法，通过逐步降低”温度”来寻找到函数的全局最优解。
基本步骤
import numpy as np
def objective_function(x):
return x**2
def simulated_annealing(objective, bounds, n_iterations, step_size, temp):
best = bounds[:, 0] + (np.random.rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0]))
best_eval = objective(best)
curr, curr_eval = best, best_eval
for i in range(n_iterations):
candidate = curr + np.random.randn(len(bounds)) * step_size
candidate_eval = objective(candidate)
if candidate_eval < best_eval or np.exp((curr_eval - candidate_eval) / temp) > np.random.rand():
curr, curr_eval = candidate, candidate_eval
if candidate_eval < best_eval:
best, best_eval = candidate, candidate_eval
temp *= 0.99
return best, best_eval
bounds = np.array([[-5.0, 5.0]])
n_iterations = 1000
step_size = 0.1
temp = 10
best, score = simulated_annealing(objective_function, bounds, n_iterations, step_size, temp)
print(f'Best solution: {best}')
print(f'Objective function value: {score}')Best solution: [-0.00021292]
Objective function value: [4.53336993e-08]粒子群优化是一种基于群体智能的优化算法，通过模拟鸟群觅食行为来寻找最优解。
基本步骤
import numpy as np
def objective_function(x):
return np.sum(x**2)
def particle_swarm_optimization(objective, bounds, n_particles, n_iterations, w=0.5, c1=1.0, c2=2.0):
dim = len(bounds)
swarm = np.random.rand(n_particles, dim)
swarm = bounds[:, 0] + swarm * (bounds[:, 1] - bounds[:, 0])
velocity = np.zeros_like(swarm)
personal_best = np.copy(swarm)
personal_best_scores = np.array([objective(p) for p in personal_best])
global_best = swarm[np.argmin(personal_best_scores)]
global_best_score = np.min(personal_best_scores)
for i in range(n_iterations):
for j in range(n_particles):
r1, r2 = np.random.rand(dim), np.random.rand(dim)
velocity[j] = w * velocity[j] + c1 * r1 * (personal_best[j] - swarm[j]) + c2 * r2 * (global_best - swarm[j])
swarm[j] += velocity[j]
score = objective(swarm[j])
if score < personal_best_scores[j]:
personal_best[j], personal_best_scores[j] = swarm[j], score
global_best = personal_best[np.argmin(personal_best_scores)]
global_best_score = np.min(personal_best_scores)
return global_best, global_best_score
bounds = np.array([[-5.0, 5.0], [-5.0, 5.0]])
n_particles = 30
n_iterations = 100
best, score = particle_swarm_optimization(objective_function, bounds, n_particles, n_iterations)
print(f'Best solution: {best}')
print(f'Objective function value: {score}')Best solution: [-6.14061071e-12 5.39887490e-12]
Objective function value: 6.68549500754692e-23上述三种高级优化算法（遗传算法、模拟退火和粒子群优化）在处理复杂的非线性、多峰和多约束优化问题时非常有效。它们通过模拟生物进化、物理退火和群体智能等自然现象的过程，能够在全局范围内搜索最优解，是解决实际优化问题的重要工具。
Gurobipy 是一个功能强大的库，能够处理更复杂的线性和非线性优化问题。我们将实现一个简单的混合整数线性规划问题。
from gurobipy import Model, GRB
# 创建模型
model = Model()
# 添加变量
x = model.addVar(vtype=GRB.CONTINUOUS, name="x")
y = model.addVar(vtype=GRB.INTEGER, name="y")
# 设置目标函数
model.setObjective(2 * x + 3 * y, GRB.MINIMIZE)
# 添加约束
model.addConstr(x + y >= 10, "c0")
model.addConstr(3 * x + 2 * y <= 25, "c1")
# 优化模型
model.optimize()
# 输出结果
for v in model.getVars():
print(f'{v.varName}: {v.x}')
print(f'Objective: {model.objVal}')Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (mac64[arm] - Darwin 24.1.0 24B83)
CPU model: Apple M3 Pro
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Model fingerprint: 0x6e3d58b8
Variable types: 1 continuous, 1 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 3e+00]
Objective range [2e+00, 3e+00]
Bounds range [0e+00, 0e+00]
RHS range [1e+01, 2e+01]
Presolve removed 2 rows and 2 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 12 available processors)
Solution count 1: 25
Optimal solution found (tolerance 1.00e-04)
Best objective 2.500000000000e+01, best bound 2.500000000000e+01, gap 0.0000%
x: 5.0
y: 5.0
Objective: 25.0