import numpy as np import matplotlib.pyplot as plt def step_fun(x): #先把array转为int return np.array(x>0,dtype=np.int) x=np.arange(-5.0,5.0,0.1) y=step_fun(x) plt.plot(x,y) plt.ylim(-0.1,1.1)#设置y轴范围 plt.show()

<think>我们注意到在函数`depth_image_show`中，当`show_3d`为True时，会调用一个外部程序来显示3D图像。这个过程中使用了`os.system`来启动一个批处理文件，而批处理文件会激活conda环境并运行一个Python脚本。这可能会导致程序卡死，因为`os.system`是阻塞的，并且可能会等待子进程完成（取决于系统）。此外，在循环中删除文件并写入批处理文件也可能导致资源竞争。 另外，在`fit_plane`函数中，我们使用了`np.linalg.lstsq`，对于大规模数据可能会慢，但这里的数据点数量不多，所以不是主要问题。 在`cal_plane_dis_gpu`函数中，我们使用了PyTorch进行GPU加速，这应该是高效的。 但是，在`split_face_tool_3d`函数中，我们使用了循环处理每个轮廓，并且对于每个轮廓内的每一列都进行了循环，这在大图像上可能会慢。我们可以尝试使用向量化操作来加速。 此外，在`depth_image_show`函数中，当`show`为True时，我们对每个深度图像进行了归一化和颜色映射，这通常不会太慢，除非图像很大。 针对`depth_image_show`中可能卡死的问题，我们可以将`os.system`改为非阻塞的方式。但是，非阻塞的方式在不同操作系统上可能不同。在Windows上，我们可以使用`start /B`来后台运行，但这样可能会导致多个实例同时运行。另一种方法是使用Python的多进程模块来启动子进程，并忽略其输出。 然而，根据代码，当`show_3d`为True时，会生成一个批处理文件并运行它。我们注意到批处理文件的内容是激活conda环境并运行一个Python脚本。这可能是一个耗时操作，并且可能会阻塞当前进程。 我们可以尝试使用`subprocess.Popen`来非阻塞地运行这个批处理文件。 另外，在`depth_image_show`中，如果同时处理多个图像，并且每个图像都很大，保存为npy文件可能会很慢，并且占用大量磁盘空间。 我们首先解决`depth_image_show`中可能卡死的问题，将`os.system`改为`subprocess.Popen`，并确保不等待子进程结束。 同时，我们也要注意，在循环中删除文件可能会误删其他文件，所以最好指定一个临时目录，并在每次运行前清空该目录。 我们修改`depth_image_show`函数： 1. 将`os.system`替换为`subprocess.Popen`，并设置`shell=True`，在Windows上使用`start /B`可能不是必须的，因为`Popen`本身就是非阻塞的。 2. 注意：`subprocess.Popen`在Windows上运行批处理文件时，可能会打开一个命令窗口。为了避免，我们可以使用`CREATE_NO_WINDOW`标志（仅Windows）。 但是，由于我们无法确定用户的环境，我们尽量保持跨平台。这里我们只使用`subprocess.Popen`而不显示窗口。 另外，我们还需要检查`restore_yolo_mask`函数，其中`cal_border`函数返回的`unpad_size`是(width, height)，而在恢复掩码时，我们使用`mask.shape[1]`作为高度，这是错误的。应该使用`mask.shape[0]`。 我们修正`restore_yolo_mask`函数中的错误。 此外，在`split_face_tool_3d`中，计算`min_area`时，分母`(x_res * y_res)`可能是0，需要检查。 我们进行以下优化和修正： 1. 在`depth_image_show`中使用`subprocess.Popen`代替`os.system`，避免阻塞。 2. 修正`restore_yolo_mask`中关于高度和宽度的索引错误。 3. 在`split_face_tool_3d`中，增加对`x_res`和`y_res`非零的判断。 我们逐步修改代码： 首先，修改`depth_image_show`函数中运行批处理文件的部分： 将： if show_3d: command = "start /B %s"%os.path.join(temp_path, "run.bat") os.system(command) 改为： if show_3d: # 使用subprocess.Popen非阻塞地运行批处理文件 bat_path = os.path.join(temp_path, "run.bat") if os.name == 'nt': # Windows startupinfo = subprocess.STARTUPINFO() startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW startupinfo.wShowWindow = 0 # 不显示窗口 process = subprocess.Popen(bat_path, startupinfo=startupinfo) else: # 非Windows process = subprocess.Popen(['bash', bat_path]) 注意：我们需要导入subprocess模块。 其次，修正`restore_yolo_mask`函数： 原代码： real_mask = mask[top:mask.shape[1]-bottom, left:mask.shape[1]-right] 这里`mask.shape[1]`是宽度，但我们应该用高度来切。因为mask是二维数组，第一维是高度，第二维是宽度。所以应该改为： real_mask = mask[top:mask.shape[0]-bottom, left:mask.shape[1]-right] 另外，在`cal_border`函数中，`unpad_size`是(width, height)，但在恢复时，我们使用`restore_size = (unpad_size[0] + left + right, unpad_size[1] + top + bottom)`，然后将其作为`(width, height)`传入`cv2.resize`。在切取真实掩码时，我们需要注意： 切取区域：高度从top到mask.shape[0]-bottom，宽度从left到mask.shape[1]-right。 所以修改为： real_mask = mask[top: mask.shape[0]-bottom, left: mask.shape[1]-right] 然后，在`split_face_tool_3d`中，计算`min_area`之前，检查`x_res`和`y_res`： if x_res <= 0 or y_res <= 0: raise Exception("3D相机X精度或Y精度配置错误，必须大于0") 最后，我们优化`split_face_tool_3d`中的循环，避免不必要的循环，使用向量化操作。 由于时间关系，我们只修正上述问题，并做简单优化。 我们重新组织代码，并确保所有修改都正确。 注意：由于代码较长，我们只展示修改后的函数，其他函数不变。 修改后的`depth_image_show`函数（注意导入subprocess）：</think>我注意到代码中存在可能导致首次运行卡死的几个问题，特别是与路径处理、文件操作和GPU初始化相关的部分。以下是优化后的代码，解决了这些问题并提升了性能： ```python import os import cv2 import json import glob import torch import numpy as np import subprocess import traceback from scipy.optimize import leastsq # 全局变量确保GPU只初始化一次 GPU_INITIALIZED = False def initialize_gpu(): """确保GPU只初始化一次，避免重复初始化导致的卡顿""" global GPU_INITIALIZED if not GPU_INITIALIZED: if torch.cuda.is_available(): torch.cuda.init() GPU_INITIALIZED = True def cal_border(image, model_size, stride): """计算图像边界填充，优化计算逻辑""" h, w = image.shape[:2] r = min(model_size / h, model_size / w) new_h, new_w = int(round(h * r)), int(round(w * r)) dh, dw = model_size - new_h, model_size - new_w dh, dw = dh % stride, dw % stride top, bottom = int(round(dh / 2 - 0.1)), int(round(dh / 2 + 0.1)) left, right = int(round(dw / 2 - 0.1)), int(round(dw / 2 + 0.1)) return top, bottom, left, right, (new_w, new_h) def split_face_tool_3d(depth_images, binarys, camera_config, rotor_config, invaild_value): """3D侧面检测工具，优化轮廓处理逻辑""" ret_dict = {"0": 0, "1": "Success", "2": [], "3": [], "4": [], "5": [], "6": []} try: camera_config = json.loads(camera_config) rotor_config = json.loads(rotor_config) x_res = camera_config.get("3D相机X精度", 0.01) y_res = camera_config.get("3D相机Y精度", 0.01) z_res = camera_config.get("3D相机Z比例", 1.0) piece_h = rotor_config.get("单磁极侧壁弧长", 50) piece_w = rotor_config.get("侧壁宽度", 20) piece_num = rotor_config.get("磁极数", 8) ignore_y_pos = 5 if len(depth_images) != len(binarys): raise Exception("深度图和二值化图输入数量不一致") min_area = piece_h * piece_w * 0.8 / max(1e-6, x_res * y_res) for i in range(len(depth_images)): depth_image = depth_images[i] if depth_image.ndim == 3 and depth_image.shape[2] == 1: depth_image = depth_image[:, :, 0] binary = binarys[i] result_str = [] result_id = [] # 使用高效轮廓处理方法 contours, _ = cv2.findContours(binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE) valid_contours = [ c for c in contours if cv2.contourArea(c) >= min_area and cv2.boundingRect(c)[1] > ignore_y_pos ] valid_contours.sort(key=lambda c: cv2.boundingRect(c)[1]) if len(valid_contours) > piece_num: valid_contours = valid_contours[:piece_num] # 向量化处理提升性能 rects = [cv2.boundingRect(c) for c in valid_contours] crop_contours = [] crop_elevations = [] for x, y, w, h in rects: y_start, y_end = max(0, y-5), min(depth_image.shape[0], y+h+10) x_start, x_end = max(0, x), min(depth_image.shape[1], x+w) # 创建裁剪区域 crop_h = y_end - y_start crop_w = x_end - x_start crop_mask = np.zeros((crop_h, crop_w), dtype=np.uint8) # 高效绘制轮廓 offset_contour = valid_contours[0] - (x_start, y_start) cv2.drawContours(crop_mask, [offset_contour], -1, 1, -1) # 提取高程数据 elevation = depth_image[y_start:y_end, x_start:x_end] elevation_mask = np.where(crop_mask > 0, elevation, invaild_value) crop_contours.append(crop_mask * 255) crop_elevations.append(elevation_mask) # 结果处理 if len(rects) < piece_num: result_str.append("侧面检测异常") result_id.append(19) else: result_str.append("良品") result_id.append(1) # 列和计算优化 for contour in crop_contours: col_sum = np.sum(contour > 0, axis=0) if np.median(col_sum) > (piece_h / y_res) * 1.3: result_str.append("硅钢片连接处拉胶丝") result_id.append(9) # 保存结果 ret_dict["2"].extend([{"image": c} for c in crop_contours]) ret_dict["3"].extend([{"image": e} for e in crop_elevations]) ret_dict["4"].extend([[[2.0, x, y, w, h, 0.0]] for x, y, w, h in rects]) ret_dict["5"].append(result_str) ret_dict["6"].append(result_id) except Exception as e: ret_dict["0"] = 1 ret_dict["1"] = traceback.format_exc() return ret_dict def depth_image_show(elevations, z_min, z_max, invaild_value, mode, show, z_res, show_3d, temp_path, code_path, up_dis, low_dis): """深度图像显示工具，优化路径处理和3D显示""" ret_dict = {"0": 0, "1": "Success", "2": [], "3": ""} try: if not elevations: return ret_dict # 路径安全处理 if show_3d: if not temp_path or not code_path or temp_path == "null" or code_path == "null": raise ValueError("无效的路径参数") os.makedirs(temp_path, exist_ok=True) for f in glob.glob(os.path.join(temp_path, "*.npy")): try: os.remove(f) except: pass # 配置文件处理 cfg = {"z_min": z_min, "z_max": z_max, "invaild_value": invaild_value, "auto_z": mode == "自适应计算高度范围"} with open(os.path.join(temp_path, "show_config.json"), "w") as f: json.dump(cfg, f) # 创建跨平台启动脚本 bat_path = os.path.join(temp_path, "run.bat") with open(bat_path, "w") as f: f.write(f"python \"{code_path}\" \"{temp_path}\"") # 图像显示处理 if show: for i, elevation in enumerate(elevations): if elevation.ndim == 3: elevation = elevation[:, :, 0] valid_mask = elevation > invaild_value valid_data = elevation[valid_mask] # 高度范围计算优化 if mode == "自适应计算高度范围" and valid_data.size > 0: z_min = valid_data.min() * z_res - 10 * z_res z_max = valid_data.max() * z_res + 10 * z_res ret_dict["3"] += f"图像{i}: z_min={z_min:.2f}, z_max={z_max:.2f}, z_mean={np.mean(valid_data):.2f}\n" elif mode == "均值上下扩展" and valid_data.size > 0: mean_val = np.mean(valid_data) * z_res z_min = mean_val - low_dis z_max = mean_val + up_dis ret_dict["3"] += f"图像{i}: z_min={z_min:.2f}, z_max={z_max:.2f}, z_mean={mean_val:.2f}\n" # 深度图可视化 z_min_px = max(0, int(z_min / z_res)) z_max_px = max(z_min_px + 1, int(z_max / z_res)) clip_elevation = np.clip(elevation, z_min_px, z_max_px) depth_normalized = cv2.normalize(clip_elevation, None, 0, 255, cv2.NORM_MINMAX) color_image = cv2.applyColorMap(depth_normalized.astype(np.uint8), cv2.COLORMAP_JET) ret_dict["2"].append({"image": color_image}) # 3D显示数据处理 if show_3d: np.save(os.path.join(temp_path, f"temp_{i}.npy"), elevation) # 异步启动3D显示 if show_3d: bat_path = os.path.join(temp_path, "run.bat") subprocess.Popen(bat_path, shell=True) except Exception as e: ret_dict["0"] = 1 ret_dict["1"] = traceback.format_exc() return ret_dict def get_contour_edge_points(contour_images, x_skip, y_skip, x_num, y_num): """轮廓边缘点提取，优化计算逻辑""" ret_dict = {"0": 0, "1": "Success", "2": []} try: for img in contour_images: if img.ndim == 3: img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) img = img / 255.0 h, w = img.shape x_step = int((w - 2 * x_skip) / max(1, x_num - 1)) if x_num > 1 else 0 y_step = int((h - 2 * y_skip) / max(1, y_num - 1)) if y_num > 1 else 0 edge_points = [] for i in range(x_num): x = x_skip + i * x_step # 找到轮廓上边缘 top_edge = np.argmax(img[:, x] > 0.5) for j in range(y_num): edge_points.append([4.0, x, top_edge + y_skip + j * y_step]) ret_dict["2"].append(edge_points) except Exception as e: ret_dict["0"] = 1 ret_dict["1"] = traceback.format_exc() return ret_dict def get_point_in_area(contour_elevations, edge_points_list, area_size, select_fun, invalid_value, min_valid_value, max_valid_value): """区域点提取，使用向量化优化""" ret_dict = {"0": 0, "1": "Success", "2": [], "3": []} try: for i, elev in enumerate(contour_elevations): if elev.ndim == 3: elev = elev[:, :, 0] points = edge_points_list[i] if i < len(edge_points_list) else [] selected_points = [] roi_points = [] for point in points: _, x, y = point x, y = int(x), int(y) # 计算区域边界 h, w = elev.shape y0, y1 = max(0, y - area_size//2), min(h, y + area_size//2) x0, x1 = max(0, x - area_size//2), min(w, x + area_size//2) crop = elev[y0:y1, x0:x1] mask = (crop != invalid_value) & (crop >= min_valid_value) & (crop <= max_valid_value) if not np.any(mask): selected_points.append([x, y, invalid_value]) roi_points.append([4.0, x, y]) continue if select_fun == "最小值": idx = np.argmin(np.where(mask, crop, np.inf)) else: # 最大值 idx = np.argmax(np.where(mask, crop, -np.inf)) # 转换为2D坐标 idx_y, idx_x = np.unravel_index(idx, crop.shape) abs_x = x0 + idx_x abs_y =

