山东省大学生数学建模竞赛代码手册

线性规划 (Linear Programming)

import numpy as np
from scipy.optimize import linprog

# 目标函数系数（最小化），例如：c = [c1, c2, ..., cn]
c = [-3, -2]  # 最大化 3x1 + 2x2 等价于最小化 -3x1 -2x2

# 不等式约束矩阵 A_ub * x <= b_ub
A_ub = [[1, 1], [2, 1], [1, 0]]   # 三个约束
b_ub = [10, 12, 4]

# 等式约束矩阵 A_eq * x == b_eq（若无等式约束，设为None）
A_eq = None
b_eq = None

# 变量的上下界
bounds = [(0, None), (0, None)]  # x1 >= 0, x2 >= 0

# 求解线性规划
res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')

print('最优解：', res.x)
print('最优值：', -res.fun)  # 因为目标函数取了负号
print('求解状态：', res.message)

参数	说明	示例
c	目标函数系数（最小化）	[-3, -2]
A_ub, b_ub	不等式约束矩阵和右侧向量	A_ub=[[1,1]], b_ub=[10]
A_eq, b_eq	等式约束矩阵和右侧向量	None
bounds	变量上下界	[(0,None), (0,None)]
method	求解算法，推荐'highs'	'highs'

整数规划 (Integer Programming)

# 安装：pip install pulp
import pulp

# 创建问题（最大化）
prob = pulp.LpProblem('Integer_Programming_Example', pulp.LpMaximize)

# 定义变量（lowBound下限，cat类型：Integer整数，Binary二进制）
x1 = pulp.LpVariable('x1', lowBound=0, cat='Integer')
x2 = pulp.LpVariable('x2', lowBound=0, cat='Integer')

# 目标函数
prob += 3*x1 + 2*x2, 'Objective'

# 约束条件
prob += x1 + x2 <= 10, 'Constraint1'
prob += 2*x1 + x2 <= 12, 'Constraint2'
prob += x1 <= 4, 'Constraint3'

# 求解
prob.solve(pulp.PULP_CBC_CMD(msg=False))

print('求解状态：', pulp.LpStatus[prob.status])
print('最优解：')
for var in prob.variables():
    print(f'{var.name} = {var.varValue}')
print('最优值：', pulp.value(prob.objective))

参数	说明	示例
cat	变量类型：'Integer'整数，'Binary'0-1，'Continuous'连续	'Integer'
lowBound, upBound	变量上下界	lowBound=0
求解器	默认使用CBC，也可用GLPK等	pulp.PULP_CBC_CMD

非线性规划 (Nonlinear Programming)

import numpy as np
from scipy.optimize import minimize

# 定义目标函数（Rosenbrock函数，经典测试函数）
def objective(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# 约束条件：x0 + x1 ≥ 1，x0² + x1² ≤ 2
def constraint1(x):
    return x[0] + x[1] - 1  # ≥0

def constraint2(x):
    return 2 - x[0]**2 - x[1]**2  # ≥0

cons = [
    {'type': 'ineq', 'fun': constraint1},
    {'type': 'ineq', 'fun': constraint2}
]

# 变量边界
bounds = [(-2, 2), (-2, 2)]

# 初始猜测
x0 = np.array([0.5, 0.5])

# 求解
res = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=cons)

print('最优解：', res.x)
print('最优值：', res.fun)
print('是否成功：', res.success)
print('迭代次数：', res.nit)

参数	说明	示例
method	求解算法：'SLSQP'（约束），'trust-constr'，'BFGS'（无约束）	'SLSQP'
constraints	约束列表，type可为'eq'或'ineq'	[{'type':'ineq', 'fun': c}]
bounds	变量边界	[(-2,2), (-2,2)]
options	算法选项，如maxiter, ftol	{'maxiter': 1000}

动态规划 (Dynamic Programming)

def knapsack_01(weights, values, capacity):
    """0-1背包问题动态规划求解"""
    n = len(weights)
    # dp[i][w] 表示前i个物品在容量w下的最大价值
    dp = [[0] * (capacity + 1) for _ in range(n + 1)]
    
    for i in range(1, n + 1):
        for w in range(capacity + 1):
            if weights[i-1] <= w:
                dp[i][w] = max(dp[i-1][w], dp[i-1][w - weights[i-1]] + values[i-1])
            else:
                dp[i][w] = dp[i-1][w]
    
    # 回溯找出选择的物品
    selected = []
    w = capacity
    for i in range(n, 0, -1):
        if dp[i][w] != dp[i-1][w]:
            selected.append(i-1)
            w -= weights[i-1]
    
    return dp[n][capacity], selected

# 示例数据
weights = [2, 3, 4, 5]   # 物品重量
values = [3, 4, 5, 6]    # 物品价值
capacity = 8             # 背包容量

max_value, selected_items = knapsack_01(weights, values, capacity)
print('最大价值：', max_value)
print('选择的物品索引：', selected_items)
print('对应重量：', [weights[i] for i in selected_items])
print('对应价值：', [values[i] for i in selected_items])

参数	说明	示例
weights	物品重量列表	[2,3,4,5]
values	物品价值列表	[3,4,5,6]
capacity	背包容量	8
dp数组	动态规划表，维度(n+1)×(capacity+1)

多目标规划 (Multi‑Objective Optimization)

import numpy as np
from scipy.optimize import minimize

# 两个目标函数
def f1(x):
    return x[0]**2 + x[1]**2

def f2(x):
    return (x[0]-1)**2 + (x[1]-1)**2

# 加权和法：将多目标转化为单目标
def weighted_sum(x, w1=0.5, w2=0.5):
    return w1 * f1(x) + w2 * f2(x)

# 约束：x0 + x1 ≥ 1
cons = ({'type': 'ineq', 'fun': lambda x: x[0] + x[1] - 1})

# 不同权重组合得到不同 Pareto 解
pareto_solutions = []
for w1 in np.linspace(0, 1, 6):
    w2 = 1 - w1
    res = minimize(weighted_sum, x0=[0, 0], args=(w1, w2), constraints=cons)
    if res.success:
        pareto_solutions.append((res.x, f1(res.x), f2(res.x)))

print('Pareto 解集（x, f1, f2）：')
for sol in pareto_solutions:
    print(f'x={sol[0]}, f1={sol[1]:.3f}, f2={sol[2]:.3f}')

# 绘制 Pareto 前沿（需 matplotlib）
import matplotlib.pyplot as plt
f1_vals = [sol[1] for sol in pareto_solutions]
f2_vals = [sol[2] for sol in pareto_solutions]
plt.scatter(f1_vals, f2_vals, c='red')
plt.xlabel('f1')
plt.ylabel('f2')
plt.title('Pareto Front (Weighted Sum)')
plt.grid(True)
plt.show()

参数	说明	示例
w1, w2	目标权重，反映偏好	w1=0.3, w2=0.7
求解方法	加权和法、ε约束法、NSGA‑II等	weighted_sum
Pareto前沿	无法改进任一目标而不损害其他目标的解集

灰色预测 GM(1,1)

import numpy as np

def gm11(x0, predict_num=5):
    """
    灰色预测 GM(1,1)
    :param x0: 原始序列，一维数组
    :param predict_num: 预测步数
    :return: 预测序列（包括历史拟合值）
    """
    # 累加生成
    x1 = np.cumsum(x0)
    # 紧邻均值生成
    z1 = (x1[:-1] + x1[1:]) / 2.0
    # 构造矩阵 B, Y
    B = np.vstack([-z1, np.ones(len(z1))]).T
    Y = x0[1:].reshape(-1, 1)
    # 最小二乘估计参数
    u = np.linalg.inv(B.T @ B) @ B.T @ Y
    a, b = u[0, 0], u[1, 0]
    # 时间响应式
    def x1_k(k):
        return (x0[0] - b / a) * np.exp(-a * k) + b / a
    # 拟合值
    fit_x1 = [x1_k(i) for i in range(len(x0))]
    fit_x0 = [fit_x1[0]] + [fit_x1[i] - fit_x1[i-1] for i in range(1, len(fit_x1))]
    # 预测值
    forecast_x1 = [x1_k(i) for i in range(len(x0), len(x0) + predict_num)]
    forecast_x0 = [forecast_x1[0] - fit_x1[-1]] + [forecast_x1[i] - forecast_x1[i-1] for i in range(1, predict_num)]
    # 精度检验
    residual = x0 - fit_x0
    S1 = np.std(x0, ddof=1)
    S2 = np.std(residual, ddof=1)
    C = S2 / S1  # 后验差比
    P = np.sum(np.abs(residual - np.mean(residual)) < 0.6745 * S1) / len(x0)  # 小误差概率
    return fit_x0, forecast_x0, a, b, C, P

# 示例
x0 = np.array([71.1, 72.4, 72.4, 72.1, 71.4, 72.0, 71.6])  # 原始数据
fit, forecast, a, b, C, P = gm11(x0, predict_num=3)
print('拟合值：', fit)
print('预测值（未来3期）：', forecast)
print('发展系数 a =', a)
print('灰色作用量 b =', b)
print('后验差比 C =', C)
print('小误差概率 P =', P)
print('精度等级：', '优' if C < 0.35 and P > 0.95 else '合格' if C < 0.5 and P > 0.8 else '勉强' if C < 0.65 and P > 0.7 else '不合格')

参数	说明	示例
x0	原始非负序列	[71.1,72.4,72.4,...]
predict_num	预测步数	3
a	发展系数，反映增长趋势	-0.01
b	灰色作用量	72.5
C	后验差比，越小越好	0.25
P	小误差概率，越大越好	0.95

