import numpy as np
import pandas as pd
import tensorflow as tf
import tensorflow_probability as tfp
import math
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats
from collections import defaultdict
import random
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
# 设置TensorFlow兼容性
tf.compat.v1.enable_eager_execution()
tf.compat.v1.disable_v2_behavior()
class MDN:
def __init__(self,
n_mixtures=-1,
dist='laplace',
input_neurons=1000,
hidden_neurons=[25],
gmm_boost=False,
optimizer='adam',
learning_rate=0.001,
early_stopping=10,
tf_mixture_family=True,
input_activation='relu',
hidden_activation='tanh',
n_outputs=1):
tf.compat.v1.reset_default_graph()
self.n_mixtures = n_mixtures
self.input_neurons = input_neurons
self.hidden_neurons = hidden_neurons
self.gmm_boost = gmm_boost
self.learning_rate = learning_rate
self.early_stopping = early_stopping
self.tf_mixture_family = tf_mixture_family
self.optimizer = optimizer
self.input_activation = input_activation
self.hidden_activation = hidden_activation
self.dist = dist
self.n_outputs = n_outputs
self.x_scaler = StandardScaler()
self.y_scaler = StandardScaler()
def fit(self, X, Y, epochs, batch_size):
# 数据归一化
X = self.x_scaler.fit_transform(X)
Y = self.y_scaler.fit_transform(Y)
EPOCHS = epochs
BATCH_SIZE = batch_size
n = len(X)
XY = np.concatenate((X, Y), axis=1)
self._X = X.copy()
hidden_neurons = self.hidden_neurons
# 确定混合分量数量
if self.n_mixtures == -1:
lowest_bic = np.infty
bic = []
n_components_range = range(1, 7)
cv_types = ['spherical', 'tied', 'diag', 'full']
for cv_type in cv_types:
for n_components in n_components_range:
gmm = GaussianMixture(n_components=n_components,
covariance_type=cv_type,
max_iter=10000)
gmm.fit(XY)
bic.append(gmm.bic(XY))
if bic[-1] < lowest_bic:
lowest_bic = bic[-1]
best_gmm = gmm
self.n_mixtures = n_components
self.n_mixtures = max(2, self.n_mixtures) # 至少2个混合分量
self._y = Y.copy()
# TensorFlow数据管道
dataset = tf.compat.v1.data.Dataset \
.from_tensor_slices((X, Y)) \
.repeat(EPOCHS).shuffle(len(X)).batch(BATCH_SIZE)
iter_ = tf.compat.v1.data.make_one_shot_iterator(dataset)
x, y = iter_.get_next()
K = self.n_mixtures
self.K = K
self.x = x
# 激活函数选择
activations = {
'crelu': tf.nn.crelu,
'relu6': tf.nn.relu6,
'elu': tf.nn.elu,
'selu': tf.nn.selu,
'leaky_relu': tf.nn.leaky_relu,
'relu': tf.nn.relu,
'swish': tf.nn.swish,
'tanh': tf.nn.tanh,
'linear': None,
'softplus': tf.nn.softplus,
'sigmoid': tf.nn.sigmoid,
'softmax': tf.nn.softmax
}
input_actv = activations.get(self.input_activation.lower(), tf.nn.relu)
h_actv = activations.get(self.hidden_activation.lower(), tf.nn.relu)
n_layer = len(hidden_neurons)
# 神经网络架构
if n_layer < 1:
self.layer_last = tf.layers.dense(x, units=self.input_neurons, activation=input_actv)
else:
self.layer_1 = tf.layers.dense(x, units=self.input_neurons, activation=input_actv)
for i in range(2, n_layer + 2):
n_neurons = hidden_neurons[i - 2]
if i == n_layer + 1:
exec(f'self.layer_last = tf.layers.dense(self.layer_{i - 1}, units=n_neurons, activation=h_actv)')
else:
exec(f'self.layer_{i} = tf.layers.dense(self.layer_{i - 1}, units=n_neurons, activation=h_actv)')
# 输出层 - 修改为适应多目标
self.mu = tf.layers.dense(self.layer_last, units=self.n_outputs * K, activation=None, name="mu")
self.var = tf.exp(tf.layers.dense(self.layer_last, units=self.n_outputs * K, activation=None, name="sigma"))
self.pi = tf.layers.dense(self.layer_last, units=K, activation=tf.nn.softmax, name="mixing")
# 损失函数
if self.tf_mixture_family:
self.mixture_distribution = tfp.distributions.Categorical(probs=self.pi)
if self.dist.lower() == 'normal':
self.distribution = tfp.distributions.Normal(
loc=tf.reshape(self.mu, [-1, K, self.n_outputs]),
scale=tf.reshape(self.var, [-1, K, self.n_outputs])
)
elif self.dist.lower() in ['laplacian', 'laplace']:
self.distribution = tfp.distributions.Laplace(
loc=tf.reshape(self.mu, [-1, K, self.n_outputs]),
scale=tf.reshape(self.var, [-1, K, self.n_outputs])
)
else:
self.distribution = tfp.distributions.Normal(
loc=tf.reshape(self.mu, [-1, K, self.n_outputs]),
scale=tf.reshape(self.var, [-1, K, self.n_outputs])
)
self.likelihood = tfp.distributions.MixtureSameFamily(
mixture_distribution=self.mixture_distribution,
components_distribution=self.distribution)
# 重塑y以匹配输出维度
y_reshaped = tf.reshape(y, [-1, 1, self.n_outputs])
self.log_likelihood = -self.likelihood.log_prob(y_reshaped)
self.mean_loss = tf.reduce_mean(self.log_likelihood)
else:
# 简化的损失函数实现
self.mean_loss = tf.reduce_mean(tf.square(y - self.mu[:, 0:self.n_outputs]))
# 优化器
self.global_step = tf.Variable(0, trainable=False)
optimizers = {
'adam': tf.compat.v1.train.AdamOptimizer,
'adadelta': tf.compat.v1.train.AdadeltaOptimizer,
'adagrad': tf.compat.v1.train.AdagradOptimizer,
'gradientdescent': tf.compat.v1.train.GradientDescentOptimizer,
'rmsprop': tf.compat.v1.train.RMSPropOptimizer
}
optimizer_class = optimizers.get(self.optimizer.lower(), tf.compat.v1.train.AdamOptimizer)
self.train_op = optimizer_class(learning_rate=self.learning_rate).minimize(self.mean_loss)
self.init = tf.compat.v1.global_variables_initializer()
# 训练模型
self.sess = tf.compat.v1.Session()
self.sess.run(self.init)
best_loss = 1e+10
self.stopping_step = 0
for i in range(EPOCHS * (n // BATCH_SIZE)):
_, loss = self.sess.run([self.train_op, self.mean_loss])
if loss < best_loss:
self.stopping_step = 0
best_loss = loss
print(f"Epoch: {i} Loss: {loss:.3f}")
else:
self.stopping_step += 1
if self.stopping_step >= self.early_stopping:
print(f"Early stopping triggered at step: {i} loss: {loss}")
break
def predict_mean(self, X_pred):
"""预测均值"""
X_pred = self.x_scaler.transform(X_pred)
mu = self.sess.run(self.mu, feed_dict={self.x: X_pred})
# 重塑mu以提取每个目标的均值
mu_reshaped = mu.reshape(-1, self.K, self.n_outputs)
# 计算加权平均
pi = self.sess.run(self.pi, feed_dict={self.x: X_pred})
weighted_mu = np.zeros((len(X_pred), self.n_outputs))
for i in range(len(X_pred)):
for j in range(self.n_outputs):
weighted_mu[i, j] = np.sum(pi[i, :] * mu_reshaped[i, :, j])
# 反归一化
weighted_mu = self.y_scaler.inverse_transform(weighted_mu)
return weighted_mu
def predict_dist(self, X_pred):
"""预测分布参数"""
X_pred = self.x_scaler.transform(X_pred)
mu, var, pi = self.sess.run([self.mu, self.var, self.pi], feed_dict={self.x: X_pred})
# 重塑参数
mu = mu.reshape(-1, self.K, self.n_outputs)
var = var.reshape(-1, self.K, self.n_outputs)
# 反归一化
mu_unscaled = np.zeros_like(mu)
for k in range(self.K):
mu_unscaled[:, k, :] = self.y_scaler.inverse_transform(mu[:, k, :])
return mu_unscaled, var, pi
class Individual:
def __init__(self):
self.solution = None
self.objective = defaultdict()
self.n = 0
self.rank = 0
self.S = []
self.distance = 0
self.feasible = True
def bound_process(self, bound_min, bound_max):
"""处理边界约束"""
for i, item in enumerate(self.solution):
if item > bound_max[i]:
self.solution[i] = bound_max[i]
elif item < bound_min[i]:
self.solution[i] = bound_min[i]
def check_constraints(self):
"""检查制造约束"""
# 壁厚约束
if self.solution[0] < 0.1 or self.solution[0] > 0.5:
self.feasible = False
# 角度约束
if self.solution[1] < 10 or self.solution[1] > 80:
self.feasible = False
# 尺寸约束
if self.solution[2] < 5 or self.solution[2] > 10:
self.feasible = False
if self.solution[3] < 10 or self.solution[3] > 20:
self.feasible = False
return self.feasible
def fast_non_dominated_sort(P):
"""非支配排序"""
F = defaultdict(list)
for p in P:
p.S = []
p.n = 0
for q in P:
if dominate(p, q):
p.S.append(q)
elif dominate(q, p):
p.n += 1
if p.n == 0:
p.rank = 1
F[1].append(p)
i = 1
while F[i]:
Q = []
for p in F[i]:
for q in p.S:
q.n -= 1
if q.n == 0:
q.rank = i + 1
Q.append(q)
i += 1
F[i] = Q
return F
def dominate(a, b):
"""判断a是否支配b"""
# 对于最小化问题，所有目标值都小于等于且至少一个严格小于
better_or_equal = True
strictly_better = False
for key in a.objective:
if a.objective[key] > b.objective[key]:
better_or_equal = False
break
elif a.objective[key] < b.objective[key]:
strictly_better = True
return better_or_equal and strictly_better
def crowding_distance_assignment(L):
"""拥挤度计算"""
l = len(L)
for i in range(l):
L[i].distance = 0
for m in L[0].objective.keys():
L.sort(key=lambda x: x.objective[m])
L[0].distance = float('inf')
L[l - 1].distance = float('inf')
f_max = L[l - 1].objective[m]
f_min = L[0].objective[m]
if f_max == f_min:
continue
for i in range(1, l - 1):
L[i].distance += (L[i + 1].objective[m] - L[i - 1].objective[m]) / (f_max - f_min)
def binary_tournament(ind1, ind2):
"""二元锦标赛选择"""
if ind1.rank != ind2.rank:
return ind1 if ind1.rank < ind2.rank else ind2
elif ind1.distance != ind2.distance:
return ind1 if ind1.distance > ind2.distance else ind2
else:
return ind1
def crossover_mutation(parent1, parent2, eta, bound_min, bound_max):
"""交叉和变异"""
poplength = len(parent1.solution)
offspring1 = Individual()
offspring2 = Individual()
offspring1.solution = np.empty(poplength)
offspring2.solution = np.empty(poplength)
# 模拟二进制交叉
for i in range(poplength):
rand = random.random()
beta = (rand * 2) ** (1 / (eta + 1)) if rand < 0.5 else (1 / (2 * (1 - rand))) ** (1.0 / (eta + 1))
offspring1.solution[i] = 0.5 * ((1 + beta) * parent1.solution[i] + (1 - beta) * parent2.solution[i])
offspring2.solution[i] = 0.5 * ((1 - beta) * parent1.solution[i] + (1 + beta) * parent2.solution[i])
# 多项式变异
for i in range(poplength):
if random.random() < 0.1: # 变异概率
mu = random.random()
delta = (2 * mu) ** (1 / (eta + 1)) if mu < 0.5 else (1 - (2 * (1 - mu)) ** (1 / (eta + 1)))
offspring1.solution[i] += delta
# 边界处理
offspring1.bound_process(bound_min, bound_max)
offspring2.bound_process(bound_min, bound_max)
return [offspring1, offspring2]
def make_new_pop(P, eta, bound_min, bound_max, objective_fun):
"""生成新种群"""
popnum = len(P)
Q = []
for _ in range(int(popnum / 2)):
# 选择父代
i, j = random.sample(range(popnum), 2)
parent1 = binary_tournament(P[i], P[j])
i, j = random.sample(range(popnum), 2)
parent2 = binary_tournament(P[i], P[j])
# 确保父代不同
while (parent1.solution == parent2.solution).all():
i, j = random.sample(range(popnum), 2)
parent2 = binary_tournament(P[i], P[j])
# 交叉变异产生子代
offspring = crossover_mutation(parent1, parent2, eta, bound_min, bound_max)
# 计算子代目标函数值
for child in offspring:
child.calculate_objective(objective_fun)
Q.append(child)
return Q
def NSGAII(mdn_ea, mdn_npr, generations=100, popnum=100, eta=1):
"""NSGA-II主算法"""
# 设计变量边界
poplength = 4 # 壁厚t, 角度θ, 水平边长L₁, 高度H
bound_min = np.array([0.1, 10, 5, 10])
bound_max = np.array([0.5, 80, 10, 20])
# 目标函数
def objective_fun(x):
individual = Individual()
individual.solution = x
if not individual.check_constraints():
return {1: float('inf'), 2: float('inf')} # 违反约束返回极大值
# 使用MDN预测
ea_pred = mdn_ea.predict_mean(x.reshape(1, -1))[0, 0]
npr_pred = mdn_npr.predict_mean(x.reshape(1, -1))[0, 0]
# 转换为最小化问题
return {1: -ea_pred, 2: npr_pred} # 最大化EA(取负)，最小化NPR
# 初始化种群
P = []
for i in range(popnum):
ind = Individual()
ind.solution = np.random.rand(poplength) * (bound_max - bound_min) + bound_min
ind.bound_process(bound_min, bound_max)
ind.calculate_objective(objective_fun)
P.append(ind)
# 非支配排序
F = fast_non_dominated_sort(P)
# 生成子代
Q = make_new_pop(P, eta, bound_min, bound_max, objective_fun)
# 主循环
for gen_cur in range(generations):
R_t = P + Q # 合并父代和子代
F = fast_non_dominated_sort(R_t)
P_n = [] # 新父代
i = 1
while len(P_n) + len(F[i]) < popnum:
crowding_distance_assignment(F[i])
P_n += F[i]
i += 1
# 按拥挤度排序并选择
F[i].sort(key=lambda x: x.distance, reverse=True)
P_n += F[i][:popnum - len(P_n)]
# 生成新子代
Q_n = make_new_pop(P_n, eta, bound_min, bound_max, objective_fun)
# 更新种群
P = P_n
Q = Q_n
# 打印进度
if gen_cur % 10 == 0:
print(f"Generation {gen_cur}: Pareto front size = {len(F[1])}")
# 提取Pareto前沿
pareto_front = [ind for ind in P if ind.rank == 1 and ind.feasible]
return pareto_front
def generate_training_data(n_samples=1000):
"""生成训练数据（替代有限元模拟）"""
# 设计变量: [壁厚t, 角度θ, 水平边长L₁, 高度H]
X = np.zeros((n_samples, 4))
X[:, 0] = np.random.uniform(0.1, 0.5, n_samples) # 壁厚: 0.1-0.5
X[:, 1] = np.random.uniform(10, 80, n_samples) # 角度: 10-80
X[:, 2] = np.random.uniform(5, 10, n_samples) # 水平边长: 5-10
X[:, 3] = np.random.uniform(10, 20, n_samples) # 高度: 10-20
# 模拟目标值（实际应用中应从有限元模拟获取）
# EA = 吸能效率，NPR = 负泊松比
Y_ea = 0.5 * X[:, 0] + 0.3 * X[:, 1] - 0.2 * X[:, 2] + 0.1 * X[:, 3] + np.random.normal(0, 0.1, n_samples)
Y_npr = -0.3 * X[:, 0] - 0.4 * X[:, 1] + 0.2 * X[:, 2] - 0.1 * X[:, 3] + np.random.normal(0, 0.05, n_samples)
# 确保EA为正，NPR为负
Y_ea = np.abs(Y_ea)
Y_npr = -np.abs(Y_npr)
return X, Y_ea, Y_npr
def plot_pareto_front(pareto_front):
"""绘制Pareto前沿"""
ea_values = [-ind.objective[1] for ind in pareto_front] # 注意：之前为了最小化取了负号
npr_values = [ind.objective[2] for ind in pareto_front]
plt.figure(figsize=(10, 6))
plt.scatter(npr_values, ea_values, alpha=0.7)
plt.xlabel('Negative Poisson Ratio (NPR)')
plt.ylabel('Energy Absorption (EA)')
plt.title('Pareto Front: EA vs NPR')
plt.grid(True)
plt.show()
def main():
"""主函数"""
print("1. 生成训练数据...")
X, Y_ea, Y_npr = generate_training_data(1000)
print("2. 训练MDN模型...")
# 划分训练测试集
X_train, X_test, Y_ea_train, Y_ea_test = train_test_split(X, Y_ea, test_size=0.2, random_state=42)
_, _, Y_npr_train, Y_npr_test = train_test_split(X, Y_npr, test_size=0.2, random_state=42)
# 训练EA模型
mdn_ea = MDN(n_mixtures=3, n_outputs=1, hidden_neurons=[20, 15])
mdn_ea.fit(X_train, Y_ea_train.reshape(-1, 1), epochs=300, batch_size=32)
# 训练NPR模型
mdn_npr = MDN(n_mixtures=3, n_outputs=1, hidden_neurons=[20, 15])
mdn_npr.fit(X_train, Y_npr_train.reshape(-1, 1), epochs=300, batch_size=32)
# 评估模型
ea_pred = mdn_ea.predict_mean(X_test)
npr_pred = mdn_npr.predict_mean(X_test)
print(f"EA模型R2: {r2_score(Y_ea_test, ea_pred):.4f}")
print(f"NPR模型R2: {r2_score(Y_npr_test, npr_pred):.4f}")
print("3. 运行NSGA-II优化...")
pareto_front = NSGAII(mdn_ea, mdn_npr, generations=50, popnum=50)
print("4. 优化结果分析...")
print(f"找到 {len(pareto_front)} 个Pareto最优解")
# 选择最优解（根据偏好）
if pareto_front:
# 选择EA最大且NPR最小的解
best_idx = 0
best_score = -pareto_front[0].objective[1] - pareto_front[0].objective[2] # 最大化EA，最小化NPR
for i, ind in enumerate(pareto_front[1:], 1):
score = -ind.objective[1] - ind.objective[2]
if score > best_score:
best_score = score
best_idx = i
best_solution = pareto_front[best_idx]
print(f"最优解参数: 壁厚={best_solution.solution[0]:.3f}, 角度={best_solution.solution[1]:.3f}, "
f"水平边长={best_solution.solution[2]:.3f}, 高度={best_solution.solution[3]:.3f}")
# 预测性能
ea_pred = mdn_ea.predict_mean(best_solution.solution.reshape(1, -1))[0, 0]
npr_pred = mdn_npr.predict_mean(best_solution.solution.reshape(1, -1))[0, 0]
print(f"预测性能: EA={ea_pred:.4f}, NPR={npr_pred:.4f}")
# 绘制Pareto前沿
plot_pareto_front(pareto_front)
else:
print("未找到可行解")
if __name__ == "__main__":
main()
File "E:/test_p/learning/main.py", line 574, in <module>
main()
File "E:/test_p/learning/main.py", line 527, in main
mdn_ea.fit(X_train, Y_ea_train.reshape(-1, 1), epochs=300, batch_size=32)
File "E:/test_p/learning/main.py", line 156, in fit
components_distribution=self.distribution)
File "E:\py\Anaconda3\envs\py37tf1_env\lib\site-packages\decorator.py", line 232, in fun
return caller(func, *(extras + args), **kw)
File "E:\py\Anaconda3\envs\py37tf1_env\lib\site-packages\tensorflow_probability\python\distributions\distribution.py", line 276, in wrapped_init
default_init(self_, *args, **kwargs)
File "E:\py\Anaconda3\envs\py37tf1_env\lib\site-packages\tensorflow_probability\python\distributions\mixture_same_family.py", line 224, in __init__
"({})".format(km, kc))
ValueError: `mixture_distribution components` (3) does not equal `components_distribution.batch_shape[-1]` (1)
本文介绍了一个使用boost::ratio_not_equal的C++测试程序，展示了如何利用Boost库比较两个比例类型的不等性。通过定义1/2和2/4的比例类型并进行比较，结果显示它们相等。Boost的比例类型提供了丰富的操作，便于处理分数和比例的计算。
订阅专栏 解锁全文
扫一扫
370
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
