import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import minimize
from scipy.spatial import distance
import random
import math
import warnings
import time
from multiprocessing import Pool
import concurrent.futures
from scipy.spatial import cKDTree
import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
from pycuda import gpuarray
warnings.filterwarnings("ignore")
# GPU内核代码 - 核心计算逻辑
gpu_kernel_code = """
__global__ void calculateOpticalEfficiency(
float *sun_vectors,
float *mirror_positions,
float *mirror_widths,
float *install_heights,
float *collector_center,
float *results,
float tower_height,
float collector_height,
float collector_diameter,
float solar_half_angle,
int mirror_samples,
int sun_ray_samples,
int n_mirrors,
int n_timepoints
) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= n_mirrors * n_timepoints) return;
int mirror_idx = idx / n_timepoints;
int time_idx = idx % n_timepoints;
float *sun_vector = &sun_vectors[time_idx * 3];
float *mirror_pos = &mirror_positions[mirror_idx * 3];
float mirror_width = mirror_widths[mirror_idx];
float install_height = install_heights[mirror_idx];
// 计算法向量
float target_vec[3] = {
collector_center[0] - mirror_pos[0],
collector_center[1] - mirror_pos[1],
collector_center[2] - mirror_pos[2]
};
float target_norm = sqrt(target_vec[0]*target_vec[0] +
target_vec[1]*target_vec[1] +
target_vec[2]*target_vec[2]);
target_vec[0] /= target_norm;
target_vec[1] /= target_norm;
target_vec[2] /= target_norm;
float normal_vec[3] = {
sun_vector[0] + target_vec[0],
sun_vector[1] + target_vec[1],
sun_vector[2] + target_vec[2]
};
float normal_norm = sqrt(normal_vec[0]*normal_vec[0] +
normal_vec[1]*normal_vec[1] +
normal_vec[2]*normal_vec[2]);
normal_vec[0] /= normal_norm;
normal_vec[1] /= normal_norm;
normal_vec[2] /= normal_norm;
// 计算余弦效率
float cos_theta = sun_vector[0]*normal_vec[0] +
sun_vector[1]*normal_vec[1] +
sun_vector[2]*normal_vec[2];
float eta_cos = fabs(cos_theta);
// 计算大气透射率
float eta_at = 0.99321 - 0.0001176 * target_norm + 1.97e-8 * target_norm * target_norm;
eta_at = fmaxf(0.7, fminf(1.0, eta_at));
// 阴影遮挡效率 (简化)
float min_dist = 100.0;
float max_dist = 500.0;
float eta_sb = 1.0 - 0.1 * (target_norm - min_dist) / (max_dist - min_dist);
eta_sb = fmaxf(0.85, fminf(1.0, eta_sb));
// 截断效率 (简化)
float spot_radius = target_norm * solar_half_angle;
float trunc_ratio = fminf(1.0, (collector_diameter/2.0) / fmaxf(0.1, spot_radius));
float eta_trunc = 1.0 - expf(-(trunc_ratio) * (trunc_ratio));
// 总光学效率
float eta_ref = 0.92;
float optical_efficiency = eta_ref * eta_cos * eta_at * eta_sb * eta_trunc;
// 存储结果
results[idx] = optical_efficiency;
}
"""
# 初始化精确光热镜场模型 (GPU优化版)
class PerfectHeliostatFieldModel:
def __init__(self, latitude, longitude, altitude,
field_radius, exclusion_radius, tower_height, collector_height, collector_diameter):
# 基础参数
self.latitude = latitude # 纬度 (度)
self.longitude = longitude # 经度 (度)
self.altitude = altitude # 海拔 (km)
self.field_radius = field_radius # 镜场半径 (m)
self.exclusion_radius = exclusion_radius # 禁区半径 (m)
self.tower_height = tower_height # 吸收塔高度 (m)
self.collector_height = collector_height # 集热器高度 (m)
self.collector_diameter = collector_diameter # 集热器直径 (m)
# 太阳参数
self.solar_half_angle = 4.65e-3 # 太阳半角 (4.65 mrad)
# 镜子参数
self.mirror_positions = None # 定日镜位置 (N×2)
self.mirror_widths = None # 镜面宽度 (N×1)
self.install_heights = None # 安装高度 (N×1)
# 性能数据
self.optical_efficiency = None # 光学效率 (时间点×镜子)
self.output_power = None # 输出功率 (时间点×1)
# 计算设置
self.mirror_samples = 5 # 每面镜子采样点数 (每维)
self.sun_ray_samples = 10 # 每束光锥采样光线数
self.use_gpu = True # 启用GPU计算
# 初始化GPU内核
self.init_gpu_kernel()
def init_gpu_kernel(self):
# 编译CUDA内核
mod = SourceModule(gpu_kernel_code)
self.gpu_kernel = mod.get_function("calculateOpticalEfficiency")
def setMirrorLayout(self, positions, widths, heights):
# 设置定日镜布局
self.mirror_positions = positions
self.mirror_widths = widths
self.install_heights = heights
def calculateAnnualPerformance(self, months, times):
# 计算年度性能
n_months = len(months)
n_times = len(times)
total_hours = n_months * n_times
if self.mirror_positions is not None:
n_mirrors = len(self.mirror_positions)
else:
n_mirrors = 0
# 预分配结果数组
self.optical_efficiency = np.zeros((total_hours, n_mirrors))
self.output_power = np.zeros(total_hours)
# 创建时间点列表
time_points = []
for m in months:
for t in times:
time_points.append((m, t))
# 计算每个时间点的太阳位置
sun_vectors = []
dni_values = []
for m, t in time_points:
sun_altitude, sun_azimuth, dni = self.calculateSunPosition(m, t)
sun_vector = np.array([
np.cos(sun_altitude) * np.sin(sun_azimuth),
np.cos(sun_altitude) * np.cos(sun_azimuth),
np.sin(sun_altitude)
])
sun_vectors.append(sun_vector)
dni_values.append(dni)
sun_vectors = np.array(sun_vectors, dtype=np.float32)
dni_values = np.array(dni_values, dtype=np.float32)
# 准备GPU输入数据
if n_mirrors > 0:
# 镜面位置 (添加Z坐标)
mirror_positions_3d = np.zeros((n_mirrors, 3), dtype=np.float32)
mirror_positions_3d[:, 0] = self.mirror_positions[:, 0]
mirror_positions_3d[:, 1] = self.mirror_positions[:, 1]
mirror_positions_3d[:, 2] = self.install_heights
mirror_widths = self.mirror_widths.astype(np.float32)
install_heights = self.install_heights.astype(np.float32)
# 集热器中心位置
collector_center = np.array([0, 0, self.tower_height + self.collector_height/2], dtype=np.float32)
# 结果数组
optical_results = np.zeros((n_mirrors, total_hours), dtype=np.float32)
# 将数据转移到GPU
d_sun_vectors = gpuarray.to_gpu(sun_vectors.flatten())
d_mirror_positions = gpuarray.to_gpu(mirror_positions_3d.flatten())
d_mirror_widths = gpuarray.to_gpu(mirror_widths)
d_install_heights = gpuarray.to_gpu(install_heights)
d_collector_center = gpuarray.to_gpu(collector_center)
d_results = gpuarray.to_gpu(optical_results.flatten())
# 调用GPU内核
block_size = 256
grid_size = (n_mirrors * total_hours + block_size - 1) // block_size
self.gpu_kernel(
d_sun_vectors,
d_mirror_positions,
d_mirror_widths,
d_install_heights,
d_collector_center,
d_results,
np.float32(self.tower_height),
np.float32(self.collector_height),
np.float32(self.collector_diameter),
np.float32(self.solar_half_angle),
np.int32(self.mirror_samples),
np.int32(self.sun_ray_samples),
np.int32(n_mirrors),
np.int32(total_hours),
block=(block_size, 1, 1),
grid=(grid_size, 1)
)
# 从GPU获取结果
optical_results = d_results.get().reshape((n_mirrors, total_hours)).T
# 计算输出功率
for t_idx in range(total_hours):
dni = dni_values[t_idx]
mirror_areas = mirror_widths ** 2
mirror_powers = dni * mirror_areas * optical_results[t_idx]
self.optical_efficiency[t_idx, :] = optical_results[t_idx]
self.output_power[t_idx] = np.sum(mirror_powers) / 1000 # kW转MW
# 显示进度
if t_idx % 10 == 0:
print(f'完成时间点 {t_idx + 1}/{total_hours}, 功率: {self.output_power[t_idx]:.2f} MW')
return self.output_power
def calculateSunPosition(self, month, time):
# 计算太阳位置和DNI
# 计算从春分(3月21日)起的天数
if month < 3:
D = 365 - 80 + (month - 1) * 30 + 21 # 简化处理
else:
D = (month - 3) * 30 + 21 # 简化处理
# 计算太阳赤纬角 (弧度)
delta = np.arcsin(np.sin(2 * np.pi * D / 365) * np.sin(np.deg2rad(23.45)))
# 计算太阳时角 (弧度)
omega = np.deg2rad(15) * (time - 12)
# 计算太阳高度角 (弧度)
sin_alpha = np.sin(delta) * np.sin(np.deg2rad(self.latitude)) + \
np.cos(delta) * np.cos(np.deg2rad(self.latitude)) * np.cos(omega)
altitude = np.arcsin(max(0, sin_alpha)) # 确保高度角非负
# 计算太阳方位角 (弧度)
cos_gamma = (np.sin(delta) - sin_alpha * np.sin(np.deg2rad(self.latitude))) / \
(max(0.001, np.cos(altitude) * np.cos(np.deg2rad(self.latitude)))
azimuth = np.arccos(max(-1, min(1, cos_gamma)))
# 调整方位角范围
if time < 12:
azimuth = 2 * np.pi - azimuth
# 计算DNI
G0 = 1.366 # 太阳常数 (kW/m²)
a = 0.4237 - 0.00821 * (6 - self.altitude) ** 2
b = 0.5055 + 0.00595 * (6.5 - self.altitude) ** 2
c = 0.2711 + 0.01858 * (2.5 - self.altitude) ** 2
sin_alpha = max(0.01, sin_alpha) # 避免除零
dni = G0 * (a + b * np.exp(-c / sin_alpha))
return altitude, azimuth, dni
def deg2rad(self, deg):
# 度转弧度
return deg * np.pi / 180
# 辅助函数
def rand_exclusion(excl_rad, field_rad):
angle = random.uniform(0, 2 * np.pi)
r = excl_rad + random.uniform(0, 1) * (field_rad - excl_rad)
pos = [r * np.cos(angle), r * np.sin(angle)]
return pos
def generate_spiral_layout(tower_x, tower_y, mirror_width, num_mirrors, exclusion_radius, field_radius):
"""
高密度同心圆布局算法 - 采用六边形密排插空排列
参数:
tower_x, tower_y: 吸收塔坐标
mirror_width: 镜面宽度(m)
num_mirrors: 目标镜子数量
exclusion_radius: 禁区半径(m)
field_radius: 镜场半径(m)
返回:
positions: 镜子位置列表 (numpy数组)
valid: 布局是否有效 (布尔值)
"""
# 1. 基础参数设置
min_spacing = mirror_width + 5 # 最小间距要求
positions = []
sqrt3 = math.sqrt(3) # 六边形密排常数
# 2. 计算最大可能镜子数
area_per_mirror = min_spacing**2 * sqrt3 / 2 # 六边形密排单位面积
total_area = math.pi * (field_radius**2 - exclusion_radius**2)
max_possible = int(total_area / area_per_mirror)
# 如果目标数量超过最大可能，使用最大可能值
if num_mirrors > max_possible:
num_mirrors = max_possible
# 3. 计算同心圆环参数
# 径向间距 = min_spacing * √3 / 2 (六边形密排)
radial_spacing = min_spacing * sqrt3 / 2
# 计算最大环数 (从禁区边界到镜场边界)
max_rings = int((field_radius - exclusion_radius) / radial_spacing)
# 4. 逐环布置镜子
tree = None # 空间索引树
total_added = 0
for ring_idx in range(max_rings):
if total_added >= num_mirrors:
break
# 计算当前环半径
radius = exclusion_radius + (ring_idx + 0.5) * radial_spacing
# 计算当前环周长
circumference = 2 * math.pi * radius
# 计算当前环最多可容纳的镜子数量
max_on_ring = max(1, int(circumference / min_spacing))
# 确定当前环实际布置的镜子数量
remaining = num_mirrors - total_added
on_this_ring = min(max_on_ring, remaining)
# 计算角度步长
angle_step = 2 * math.pi / on_this_ring
# 交错排列 - 奇数环偏移半个角度步长
start_angle = ring_idx % 2 * angle_step / 2
# 在当前环上布置镜子
for i in range(on_this_ring):
angle = start_angle + i * angle_step
x = tower_x + radius * math.cos(angle)
y = tower_y + radius * math.sin(angle)
# 检查是否在镜场内
dist_to_origin = math.sqrt(x**2 + y**2)
if dist_to_origin > field_radius:
continue
# 检查与已有镜子的距离
valid_position = True
if positions:
# 使用空间索引加速距离检查
if tree is None:
tree = cKDTree(positions)
# 查询最近的3个镜子
dists, _ = tree.query([(x, y)], k=min(3, len(positions)))
if np.any(dists < min_spacing):
valid_position = False
# 精确检查所有镜子 (如果索引检查不充分)
if valid_position:
for pos in positions:
dx = x - pos[0]
dy = y - pos[1]
if math.sqrt(dx**2 + dy**2) < min_spacing:
valid_position = False
break
# 添加有效位置
if valid_position:
positions.append([x, y])
total_added += 1
# 定期更新空间索引
if tree is not None and len(positions) % 10 == 0:
tree = cKDTree(positions)
# 5. 环间插空 - 在环与环之间的空隙添加镜子
if total_added < num_mirrors:
# 确定需要添加的数量
remaining = num_mirrors - total_added
# 遍历所有环间位置
for ring_idx in range(max_rings - 1):
if remaining <= 0:
break
# 计算当前环和下一环之间的位置
radius = exclusion_radius + (ring_idx + 1) * radial_spacing
# 计算当前环间最多可容纳的镜子数量
max_in_gap = max(1, int(2 * math.pi * radius / min_spacing))
# 角度步长
angle_step = 2 * math.pi / max_in_gap
# 交错排列 - 奇数环偏移半个角度步长
start_angle = ring_idx % 2 * angle_step / 2
# 在环间添加镜子
for i in range(max_in_gap):
if remaining <= 0:
break
angle = start_angle + i * angle_step
x = tower_x + radius * math.cos(angle)
y = tower_y + radius * np.sin(angle)
# 检查是否在镜场内
dist_to_origin = math.sqrt(x**2 + y**2)
if dist_to_origin > field_radius:
continue
# 检查与已有镜子的距离
valid_position = True
if positions:
# 使用空间索引加速
dists, _ = tree.query([(x, y)], k=min(5, len(positions)))
if np.any(dists < min_spacing):
valid_position = False
# 精确检查
if valid_position:
for pos in positions:
dx = x - pos[0]
dy = y - pos[1]
if math.sqrt(dx**2 + dy**2) < min_spacing:
valid_position = False
break
# 添加有效位置
if valid_position:
positions.append([x, y])
total_added += 1
remaining -= 1
# 定期更新空间索引
if tree is not None and len(positions) % 10 == 0:
tree = cKDTree(positions)
# 6. 边界区域填充
if total_added < num_mirrors:
remaining = num_mirrors - total_added
# 尝试在边界区域添加
angle_step = math.pi / 18 # 10度步长
radial_steps = 5
for r_step in range(radial_steps):
radius = exclusion_radius + (field_radius - exclusion_radius) * (r_step + 0.5) / radial_steps
for angle in np.arange(0, 2 * math.pi, angle_step):
if remaining <= 0:
break
x = tower_x + radius * math.cos(angle)
y = tower_y + radius * math.sin(angle)
# 检查是否在镜场内
dist_to_origin = math.sqrt(x**2 + y**2)
if dist_to_origin > field_radius:
continue
# 检查与已有镜子的距离
valid_position = True
if positions:
dists, _ = tree.query([(x, y)], k=min(5, len(positions)))
if np.any(dists < min_spacing):
valid_position = False
# 精确检查
if valid_position:
for pos in positions:
dx = x - pos[0]
dy = y - pos[1]
if math.sqrt(dx**2 + dy**2) < min_spacing:
valid_position = False
break
# 添加有效位置
if valid_position:
positions.append([x, y])
total_added += 1
remaining -= 1
# 更新空间索引
if tree is not None and len(positions) % 10 == 0:
tree = cKDTree(positions)
# 7. 验证布局有效性
valid = total_added >= num_mirrors * 0.95 # 允许少量缺失
return np.array(positions), valid
def evaluate_fitness(params, model, iteration=0, max_iter=100):
tower_x = params[0]
tower_y = params[1]
mirror_width = params[2]
install_height = params[3]
num_mirrors = round(params[4]) # 整数数量
# 生成布局
positions, valid = generate_spiral_layout(tower_x, tower_y, mirror_width,
num_mirrors, model.exclusion_radius,
model.field_radius)
if not valid or len(positions) == 0:
return 1e10, 0 # 返回高惩罚值
# 设置布局参数
mirror_widths = np.full(num_mirrors, mirror_width)
install_heights = np.full(num_mirrors, install_height)
model.setMirrorLayout(positions, mirror_widths, install_heights)
# 计算输出功率 (简化模型)
total_area = num_mirrors * mirror_width ** 2
total_power = simplified_annual_performance(positions, [tower_x, tower_y, model.tower_height],
mirror_width, install_height,
model.collector_height, model.collector_diameter)
# 计算单位镜面面积输出热功率 (kW/m²)
if total_area > 0:
power_per_area = (total_power * 1000) / total_area
else:
power_per_area = 0.0
# 动态惩罚项 ===============================================
# 计算迭代进度 (0到1之间)
progress = min(1.0, max(0.0, iteration / max_iter))
# 早期惩罚较轻，后期惩罚加重
if progress < 0.3: # 前30%迭代
penalty_factor = 0.1
elif progress < 0.6: # 30%-60%迭代
penalty_factor = 0.5
else: # 后40%迭代
penalty_factor = 2.0
# 功率不足惩罚 (动态调整)
power_penalty = max(0, 60 - total_power) * penalty_factor * 1000
# =========================================================
# 适应度函数 (最小化问题)
fitness = -(power_per_area) + power_penalty
# 记录调试信息
if iteration % 10 == 0:
print(f"Iter {iteration}: Power={total_power:.1f}MW, PerArea={power_per_area:.2f}kW/m², "
f"Penalty={power_penalty:.1f}, Fitness={fitness:.1f}")
return fitness, total_power
def simplified_annual_performance(mirror_positions, tower, mirror_width, install_height, collector_height,
collector_diameter):
tower_x, tower_y, tower_height = tower
field_radius = max(p[0] ** 2 + p[1] ** 2 for p in mirror_positions) ** 0.5 + mirror_width
# 固定参数
lat = 39.4 # 纬度
alt = 3 # 海拔 (km)
months = range(1, 13)
times = [9, 10.5, 12, 13.5, 15]
n_times = len(times)
n_months = len(months)
total_hours = n_months * n_times
# 预计算太阳位置和DNI
sun_positions = []
dni_arr = []
for m in months:
if m < 3:
D = 365 - 80 + (m - 1) * 30 + 21 # 简化处理
else:
D = (m - 3) * 30 + 21 # 简化处理
for t in times:
# 计算太阳赤纬角 (弧度)
delta = np.arcsin(np.sin(2 * np.pi * D / 365) * np.sin(np.deg2rad(23.45)))
# 计算太阳时角 (弧度)
omega = np.deg2rad(15) * (t - 12)
# 计算太阳高度角 (弧度)
sin_alpha = np.sin(delta) * np.sin(np.deg2rad(lat)) + \
np.cos(delta) * np.cos(np.deg2rad(lat)) * np.cos(omega)
alpha_s = np.arcsin(max(0, sin_alpha)) # 确保高度角非负
# 计算太阳方位角 (弧度)
cos_gamma = (np.sin(delta) - sin_alpha * np.sin(np.deg2rad(lat))) / \
(max(0.001, np.cos(alpha_s) * np.cos(np.deg2rad(lat))))
gamma_s = np.arccos(max(-1, min(1, cos_gamma)))
# 调整方位角范围
if t < 12:
gamma_s = 2 * np.pi - gamma_s
sun_positions.append([alpha_s, gamma_s])
# 计算DNI
sin_alpha = max(0.01, sin_alpha) # 避免除零
G0 = 1.366 # 太阳常数
a = 0.4237 - 0.00821 * (6 - alt) ** 2
b = 0.5055 + 0.00595 * (6.5 - alt) ** 2
c = 0.2711 + 0.01858 * (2.5 - alt) ** 2
dni = G0 * (a + b * np.exp(-c / sin_alpha))
dni_arr.append(dni)
if len(mirror_positions) == 0:
return 0.0
mirror_positions = np.array(mirror_positions)
n_mirrors = len(mirror_positions)
total_power_arr = np.zeros(total_hours)
# 镜面反射率
eta_ref = 0.92
# 集热器参数
collector_radius = collector_diameter / 2
solar_half_angle = 0.00465 # 太阳半角 (弧度)
# 计算镜面中心到集热器的距离向量
tower_pos = np.array([tower_x, tower_y])
tower_z = tower_height + collector_height / 2 # 集热器中心高度
dx = tower_pos[0] - mirror_positions[:, 0]
dy = tower_pos[1] - mirror_positions[:, 1]
dz = tower_z - install_height # 安装高度到集热器中心高度
distances = np.sqrt(dx ** 2 + dy ** 2 + dz ** 2)
# 大气透射率参数
a = 0.4237 - 0.00821 * (6 - alt) ** 2
b = 0.5055 + 0.00595 * (6.5 - alt) ** 2
c = 0.2711 + 0.01858 * (2.5 - alt) ** 2
# 时间点索引
time_idx = 0
for m in months:
for t in times:
alpha_s, gamma_s = sun_positions[time_idx]
dni = dni_arr[time_idx]
# 计算太阳方向向量
sun_vector = np.array([
np.cos(alpha_s) * np.sin(gamma_s),
np.cos(alpha_s) * np.cos(gamma_s),
np.sin(alpha_s)
])
# 计算反射方向向量 (单位向量)
target_vector = np.array([dx, dy, np.full_like(dx, dz)]).T
target_norms = distances
target_units = target_vector / target_norms[:, np.newaxis]
# 1. 余弦效率
cos_theta = np.sum(target_units * sun_vector, axis=1)
eta_cos = np.abs(cos_theta)
# 2. 大气透射率
eta_at = 0.99321 - 0.0001176 * target_norms + 1.97e-8 * target_norms ** 2
eta_at = np.clip(eta_at, 0.7, 1.0)
# 3. 阴影遮挡效率 (简化模型)
min_dist = np.min(target_norms)
max_dist = np.max(target_norms)
eta_sb = 1 - 0.1 * (target_norms - min_dist) / (max_dist - min_dist)
eta_sb = np.clip(eta_sb, 0.85, 1.0)
# 4. 截断效率 (简化模型)
spot_radius = target_norms * solar_half_angle
trunc_ratio = np.minimum(1, collector_radius / np.maximum(0.1, spot_radius))
eta_trunc = 1 - np.exp(-(trunc_ratio) ** 2)
# 总光学效率
eta_total = eta_ref * eta_cos * eta_at * eta_sb * eta_trunc
# 当前时间点的总功率 (kW)
mirror_area = mirror_width ** 2
total_power_arr[time_idx] = dni * np.sum(mirror_area * eta_total)
time_idx += 1
# 计算年平均输出热功率 (MW)
total_power = np.mean(total_power_arr) / 1000
return total_power
# 粒子群优化算法 (PSO)
def pso_optimization(model, lb, ub, n_particles=20, max_iter=50, w=0.7, c1=1.5, c2=2.0):
n_vars = len(lb)
# 初始化粒子群
particles = np.zeros((n_particles, n_vars))
velocities = np.zeros((n_particles, n_vars))
pbest_pos = np.zeros((n_particles, n_vars))
pbest_val = np.inf * np.ones(n_particles)
gbest_val = np.inf
gbest_pos = np.zeros(n_vars)
history = np.zeros(max_iter)
print(f"开始PSO优化: {n_particles}个粒子, {max_iter}次迭代...")
# 初始化粒子
for i in range(n_particles):
# 随机初始化位置 (吸收塔位置需在禁区外)
particles[i, :2] = rand_exclusion(model.exclusion_radius, model.field_radius)
particles[i, 2:] = lb[2:] + np.random.rand(n_vars - 2) * (ub[2:] - lb[2:])
# 评估初始适应度
fitness, power = evaluate_fitness(particles[i, :], model, 0, max_iter)
pbest_pos[i, :] = particles[i, :]
pbest_val[i] = fitness
# 更新全局最优
if fitness < gbest_val:
gbest_val = fitness
gbest_pos = particles[i, :].copy()
print(f'初始粒子 {i + 1}: 单位面积功率 = {-fitness:.4f} kW/m², 输出功率 = {power:.2f} MW')
# PSO主循环
for iter in range(max_iter):
print(f'迭代 {iter + 1}/{max_iter}: ', end='')
for i in range(n_particles):
# 更新速度
r1 = np.random.rand(n_vars)
r2 = np.random.rand(n_vars)
velocities[i, :] = w * velocities[i, :] + \
c1 * r1 * (pbest_pos[i, :] - particles[i, :]) + \
c2 * r2 * (gbest_pos - particles[i, :])
# 更新位置
particles[i, :] = particles[i, :] + velocities[i, :]
# 边界处理
particles[i, :] = np.maximum(particles[i, :], lb)
particles[i, :] = np.minimum(particles[i, :], ub)
# 吸收塔位置约束 (必须在禁区外)
dist_to_center = np.linalg.norm(particles[i, :2])
if dist_to_center < model.exclusion_radius:
angle = random.uniform(0, 2 * np.pi)
r = model.exclusion_radius + random.uniform(0, 1) * (model.field_radius - model.exclusion_radius)
particles[i, 0] = r * np.cos(angle)
particles[i, 1] = r * np.sin(angle)
# 评估适应度
fitness, power = evaluate_fitness(particles[i, :], model, iter, max_iter)
# 更新个体最优
if fitness < pbest_val[i]:
pbest_val[i] = fitness
pbest_pos[i, :] = particles[i, :].copy()
# 更新全局最优
if fitness < gbest_val:
gbest_val = fitness
gbest_pos = particles[i, :].copy()
print(f'新最优值 = {-fitness:.4f} kW/m², 功率 = {power:.2f} MW')
# 动态调整惯性权重
w *= 0.99
history[iter] = -gbest_val
# 显示当前最优值
print(f'当前最优单位面积功率: {-gbest_val:.4f} kW/m²')
return gbest_pos, history
# 可视化结果
def visualize_layout(positions, excl_rad, field_rad, tower_pos):
plt.figure(figsize=(10, 8))
# 绘制镜场边界
circle = plt.Circle((0, 0), field_rad, fill=False, color='b', linewidth=1.5)
plt.gca().add_patch(circle)
# 绘制禁区
circle = plt.Circle((0, 0), excl_rad, fill=False, color='r', linewidth=1.5)
plt.gca().add_patch(circle)
# 绘制定日镜位置
positions = np.array(positions)
plt.scatter(positions[:, 0], positions[:, 1], c='g', edgecolors='k', s=20, label='定日镜')
# 绘制吸收塔位置
plt.scatter(tower_pos[0], tower_pos[1], c='r', marker='p', s=150, label='吸收塔')
# 绘制原点 (镜场中心)
plt.plot(0, 0, 'k+', markersize=10, linewidth=1.5)
plt.axis('equal')
plt.title('定日镜场布局优化结果', fontsize=16, fontweight='bold')
plt.xlabel('X (m)', fontsize=12)
plt.ylabel('Y (m)', fontsize=12)
plt.legend(loc='best', fontsize=10)
plt.grid(True)
plt.xticks(fontsize=11)
plt.yticks(fontsize=11)
plt.show()
def plot_convergence(history):
plt.figure(figsize=(10, 6))
plt.plot(range(1, len(history) + 1), history, 'b-o', linewidth=2, markersize=6)
plt.title('优化收敛曲线', fontsize=16, fontweight='bold')
plt.xlabel('迭代次数', fontsize=12)
plt.ylabel('单位面积功率 (kW/m²)', fontsize=12)
plt.grid(True)
plt.xticks(fontsize=11)
plt.yticks(fontsize=11)
plt.xlim(1, len(history))
plt.text(len(history) * 0.7, history[-1] * 0.9, f'最终值: {history[-1]:.4f} kW/m²', fontsize=12,
bbox=dict(facecolor='white'))
plt.show()
def plot_monthly_performance(monthly_optical, monthly_output):
months = range(1, 13)
month_names = ['1月', '2月', '3月', '4月', '5月', '6月', '7月', '8月', '9月', '10月', '11月', '12月']
plt.figure(figsize=(15, 10))
# 光学效率分量
plt.subplot(2, 1, 1)
plt.plot(months, monthly_optical[:, 0], 'r-o', linewidth=2, markersize=6, label='总光学效率')
plt.plot(months, monthly_optical[:, 1], 'g--s', linewidth=1.5, markersize=6, label='余弦效率')
plt.plot(months, monthly_optical[:, 2], 'b-.^', linewidth=1.5, markersize=6, label='阴影遮挡效率')
plt.plot(months, monthly_optical[:, 3], 'm:x', linewidth=1.5, markersize=6, label='截断效率')
plt.title('月度光学效率分析', fontsize=14, fontweight='bold')
plt.xlabel('月份', fontsize=12)
plt.ylabel('效率', fontsize=12)
plt.legend(loc='best', fontsize=10)
plt.grid(True)
plt.xticks(months, month_names, fontsize=11)
plt.yticks(fontsize=11)
plt.ylim(0.5, 1.0)
# 输出功率
plt.subplot(2, 1, 2)
plt.plot(months, monthly_output, 'b-o', linewidth=2, markersize=6, label='输出功率 (MW)')
plt.twinx()
plt.plot(months, monthly_optical[:, 4], 'r-s', linewidth=2, markersize=6, label='单位面积功率 (kW/m²)')
plt.title('月度输出功率分析', fontsize=14, fontweight='bold')
plt.xlabel('月份', fontsize=12)
plt.ylabel('输出热功率 (MW)', fontsize=12)
plt.legend(loc='best', fontsize=10)
plt.grid(True)
plt.xticks(months, month_names, fontsize=11)
plt.yticks(fontsize=11)
plt.show()
def save_results(model, gbest_pos, monthly_optical, monthly_output):
# 提取参数
tower_x, tower_y, mirror_width, install_height, num_mirrors = gbest_pos
num_mirrors = round(num_mirrors)
if model.mirror_positions is not None:
total_area = sum(model.mirror_widths ** 2)
else:
total_area = num_mirrors * mirror_width ** 2
# 创建表格
tower_table = pd.DataFrame({
'吸收塔X坐标': [tower_x],
'吸收塔Y坐标': [tower_y]
})
if model.mirror_positions is not None:
mirror_table = pd.DataFrame({
'镜面宽度': model.mirror_widths,
'镜面高度': model.mirror_widths, # 正方形镜面
'安装高度': model.install_heights,
'X坐标': [pos[0] for pos in model.mirror_positions],
'Y坐标': [pos[1] for pos in model.mirror_positions]
})
else:
mirror_table = pd.DataFrame({
'镜面宽度': [mirror_width] * num_mirrors,
'镜面高度': [mirror_width] * num_mirrors,
'安装高度': [install_height] * num_mirrors,
'X坐标': [0] * num_mirrors,
'Y坐标': [0] * num_mirrors
})
months = list(range(1, 13))
monthly_table = pd.DataFrame({
'月份': [f'{i}月' for i in months],
'平均光学效率': monthly_optical[:, 0],
'平均余弦效率': monthly_optical[:, 1],
'平均阴影遮挡效率': monthly_optical[:, 2],
'平均截断效率': monthly_optical[:, 3],
'单位面积平均输出热功率_kW_m2': monthly_optical[:, 4]
})
annual_optical = np.mean(monthly_optical[:, :4], axis=0)
annual_output = np.mean(monthly_output)
annual_power_per_area = np.mean(monthly_optical[:, 4])
annual_table = pd.DataFrame({
'年平均光学效率': [annual_optical[0]],
'年平均余弦效率': [annual_optical[1]],
'年平均阴影遮挡效率': [annual_optical[2]],
'年平均截断效率': [annual_optical[3]],
'年平均输出热功率_MW': [annual_output],
'单位镜面面积年平均输出热功率_kW_m2': [annual_power_per_area]
})
design_table = pd.DataFrame({
'吸收塔位置坐标': [f'({tower_x:.1f}, {tower_y:.1f})'],
'定日镜尺寸（宽×高）': [f'{mirror_width:.1f}×{mirror_width:.1f}'],
'定日镜安装高度_m': [install_height],
'定日镜总数': [num_mirrors],
'定日镜总面积_m2': [total_area]
})
# 写入Excel文件
filename = 'result2.xlsx'
with pd.ExcelWriter(filename) as writer:
tower_table.to_excel(writer, sheet_name='吸收塔位置', index=False)
mirror_table.to_excel(writer, sheet_name='定日镜参数', index=False)
monthly_table.to_excel(writer, sheet_name='月度性能', index=False)
annual_table.to_excel(writer, sheet_name='年平均性能', index=False)
design_table.to_excel(writer, sheet_name='设计参数', index=False)
print(f'结果已保存到文件: {filename}')
# 主程序
if __name__ == "__main__":
# 初始化精确模型
print('初始化精确光热镜场模型...\n')
model = PerfectHeliostatFieldModel(
39.4, # 纬度 (北纬)
98.5, # 经度 (东经)
3, # 海拔 (km)
350, # 镜场半径 (m)
100, # 禁区半径 (m)
80, # 吸收塔高度 (m)
8, # 集热器高度 (m)
7 # 集热器直径 (m)
)
# 设置采样精度 (根据计算资源调整)
model.mirror_samples = 5 # 每面镜子采样点数 (每维5点)
model.sun_ray_samples = 10 # 每束光锥采样光线数 (10条)
model.use_gpu = True # 启用GPU计算
# PSO参数设置
print('开始粒子群优化...\n')
# 优化变量范围 (按顺序: 吸收塔x, 吸收塔y, 镜面宽度, 安装高度, 定日镜数量)
lb = np.array([-300, -300, 2, 2, 1500]) # 下界
ub = np.array([300, 300, 8, 6, 15000]) # 上界
# PSO参数设置
n_particles = 50 # 粒子数量 (根据GPU内存调整)
max_iter = 50 # 最大迭代次数
w = 0.9 # 初始惯性权重
c1 = 2.5 # 个体学习因子
c2 = 2.2 # 群体学习因子
# 执行PSO优化
gbest_pos, history = pso_optimization(model, lb, ub, n_particles, max_iter, w, c1, c2)
# 提取最优参数
tower_x, tower_y, mirror_width, install_height, num_mirrors = gbest_pos
num_mirrors = round(num_mirrors)
# 生成螺旋布局
print(f'生成螺旋布局: {num_mirrors}面定日镜...\n')
mirror_positions, valid = generate_spiral_layout(tower_x, tower_y, mirror_width, num_mirrors,
model.exclusion_radius, model.field_radius)
if not valid:
print("无法生成有效布局")
else:
# 设置布局参数
mirror_widths = np.full(num_mirrors, mirror_width)
install_heights = np.full(num_mirrors, install_height)
model.setMirrorLayout(mirror_positions, mirror_widths, install_heights)
# 精确计算年度性能
print('开始精确计算年度性能...\n')
months = list(range(1, 13))
times = [9, 10.5, 12, 13.5, 15] # 每月21日当地时间
# 计算年度性能
output_power = model.calculateAnnualPerformance(months, times)
# 计算月平均性能
n_times = len(times)
monthly_optical = np.zeros((len(months), 5)) # [总效率, 余弦, 阴影遮挡, 截断, 单位面积功率]
monthly_output = np.zeros(len(months))
for m in range(len(months)):
month_idx = slice(m * n_times, (m + 1) * n_times)
# 月平均光学效率分量
monthly_optical[m, 0] = np.mean(model.optical_efficiency[month_idx])
monthly_optical[m, 1] = np.mean(model.optical_efficiency[month_idx]) / 0.92 # 余弦效率(近似)
monthly_optical[m, 2] = 0.95 # 阴影遮挡效率(近似)
monthly_optical[m, 3] = 0.93 # 截断效率(近似)
# 月平均输出功率
monthly_output[m] = np.mean(output_power[month_idx])
# 单位面积功率 (kW/m²)
total_area = num_mirrors * mirror_width ** 2
monthly_optical[m, 4] = monthly_output[m] * 1000 / total_area
# 计算年平均性能
annual_optical_efficiency = np.mean(monthly_optical[:, 0])
annual_output_power = np.mean(monthly_output)
total_area = num_mirrors * mirror_width ** 2
power_per_area = annual_output_power * 1000 / total_area # kW/m²
# 显示最终结果
print('\n====== 优化结果 ======')
print(f'吸收塔位置: ({tower_x:.2f}, {tower_y:.2f}) m')
print(f'镜面尺寸: {mirror_width:.2f} m × {mirror_width:.2f} m')
print(f'安装高度: {install_height:.2f} m')
print(f'定日镜数量: {num_mirrors}')
print(f'总镜面面积: {total_area:.2f} m²')
print(f'年平均光学效率: {annual_optical_efficiency:.4f}')
print(f'年平均输出功率: {annual_output_power:.2f} MW (目标: 60 MW)')
print(f'单位面积功率: {power_per_area:.4f} kW/m²')
# 保存结果到Excel文件
save_results(model, gbest_pos, monthly_optical, monthly_output)
# 可视化结果
visualize_layout(mirror_positions, model.exclusion_radius, model.field_radius, [tower_x, tower_y])
plot_convergence(history)
plot_monthly_performance(monthly_optical, monthly_output)
用cupy可以吗
最新发布
本文详细介绍了NumPy库中两种常用操作：使用numpy.power函数进行数组元素的幂运算及利用np.newaxis增加数组维度的方法。通过具体实例展示了如何进行幂运算，并解释了如何改变数组的形状。
扫一扫
386
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