时间序列 ARIMA

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 示例数据（模拟）
np.random.seed(42)
n = 100
time = pd.date_range('2020-01-01', periods=n, freq='M')
# 生成含有趋势和季节性的序列
trend = 0.05 * np.arange(n)
seasonal = 10 * np.sin(2 * np.pi * np.arange(n) / 12)
noise = np.random.normal(0, 2, n)
series = trend + seasonal + noise
ts = pd.Series(series, index=time)

# 1. 平稳性检验
def adf_test(series):
    result = adfuller(series)
    print('ADF统计量：', result[0])
    print('p值：', result[1])
    print('临界值：', result[4])
    if result[1] < 0.05:
        print('序列平稳')
    else:
        print('序列非平稳，需要差分')
    return result[1]

p_value = adf_test(ts)
if p_value >= 0.05:
    # 差分
    ts_diff = ts.diff().dropna()
    print('一阶差分后ADF检验：')
    adf_test(ts_diff)
else:
    ts_diff = ts

# 2. 确定 p, q（通过ACF/PACF）
fig, axes = plt.subplots(1, 2, figsize=(12,4))
plot_acf(ts_diff, lags=20, ax=axes[0])
plot_pacf(ts_diff, lags=20, ax=axes[1])
plt.show()

# 3. 拟合 ARIMA 模型（示例 p=1,d=1,q=1）
model = ARIMA(ts, order=(1,1,1))
model_fit = model.fit()
print(model_fit.summary())

# 4. 残差诊断
residuals = model_fit.resid
fig, axes = plt.subplots(1, 2, figsize=(12,4))
axes[0].plot(residuals)
axes[0].set_title('残差序列')
plot_acf(residuals, lags=20, ax=axes[1])
axes[1].set_title('残差ACF')
plt.show()

# 5. 预测未来10期
forecast_steps = 10
forecast = model_fit.forecast(steps=forecast_steps)
print('未来10期预测值：', forecast.values)
print('预测标准差：', np.sqrt(model_fit.forecast_variance))

参数	说明	示例
order (p,d,q)	自回归阶数、差分阶数、移动平均阶数	(1,1,1)
seasonal_order	季节性参数 (P,D,Q,s)	(1,1,1,12)
ADF检验	p<0.05表示平稳	p=0.01
ACF/PACF	确定p,q：ACF截尾定q，PACF截尾定p

线性回归

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
import matplotlib.pyplot as plt

# 生成示例数据
np.random.seed(42)
n = 100
X = np.random.randn(n, 3)  # 3个特征
true_coef = np.array([2.5, -1.3, 0.8])
y = X @ true_coef + np.random.randn(n) * 0.5  # 线性关系加噪声

# 划分训练集/测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 创建线性回归模型
model = LinearRegression()
model.fit(X_train, y_train)

# 输出系数
print('回归系数：', model.coef_)
print('截距：', model.intercept_)
print('R²（训练集）：', model.score(X_train, y_train))

# 预测测试集
y_pred = model.predict(X_test)
print('R²（测试集）：', r2_score(y_test, y_pred))
print('RMSE：', np.sqrt(mean_squared_error(y_test, y_pred)))
print('MAE：', mean_absolute_error(y_test, y_pred))

# 可视化预测 vs 真实值
plt.scatter(y_test, y_pred, alpha=0.7)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2)
plt.xlabel('真实值')
plt.ylabel('预测值')
plt.title('线性回归预测效果')
plt.grid(True)
plt.show()

# 使用模型进行新预测
X_new = np.array([[0.5, -0.2, 1.1]])
y_new = model.predict(X_new)
print('新样本预测值：', y_new)

参数	说明	示例
特征矩阵 X	自变量，形状 (n_samples, n_features)	100×3
目标向量 y	因变量，形状 (n_samples,)	100
test_size	测试集比例	0.2
评估指标	R²（越近1越好）、RMSE、MAE	R²=0.85

BP神经网络

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
import matplotlib.pyplot as plt

# 生成非线性数据
np.random.seed(42)
n = 500
X = np.random.randn(n, 5)
# 非线性关系：y = x1^2 + sin(x2) + x3*x4 + noise
y = X[:,0]**2 + np.sin(X[:,1]) + X[:,2]*X[:,3] + np.random.randn(n)*0.1

# 划分数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 标准化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 构建神经网络模型
model = keras.Sequential([
    layers.Dense(64, activation='relu', input_shape=(X_train.shape[1],)),
    layers.Dropout(0.2),
    layers.Dense(32, activation='relu'),
    layers.Dense(16, activation='relu'),
    layers.Dense(1)  # 输出层，回归问题无需激活函数
])

# 编译模型
model.compile(optimizer='adam', loss='mse', metrics=['mae'])

# 训练模型
history = model.fit(X_train_scaled, y_train, 
                    validation_split=0.2, 
                    epochs=100, 
                    batch_size=32, 
                    verbose=0)

# 绘制训练曲线
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(history.history['loss'], label='训练损失')
plt.plot(history.history['val_loss'], label='验证损失')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.legend()
plt.subplot(1,2,2)
plt.plot(history.history['mae'], label='训练MAE')
plt.plot(history.history['val_mae'], label='验证MAE')
plt.xlabel('Epoch')
plt.ylabel('MAE')
plt.legend()
plt.show()

# 评估测试集
test_loss, test_mae = model.evaluate(X_test_scaled, y_test, verbose=0)
print('测试集 MSE：', test_loss)
print('测试集 MAE：', test_mae)

# 预测
y_pred = model.predict(X_test_scaled).flatten()
print('前5个预测值：', y_pred[:5])
print('对应真实值：', y_test[:5])

# 保存模型（可选）
# model.save('bp_model.h5')

参数	说明	示例
网络结构	各层神经元数量	[64,32,16,1]
激活函数	隐藏层常用ReLU，输出层线性	'relu'
损失函数	回归常用MSE，分类用交叉熵	'mse'
优化器	Adam, SGD, RMSprop	'adam'
epochs, batch_size	训练轮数、批大小	100, 32

Logistic回归

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score, confusion_matrix, classification_report
import matplotlib.pyplot as plt
import seaborn as sns

# 生成二分类数据
np.random.seed(42)
n = 300
X = np.random.randn(n, 4)
# 生成概率：使用线性组合 + sigmoid
z = 0.8*X[:,0] - 1.2*X[:,1] + 0.5*X[:,2] + 0.3*X[:,3]
prob = 1 / (1 + np.exp(-z))
y = (prob > 0.5).astype(int)  # 根据概率生成标签

# 划分数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 标准化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 创建Logistic回归模型
model = LogisticRegression(penalty='l2', C=1.0, solver='lbfgs', max_iter=1000)
model.fit(X_train_scaled, y_train)

# 预测
y_pred = model.predict(X_test_scaled)
y_pred_proba = model.predict_proba(X_test_scaled)[:, 1]  # 正类概率

# 评估
print('准确率：', accuracy_score(y_test, y_pred))
print('精确率：', precision_score(y_test, y_pred))
print('召回率：', recall_score(y_test, y_pred))
print('F1分数：', f1_score(y_test, y_pred))
print('AUC：', roc_auc_score(y_test, y_pred_proba))
print('\n分类报告：')
print(classification_report(y_test, y_pred))

# 混淆矩阵
cm = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(6,5))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues')
plt.xlabel('预测标签')
plt.ylabel('真实标签')
plt.title('混淆矩阵')
plt.show()

# 新样本预测
X_new = np.array([[0.5, -0.2, 1.1, 0.3]])
X_new_scaled = scaler.transform(X_new)
pred_class = model.predict(X_new_scaled)
pred_prob = model.predict_proba(X_new_scaled)
print('新样本预测类别：', pred_class[0])
print('新样本预测概率（类0, 类1）：', pred_prob[0])

参数	说明	示例
penalty	正则化类型：'l1', 'l2', 'elasticnet', None	'l2'
C	正则化强度，越小正则化越强	1.0
solver	优化算法：'lbfgs', 'liblinear', 'saga'等	'lbfgs'
max_iter	最大迭代次数	1000
multi_class	多分类策略：'ovr', 'multinomial'	'ovr'

层次分析法 (AHP)

import numpy as np

def ahp_weights(judgment_matrix):
    """
    计算判断矩阵的权重向量并进行一致性检验
    :param judgment_matrix: n×n 判断矩阵
    :return: 权重向量, CI, RI, CR
    """
    n = judgment_matrix.shape[0]
    # 计算特征值和特征向量
    eigenvalues, eigenvectors = np.linalg.eig(judgment_matrix)
    max_eig = np.max(eigenvalues.real)
    max_idx = np.argmax(eigenvalues.real)
    eig_vec = eigenvectors[:, max_idx].real
    # 权重向量（归一化）
    weights = eig_vec / np.sum(eig_vec)
    # 一致性指标 CI
    CI = (max_eig - n) / (n - 1)
    # 随机一致性指标 RI（n≤15）
    RI_dict = {1:0, 2:0, 3:0.58, 4:0.90, 5:1.12, 6:1.24, 7:1.32, 8:1.41, 9:1.45, 10:1.49, 11:1.51, 12:1.48, 13:1.56, 14:1.57, 15:1.59}
    RI = RI_dict.get(n, 1.6)
    # 一致性比率 CR
    CR = CI / RI
    return weights, CI, RI, CR

