import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import matplotlib as mpl
# 设置中文字体
mpl.rcParams["font.family"] = ["Microsoft YaHei", "SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
mpl.rcParams["text.usetex"] = False # 禁用LaTeX避免额外依赖
# --------------------------
# 1. 系统参数与场景2数据读取
# --------------------------
# 系统参数（附录1）
M = 2800 # 总质量：车身2000kg + 8个偏心块(100kg×8)
m0 = 100 # 单个偏心块质量(kg)
r = 0.2 # 旋转半径(m)
k = 7225200 # 等效刚度(N/m)
c = 3600 # 等效阻尼(N·s/m)
omega_max = 100 # 最大角速度(rad/s)
alpha_max = 5000 # 最大角加速度(rad/s²)
T = 10 # 总时间(s)
t_eval = np.linspace(0, T, 1000) # 时间轴
dt = t_eval[1] - t_eval[0] # 时间步长
# 读取场景2数据
data = pd.read_excel("场景2.xlsx", skiprows=[1])
t_sam = data["采样时刻"].values
F_dis_sam = data["横向扰动力(单位：N)"].values
def F_dis(t):
return np.interp(t, t_sam, F_dis_sam)
# --------------------------
# 2. 粒子群优化（PSO）算法优化PID参数
# --------------------------
class PSO:
def __init__(self, func, dim, x_min, x_max, pop_size=30, max_iter=50):
self.func = func # 适应度函数（目标：最小化Ih）
self.dim = dim # 参数维度（Kp, Kd, Ki → 3）
self.x_min = x_min # 参数最小值
self.x_max = x_max # 参数最大值
self.pop_size = pop_size # 粒子数量
self.max_iter = max_iter # 迭代次数
# 初始化粒子位置和速度
self.x = np.random.uniform(x_min, x_max, (pop_size, dim)) # 位置
self.v = np.random.uniform(-1, 1, (pop_size, dim)) # 速度
self.pbest_x = self.x.copy() # 个体最优位置
self.pbest_f = np.array([np.inf] * pop_size) # 个体最优适应度
self.gbest_x = np.zeros(dim) # 全局最优位置
self.gbest_f = np.inf # 全局最优适应度
self.iter_history = [] # 记录迭代过程
def update(self):
"""执行PSO优化，返回最优参数和最优适应度"""
w = 0.8 - 0.5 * (np.arange(self.max_iter) / self.max_iter) # 惯性权重衰减
c1, c2 = 2, 2 # 学习因子
for iter in range(self.max_iter):
# 计算适应度
f_values = np.array([self.func(x) for x in self.x])
# 更新个体最优和全局最优
for i in range(self.pop_size):
if f_values[i] < self.pbest_f[i]:
self.pbest_f[i] = f_values[i]
self.pbest_x[i] = self.x[i].copy()
if f_values[i] < self.gbest_f:
self.gbest_f = f_values[i]
self.gbest_x = self.x[i].copy()
# 记录历史
self.iter_history.append(self.gbest_f)
# 更新速度和位置
r1, r2 = np.random.rand(self.pop_size, self.dim), np.random.rand(self.pop_size, self.dim)
self.v = w[iter] * self.v + c1 * r1 * (self.pbest_x - self.x) + c2 * r2 * (self.gbest_x - self.x)
self.x = np.clip(self.x + self.v, self.x_min, self.x_max) # 限制位置在范围内
# 打印迭代信息
if (iter + 1) % 10 == 0:
print(f"PSO迭代 {iter + 1}/{self.max_iter}：最优Ih = {self.gbest_f:.6f}")
return self.gbest_x, self.gbest_f
# --------------------------
# 3. 适应度函数（计算给定PID参数的振动指数Ih）
# --------------------------
def pid_fitness(params):
Kp, Kd, Ki = params # 待优化的PID参数
# 初始化状态变量
integral_y = 0.0
y_prev, dot_y_prev = 0.01, 0.0 # 初始位移和速度
theta = [0.0, 0.0, 0.0, 0.0] # 4组角位移
omega = [0.0, 0.0, 0.0, 0.0] # 4组角速度
# 存储系统状态
y_list, dot_y_list = [], []
for t in t_eval:
# 计算当前无控制时的加速度
y = y_prev
dot_y = dot_y_prev
ddot_y_unctrl = (F_dis(t) - k * y - c * dot_y) / M
# PID控制计算
integral_y += y * dt
integral_y = np.clip(integral_y, -1e-3, 1e-3) # 抗积分饱和
F_act = - (Kp * y + Kd * dot_y + Ki * integral_y)
# 力转换为角位移（含约束）
F_per_group = [F_act / 4] * 4 # 4组力的列表（修正：确保可迭代）
omega_target = []
for f in F_per_group:
if f <= 0:
omega_t = np.sqrt(-f / (2 * m0 * r)) if f != 0 else 0.0
else:
omega_t = np.sqrt(f / (2 * m0 * r))
omega_t = np.clip(omega_t, -omega_max, omega_max)
omega_target.append(omega_t)
# 角加速度约束
for i in range(4):
alpha = (omega_target[i] - omega[i]) / dt
alpha = np.clip(alpha, -alpha_max, alpha_max)
omega[i] += alpha * dt
theta[i] += omega[i] * dt
theta[i] %= (2 * np.pi) # 归一化
# 计算实际控制力
F_act_real = 0.0
for i in range(4):
F_act_real += 2 * m0 * r * (omega[i] ** 2) * np.sin(theta[i])
# 更新系统状态
ddot_y = (F_dis(t) - k * y - c * dot_y + F_act_real) / M
dot_y_new = dot_y + ddot_y * dt
y_new = y + dot_y_new * dt
# 保存状态
y_list.append(y_new)
dot_y_list.append(dot_y_new)
# 更新历史值
y_prev, dot_y_prev = y_new, dot_y_new
# 计算振动指数Ih
y_dot = np.array(dot_y_list)
y_ddot = np.gradient(y_dot, t_eval)
Ih = (1 / T) * np.trapz(y_ddot ** 2, t_eval)
return Ih
# --------------------------
# 4. 执行PSO优化PID参数
# --------------------------
# PID参数搜索范围（根据经验设定）
x_min = [1e4, 1e3, 1e2] # Kp, Kd, Ki最小值
x_max = [1e6, 1e5, 1e4] # Kp, Kd, Ki最大值
# 初始化PSO并优化（修正：调用update方法）
pso = PSO(
func=pid_fitness,
dim=3,
x_min=x_min,
x_max=x_max,
pop_size=30,
max_iter=50
)
best_params, best_Ih = pso.update() # 关键修正：方法名从optimize改为update
Kp_opt, Kd_opt, Ki_opt = best_params
print("\n===== PSO优化结果 =====")
print(f"最优PID参数：Kp = {Kp_opt:.0f}, Kd = {Kd_opt:.0f}, Ki = {Ki_opt:.0f}")
print(f"优化后的最小振动指数Ih = {best_Ih:.6f} (m/s²)²")
# --------------------------
# 5. 用优化后的参数重新仿真并输出结果
# --------------------------
# 重置全局变量，用最优参数仿真
integral_y = 0.0
y_prev, dot_y_prev = 0.01, 0.0
theta = [0.0, 0.0, 0.0, 0.0]
omega = [0.0, 0.0, 0.0, 0.0]
theta_history = []
y_list, dot_y_list = [], []
for t in t_eval:
# 系统状态
y = y_prev
dot_y = dot_y_prev
# PID控制（用优化参数）
integral_y += y * dt
integral_y = np.clip(integral_y, -1e-3, 1e-3)
F_act = - (Kp_opt * y + Kd_opt * dot_y + Ki_opt * integral_y)
# 力转角位移
F_per_group = [F_act / 4] * 4
omega_target = []
for f in F_per_group:
if f <= 0:
omega_t = np.sqrt(-f / (2 * m0 * r)) if f != 0 else 0.0
else:
omega_t = np.sqrt(f / (2 * m0 * r))
omega_t = np.clip(omega_t, -omega_max, omega_max)
omega_target.append(omega_t)
# 角加速度约束
for i in range(4):
alpha = (omega_target[i] - omega[i]) / dt
alpha = np.clip(alpha, -alpha_max, alpha_max)
omega[i] += alpha * dt
theta[i] += omega[i] * dt
theta[i] %= (2 * np.pi)
theta_history.append(theta.copy())
# 实际控制力与状态更新
F_act_real = 0.0
for i in range(4):
F_act_real += 2 * m0 * r * (omega[i] ** 2) * np.sin(theta[i])
ddot_y = (F_dis(t) - k * y - c * dot_y + F_act_real) / M
dot_y_new = dot_y + ddot_y * dt
y_new = y + dot_y_new * dt
y_list.append(y_new)
dot_y_list.append(dot_y_new)
y_prev, dot_y_prev = y_new, dot_y_new
# 计算振动指数
y_dot = np.array(dot_y_list)
y_ddot_with_act = np.gradient(y_dot, t_eval)
Ih_with_act = (1 / T) * np.trapz(y_ddot_with_act ** 2, t_eval)
# 无控制时的振动指数（对比用）
def no_control_simulation():
y, dot_y = 0.01, 0.0
y_dot_list = []
for t in t_eval:
ddot_y = (F_dis(t) - k * y - c * dot_y) / M
dot_y += ddot_y * dt
y += dot_y * dt
y_dot_list.append(dot_y)
y_ddot = np.gradient(y_dot_list, t_eval)
return (1 / T) * np.trapz(y_ddot ** 2, t_eval)
Ih_no_act = no_control_simulation()
# --------------------------
# 6. 结果输出与绘图
# --------------------------
print("\n===== 问题3最终结果 =====")
print(f"1. 无控制时振动指数 Ih = {Ih_no_act:.6f} (m/s²)²")
print(f"2. 优化控制后振动指数 Ih = {Ih_with_act:.6f} (m/s²)²")
print(f"3. 振动抑制率 = {(Ih_no_act - Ih_with_act) / Ih_no_act * 100:.2f}%")
# 图1：PSO优化迭代曲线
plt.figure(figsize=(10, 6))
plt.plot(range(1, pso.max_iter + 1), pso.iter_history, color='blue', linewidth=2)
plt.xlabel('迭代次数', fontsize=12)
plt.ylabel('最优振动指数 Ih', fontsize=12)
plt.title('PSO优化PID参数的迭代过程', fontsize=14)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('PSO优化过程.png', dpi=300)
plt.show()
# 图2：横向位移对比
plt.figure(figsize=(10, 6))
# 无控制位移
y_no_act = []
y, dot_y = 0.01, 0.0
for t in t_eval:
ddot_y = (F_dis(t) - k * y - c * dot_y) / M
dot_y += ddot_y * dt
y += dot_y * dt
y_no_act.append(y)
plt.plot(t_eval, y_no_act, label='无控制', linestyle='--', color='gray', linewidth=2)
plt.plot(t_eval, y_list, label='PSO优化PID控制', color='blue', linewidth=2)
plt.xlabel('时间 t (s)', fontsize=12)
plt.ylabel('横向位移 y (m)', fontsize=12)
plt.title('车身横向位移对比（场景2）', fontsize=14)
plt.legend(fontsize=12)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('场景2_位移曲线.png', dpi=300)
plt.show()
# 图3：横向加速度对比
y_ddot_no_act = np.gradient(np.array(y_no_act), t_eval)
plt.figure(figsize=(10, 6))
plt.plot(t_eval, y_ddot_no_act, label='无控制', linestyle='--', color='gray', linewidth=2)
plt.plot(t_eval, y_ddot_with_act, label='PSO优化PID控制', color='red', linewidth=2)
plt.xlabel('时间 t (s)', fontsize=12)
plt.ylabel(r'横向加速度 $\ddot{y}$ (m/s²)', fontsize=12)
plt.title('车身横向加速度对比（场景2）', fontsize=14)
plt.legend(fontsize=12)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('场景2_加速度曲线.png', dpi=300)
plt.show()
# 图4：作动器角位移曲线
theta_history = np.array(theta_history)
plt.figure(figsize=(10, 6))
labels = ['作动器1第1组', '作动器1第2组', '作动器2第1组', '作动器2第2组']
colors = ['red', 'blue', 'green', 'orange']
for i in range(4):
plt.plot(t_eval, theta_history[:, i], label=labels[i], color=colors[i], linewidth=1.5)
plt.xlabel('时间 t (s)', fontsize=12)
plt.ylabel(r'角位移 $\theta$ (rad)', fontsize=12)
plt.title('离心式作动器角位移变化（场景2）', fontsize=14)
plt.legend(fontsize=10)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('场景2_角位移曲线.png', dpi=300)
plt.show()
解读该程序，现在我的输出结果如下：
E:\Anaconda\python.exe "E:\Users\Master\Desktop\SHUWEICUP2534088\problem 2\solution.py"
PSO迭代 10/50：最优Ih = 5130.602517
PSO迭代 20/50：最优Ih = 5130.602517
PSO迭代 30/50：最优Ih = 5130.602517
PSO迭代 40/50：最优Ih = 5130.602517
PSO迭代 50/50：最优Ih = 5130.602517
===== PSO优化结果 =====
最优PID参数：Kp = 10000, Kd = 1000, Ki = 100
优化后的最小振动指数Ih = 5130.602517 (m/s²)²
===== 问题3最终结果 =====
1. 无控制时振动指数 Ih = 5129.593284 (m/s²)²
2. 优化控制后振动指数 Ih = 5130.602517 (m/s²)²
3. 振动抑制率 = -0.02%
进程已结束，退出代码为 0
目前出现一下问题：1.PSO优化算法过于僵硬，输出的最优PID参数过于刻意 2.且优化后优化控制后振动指数Ih反而增大了，这与我要解决的问题背道而驰
最新发布
本文介绍了一种使用迭代法计算平方根的C语言程序实现。该程序通过不断逼近的方法求解给定数值的平方根，并确保计算精度达到指定阈值。输入一个正数X，程序将输出其平方根，结果精确到小数点后三位。
扫一扫
1986
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
