clc; clear;
% ========== 初始化COMSOL模型 ==========
try
model = mphopen('double-layer3.mph'); % 替换为实际模型路径
model.param.set('lambda', '1.064[um]', '工作波长');
catch ME
error('无法打开COMSOL模型文件: %s', ME.message);
end
% ========== 光纤参数定义 ==========
dco_value = 80; % μm (纤芯直径)
core_radius = dco_value / 2 * 1e-6; % 转换为米
n_high = 1.4738; % 反谐振单元折射率(As2S3)
n_low = 1.4496; % 背景材料折射率(SiO2)
n_air = 1.0; % 环境折射率
% ========== PSO参数设置 ==========
lb = [0.6, 0.3]; % [dcl2/dco, dcl1/dcl2] 下限
ub = [0.8, 0.4]; % 上限
nPop = 20; % 粒子数
maxIter = 15; % 迭代次数
w = 0.7; % 惯性权重
c1 = 1.5; % 认知系数
c2 = 1.5; % 社会系数
% ========== 初始化优化算法 ==========
particles = lb + rand(nPop, 2) .* (ub - lb);
velocities = 0.1 * rand(nPop, 2);
particles(1, :) = [0.61, 0.35]; % 种子粒子（可自定义初始值）
% ========== 结果存储初始化 ==========
nx = 50; ny = 50;
[dcl2_grid, dcl1_ratio_grid] = meshgrid(linspace(lb(1), ub(1), nx), linspace(lb(2), ub(2), ny));
fund_loss_grid = nan(nx, ny);
ho_loss_grid = nan(nx, ny);
mfa_grid = nan(nx, ny);
homsr_grid = nan(nx, ny);
convergence_MFA = zeros(maxIter, 1);
% ========== 最优解初始化 ==========
pBest = particles;
pBestScore = inf(nPop, 1);
gBest = particles(1, :);
gBestScore = inf;
best_MFA = 0;
best_neff = 0;
best_loss = 0;
best_ho_loss = NaN;
best_ho_neff = NaN;
best_HOMSR = Inf;
% ========== 主优化循环 ==========
for iter = 1:maxIter
fprintf('========== 迭代 %d/%d ==========\n', iter, maxIter);
for i = 1:nPop
% 更新结构参数
dcl2_ratio = particles(i, 1);
dcl1_ratio = particles(i, 2);
model.param.set('dco', sprintf('%.2f', dco_value));
model.param.set('k1', sprintf('%.4f', dcl2_ratio)); % dcl2/dco
model.param.set('k2', sprintf('%.4f', dcl1_ratio)); % dcl1/dcl2
try
% ===== 生成计算坐标 =====
x = linspace(-1.5*core_radius, 1.5*core_radius, 300);
y = linspace(-1.5*core_radius, 1.5*core_radius, 300);
dx = x(2) - x(1);
dy = y(2) - y(1);
area_element = dx * dy;
[X, Y] = meshgrid(x, y);
coord = [X(:)'; Y(:)']; % 2xN 坐标矩阵
num_points = size(coord, 2); % 记录坐标点数量（关键参考值）
fprintf('坐标点数量: %d\n', num_points);
% ===== 第一次求解：基模搜索（基准值1.44957） =====
eigen_solver = 'sol1';
dataset_name = 'dset1';
% 清除现有特征值设置
tags_java = model.sol(eigen_solver).feature.tags;
tags_cell = cellfun(@(x) char(x), cell(tags_java), 'UniformOutput', false);
if any(strcmp('e1', tags_cell))
model.sol(eigen_solver).feature.remove('e1');
end
% 设置基模搜索参数
e1 = model.sol(eigen_solver).feature.create('e1', 'Eigenvalue');
e1.set('shift', '1.44957'); % 基模基准值
e1.set('neigs', 20); % 求解20个模式
model.sol(eigen_solver).runAll;
% 获取基模候选数据
neff_fund_candidates = mphinterp(model, 'emw.neff', 'coord', [0;0], 'dataset', dataset_name, 'solnum', 'all');
loss_fund_candidates = mphinterp(model, 'emw.dampzdB', 'coord', [0;0], 'dataset', dataset_name, 'solnum', 'all');
mfa_fund_candidates = mphinterp(model, 'MFA', 'coord', [0;0], 'dataset', dataset_name, 'solnum', 'all') * 1e12;
% 获取电场分布并强制维度匹配
n_modes_fund = length(neff_fund_candidates);
E_field_fund = cell(1, n_modes_fund);
for k = 1:n_modes_fund
Ex = mphinterp(model, 'emw.Ex', 'coord', coord, 'dataset', dataset_name, 'solnum', k);
Ey = mphinterp(model, 'emw.Ey', 'coord', coord, 'dataset', dataset_name, 'solnum', k);
% 强制调整为与坐标点数量一致（关键修复）
if length(Ex) ~= num_points
warning('粒子 %d 模式 %d 电场Ex维度不匹配，强制调整: %d -> %d', ...
i, k, length(Ex), num_points);
% 使用线性插值强制匹配维度
Ex = interp1(linspace(1, length(Ex), length(Ex)), Ex, ...
linspace(1, length(Ex), num_points), 'linear', 0);
Ey = interp1(linspace(1, length(Ey), length(Ey)), Ey, ...
linspace(1, length(Ey), num_points), 'linear', 0);
end
% 确保是列向量
Ex = Ex(:);
Ey = Ey(:);
E_field_fund{k} = [Ex, Ey];
end
% ===== 基模识别（LP01模式特征） =====
fund_mode_idx = NaN;
for k = 1:n_modes_fund
% 基模特征：1. 轴对称 2. 纤芯能量占比高 3. 无径向节点
is_axsym = is_axially_symmetric(E_field_fund{k}, X, Y);
core_ratio = compute_core_energy_ratio(E_field_fund{k}, coord, core_radius, n_high, n_low, area_element);
radial_nodes = count_radial_nodes(E_field_fund{k}, coord, core_radius);
if is_axsym && (core_ratio > 0.7) && (radial_nodes == 0)
fund_mode_idx = k;
break;
end
end
if isnan(fund_mode_idx)
error('未在第一次求解中识别到基模！');
end
% 基模特性
fundamental_neff = real(neff_fund_candidates(fund_mode_idx));
fundamental_loss = loss_fund_candidates(fund_mode_idx);
fundamental_MFA = mfa_fund_candidates(fund_mode_idx);
fprintf('--> 基模识别: 模式 %d (neff=%.6f, 损耗=%.4f dB/m)\n', ...
fund_mode_idx, fundamental_neff, fundamental_loss);
% ===== 第二次求解：高阶模搜索（基准值1.44953） =====
% 清除现有特征值设置
if any(strcmp('e1', tags_cell))
model.sol(eigen_solver).feature.remove('e1');
end
% 设置高阶模搜索参数
e1 = model.sol(eigen_solver).feature.create('e1', 'Eigenvalue');
e1.set('shift', '1.44953'); % 高阶模基准值
e1.set('neigs', 60); % 求解更多模式
model.sol(eigen_solver).runAll;
% 获取高阶模候选数据
neff_ho_candidates = mphinterp(model, 'emw.neff', 'coord', [0;0], 'dataset', dataset_name, 'solnum', 'all');
loss_ho_candidates = mphinterp(model, 'emw.dampzdB', 'coord', [0;0], 'dataset', dataset_name, 'solnum', 'all');
% 获取电场分布并强制维度匹配
n_modes_ho = length(neff_ho_candidates);
E_field_ho = cell(1, n_modes_ho);
for k = 1:n_modes_ho
Ex = mphinterp(model, 'emw.Ex', 'coord', coord, 'dataset', dataset_name, 'solnum', k);
Ey = mphinterp(model, 'emw.Ey', 'coord', coord, 'dataset', dataset_name, 'solnum', k);
% 强制调整为与坐标点数量一致（关键修复）
if length(Ex) ~= num_points
warning('粒子 %d 高阶模式 %d 电场Ex维度不匹配，强制调整: %d -> %d', ...
i, k, length(Ex), num_points);
Ex = interp1(linspace(1, length(Ex), length(Ex)), Ex, ...
linspace(1, length(Ex), num_points), 'linear', 0);
Ey = interp1(linspace(1, length(Ey), length(Ey)), Ey, ...
linspace(1, length(Ey), num_points), 'linear', 0);
end
% 确保是列向量
Ex = Ex(:);
Ey = Ey(:);
E_field_ho{k} = [Ex, Ey];
end
% ===== 高阶模识别与最小损耗筛选 =====
ho_losses = [];
ho_indices = [];
for k = 1:n_modes_ho
% 高阶模特征：1. 非轴对称 2. 有径向节点 3. 排除基模特征
is_high_order = ~is_axially_symmetric(E_field_ho{k}, X, Y);
core_ratio = compute_core_energy_ratio(E_field_ho{k}, coord, core_radius, n_high, n_low, area_element);
radial_nodes = count_radial_nodes(E_field_ho{k}, coord, core_radius);
% 判定为高阶模的条件
if (is_high_order || radial_nodes >= 1) && (core_ratio < 0.7)
ho_indices(end+1) = k;
ho_losses(end+1) = loss_ho_candidates(k);
end
end
% 寻找最小损耗高阶模
if ~isempty(ho_losses)
[min_ho_loss, min_idx] = min(ho_losses);
min_ho_idx = ho_indices(min_idx);
best_ho_neff = real(neff_ho_candidates(min_ho_idx));
HOMSR = min_ho_loss / fundamental_loss;
fprintf('--> 高阶模识别: 模式 %d (损耗=%.4f dB/m, HOMSR=%.2f)\n', ...
min_ho_idx, min_ho_loss, HOMSR);
else
min_ho_loss = NaN;
best_ho_neff = NaN;
HOMSR = Inf;
fprintf('警告: 未识别到有效高阶模！\n');
end
% ===== 适应度计算 =====
if ~isnan(min_ho_loss) && (fundamental_loss < 0.1) && (HOMSR > 100)
fitness = -fundamental_MFA; % 最大化模场面积
else
fitness = inf;
end
% 存储绘图数据
[~, xi] = min(abs(dcl2_grid(1,:) - dcl2_ratio));
[~, yi] = min(abs(dcl1_ratio_grid(:,1) - dcl1_ratio));
fund_loss_grid(xi, yi) = fundamental_loss;
ho_loss_grid(xi, yi) = min_ho_loss;
mfa_grid(xi, yi) = fundamental_MFA;
homsr_grid(xi, yi) = HOMSR;
catch ME
warning('粒子 %d 迭代失败: %s (文件: %s, 行号: %d)', ...
i, ME.message, ME.stack(1).name, ME.stack(1).line);
fitness = inf;
continue;
end
% 更新个体最优
if fitness < pBestScore(i)
pBest(i, :) = particles(i, :);
pBestScore(i) = fitness;
end
% 更新全局最优
if fitness < gBestScore
gBest = particles(i, :);
gBestScore = fitness;
best_MFA = fundamental_MFA;
best_neff = fundamental_neff;
best_loss = fundamental_loss;
best_ho_loss = min_ho_loss;
best_ho_neff = best_ho_neff;
best_HOMSR = HOMSR;
end
end
% 记录收敛历史
convergence_MFA(iter) = best_MFA;
% 输出迭代信息
fprintf('迭代 %d 最优: MFA=%.4f μm², 基模损耗=%.4f dB/m, HOMSR=%.2f\n', ...
iter, best_MFA, best_loss, best_HOMSR);
% 更新粒子位置
r1 = rand(nPop, 1); r2 = rand(nPop, 1);
velocities = w*velocities + c1*r1.*(pBest - particles) + c2*r2.*(gBest - particles);
particles = max(min(particles + velocities, ub), lb);
end
% ========== 结果可视化 ==========
% 1. 结构可视化
try
figure;
mphplot(model, 'pg2', 'surface', 'on');
title('反谐振光纤结构');
catch
warning('结构可视化失败');
end
% 2. 基模场分布
try
figure;
mphplot(model, 'pg1', 'dataset', dataset_name, 'solnum', fund_mode_idx, 'data', 'emw.normE');
title(sprintf('基模场分布 (neff=%.4f)', best_neff));
catch
warning('基模场分布可视化失败');
end
% 3. 优化结果热力图
figure('Position', [100, 100, 1200, 800]);
subplot(2,2,1);
if any(~isnan(fund_loss_grid(:)))
contourf(dcl2_grid, dcl1_ratio_grid, fund_loss_grid', 20, 'LineColor', 'none');
colorbar; title('基模损耗 (dB/m)'); xlabel('dcl2/dco'); ylabel('dcl1/dcl2');
else
title('基模损耗 (无有效数据)');
end
subplot(2,2,2);
if any(~isnan(ho_loss_grid(:)))
contourf(dcl2_grid, dcl1_ratio_grid, ho_loss_grid', 20, 'LineColor', 'none');
colorbar; title('高阶模最小损耗 (dB/m)'); xlabel('dcl2/dco'); ylabel('dcl1/dcl2');
else
title('高阶模最小损耗 (无有效数据)');
end
subplot(2,2,3);
if any(~isnan(mfa_grid(:)))
contourf(dcl2_grid, dcl1_ratio_grid, mfa_grid', 20, 'LineColor', 'none');
colorbar; title('模场面积 (μm²)'); xlabel('dcl2/dco'); ylabel('dcl1/dcl2');
else
title('模场面积 (无有效数据)');
end
subplot(2,2,4);
if any(~isnan(homsr_grid(:)))
contourf(dcl2_grid, dcl1_ratio_grid, homsr_grid', 20, 'LineColor', 'none');
colorbar; title('高阶模抑制比'); xlabel('dcl2/dco'); ylabel('dcl1/dcl2');
else
title('高阶模抑制比 (无有效数据)');
end
% 4. 收敛曲线
if any(convergence_MFA > 0)
figure;
plot(1:maxIter, convergence_MFA, 'r-s', 'LineWidth', 1.5);
xlabel('迭代次数'); ylabel('模场面积 (μm²)');
title('MFA收敛曲线'); grid on;
else
warning('无有效收敛数据');
end
% ========== 输出最终结果 ==========
fprintf('\n====== 最优结构参数 ======\n');
fprintf('dcl2/dco 比例: %.4f\n', gBest(1));
fprintf('dcl1/dcl2 比例: %.4f\n', gBest(2));
fprintf('实际 dcl2: %.2f μm\n', dco_value * gBest(1));
fprintf('实际 dcl1: %.2f μm\n', dco_value * gBest(1) * gBest(2));
fprintf('\n====== 模式特性 ======\n');
fprintf('基模有效折射率(neff): %.6f\n', best_neff);
fprintf('基模损耗: %.4f dB/m\n', best_loss);
fprintf('模场面积(MFA): %.2f μm²\n', best_MFA);
fprintf('高阶模最小损耗: %.4f dB/m\n', best_ho_loss);
fprintf('高阶模抑制比(HOMSR): %.2f\n', best_HOMSR);
% ========== 辅助函数 ==========
function is_sym = is_axially_symmetric(E_field, X, Y)
% 判断场分布是否轴对称（旋转不变性）
E_mag = sqrt(sum(abs(E_field).^2, 2)); % 电场模值
E_mag = reshape(E_mag, size(X)); % 恢复2D场分布
% 测试多个角度的旋转对称性
angles = 0:30:150; % 测试角度
corr_scores = zeros(size(angles));
for j = 1:length(angles)
theta = deg2rad(angles(j));
X_rot = X*cos(theta) - Y*sin(theta);
Y_rot = X*sin(theta) + Y*cos(theta);
% 插值旋转后的场分布
E_rot = griddata(X(:), Y(:), E_mag(:), X_rot, Y_rot, 'cubic');
E_rot(isnan(E_rot)) = 0; % 处理NaN值
corr_scores(j) = corr2(E_mag, E_rot);
end
is_sym = all(corr_scores > 0.9); % 高相关性判定为轴对称
end
function ratio = compute_core_energy_ratio(E_field, coord, core_radius, n_core, n_clad, area_element)
% 计算纤芯能量占比，增强维度一致性保障
r = sqrt(coord(1,:).^2 + coord(2,:).^2);
n_map = n_clad * ones(size(r));
n_map(r <= core_radius) = n_core;
% 提取电场分量并确保为列向量
Ex = E_field(:,1);
Ey = E_field(:,2);
% 最终维度检查与修复（最后一道防线）
if length(Ex) ~= length(n_map)
warning('强制修复能量计算维度不匹配: Ex长度=%d, n_map长度=%d', length(Ex), length(n_map));
% 强制调整到相同长度
target_len = max(length(Ex), length(n_map));
Ex = interp1(1:length(Ex), Ex, linspace(1, length(Ex), target_len));
Ey = interp1(1:length(Ey), Ey, linspace(1, length(Ey), target_len));
n_map = interp1(1:length(n_map), n_map, linspace(1, length(n_map), target_len));
r = interp1(1:length(r), r, linspace(1, length(r), target_len));
end
energy_density = (n_map.^2) .* (abs(Ex).^2 + abs(Ey).^2);
core_energy = sum(energy_density(r <= core_radius)) * area_element;
total_energy = sum(energy_density) * area_element;
% 避免除零错误
if total_energy == 0 || isnan(total_energy)
ratio = 0;
else
ratio = core_energy / total_energy;
end
end
function count = count_radial_nodes(E_field, coord, core_radius)
% 计算径向节点数（场强为零的环数）
r = sqrt(coord(1,:).^2 + coord(2,:).^2);
E_mag = sqrt(sum(abs(E_field).^2, 2));
% 确保维度一致
if length(E_mag) ~= length(r)
warning('强制修复径向节点计算维度不匹配: E_mag长度=%d, r长度=%d', length(E_mag), length(r));
target_len = max(length(E_mag), length(r));
E_mag = interp1(1:length(E_mag), E_mag, linspace(1, length(E_mag), target_len));
r = interp1(1:length(r), r, linspace(1, length(r), target_len));
end
% 筛选纤芯内的点并按半径排序
idx = r <= core_radius & r > 0;
[r_sorted, idx_sort] = sort(r(idx));
E_sorted = E_mag(idx(idx_sort));
% 平滑处理
E_smoothed = smooth(E_sorted, 10);
% 寻找过零点（节点）
crossings = 0;
for j = 2:length(E_smoothed)
if sign(E_smoothed(j)) ~= sign(E_smoothed(j-1)) && ...
abs(E_smoothed(j)) < 0.3*max(E_smoothed) % 排除边缘效应
crossings = crossings + 1;
end
end
count = crossings;
end
========== 迭代 1/15 ==========
坐标点数量: 90000
警告: 粒子 1 迭代失败: 逻辑 AND (&&)和 OR (||)运算符的操作数必须可转换为标量逻辑值。请使用 ANY 或 ALL 函数将操作数简化为标量逻辑值。 (文件:
compute_core_energy_ratio, 行号: 28)
> 位置：demo2 (第 233 行)
坐标点数量: 90000
根据报错修改代码
本文详细解析了MATLAB中size函数的使用方法及应用案例，包括如何获取矩阵的行数、列数以及理解占位符~的作用，特别强调了自2009版之后的版本变化。
扫一扫
336
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