# 示例：选择旅游地（目标层：景色、费用、居住、饮食、交通）
# 准则层对目标层的判断矩阵（5个准则）
judgment = np.array([
    [1, 1/2, 4, 3, 3],
    [2, 1, 7, 5, 5],
    [1/4, 1/7, 1, 1/2, 1/3],
    [1/3, 1/5, 2, 1, 1],
    [1/3, 1/5, 3, 1, 1]
])

weights, CI, RI, CR = ahp_weights(judgment)
print('准则层权重：', weights)
print('CI =', CI)
print('RI =', RI)
print('CR =', CR)
if CR < 0.1:
    print('一致性检验通过（CR < 0.1）')
else:
    print('一致性检验未通过，需调整判断矩阵')

# 假设有三个方案（A,B,C）对每个准则的判断矩阵（此处简化，实际需分别构造）
# 通常需要构建5个3×3矩阵，计算每个准则下方案的权重，然后合成总排序
print('\n假设方案层权重矩阵已计算，总排序 = 准则权重 × 方案权重矩阵')

参数	说明	示例
判断矩阵	正互反矩阵，aij表示i相对于j的重要性（1‑9标度）	见示例
权重向量	各准则的相对重要性	[0.3,0.5,0.1,0.05,0.05]
CI	一致性指标，CI=(λ_max‑n)/(n‑1)	0.03
RI	随机一致性指标（查表）	1.12 (n=5)
CR	一致性比率，CR=CI/RI，应<0.1	0.027

熵权法

import numpy as np

def entropy_weight(data):
    """
    熵权法计算指标权重
    :param data: m×n 矩阵，m个样本，n个指标（要求指标为正向指标）
    :return: 权重向量 (n,)
    """
    # 数据标准化（极差法，转化为[0,1]）
    min_val = np.min(data, axis=0)
    max_val = np.max(data, axis=0)
    norm_data = (data - min_val) / (max_val - min_val + 1e-8)  # 防止除零
    # 计算比重
    p = norm_data / np.sum(norm_data, axis=0, keepdims=True)
    # 计算熵值
    k = 1 / np.log(data.shape[0])
    entropy = -k * np.sum(p * np.log(p + 1e-8), axis=0)  # 防止log(0)
    # 差异系数
    d = 1 - entropy
    # 权重
    weights = d / np.sum(d)
    return weights

# 示例：5个样本，4个指标（均为正向指标，越大越好）
data = np.array([
    [100, 80, 90, 70],
    [90, 85, 88, 75],
    [95, 78, 92, 80],
    [85, 90, 85, 85],
    [80, 82, 95, 90]
])

weights = entropy_weight(data)
print('各指标权重：', weights)
print('权重和：', np.sum(weights))

# 若存在负向指标（越小越好），需先正向化
# 方法：取倒数或 max‑x 变换
# 例如，对于负向指标列 j：data[:, j] = max(data[:, j]) - data[:, j]

参数	说明	示例
data	原始数据矩阵，m×n	5×4
正向指标	越大越好，需统一	成绩、收益
负向指标	越小越好，需正向化	成本、污染
熵值	反映指标变异程度，熵越大变异越小	0.92
差异系数	d=1‑e，越大权重越大	0.08

TOPSIS（逼近理想解排序法）

import numpy as np

def topsis(data, weights=None, positive_idx=None):
    """
    TOPSIS 综合评价
    :param data: m×n 矩阵，m个方案，n个指标
    :param weights: 权重向量，长度n（若None则等权）
    :param positive_idx: 正向指标列表（默认全为正向），负向指标需在列表中剔除
    :return: 相对接近度，排序索引
    """
    m, n = data.shape
    # 默认等权
    if weights is None:
        weights = np.ones(n) / n
    # 默认全为正向指标
    if positive_idx is None:
        positive_idx = list(range(n))
    # 1. 数据归一化（向量归一化）
    norm_data = data / np.sqrt(np.sum(data**2, axis=0))
    # 2. 加权规范化矩阵
    weighted = norm_data * weights
    # 3. 正负理想解
    positive_ideal = np.zeros(n)
    negative_ideal = np.zeros(n)
    for j in range(n):
        if j in positive_idx:
            positive_ideal[j] = np.max(weighted[:, j])
            negative_ideal[j] = np.min(weighted[:, j])
        else:  # 负向指标
            positive_ideal[j] = np.min(weighted[:, j])
            negative_ideal[j] = np.max(weighted[:, j])
    # 4. 距离计算
    d_positive = np.sqrt(np.sum((weighted - positive_ideal)**2, axis=1))
    d_negative = np.sqrt(np.sum((weighted - negative_ideal)**2, axis=1))
    # 5. 相对接近度
    closeness = d_negative / (d_positive + d_negative + 1e-8)
    # 排序（从大到小）
    rank = np.argsort(-closeness)
    return closeness, rank

# 示例：5个方案，4个指标（前3个正向，第4个负向）
data = np.array([
    [100, 80, 90, 30],   # 方案1
    [90, 85, 88, 25],    # 方案2
    [95, 78, 92, 20],    # 方案3
    [85, 90, 85, 35],    # 方案4
    [80, 82, 95, 40]     # 方案5
])
weights = np.array([0.3, 0.2, 0.3, 0.2])  # 权重
positive_idx = [0, 1, 2]  # 前三个为正向，第四个为负向

closeness, rank = topsis(data, weights, positive_idx)
print('相对接近度：', closeness)
print('排序（方案索引，从优到劣）：', rank)
print('最优方案：', rank[0] + 1)
print('最劣方案：', rank[-1] + 1)

# 可视化
import matplotlib.pyplot as plt
plt.bar(range(1, len(closeness)+1), closeness[rank])
plt.xlabel('方案排名')
plt.ylabel('相对接近度')
plt.title('TOPSIS 综合评价结果')
plt.grid(axis='y')
plt.show()

参数	说明	示例
data	原始数据矩阵，m×n	5×4
weights	指标权重，长度n	[0.3,0.2,0.3,0.2]
positive_idx	正向指标索引列表	[0,1,2]
相对接近度	值越大方案越优，范围[0,1]	0.75

模糊综合评价

import numpy as np

def fuzzy_evaluation(weight, membership_matrix, score_vector=None):
    """
    模糊综合评价
    :param weight: 权重向量 (n,)
    :param membership_matrix: 隶属度矩阵 n×k，n个因素，k个评语等级
    :param score_vector: 各评语等级的分值向量 (k,)，用于加权得分（可选）
    :return: 综合评价向量 (k,)，加权得分（若提供score_vector）
    """
    # 模糊合成（加权平均型）
    result = np.dot(weight, membership_matrix)
    # 归一化（使和为1）
    result = result / np.sum(result)
    # 若提供分值向量，计算加权得分
    if score_vector is not None:
        weighted_score = np.dot(result, score_vector)
        return result, weighted_score
    return result

# 示例：评价某产品的“质量”、“价格”、“服务”
# 因素集 U = {质量, 价格, 服务}
# 评语集 V = {很好, 好, 一般, 差}
# 权重向量（质量0.5, 价格0.3, 服务0.2）
weight = np.array([0.5, 0.3, 0.2])

# 隶属度矩阵（每一行表示一个因素对各个评语的隶属度）
membership = np.array([
    [0.6, 0.3, 0.1, 0.0],   # 质量：60%很好，30%好，10%一般，0%差
    [0.1, 0.4, 0.4, 0.1],   # 价格
    [0.2, 0.5, 0.2, 0.1]    # 服务
])

# 评语分值（很好=100，好=80，一般=60，差=40）
score = np.array([100, 80, 60, 40])

result, total_score = fuzzy_evaluation(weight, membership, score)
print('模糊综合评价向量：', result)
print('加权总分：', total_score)
print('最大隶属度评语：', ['很好','好','一般','差'][np.argmax(result)])

# 多方案比较
membership2 = np.array([
    [0.4, 0.4, 0.2, 0.0],
    [0.2, 0.5, 0.2, 0.1],
    [0.3, 0.4, 0.2, 0.1]
])
result2, score2 = fuzzy_evaluation(weight, membership2, score)
print('\n方案二综合评价向量：', result2)
print('方案二加权总分：', score2)
print('方案一 vs 方案二：', '方案一更优' if total_score > score2 else '方案二更优')

参数	说明	示例
weight	因素权重向量	[0.5,0.3,0.2]
membership_matrix	隶属度矩阵，n×k	3×4
score_vector	评语分值向量	[100,80,60,40]
合成算子	常用加权平均型，也可用取大取小型	加权平均

主成分分析 (PCA)

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# 生成示例数据：10个样本，5个指标
np.random.seed(42)
data = np.random.randn(10, 5) * 10 + np.array([50, 30, 20, 40, 60])

# 1. 标准化
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)

# 2. PCA
pca = PCA()
pca.fit(data_scaled)

# 3. 特征值、贡献率
print('各主成分方差（特征值）：', pca.explained_variance_)
print('各主成分贡献率：', pca.explained_variance_ratio_)
print('累计贡献率：', np.cumsum(pca.explained_variance_ratio_))

# 选择主成分数（累计贡献率≥85%）
n_components = np.argmax(np.cumsum(pca.explained_variance_ratio_) >= 0.85) + 1
print('建议主成分数：', n_components)

# 重新PCA，保留n个主成分
pca2 = PCA(n_components=n_components)
principal_components = pca2.fit_transform(data_scaled)
print('主成分得分矩阵形状：', principal_components.shape)