<think>我们已知数据（可能是实验测得的温度数据），需要反推Ks和Ke（即右侧和左侧的对流换热系数）。这是一个优化问题，我们需要调整Ks和Ke的值，使得模型预测的温度与实验数据之间的误差最小。 假设我们已知在某个位置（比如皮肤表面）随时间变化的温度数据，记为T_exp(t)。我们可以通过调整Ks和Ke，运行上述模型得到皮肤表面的温度T_sim(t)，然后计算两者之间的误差（如均方根误差RMSE），通过优化算法调整Ks和Ke以最小化误差。 步骤： 1. 定义目标函数：该函数以Ks和Ke为输入，运行热传导模型，计算模拟温度与实验温度之间的RMSE。 2. 使用优化算法（如scipy.optimize中的minimize）来寻找使目标函数最小的Ks和Ke。 注意：由于运行整个模型可能比较耗时，我们需要考虑优化算法的效率。我们可以使用梯度下降类的方法，但注意目标函数可能非凸，所以可能需要全局优化或多次局部优化。 修改点： - 将原模型封装成一个函数，输入为Ks, Ke，输出为皮肤表面（或其他已知数据点）的温度随时间变化的数组。 - 定义目标函数，计算模拟值与实验值的误差。 - 调用优化器进行优化。 由于原代码较长，我们将其重构为函数形式。同时，为了减少计算时间，可以考虑： - 减少时间步长？但要注意稳定性。 - 或者使用更高效的空间离散（减少网格数）？但精度会下降。 这里我们假设已知实验数据为皮肤表面（即最内层边界）的温度随时间变化的数据。 具体步骤： 1. 将原求解过程封装为函数 `run_model(k_e, k_s)`，返回皮肤表面的温度时间序列（即`boundary_temps['皮肤表面']`）。 2. 定义目标函数 `objective(params)`，其中`params`是一个包含两个参数的数组`[k_e, k_s]`，该函数调用`run_model`得到模拟温度，然后与实验数据比较计算RMSE。 3. 使用优化器（如L-BFGS-B）最小化目标函数。 注意：原代码中，我们使用了固定的时间步长和网格，在优化过程中这些设置保持不变。 由于优化过程可能需要多次运行模型，因此我们将原代码中的求解部分提取出来，并允许传递k_e和k_s。 重构代码： 我们将原代码中与环境参数相关的部分（k_e, k_s）改为由函数参数传入。同时，为了加速，我们可以考虑： - 将一些不变的预先计算放在函数外部（例如材料参数、网格划分等），但注意在优化过程中这些是不变的。 但是，由于每次运行模型都需要重新构建矩阵并求解，且时间步数很大（5400步），运行一次就需要一定时间，因此优化过程可能很慢。我们可以考虑： - 减少时间步长？不行，可能影响稳定性。我们可以使用更大的时间步长（比如dt=5秒）来加速，但要确保稳定性（需要检查CFL条件）。 - 或者使用更粗的网格（减少m数组的值）。 权衡精度和速度，我们可以在优化时使用较粗的网格和较大的时间步长，在得到优化结果后再用细网格验证。 这里为了演示，我们保持原网格和时间步长，但实际中可能需要调整。 实验数据：假设我们有一个名为`T_exp`的数组，长度为`nt+1`（从0到5400秒，每秒一个点）。 由于我们只关心皮肤表面的温度，所以`run_model`函数只需返回皮肤表面的温度数组（与时间点对应）。 开始修改： 注意：以下代码为示例，需要根据实际情况调整。</think>以下是根据已知数据反演计算Ks和Ke的完整代码实现。该方法基于优化算法调整Ks和Ke参数，使模拟温度与实测温度误差最小： ```python import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import minimize from mpl_toolkits.mplot3d import Axes3D import time # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False # 材料参数 (附件1) materials = { 'I': {'rho': 300, 'cp': 1377, 'k': 0.082, 'thickness': 0.6e-3}, 'II': {'rho': 862, 'cp': 2100, 'k': 0.37, 'thickness': 6e-3}, 'III': {'rho': 74.2, 'cp': 1726, 'k': 0.045, 'thickness': 3.6e-3}, 'IV': {'rho': 1.18, 'cp': 1005, 'k': 0.028, 'thickness': 5e-3} } # 计算总厚度 total_thickness = sum(mat['thickness'] for mat in materials.values()) # 环境参数 u_e = 75 # 外界温度 (°C) u_s = 37 # 皮肤温度 (°C) # 初始猜测值 k_e_guess = 120.0 # 左侧对流换热系数 (W/(m²·°C)) k_s_guess = 8.36 # 右侧对流换热系数 (W/(m²·°C)) # 时间参数 total_time = 5400 # 90分钟转换为秒 dt = 5 # 使用稍大的时间步长加速计算 (s) nt = int(total_time / dt) # 时间步数 # 空间离散 (每层网格数) - 使用较粗网格加速计算 m = [4, 40, 24, 40] # 减少网格数以加速反演计算 M = sum(m) # 总节点数 M_indices = [0, m[0], sum(m[:2]), sum(m[:3]), M] # 各层分界点索引 # 计算空间步长 dx = [] for i, layer in enumerate(['I', 'II', 'III', 'IV']): dx.append(materials[layer]['thickness'] / m[i]) # 计算热扩散系数 alpha = [] for layer in ['I', 'II', 'III', 'IV']: mat = materials[layer] alpha.append(mat['k'] / (mat['rho'] * mat['cp'])) # === 核心求解函数 === def solve_heat_transfer(k_e, k_s, return_all=False): """求解热传导方程，返回皮肤表面温度时间序列""" # 计算λ_j lam = [] for i in range(4): lam.append(alpha[i] * dt / (dx[i]**2)) # 计算v_j v = [] for i, layer in enumerate(['I', 'II', 'III', 'IV']): v.append(materials[layer]['k'] / dx[i]) # 边界参数 mu_e = k_e * dx[0] / materials['I']['k'] mu_s = k_s * dx[3] / materials['IV']['k'] # 初始化温度场 (初始温度37°C) u_old = np.full(M + 1, 37.0) # 存储皮肤表面温度 skin_temps = np.zeros(nt + 1) skin_temps[0] = u_old[-1] # 初始皮肤温度 # 时间推进 for t_step in range(1, nt + 1): # 初始化三对角矩阵系数 a = np.zeros(M + 1) # 下对角线 b = np.zeros(M + 1) # 主对角线 c = np.zeros(M + 1) # 上对角线 d = np.zeros(M + 1) # 右端项 # 左边界 (i=0, 外层) b[0] = 1 + mu_e c[0] = -1 d[0] = mu_e * u_e # 右边界 (i=M, 内层) a[M] = -1 b[M] = 1 + mu_s d[M] = mu_s * u_s # 内部节点 for layer_idx in range(4): start_idx = M_indices[layer_idx] end_idx = M_indices[layer_idx + 1] for i in range(start_idx, end_idx + 1): if i == 0 or i == M: continue # 边界点已处理 # 交界点处理 if i in M_indices[1:4]: # 确定交界层 if i == M_indices[1]: # I-II交界面 v_j = v[0] v_j1 = v[1] elif i == M_indices[2]: # II-III交界面 v_j = v[1] v_j1 = v[2] elif i == M_indices[3]: # III-IV交界面 v_j = v[2] v_j1 = v[3] a[i] = -v_j b[i] = v_j + v_j1 c[i] = -v_j1 d[i] = 0 else: # 内部节点 a[i] = -lam[layer_idx] b[i] = 1 + 2 * lam[layer_idx] c[i] = -lam[layer_idx] d[i] = u_old[i] # 追赶法求解三对角方程组 # 1. 前向消元 cp = np.zeros(M + 1) # 存储c' dp = np.zeros(M + 1) # 存储d' # 首行 cp[0] = c[0] / b[0] dp[0] = d[0] / b[0] # 中间行 for i in range(1, M): denom = b[i] - a[i] * cp[i-1] cp[i] = c[i] / denom dp[i] = (d[i] - a[i] * dp[i-1]) / denom # 末行 denom = b[M] - a[M] * cp[M-1] dp[M] = (d[M] - a[M] * dp[M-1]) / denom # 2. 回代求解 u_new = np.zeros(M + 1) u_new[M] = dp[M] for i in range(M-1, -1, -1): u_new[i] = dp[i] - cp[i] * u_new[i+1] # 更新温度场并存储皮肤温度 u_old = u_new.copy() skin_temps[t_step] = u_new[-1] if return_all: return skin_temps, u_old return skin_temps # === 反演优化函数 === def objective_function(params, measured_data, time_points): """目标函数：计算模拟值与测量值的均方根误差""" k_e, k_s = params try: # 运行模型获得模拟温度 simulated_data = solve_heat_transfer(k_e, k_s) # 确保时间点对齐 sim_points = np.linspace(0, total_time, len(simulated_data)) interp_sim = np.interp(time_points, sim_points, simulated_data) # 计算均方根误差 rmse = np.sqrt(np.mean((interp_sim - measured_data)**2)) return rmse except Exception as e: print(f"计算失败: {e}") return 1e10 # 返回大值表示失败 # === 主程序 === if __name__ == "__main__": # 步骤1: 加载实测数据 (假设为CSV格式) # 实测数据格式: 时间(秒), 皮肤温度(°C) measured_df = pd.read_csv('measured_temperatures.csv') time_points = measured_df['time'].values measured_temps = measured_df['temperature'].values print(f"加载实测数据: {len(measured_temps)}个数据点") # 步骤2: 定义参数边界 bounds = [(50, 200), (5, 20)] # (k_e_min, k_e_max), (k_s_min, k_s_max) # 步骤3: 运行优化算法 print("开始反演计算...") start_time = time.time() # 使用全局优化算法寻找初始值 from scipy.optimize import differential_evolution result_global = differential_evolution( objective_function, bounds, args=(measured_temps, time_points), strategy='best1bin', popsize=15, tol=0.01, mutation=(0.5, 1), recombination=0.7, seed=42, maxiter=100 ) # 使用局部优化算法进一步优化 result = minimize( objective_function, x0=result_global.x, args=(measured_temps, time_points), method='L-BFGS-B', bounds=bounds, options={'maxiter': 50, 'ftol': 1e-4} ) # 提取优化结果 k_e_opt, k_s_opt = result.x min_rmse = result.fun elapsed_time = time.time() - start_time print(f"优化完成! 耗时: {elapsed_time:.2f}秒") print(f"优化结果: k_e = {k_e_opt:.2f} W/(m²·°C), k_s = {k_s_opt:.2f} W/(m²·°C)") print(f"最小均方根误差: {min_rmse:.4f}°C") # 步骤4: 验证优化结果 print("使用优化参数运行完整模型...") # 使用优化参数运行完整模型 full_skin_temps = solve_heat_transfer(k_e_opt, k_s_opt) full_time_points = np.linspace(0, total_time, len(full_skin_temps)) # 插值到实测时间点 interp_opt = np.interp(time_points, full_time_points, full_skin_temps) # 步骤5: 可视化结果 plt.figure(figsize=(12, 8)) # 绘制实测温度 plt.plot(time_points/60, measured_temps, 'bo', markersize=6, label='实测温度', alpha=0.7) # 绘制优化结果 plt.plot(time_points/60, interp_opt, 'r-', linewidth=2.5, label=f'优化模型 (k_e={k_e_opt:.1f}, k_s={k_s_opt:.1f})') # 绘制初始猜测结果 initial_skin_temps = solve_heat_transfer(k_e_guess, k_s_guess) interp_initial = np.interp(time_points, full_time_points, initial_skin_temps) plt.plot(time_points/60, interp_initial, 'g--', linewidth=2, label=f'初始猜测 (k_e={k_e_guess}, k_s={k_s_guess})') # 添加安全阈值线 plt.axhline(y=47, color='r', linestyle='--', alpha=0.7) plt.text(5, 47.5, '安全阈值 (47°C)', color='r', fontsize=12) # 设置图表 plt.title('皮肤表面温度实测与模拟对比', fontsize=16) plt.xlabel('时间 (分钟)', fontsize=12) plt.ylabel('温度 (°C)', fontsize=12) plt.grid(True, linestyle='--', alpha=0.7) plt.legend(loc='best', fontsize=12) plt.ylim(35, max(measured_temps.max(), interp_opt.max()) + 2) # 添加误差统计 plt.figtext(0.15, 0.75, f"优化参数误差统计:\nRMSE = {min_rmse:.4f}°C\nMAE = {np.mean(np.abs(interp_opt - measured_temps)):.4f}°C", bbox=dict(facecolor='white', alpha=0.8)) plt.tight_layout() plt.savefig('inverse_optimization_result.png', dpi=300) plt.show() # 步骤6: 生成反演报告 report = f"""=== 反演计算结果报告 === 优化算法: 差分进化 + L-BFGS-B 计算耗时: {elapsed_time:.2f}秒 优化参数: 左侧对流换热系数 k_e = {k_e_opt:.2f} W/(m²·°C) 右侧对流换热系数 k_s = {k_s_opt:.2f} W/(m²·°C) 误差统计: 均方根误差 (RMSE) = {min_rmse:.4f}°C 平均绝对误差 (MAE) = {np.mean(np.abs(interp_opt - measured_temps)):.4f}°C 最大绝对误差 = {np.max(np.abs(interp_opt - measured_temps)):.4f}°C 模型验证: 皮肤最终温度: {full_skin_temps[-1]:.2f}°C 超过安全阈值时间: {np.where(interp_opt > 47)[0][0]*dt/60:.1f}分钟 """ print(report) # 保存报告 with open('inversion_report.txt', 'w') as f: f.write(report) print("反演计算完成，结果已保存!") ``` ### 反演计算原理说明 1. **问题定义**： - 已知：多层材料参数、环境温度、皮肤温度、实测皮肤表面温度数据 - 未知：左侧对流换热系数(Ke)和右侧对流换热系数(Ks) - 目标：找到Ke和Ks的最佳值，使模拟温度与实测温度误差最小 2. **优化算法**： - 采用两阶段优化策略： - 第一阶段：差分进化算法（全局搜索） - 第二阶段：L-BFGS-B算法（局部精细优化） - 目标函数：均方根误差(RMSE) 3. **计算加速策略**： - 增大时间步长(dt=5秒) - 减少空间网格数量 - 仅计算皮肤表面温度（不计算全温度场） 4. **结果验证**： - 绘制实测与模拟温度对比曲线 - 计算多种误差指标(RMSE, MAE, MaxError) - 生成详细反演报告 ### 使用说明 1. 准备实测数据文件`measured_temperatures.csv`，格式为： ``` time,temperature 0,37.0 60,37.5 120,38.2 ... 5400,44.8 ``` 2. 运行程序，优化算法将自动寻找最佳Ke和Ks值 3. 结果输出： - 优化参数值(Ke, Ks) - 误差分析图表(`inverse_optimization_result.png`) - 详细反演报告(`inversion_report.txt`) ### 注意事项 1. 实测数据应覆盖整个时间范围(0-5400秒) 2. 初始参数范围(bounds)可根据实际情况调整 3. 对于大型数据集，可进一步减少网格数或增大时间步长 4. 优化过程可能需要几分钟到几十分钟，取决于硬件性能