# 4. 综合得分（以第一主成分贡献率为权重加权）
weights = pca2.explained_variance_ratio_ / np.sum(pca2.explained_variance_ratio_)
comprehensive_score = np.dot(principal_components, weights)
print('综合得分：', comprehensive_score)

# 排序
rank = np.argsort(-comprehensive_score)
print('样本排序（从优到劣）：', rank)

# 可视化
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(range(1, len(pca.explained_variance_ratio_)+1), pca.explained_variance_ratio_, 'bo-')
plt.plot(range(1, len(pca.explained_variance_ratio_)+1), np.cumsum(pca.explained_variance_ratio_), 'ro-')
plt.xlabel('主成分')
plt.ylabel('贡献率')
plt.legend(['单个贡献','累计贡献'])
plt.grid(True)
plt.subplot(1,2,2)
plt.scatter(principal_components[:,0], principal_components[:,1])
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('主成分散点图')
plt.grid(True)
plt.show()

参数	说明	示例
n_components	保留的主成分数	2
explained_variance_ratio_	各主成分贡献率	[0.6,0.3,...]
累计贡献率	通常要求≥85%	0.9
主成分得分	样本在新坐标系下的坐标	矩阵

数据包络分析 (DEA)

# 安装：pip install pyDEA
import numpy as np
from pyDEA.dea_models import DEA

# 示例：5个DMU，2个输入（人力、资金），2个输出（产量、利润）
inputs = np.array([
    [10, 200],  # DMU1
    [15, 300],  # DMU2
    [12, 250],  # DMU3
    [8, 180],   # DMU4
    [20, 400]   # DMU5
])
outputs = np.array([
    [100, 50],  # DMU1
    [120, 60],  # DMU2
    [110, 55],  # DMU3
    [90, 45],   # DMU4
    [130, 65]   # DMU5
])

# 创建DEA模型（默认CCR，规模报酬不变）
dea = DEA(inputs, outputs)
# 求解
efficiencies = dea.fit()
print('各DMU的效率值（CCR）：', efficiencies)

# 若需BCC模型（规模报酬可变）
from pyDEA.dea_models import BCC
bcc = BCC(inputs, outputs)
bcc_efficiencies = bcc.fit()
print('各DMU的效率值（BCC）：', bcc_efficiencies)

# 判断规模报酬（CCR效率 vs BCC效率）
scale_efficiency = efficiencies / bcc_efficiencies
print('规模效率：', scale_efficiency)
print('规模报酬：', ['递增' if se < 1 else '不变' if se == 1 else '递减' for se in scale_efficiency])

# 无效DMU的改进目标（松弛变量）
for i in range(len(inputs)):
    if efficiencies[i] < 0.999:  # 非有效
        print(f'\nDMU{i+1} 无效，效率={efficiencies[i]:.3f}')
        # 此处可调用dea的松弛变量结果，但pyDEA未直接提供，可手动计算
        # 简化为输出改进方向
        print('  应减少输入或增加输出')

参数	说明	示例
inputs	输入矩阵，m×p	5×2
outputs	输出矩阵，m×q	5×2
模型类型	CCR（CRS）、BCC（VRS）	CCR
效率值	0‑1，1表示DEA有效	0.95
规模报酬	递增、不变、递减	递增

SIR传染病模型

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def sir_model(y, t, beta, gamma):
    """SIR模型微分方程组"""
    S, I, R = y
    dSdt = -beta * S * I
    dIdt = beta * S * I - gamma * I
    dRdt = gamma * I
    return [dSdt, dIdt, dRdt]

# 参数设置
beta = 0.3   # 感染率
gamma = 0.1  # 恢复率
R0 = beta / gamma
print('基本再生数 R0 =', R0)

# 初始条件（总人口归一化为1）
S0 = 0.99    # 易感者比例
I0 = 0.01    # 感染者比例
R0_init = 0.0  # 康复者比例
y0 = [S0, I0, R0_init]

# 时间点
t = np.linspace(0, 200, 200)

# 求解微分方程
solution = odeint(sir_model, y0, t, args=(beta, gamma))
S, I, R = solution.T

# 绘图
plt.figure(figsize=(10,6))
plt.plot(t, S, label='易感者 S')
plt.plot(t, I, label='感染者 I')
plt.plot(t, R, label='康复者 R')
plt.xlabel('时间（天）')
plt.ylabel('人口比例')
plt.title('SIR传染病模型模拟（β=0.3, γ=0.1）')
plt.legend()
plt.grid(True)
plt.show()

# 计算峰值
peak_idx = np.argmax(I)
peak_time = t[peak_idx]
peak_value = I[peak_idx]
print('感染峰值时间：', peak_time)
print('感染峰值比例：', peak_value)
print('最终易感者比例：', S[-1])

参数	说明	示例
β (beta)	感染率，每个感染者每天接触并感染易感者的概率	0.3
γ (gamma)	恢复率，感染者每天恢复的概率	0.1
R0	基本再生数，R0=β/γ，大于1则疫情扩散	3.0
S0, I0, R0_init	初始易感、感染、康复者比例	0.99,0.01,0

SEIR模型

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def seir_model(y, t, beta, sigma, gamma):
    """SEIR模型微分方程组"""
    S, E, I, R = y
    dSdt = -beta * S * I
    dEdt = beta * S * I - sigma * E
    dIdt = sigma * E - gamma * I
    dRdt = gamma * I
    return [dSdt, dEdt, dIdt, dRdt]

# 参数设置
beta = 0.4    # 感染率
sigma = 1/5   # 潜伏期转染率（潜伏期平均5天）
gamma = 1/10  # 恢复率（感染期平均10天）
R0 = beta / gamma
print('基本再生数 R0 =', R0)

# 初始条件
S0 = 0.99
E0 = 0.0
I0 = 0.01
R0_init = 0.0
y0 = [S0, E0, I0, R0_init]

# 时间点
t = np.linspace(0, 300, 300)

# 求解
solution = odeint(seir_model, y0, t, args=(beta, sigma, gamma))
S, E, I, R = solution.T

# 绘图
plt.figure(figsize=(10,6))
plt.plot(t, S, label='易感者 S')
plt.plot(t, E, label='潜伏者 E')
plt.plot(t, I, label='感染者 I')
plt.plot(t, R, label='康复者 R')
plt.xlabel('时间（天）')
plt.ylabel('人口比例')
plt.title('SEIR模型模拟（β=0.4, σ=0.2, γ=0.1）')
plt.legend()
plt.grid(True)
plt.show()

# 计算峰值
peak_idx = np.argmax(I)
peak_time = t[peak_idx]
peak_value = I[peak_idx]
print('感染峰值时间：', peak_time)
print('感染峰值比例：', peak_value)
print('最终易感者比例：', S[-1])

参数	说明	示例
β (beta)	感染率	0.4
σ (sigma)	潜伏期转染率，σ=1/潜伏期天数	0.2
γ (gamma)	恢复率	0.1
R0	基本再生数	4.0

Logistic增长模型

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 1. 解析解模拟
def logistic_analytic(t, N0, r, K):
    """Logistic方程的解析解"""
    return K / (1 + (K/N0 - 1) * np.exp(-r * t))

# 参数
N0 = 10      # 初始数量
r = 0.2      # 内禀增长率
K = 1000     # 环境容量

# 时间点
t = np.linspace(0, 50, 100)
N = logistic_analytic(t, N0, r, K)

plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.plot(t, N, 'b-', linewidth=2)
plt.xlabel('时间')
plt.ylabel('种群数量')
plt.title('Logistic增长曲线（解析解）')
plt.grid(True)

# 2. 微分方程数值解
def logistic_ode(N, t, r, K):
    """Logistic微分方程"""
    dNdt = r * N * (1 - N/K)
    return dNdt

sol = odeint(logistic_ode, N0, t, args=(r, K))
plt.subplot(1,2,2)
plt.plot(t, sol, 'r-', linewidth=2)
plt.xlabel('时间')
plt.ylabel('种群数量')
plt.title('Logistic增长曲线（数值解）')
plt.grid(True)
plt.tight_layout()
plt.show()

# 3. 参数估计（拟合实际数据）
# 假设有观测数据
t_data = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40])
N_data = np.array([10, 30, 80, 200, 500, 750, 900, 950, 980])

# 定义拟合函数
def logistic_fit(t, r, K):
    N0_fixed = 10  # 已知初始值
    return K / (1 + (K/N0_fixed - 1) * np.exp(-r * t))

# 拟合
popt, pcov = curve_fit(logistic_fit, t_data, N_data, bounds=([0.01, 500], [1, 2000]))
r_fit, K_fit = popt
print('拟合参数：r =', r_fit, ', K =', K_fit)
print('参数标准差：', np.sqrt(np.diag(pcov)))

# 绘制拟合效果
t_fine = np.linspace(0, 50, 100)
N_fit = logistic_fit(t_fine, r_fit, K_fit)
plt.figure(figsize=(8,5))
plt.scatter(t_data, N_data, color='red', label='观测数据')
plt.plot(t_fine, N_fit, 'b-', label='拟合曲线')
plt.xlabel('时间')
plt.ylabel('种群数量')
plt.title('Logistic模型拟合实际数据')
plt.legend()
plt.grid(True)
plt.show()

参数	说明	示例
r	内禀增长率，单位时间内每个个体的增长率	0.2
K	环境容量，资源限制下的最大种群数量	1000
N0	初始种群数量	10
拐点	增长速度最快的点，出现在 N=K/2 时	t=17.3

Lotka‑Volterra捕食模型

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def lotka_volterra(z, t, alpha, beta, delta, gamma):
    """Lotka‑Volterra模型"""
    x, y = z  # x:猎物, y:捕食者
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# 参数（典型值）
alpha = 0.1   # 猎物自然增长率
beta = 0.02   # 捕食者对猎物的捕食率
delta = 0.01  # 捕食者增长率（取决于捕食效率）
gamma = 0.1   # 捕食者自然死亡率

# 初始数量
x0 = 40  # 猎物
y0 = 9   # 捕食者
z0 = [x0, y0]

# 时间点
t = np.linspace(0, 200, 1000)

# 求解
sol = odeint(lotka_volterra, z0, t, args=(alpha, beta, delta, gamma))
x, y = sol.T

# 绘制时间序列
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(t, x, 'g-', label='猎物')
plt.plot(t, y, 'r-', label='捕食者')
plt.xlabel('时间')
plt.ylabel('种群数量')
plt.title('Lotka‑Volterra模型时间序列')
plt.legend()
plt.grid(True)

# 绘制相图
plt.subplot(1,2,2)
plt.plot(x, y, 'b-', alpha=0.7)
plt.scatter(x[0], y[0], color='red', label='起点')
plt.xlabel('猎物数量')
plt.ylabel('捕食者数量')
plt.title('相图（极限环）')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 平衡点分析
# 平衡点1：x=0, y=0（灭绝）
# 平衡点2：x*=gamma/delta, y*=alpha/beta
x_star = gamma / delta
y_star = alpha / beta
print('非零平衡点：猎物 =', x_star, '捕食者 =', y_star)

# 计算周期振荡的振幅和周期（近似）
peak_indices, _ = find_peaks(x)
if len(peak_indices) > 1:
    period = np.mean(np.diff(t[peak_indices]))
    print('猎物振荡周期（近似）：', period, '时间单位')

参数	说明	示例
α (alpha)	猎物自然增长率（无捕食者时）	0.1
β (beta)	捕食率，反映捕食者捕食效率	0.02
δ (delta)	捕食者增长率系数	0.01
γ (gamma)	捕食者自然死亡率	0.1
平衡点	(γ/δ, α/β)，系统围绕该点振荡	(10,5)

K‑means聚类

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score
from sklearn.datasets import make_blobs

# 生成示例数据
X, y_true = make_blobs(n_samples=300, centers=4, cluster_std=0.8, random_state=42)

# 标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 1. 肘部法选择K
inertia = []
K_range = range(1, 11)
for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init='auto')
    kmeans.fit(X_scaled)
    inertia.append(kmeans.inertia_)

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(K_range, inertia, 'bo-')
plt.xlabel('簇数 K')
plt.ylabel('误差平方和 (Inertia)')
plt.title('肘部法')
plt.grid(True)

# 2. 轮廓系数
silhouette = []
for k in K_range[1:]:  # 轮廓系数要求至少2个簇
    kmeans = KMeans(n_clusters=k, random_state=42, n_init='auto')
    labels = kmeans.fit_predict(X_scaled)
    silhouette.append(silhouette_score(X_scaled, labels))

plt.subplot(1,2,2)
plt.plot(K_range[1:], silhouette, 'ro-')
plt.xlabel('簇数 K')
plt.ylabel('轮廓系数')
plt.title('轮廓系数法')
plt.grid(True)
plt.tight_layout()
plt.show()

# 根据轮廓系数选择最佳K（越大越好）
best_k = K_range[1:][np.argmax(silhouette)]
print('建议簇数：', best_k)

# 3. 使用最佳K进行聚类
kmeans = KMeans(n_clusters=best_k, random_state=42, n_init='auto')
labels = kmeans.fit_predict(X_scaled)
centers = scaler.inverse_transform(kmeans.cluster_centers_)  # 还原到原始尺度

# 可视化
plt.figure(figsize=(8,6))
scatter = plt.scatter(X[:,0], X[:,1], c=labels, cmap='viridis', alpha=0.7)
plt.scatter(centers[:,0], centers[:,1], c='red', marker='X', s=200, label='聚类中心')
plt.xlabel('特征1')
plt.ylabel('特征2')
plt.title('K‑means聚类结果')
plt.legend()
plt.colorbar(scatter, label='簇标签')
plt.grid(True)
plt.show()

print('聚类中心（原始尺度）：')
print(centers)
print('每个簇的样本数：', np.bincount(labels))

参数	说明	示例
n_clusters	簇数K	4
init	初始化方法：'random', 'k‑means++'	'k‑means++'
max_iter	最大迭代次数	300
inertia	误差平方和，越小聚类越紧	150.2
轮廓系数	范围[-1,1]，越大聚类越好	0.65

层次聚类

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_blobs

# 生成数据
X, _ = make_blobs(n_samples=50, centers=3, cluster_std=0.8, random_state=42)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 1. 计算 linkage 矩阵
Z = linkage(X_scaled, method='ward')  # ward方法最小化簇内方差
print('Linkage矩阵形状：', Z.shape)
print('前5行：')
print(Z[:5])

# 2. 绘制树状图
plt.figure(figsize=(12,5))
dendrogram(Z, truncate_mode='lastp', p=12, show_leaf_counts=True, leaf_rotation=90)
plt.xlabel('样本索引')
plt.ylabel('距离')
plt.title('层次聚类树状图（Ward linkage）')
plt.grid(axis='y')
plt.show()

# 3. 根据距离阈值切割
threshold = 5  # 通过树状图观察选择合适的阈值
labels = fcluster(Z, t=threshold, criterion='distance')
print('聚类结果（阈值={}）：'.format(threshold))
print('簇标签：', labels)
print('簇数：', len(np.unique(labels)))

# 4. 可视化聚类结果
plt.figure(figsize=(8,6))
scatter = plt.scatter(X[:,0], X[:,1], c=labels, cmap='tab10', s=80)
plt.xlabel('特征1')
plt.ylabel('特征2')
plt.title('层次聚类结果（阈值切割）')
plt.colorbar(scatter, label='簇标签')
plt.grid(True)
plt.show()

# 5. 尝试不同 linkage 方法
methods = ['ward', 'complete', 'average', 'single']
fig, axes = plt.subplots(2, 2, figsize=(12,10))
axes = axes.flatten()
for idx, method in enumerate(methods):
    Z_m = linkage(X_scaled, method=method)
    dendrogram(Z_m, truncate_mode='lastp', p=12, ax=axes[idx])
    axes[idx].set_title(f'Linkage: {method}')
    axes[idx].set_xlabel('样本')
    axes[idx].set_ylabel('距离')
    axes[idx].grid(axis='y')
plt.tight_layout()
plt.show()

参数	说明	示例
linkage	簇间距离计算方式：ward, complete, average, single	'ward'
metric	样本间距离度量：euclidean, cityblock, cosine等	'euclidean'
criterion	切割准则：'distance', 'maxclust'	'distance'
threshold	切割阈值（距离）或簇数	5.0

蒙特卡洛模拟

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

# 示例1：估计圆周率π（投点法）
np.random.seed(42)
n_samples = 100000
# 在[0,1]×[0,1]正方形内均匀投点
x = np.random.rand(n_samples)
y = np.random.rand(n_samples)
# 计算落在单位圆内的点数
inside = (x**2 + y**2) <= 1
pi_estimate = 4 * np.sum(inside) / n_samples
print('蒙特卡洛估计 π =', pi_estimate)
print('相对误差：', abs(pi_estimate - np.pi) / np.pi * 100, '%')

# 可视化
plt.figure(figsize=(6,6))
plt.scatter(x[inside], y[inside], color='blue', s=1, alpha=0.5, label='圆内')
plt.scatter(x[~inside], y[~inside], color='red', s=1, alpha=0.5, label='圆外')
circle = plt.Circle((0,0), 1, color='black', fill=False, linewidth=2)
plt.gca().add_artist(circle)
plt.xlim(0,1)
plt.ylim(0,1)
plt.xlabel('x')
plt.ylabel('y')
plt.title('蒙特卡洛法估计π（{}个样本）'.format(n_samples))
plt.legend(loc='upper right', markerscale=10)
plt.grid(True)
plt.show()

# 示例2：风险分析（项目利润模拟）
np.random.seed(42)
n_trials = 50000
# 假设利润 = 收入 - 成本，收入服从正态分布，成本服从均匀分布
mean_revenue = 100
std_revenue = 20
cost_low = 40
cost_high = 70

revenue = np.random.normal(mean_revenue, std_revenue, n_trials)
cost = np.random.uniform(cost_low, cost_high, n_trials)
profit = revenue - cost

# 统计结果
print('\n项目利润模拟（{}次试验）'.format(n_trials))
print('平均利润：', np.mean(profit))
print('利润标准差：', np.std(profit))
print('利润小于0的概率（亏损概率）：', np.sum(profit < 0) / n_trials)
print('5%分位数（VaR）：', np.percentile(profit, 5))
print('95%分位数：', np.percentile(profit, 95))

# 绘制利润分布直方图
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.hist(profit, bins=50, density=True, alpha=0.7, color='skyblue')
plt.axvline(np.mean(profit), color='red', linestyle='--', label='均值')
plt.xlabel('利润')
plt.ylabel('密度')
plt.title('利润分布直方图')
plt.legend()
plt.grid(True)

# 累积分布函数
plt.subplot(1,2,2)
sorted_profit = np.sort(profit)
cdf = np.arange(1, n_trials+1) / n_trials
plt.plot(sorted_profit, cdf, 'b-', linewidth=2)
plt.axhline(0.05, color='red', linestyle='--', label='5%分位')
plt.axvline(np.percentile(profit, 5), color='red', linestyle='--')
plt.xlabel('利润')
plt.ylabel('累积概率')
plt.title('利润累积分布函数（CDF）')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

参数	说明	示例
n_samples	随机样本数，越多精度越高	100000
随机分布	根据问题选择：均匀、正态、指数等	normal(100,20)
输出统计量	均值、标准差、分位数、概率	均值=60
标准误差	估计精度，≈标准差/√n	0.1

支持向量机 (SVM)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score
from sklearn.datasets import make_moons, make_circles

# 生成非线性可分数据（月亮形）
X, y = make_moons(n_samples=300, noise=0.2, random_state=42)

# 划分训练集/测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 标准化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 1. 线性SVM（对比）
svm_linear = SVC(kernel='linear', C=1.0)
svm_linear.fit(X_train_scaled, y_train)
y_pred_linear = svm_linear.predict(X_test_scaled)
print('线性SVM准确率：', accuracy_score(y_test, y_pred_linear))

# 2. RBF核SVM（默认）
svm_rbf = SVC(kernel='rbf', C=1.0, gamma='scale')
svm_rbf.fit(X_train_scaled, y_train)
y_pred_rbf = svm_rbf.predict(X_test_scaled)
print('RBF核SVM准确率：', accuracy_score(y_test, y_pred_rbf))

# 3. 网格搜索调参
param_grid = {
    'C': [0.1, 1, 10, 100],
    'gamma': ['scale', 'auto', 0.1, 1, 10],
    'kernel': ['rbf', 'poly']
}
grid = GridSearchCV(SVC(), param_grid, cv=5, scoring='accuracy', n_jobs=-1)
grid.fit(X_train_scaled, y_train)
print('最佳参数：', grid.best_params_)
print('最佳交叉验证准确率：', grid.best_score_)

# 用最佳模型预测
best_svm = grid.best_estimator_
y_pred_best = best_svm.predict(X_test_scaled)
print('测试集准确率：', accuracy_score(y_test, y_pred_best))
print('\n分类报告：')
print(classification_report(y_test, y_pred_best))

# 可视化决策边界
def plot_decision_boundary(clf, X, y, title):
    h = 0.02  # 网格步长
    x_min, x_max = X[:,0].min()-0.5, X[:,0].max()+0.5
    y_min, y_max = X[:,1].min()-0.5, X[:,1].max()+0.5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                         np.arange(y_min, y_max, h))
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    
    plt.figure(figsize=(8,6))
    plt.contourf(xx, yy, Z, alpha=0.3, cmap='coolwarm')
    plt.scatter(X[:,0], X[:,1], c=y, edgecolors='k', cmap='coolwarm', s=50)
    plt.xlabel('特征1')
    plt.ylabel('特征2')
    plt.title(title)
    plt.grid(True)
    plt.show()

plot_decision_boundary(svm_linear, X_train_scaled, y_train, '线性SVM决策边界')
plot_decision_boundary(svm_rbf, X_train_scaled, y_train, 'RBF核SVM决策边界')
plot_decision_boundary(best_svm, X_train_scaled, y_train, '调参后SVM决策边界')

参数	说明	示例
C	正则化参数，越小容忍错误越多，泛化越强	1.0
kernel	核函数：'linear', 'poly', 'rbf', 'sigmoid'	'rbf'
gamma	RBF核参数，越大模型越复杂	'scale'
degree	多项式核的次数	3
支持向量	决定超平面的关键样本	数量

最短路径（Dijkstra算法）

import heapq

def dijkstra(graph, start):
    """
    使用优先队列实现Dijkstra算法
    :param graph: 邻接表，graph[u] = [(v, weight), ...]
    :param start: 源点
    :return: 最短距离字典，前驱节点字典
    """
    n = len(graph)
    dist = {i: float('inf') for i in range(n)}
    prev = {i: None for i in range(n)}
    dist[start] = 0
    # 优先队列（距离，顶点）
    pq = [(0, start)]
    
    while pq:
        current_dist, u = heapq.heappop(pq)
        if current_dist > dist[u]:
            continue
        for v, w in graph[u]:
            new_dist = dist[u] + w
            if new_dist < dist[v]:
                dist[v] = new_dist
                prev[v] = u
                heapq.heappush(pq, (new_dist, v))
    return dist, prev

def reconstruct_path(prev, target):
    """从prev字典重建路径"""
    path = []
    while target is not None:
        path.append(target)
        target = prev[target]
    return path[::-1]

# 示例图（有向图，顶点0‑5）
graph = {
    0: [(1, 4), (2, 2)],
    1: [(2, 1), (3, 5)],
    2: [(1, 1), (3, 8), (4, 10)],
    3: [(4, 2), (5, 6)],
    4: [(5, 3)],
    5: []
}

start = 0
dist, prev = dijkstra(graph, start)
print('从顶点 {} 出发的最短距离：'.format(start))
for v in range(len(graph)):
    print('到顶点 {}: 距离 = {}，路径 = {}'.format(v, dist[v], reconstruct_path(prev, v)))

# 可视化（简单文本）
print('\n邻接表：')
for u in graph:
    print('{} -> {}'.format(u, graph[u]))

参数	说明	示例
graph	邻接表，字典或列表	{0:[(1,4),(2,2)]}
start	源点	0
dist	最短距离字典	{0:0,1:3,...}
prev	前驱节点字典，用于重建路径	{1:0,2:0,...}

最小生成树（Prim算法）

import heapq

def prim(graph, start=0):
    """
    Prim算法求最小生成树（邻接表，无向图）
    :param graph: 邻接表，graph[u] = [(v, weight), ...]
    :param start: 起始顶点
    :return: 生成树边列表 [(u,v,weight), ...]，总权重
    """
    n = len(graph)
    visited = [False] * n
    mst_edges = []
    total_weight = 0
    # 优先队列存储（权重, u, v），v尚未访问
    pq = []
    visited[start] = True
    for v, w in graph[start]:
        heapq.heappush(pq, (w, start, v))
    
    while pq and len(mst_edges) < n-1:
        weight, u, v = heapq.heappop(pq)
        if visited[v]:
            continue
        visited[v] = True
        mst_edges.append((u, v, weight))
        total_weight += weight
        # 将v的未访问邻居加入队列
        for neighbor, w in graph[v]:
            if not visited[neighbor]:
                heapq.heappush(pq, (w, v, neighbor))
    
    return mst_edges, total_weight

# 示例无向图（邻接表）
graph = {
    0: [(1, 2), (3, 6)],
    1: [(0, 2), (2, 3), (3, 8), (4, 5)],
    2: [(1, 3), (4, 7)],
    3: [(0, 6), (1, 8), (4, 9)],
    4: [(1, 5), (2, 7), (3, 9)]
}

edges, total = prim(graph, start=0)
print('最小生成树边：')
for u, v, w in edges:
    print('{} -- {} (权重 {})'.format(u, v, w))
print('总权重：', total)

# 也可使用Kruskal算法（并查集）
print('\nKruskal算法实现（略）')

参数	说明	示例
graph	无向图邻接表	{0:[(1,2),(3,6)]}
start	起始顶点	0
mst_edges	生成树的边列表	[(0,1,2),...]
total_weight	最小生成树总权重	16

最大流（Ford‑Fulkerson）

from collections import deque

def bfs(capacity, source, sink, parent):
    """BFS寻找增广路径"""
    n = len(capacity)
    visited = [False] * n
    queue = deque([source])
    visited[source] = True
    while queue:
        u = queue.popleft()
        for v in range(n):
            if not visited[v] and capacity[u][v] > 0:
                visited[v] = True
                parent[v] = u
                if v == sink:
                    return True
                queue.append(v)
    return False

def edmonds_karp(capacity, source, sink):
    """
    Edmonds‑Karp算法求最大流
    :param capacity: 容量矩阵，capacity[u][v] = 边(u→v)的容量
    :param source: 源点
    :param sink: 汇点
    :return: 最大流量，流矩阵
    """
    n = len(capacity)
    # 流矩阵，初始为0
    flow = [[0] * n for _ in range(n)]
    max_flow = 0
    parent = [-1] * n
    
    while bfs(capacity, source, sink, parent):
        # 计算增广路径上的最小剩余容量
        path_flow = float('inf')
        v = sink
        while v != source:
            u = parent[v]
            path_flow = min(path_flow, capacity[u][v])
            v = u
        # 更新残量图和流矩阵
        v = sink
        while v != source:
            u = parent[v]
            capacity[u][v] -= path_flow
            capacity[v][u] += path_flow  # 反向边
            flow[u][v] += path_flow
            v = u
        max_flow += path_flow
        parent = [-1] * n  # 重置parent
    return max_flow, flow

# 示例网络（容量矩阵）
# 顶点：0(source),1,2,3(sink)
capacity = [
    [0, 10, 5, 0],
    [0, 0, 15, 10],
    [0, 0, 0, 10],
    [0, 0, 0, 0]
]
source, sink = 0, 3
max_flow, flow_matrix = edmonds_karp([row[:] for row in capacity], source, sink)
print('最大流量：', max_flow)
print('流矩阵：')
for row in flow_matrix:
    print(row)

# 最小割：在残量图中从源点可达的顶点集合
print('\n最小割（源点侧）：', [i for i in range(len(capacity)) if capacity[i][i]==0])  # 简化表示

参数	说明	示例
capacity	容量矩阵，n×n	见代码
source, sink	源点、汇点	0,3
flow matrix	流矩阵，记录每条边的流量	[[0,8,2,0],...]
最大流	源点到汇点的最大流量	15

旅行商问题（TSP）

import numpy as np
import random
import math

def distance_matrix(coords):
    """根据坐标计算欧氏距离矩阵"""
    n = len(coords)
    dist = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            if i != j:
                dist[i][j] = np.linalg.norm(coords[i] - coords[j])
    return dist

def total_distance(path, dist):
    """计算路径总距离"""
    total = 0
    for i in range(len(path)-1):
        total += dist[path[i]][path[i+1]]
    total += dist[path[-1]][path[0]]  # 回到起点
    return total

def simulated_annealing_tsp(dist, start_temp=1000, end_temp=1, alpha=0.99, max_iter=5000):
    """
    模拟退火求解TSP
    :param dist: 距离矩阵
    :param start_temp: 初始温度
    :param end_temp: 终止温度
    :param alpha: 降温系数
    :param max_iter: 每个温度下的迭代次数
    :return: 最优路径，最优距离
    """
    n = len(dist)
    # 初始解（随机排列）
    current_path = list(range(n))
    random.shuffle(current_path)
    current_dist = total_distance(current_path, dist)
    best_path = current_path[:]
    best_dist = current_dist
    
    temp = start_temp
    while temp > end_temp:
        for _ in range(max_iter):
            # 生成邻域解：随机交换两个城市
            i, j = random.sample(range(n), 2)
            new_path = current_path[:]
            new_path[i], new_path[j] = new_path[j], new_path[i]
            new_dist = total_distance(new_path, dist)
            delta = new_dist - current_dist
            # 接受准则
            if delta < 0 or random.random() < math.exp(-delta / temp):
                current_path = new_path
                current_dist = new_dist
                if current_dist < best_dist:
                    best_path = current_path[:]
                    best_dist = current_dist
        temp *= alpha  # 降温
    return best_path, best_dist

# 生成随机城市坐标
np.random.seed(42)
n_cities = 15
coords = np.random.rand(n_cities, 2) * 100
dist = distance_matrix(coords)

# 求解
best_path, best_dist = simulated_annealing_tsp(dist, start_temp=1000, end_temp=1, alpha=0.995, max_iter=1000)
print('最优路径（城市索引）：', best_path)
print('最优距离：', best_dist)

# 可视化
import matplotlib.pyplot as plt
plt.figure(figsize=(8,6))
plt.scatter(coords[:,0], coords[:,1], c='red', s=100)
for i, (x, y) in enumerate(coords):
    plt.text(x+1, y+1, str(i), fontsize=10)
# 绘制路径
tour = coords[best_path + [best_path[0]]]
plt.plot(tour[:,0], tour[:,1], 'b-', linewidth=1, alpha=0.7)
plt.xlabel('X')
plt.ylabel('Y')
plt.title('旅行商问题近似解（模拟退火）')
plt.grid(True)
plt.show()

参数	说明	示例
dist	距离矩阵，n×n	欧氏距离
start_temp	初始温度，影响接受劣解的概率	1000
alpha	降温系数，每次温度乘以alpha	0.995
max_iter	每个温度的迭代次数	1000
邻域操作	交换、逆转、插入等	交换两城市

排队论（M/M/1模型）

import numpy as np
import matplotlib.pyplot as plt

def mm1_queue(lambd, mu, n=1000):
    """
    M/M/1排队模型理论计算
    :param lambd: 到达率（单位时间平均到达数）
    :param mu: 服务率（单位时间平均服务数）
    :param n: 模拟顾客数（用于仿真）
    :return: 理论指标字典，仿真结果字典
    """
    # 理论计算
    rho = lambd / mu
    if rho >= 1:
        print('警告：系统不稳定（ρ ≥ 1）')
        return None, None
    L = rho / (1 - rho)  # 平均队长（系统中顾客数）
    Lq = rho**2 / (1 - rho)  # 平均排队长（等待服务顾客数）
    W = 1 / (mu - lambd)  # 平均逗留时间
    Wq = rho / (mu - lambd)  # 平均等待时间
    P0 = 1 - rho  # 系统空闲概率
    
    theory = {
        'ρ': rho, 'L': L, 'Lq': Lq, 'W': W, 'Wq': Wq, 'P0': P0,
        'Pn公式': f'P_n = (1-ρ)ρ^n'
    }
    
    # 离散事件仿真（简单版本）
    np.random.seed(42)
    inter_arrival = np.random.exponential(1/lambd, n)  # 到达间隔
    service_time = np.random.exponential(1/mu, n)      # 服务时间
    
    arrival_time = np.cumsum(inter_arrival)
    departure_time = np.zeros(n)
    departure_time[0] = arrival_time[0] + service_time[0]
    for i in range(1, n):
        departure_time[i] = max(arrival_time[i], departure_time[i-1]) + service_time[i]
    
    wait_time = departure_time - arrival_time - service_time
    system_time = departure_time - arrival_time
    queue_length = np.zeros(n)
    for i in range(n):
        # 计算时刻t的队长（粗略）
        t = arrival_time[i]
        queue_length[i] = np.sum((arrival_time <= t) & (departure_time > t)) - 1  # 减去正在服务的
    
    simulation = {
        '平均队长（仿真）': np.mean(queue_length),
        '平均等待时间（仿真）': np.mean(wait_time),
        '平均逗留时间（仿真）': np.mean(system_time)
    }
    return theory, simulation

# 示例：λ=4人/小时，μ=5人/小时
lambd = 4
mu = 5
theory, sim = mm1_queue(lambd, mu, n=5000)
print('M/M/1排队模型理论指标：')
for key, val in theory.items():
    print(f'{key}: {val}')
print('\n仿真结果（{}个顾客）：'.format(5000))
for key, val in sim.items():
    print(f'{key}: {val:.3f}')

# 绘制队长概率分布
rho = lambd / mu
max_n = 10
probs = [(1-rho) * rho**n for n in range(max_n+1)]
plt.figure(figsize=(8,5))
plt.bar(range(max_n+1), probs, alpha=0.7)
plt.xlabel('系统中顾客数 n')
plt.ylabel('概率 P_n')
plt.title('M/M/1系统队长概率分布（ρ={:.2f}）'.format(rho))
plt.grid(axis='y')
plt.show()

参数	说明	示例
λ (lambd)	到达率，单位时间平均到达顾客数	4
μ (mu)	服务率，单位时间平均服务顾客数	5
ρ (rho)	系统利用率，ρ=λ/μ，必须<1	0.8
L, Lq	平均队长、平均排队长	4,3.2
W, Wq	平均逗留时间、平均等待时间	1,0.8

遗传算法 (Genetic Algorithm)

import numpy as np
import matplotlib.pyplot as plt

def fitness(x):
    """目标函数（求最大值）"""
    return np.sin(x) + np.sin(2*x) + np.cos(3*x)  # 在[0, 2π]内

def genetic_algorithm(pop_size=50, n_generations=100, crossover_rate=0.8, mutation_rate=0.1):
    # 参数
    x_min, x_max = 0, 2*np.pi
    chromosome_length = 20  # 二进制编码长度
    
    # 初始化种群（二进制串）
    population = np.random.randint(2, size=(pop_size, chromosome_length))
    
    best_fitness_history = []
    avg_fitness_history = []
    
    for gen in range(n_generations):
        # 解码
        decoded = population.dot(2**np.arange(chromosome_length)[::-1]) / (2**chromosome_length - 1) * (x_max - x_min) + x_min
        # 计算适应度
        fitness_values = fitness(decoded)
        best_fitness = np.max(fitness_values)
        avg_fitness = np.mean(fitness_values)
        best_fitness_history.append(best_fitness)
        avg_fitness_history.append(avg_fitness)
        
        # 选择（轮盘赌）
        probs = fitness_values / np.sum(fitness_values)
        selected_indices = np.random.choice(pop_size, size=pop_size, p=probs)
        selected_population = population[selected_indices]
        
        # 交叉（单点交叉）
        crossover_mask = np.random.rand(pop_size) < crossover_rate
        crossover_points = np.random.randint(1, chromosome_length, size=np.sum(crossover_mask))
        new_population = selected_population.copy()
        idx = 0
        for i in range(pop_size):
            if crossover_mask[i]:
                partner = np.random.randint(pop_size)
                point = crossover_points[idx]
                new_population[i, point:] = selected_population[partner, point:]
                idx += 1
        
        # 变异（位翻转）
        mutation_mask = np.random.rand(pop_size, chromosome_length) < mutation_rate
        new_population = np.where(mutation_mask, 1 - new_population, new_population)
        
        population = new_population
    
    # 最后一代的最佳个体
    decoded = population.dot(2**np.arange(chromosome_length)[::-1]) / (2**chromosome_length - 1) * (x_max - x_min) + x_min
    fitness_values = fitness(decoded)
    best_idx = np.argmax(fitness_values)
    best_x = decoded[best_idx]
    best_f = fitness_values[best_idx]
    
    return best_x, best_f, best_fitness_history, avg_fitness_history

best_x, best_f, best_hist, avg_hist = genetic_algorithm(pop_size=50, n_generations=100)
print('最优解 x =', best_x)
print('最优值 f =', best_f)

# 绘制适应度进化曲线
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
x_plot = np.linspace(0, 2*np.pi, 400)
y_plot = fitness(x_plot)
plt.plot(x_plot, y_plot, 'b-', label='目标函数')
plt.axvline(best_x, color='red', linestyle='--', label='GA最优解')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('目标函数与GA解')
plt.legend()
plt.grid(True)

plt.subplot(1,2,2)
plt.plot(best_hist, 'r-', label='最佳适应度')
plt.plot(avg_hist, 'g-', label='平均适应度')
plt.xlabel('迭代次数')
plt.ylabel('适应度')
plt.title('遗传算法进化过程')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

参数	说明	示例
pop_size	种群大小	50
n_generations	迭代代数	100
crossover_rate	交叉概率	0.8
mutation_rate	变异概率	0.1
编码方式	二进制、实数、排列等	二进制

模拟退火 (Simulated Annealing)

import numpy as np
import math
import matplotlib.pyplot as plt

def target_function(x):
    """多峰函数，求最小值"""
    return np.sin(x) + 0.5 * np.sin(3*x) + 0.2 * np.cos(5*x)

def simulated_annealing(initial_temp=100, final_temp=0.1, alpha=0.95, max_iter=1000):
    x_min, x_max = -10, 10
    # 初始解
    current_x = np.random.uniform(x_min, x_max)
    current_value = target_function(current_x)
    best_x = current_x
    best_value = current_value
    
    temp = initial_temp
    history = []
    
    while temp > final_temp:
        for _ in range(max_iter):
            # 生成邻域解（随机扰动）
            neighbor_x = current_x + np.random.normal(0, 1)  # 高斯扰动
            # 边界处理
            neighbor_x = np.clip(neighbor_x, x_min, x_max)
            neighbor_value = target_function(neighbor_x)
            
            delta = neighbor_value - current_value
            # 接受准则
            if delta < 0 or np.random.rand() < math.exp(-delta / temp):
                current_x = neighbor_x
                current_value = neighbor_value
                if current_value < best_value:
                    best_x = current_x
                    best_value = current_value
            history.append((current_x, current_value, temp))
        temp *= alpha  # 降温
    
    return best_x, best_value, history

best_x, best_value, history = simulated_annealing(initial_temp=100, final_temp=0.1, alpha=0.95, max_iter=100)
print('最优解 x =', best_x)
print('最优值 f =', best_value)

# 可视化
x_plot = np.linspace(-10, 10, 1000)
y_plot = target_function(x_plot)
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(x_plot, y_plot, 'b-', label='目标函数')
plt.axvline(best_x, color='red', linestyle='--', label='SA最优解')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('模拟退火求解多峰函数')
plt.legend()
plt.grid(True)

plt.subplot(1,2,2)
history_x = [h[0] for h in history]
history_f = [h[1] for h in history]
history_t = [h[2] for h in history]
plt.scatter(history_x, history_f, c=history_t, cmap='coolwarm', s=5, alpha=0.5)
plt.colorbar(label='温度')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('搜索轨迹（颜色表示温度）')
plt.grid(True)
plt.tight_layout()
plt.show()

参数	说明	示例
initial_temp	初始温度，影响接受劣解的概率	100
final_temp	终止温度	0.1
alpha	降温系数	0.95
max_iter	每个温度的迭代次数	100
邻域生成	随机扰动方式	高斯扰动

粒子群优化 (Particle Swarm Optimization)

import numpy as np
import matplotlib.pyplot as plt

def sphere_function(x):
    """Sphere函数，最小值在原点"""
    return np.sum(x**2)

def pso(n_particles=30, n_dimensions=2, max_iter=100, w=0.7, c1=1.5, c2=1.5):
    # 搜索范围
    x_min, x_max = -5, 5
    v_max = 0.2 * (x_max - x_min)
    
    # 初始化粒子位置和速度
    positions = np.random.uniform(x_min, x_max, (n_particles, n_dimensions))
    velocities = np.random.uniform(-v_max, v_max, (n_particles, n_dimensions))
    
    # 个体最优位置和适应度
    pbest_positions = positions.copy()
    pbest_values = np.array([sphere_function(p) for p in positions])
    
    # 全局最优
    gbest_idx = np.argmin(pbest_values)
    gbest_position = pbest_positions[gbest_idx].copy()
    gbest_value = pbest_values[gbest_idx]
    
    # 记录历史
    gbest_history = [gbest_value]
    positions_history = [positions.copy()]
    
    for iter in range(max_iter):
        for i in range(n_particles):
            # 更新速度
            r1, r2 = np.random.rand(n_dimensions), np.random.rand(n_dimensions)
            velocities[i] = w * velocities[i] + \
                           c1 * r1 * (pbest_positions[i] - positions[i]) + \
                           c2 * r2 * (gbest_position - positions[i])
            # 速度边界
            velocities[i] = np.clip(velocities[i], -v_max, v_max)
            # 更新位置
            positions[i] += velocities[i]
            positions[i] = np.clip(positions[i], x_min, x_max)
            
            # 计算适应度
            fitness = sphere_function(positions[i])
            # 更新个体最优
            if fitness < pbest_values[i]:
                pbest_positions[i] = positions[i].copy()
                pbest_values[i] = fitness
                # 更新全局最优
                if fitness < gbest_value:
                    gbest_position = positions[i].copy()
                    gbest_value = fitness
        
        gbest_history.append(gbest_value)
        positions_history.append(positions.copy())
    
    return gbest_position, gbest_value, gbest_history, positions_history

gbest_pos, gbest_val, gbest_hist, pos_hist = pso(n_particles=30, n_dimensions=2, max_iter=100)
print('全局最优位置：', gbest_pos)
print('全局最优值：', gbest_val)

# 可视化
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(gbest_hist, 'r-', linewidth=2)
plt.xlabel('迭代次数')
plt.ylabel('全局最优适应度')
plt.title('PSO收敛曲线')
plt.grid(True)

# 绘制粒子运动轨迹（二维时）
if len(gbest_pos) == 2:
    plt.subplot(1,2,2)
    # 绘制等高线
    x = np.linspace(-5, 5, 100)
    y = np.linspace(-5, 5, 100)
    X, Y = np.meshgrid(x, y)
    Z = X**2 + Y**2
    plt.contour(X, Y, Z, levels=20, cmap='gray')
    # 绘制粒子轨迹（取前5个粒子）
    for i in range(min(5, len(pos_hist[0]))):
        traj_x = [pos[i][0] for pos in pos_hist]
        traj_y = [pos[i][1] for pos in pos_hist]
        plt.plot(traj_x, traj_y, '.-', alpha=0.5, markersize=3)
    plt.scatter(gbest_pos[0], gbest_pos[1], c='red', s=100, marker='*', label='全局最优')
    plt.xlabel('x1')
    plt.ylabel('x2')
    plt.title('粒子运动轨迹（二维）')
    plt.legend()
    plt.grid(True)
plt.tight_layout()
plt.show()

参数	说明	示例
n_particles	粒子数	30
n_dimensions	问题维度	2
w	惯性权重，平衡全局与局部搜索	0.7
c1, c2	加速常数，分别控制个体和群体认知	1.5,1.5
v_max	最大速度，防止粒子飞出搜索空间	2.0

元胞自动机 (Cellular Automaton)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def game_of_life(grid):
    """应用生命游戏规则更新网格"""
    new_grid = grid.copy()
    rows, cols = grid.shape
    for i in range(rows):
        for j in range(cols):
            # 计算 Moore 邻居数（8邻域）
            total = np.sum(grid[i-1:i+2, j-1:j+2]) - grid[i, j]
            # 规则
            if grid[i, j] == 1:
                if total < 2 or total > 3:
                    new_grid[i, j] = 0  # 死亡
            else:
                if total == 3:
                    new_grid[i, j] = 1  # 繁殖
    return new_grid

# 初始化网格（随机）
size = 50
grid = np.random.choice([0, 1], size=(size, size), p=[0.7, 0.3])

# 模拟多步并记录
steps = 100
history = [grid.copy()]
for _ in range(steps):
    grid = game_of_life(grid)
    history.append(grid.copy())

# 绘制演化过程
fig, axes = plt.subplots(2, 3, figsize=(12,8))
axes = axes.flatten()
for idx, step in enumerate([0, 1, 10, 30, 60, 99]):
    ax = axes[idx]
    ax.imshow(history[step], cmap='binary', interpolation='nearest')
    ax.set_title('Step {}'.format(step))
    ax.set_xticks([])
    ax.set_yticks([])
plt.suptitle('生命游戏演化（随机初始）', fontsize=16)
plt.tight_layout()
plt.show()

# 动画（可选）
def update(frame):
    plt.clf()
    plt.imshow(history[frame], cmap='binary', interpolation='nearest')
    plt.title('Step {}'.format(frame))
    plt.axis('off')

# 如需保存动画，取消注释以下代码
# ani = FuncAnimation(plt.figure(), update, frames=len(history), interval=200, repeat=False)
# plt.show()

# 计算统计量
live_cells = [np.sum(h) for h in history]
plt.figure(figsize=(8,4))
plt.plot(live_cells, 'b-')
plt.xlabel('时间步')
plt.ylabel('活细胞数')
plt.title('活细胞数量变化')
plt.grid(True)
plt.show()

print('初始活细胞数：', live_cells[0])
print('最终活细胞数：', live_cells[-1])
print('最大活细胞数：', np.max(live_cells))

参数	说明	示例
网格大小	rows × cols	50×50
邻居类型	Moore（8邻域）、Von Neumann（4邻域）	Moore
状态	通常为0（死）和1（活）	0/1
规则	生命游戏：B3/S23（出生3/存活2或3）	自定义