Commit 52a15097 by wzc-a

增加批量测试程序,与matpower对比

parent bc8fe642
// case18 Power flow data for 18 bus distribution system
// Data from ...
// W. M. Grady, M. J. Samotyj and A. H. Noyola, "The application of
// network objective functions for actively minimizing the impact of
// voltage harmonics in power systems," IEEE Transactions on Power
// Delivery, vol. 7, no. 3, pp. 1379-1386, Jul 1992.
// https://doi.org/10.1109/61.141855
//
// Modifications (v2 - 2020-09-30):
// - BaseMVA changed to 10 MVA
// - Original non-consecutive bus numbers restored
// - Slack bus Vmin = Vmax = 1.05
// - Gen limits set to 100 (instead of 999)
// - Branch flow limits disabled (set to 0)
// ----- Power Flow Data -----
// system MVA base
baseMVA = 10;
// bus data
// bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
bus = [
[1, 1, 0, 0, 0, 0, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[2, 1, 0.2, 0.12, 0, 1.05, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[3, 1, 0.4, 0.25, 0, 0.6, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[4, 1, 1.5, 0.93, 0, 0.6, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[5, 1, 3, 2.26, 0, 1.8, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[6, 1, 0.8, 0.5, 0, 0, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[7, 1, 0.2, 0.12, 0, 0.6, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[8, 1, 1, 0.62, 0, 0, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[9, 1, 0.5, 0.31, 0, 0, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[20, 1, 1, 0.62, 0, 0.6, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[21, 1, 0.3, 0.19, 0, 1.2, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[22, 1, 0.2, 0.12, 0, 0, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[23, 1, 0.8, 0.5, 0, 0, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[24, 1, 0.5, 0.31, 0, 1.5, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[25, 1, 1, 0.62, 0, 0.9, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[26, 1, 0.2, 0.12, 0, 0, 1, 1, 0, 12.5, 1, 1.1, 0.9],
[50, 1, 0, 0, 0, 1.2, 1, 1, 0, 138, 1, 1.1, 0.9],
[51, 3, 0, 0, 0, 0, 1, 1, 0, 138, 1, 1.05, 1.05]
];
// generator data
// bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
gen = [
[51, 0, 0, 100, -100, 1.05, 100, 1, 100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
];
// branch data
// fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
branch = [
[1, 2, 0.00431, 0.01204, 0.000035, 0, 0, 0, 0, 0, 1, -360, 360],
[2, 3, 0.00601, 0.01677, 0.000049, 0, 0, 0, 0, 0, 1, -360, 360],
[3, 4, 0.00316, 0.00882, 0.000026, 0, 0, 0, 0, 0, 1, -360, 360],
[4, 5, 0.00896, 0.02502, 0.000073, 0, 0, 0, 0, 0, 1, -360, 360],
[5, 6, 0.00295, 0.00824, 0.000024, 0, 0, 0, 0, 0, 1, -360, 360],
[6, 7, 0.0172, 0.0212, 0.000046, 0, 0, 0, 0, 0, 1, -360, 360],
[7, 8, 0.0407, 0.03053, 0.000051, 0, 0, 0, 0, 0, 1, -360, 360],
[2, 9, 0.01706, 0.02209, 0.000043, 0, 0, 0, 0, 0, 1, -360, 360],
[1, 20, 0.0291, 0.03768, 0.000074, 0, 0, 0, 0, 0, 1, -360, 360],
[20, 21, 0.02222, 0.02877, 0.000056, 0, 0, 0, 0, 0, 1, -360, 360],
[21, 22, 0.04803, 0.06218, 0.000122, 0, 0, 0, 0, 0, 1, -360, 360],
[21, 23, 0.03985, 0.0516, 0.000101, 0, 0, 0, 0, 0, 1, -360, 360],
[23, 24, 0.02910, 0.03768, 0.000074, 0, 0, 0, 0, 0, 1, -360, 360],
[23, 25, 0.03727, 0.04593, 0.0001, 0, 0, 0, 0, 0, 1, -360, 360],
[25, 26, 0.01104, 0.0136, 0.000118, 0, 0, 0, 0, 0, 1, -360, 360],
[50, 1, 0.00312, 0.06753, 0, 0, 0, 0, 0, 0, 1, -360, 360],
[50, 51, 0.0005, 0.00344, 0, 0, 0, 0, 0, 0, 1, -360, 360]
];
//----- OPF Data -----
// generator cost data
// 1 startup shutdown n x1 y1 ... xn yn
// 2 startup shutdown n c(n-1) ... c0
gencost = [
[2, 0, 0, 3, 0, 20, 0]
];
\ No newline at end of file
// case5 Power flow data for modified 5 bus, 5 gen case based on PJM 5-bus system
// Created by Rui Bo in 2006, modified in 2010, 2014.
// Distributed with permission.
// ----- Power Flow Data -----
// system MVA base
baseMVA = 100;
// bus data
// bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
bus = [
[1, 2, 0, 0, 0, 0, 1, 1, 0, 230, 1, 1.1, 0.9],
[2, 1, 300, 98.61, 0, 0, 1, 1, 0, 230, 1, 1.1, 0.9],
[3, 2, 300, 98.61, 0, 0, 1, 1, 0, 230, 1, 1.1, 0.9],
[4, 3, 400, 131.47, 0, 0, 1, 1, 0, 230, 1, 1.1, 0.9],
[5, 2, 0, 0, 0, 0, 1, 1, 0, 230, 1, 1.1, 0.9]
];
// generator data
// bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
gen = [
[1, 40, 0, 30, -30, 1, 100, 1, 40, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 170, 0, 127.5, -127.5, 1, 100, 1, 170, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[3, 323.49, 0, 390, -390, 1, 100, 1, 520, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[4, 0, 0, 150, -150, 1, 100, 1, 200, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[5, 466.51, 0, 450, -450, 1, 100, 1, 600, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
];
// branch data
// fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
branch = [
[1, 2, 0.00281, 0.0281, 0.00712, 400, 400, 400, 0, 0, 1, -360, 360],
[1, 4, 0.00304, 0.0304, 0.00658, 0, 0, 0, 0, 0, 1, -360, 360],
[1, 5, 0.00064, 0.0064, 0.03126, 0, 0, 0, 0, 0, 1, -360, 360],
[2, 3, 0.00108, 0.0108, 0.01852, 0, 0, 0, 0, 0, 1, -360, 360],
[3, 4, 0.00297, 0.0297, 0.00674, 0, 0, 0, 0, 0, 1, -360, 360],
[4, 5, 0.00297, 0.0297, 0.00674, 240, 240, 240, 0, 0, 1, -360, 360]
];
//----- OPF Data -----
// generator cost data
// 1 startup shutdown n x1 y1 ... xn yn
// 2 startup shutdown n c(n-1) ... c0
gencost = [
[2, 0, 0, 2, 14, 0],
[2, 0, 0, 2, 15, 0],
[2, 0, 0, 2, 30, 0],
[2, 0, 0, 2, 40, 0],
[2, 0, 0, 2, 10, 0]
];
\ No newline at end of file
// case9 Power flow data for 9 bus, 3 generator case.
// Based on data from p. 70 of:
// Chow, J. H., editor. Time-Scale Modeling of Dynamic Networks with
// Applications to Power Systems. Springer-Verlag, 1982.
//
// which in turn appears to come from:
// R.P. Schulz, A.E. Turner and D.N. Ewart, "Long Term Power System
// Dynamics," EPRI Report 90-7-0, Palo Alto, California, 1974.
// ----- Power Flow Data -----
// system MVA base
baseMVA = 100;
// bus data
// bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
bus = [
[1, 3, 0, 0, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[2, 2, 0, 0, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[3, 2, 0, 0, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[4, 1, 0, 0, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[5, 1, 90, 30, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[6, 1, 0, 0, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[7, 1, 100, 35, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[8, 1, 0, 0, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9],
[9, 1, 125, 50, 0, 0, 1, 1, 0, 345, 1, 1.1, 0.9]
];
// generator data
// bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
gen = [
[1, 72.3, 27.03, 300, -300, 1.04, 100, 1, 250, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[2, 163, 6.54, 300, -300, 1.025, 100, 1, 300, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[3, 85, -10.95, 300, -300, 1.025, 100, 1, 270, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
];
// branch data
// fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
branch = [
[1, 4, 0, 0.0576, 0, 250, 250, 250, 0, 0, 1, -360, 360],
[4, 5, 0.017, 0.092, 0.158, 250, 250, 250, 0, 0, 1, -360, 360],
[5, 6, 0.039, 0.17, 0.358, 150, 150, 150, 0, 0, 1, -360, 360],
[3, 6, 0, 0.0586, 0, 300, 300, 300, 0, 0, 1, -360, 360],
[6, 7, 0.0119, 0.1008, 0.209, 150, 150, 150, 0, 0, 1, -360, 360],
[7, 8, 0.0085, 0.072, 0.149, 250, 250, 250, 0, 0, 1, -360, 360],
[8, 2, 0, 0.0625, 0, 250, 250, 250, 0, 0, 1, -360, 360],
[8, 9, 0.032, 0.161, 0.306, 250, 250, 250, 0, 0, 1, -360, 360],
[9, 4, 0.01, 0.085, 0.176, 250, 250, 250, 0, 0, 1, -360, 360]
];
//----- OPF Data -----
// generator cost data
// 1 startup shutdown n x1 y1 ... xn yn
// 2 startup shutdown n c(n-1) ... c0
gencost = [
[2, 1500, 0, 3, 0.11, 5, 150],
[2, 2000, 0, 3, 0.085, 1.2, 600],
[2, 3000, 0, 3, 0.1225, 1, 335]
];
\ No newline at end of file
function compare_results(case_name)
% 比较 execute_and_parse 结果与 matpower 计算结果
% case_name: 算例名称,如 'case14', 'case30', 'case57' 等
if nargin < 1
case_name = 'case14'; % 默认使用 case14
end
fprintf('==============开始为算例 %s 执行测试...==============\n', case_name);
% 更新测试文件以使用指定算例
update_test_files(case_name);
fprintf('\n==============开始执行并解析测试文件...==============\n');
% 执行并解析所有测试文件
parsed_results = execute_and_parse();
fprintf('\n==============开始与 matpower 结果比较...==============\n');
% 比较 Jacobian 矩阵
compare_jacobian(parsed_results, case_name);
% 比较 Ybus 矩阵
compare_ybus(parsed_results, case_name);
% 比较 Sbus 矩阵
compare_sbus(parsed_results, case_name);
fprintf('\n所有比较完成!\n\n');
end
function compare_jacobian(parsed_results, case_name)
% 比较 Jacobian 矩阵
fprintf('\n=== 比较 Jacobian 矩阵 ===\n');
% 获取解析的 Jacobian 矩阵
jac_field = find_field_by_name(parsed_results, 'jac');
if isempty(jac_field)
fprintf('未找到 Jacobian 矩阵结果\n');
return;
end
tensor_jac = parsed_results.(jac_field);
if isempty(tensor_jac)
fprintf('Jacobian 矩阵为空\n');
return;
end
% 计算 matpower 的 Jacobian 矩阵
matpower_jac = cal_makejac(case_name);
% 比较矩阵
compare_matrices(tensor_jac, matpower_jac, 'Jacobian');
end
function compare_ybus(parsed_results, case_name)
% 比较 Ybus 矩阵
fprintf('\n=== 比较 Ybus 矩阵 ===\n');
% 获取解析的 Ybus 矩阵
ybus_field = find_field_by_name(parsed_results, 'ybus');
if isempty(ybus_field)
fprintf('未找到 Ybus 矩阵结果\n');
return;
end
tensor_ybus = parsed_results.(ybus_field);
if isempty(tensor_ybus)
fprintf('Ybus 矩阵为空\n');
return;
end
% 计算 matpower 的 Ybus 矩阵
matpower_ybus = cal_makeybus(case_name);
% 比较矩阵
compare_matrices(tensor_ybus, matpower_ybus, 'Ybus');
end
function compare_sbus(parsed_results, case_name)
% 比较 Sbus 矩阵/向量
fprintf('\n=== 比较 Sbus 矩阵 ===\n');
% 获取解析的 Sbus 矩阵
sbus_field = find_field_by_name(parsed_results, 'sbus');
if isempty(sbus_field)
fprintf('未找到 Sbus 矩阵结果\n');
return;
end
tensor_sbus = parsed_results.(sbus_field);
if isempty(tensor_sbus)
fprintf('Sbus 矩阵为空\n');
return;
end
% 这里可以添加 matpower 计算 Sbus 的代码
% matpower_sbus = cal_makesbus(case_name); % 需要实现这个函数
fprintf('Sbus 解析结果大小: %dx%d\n', size(tensor_sbus, 1), size(tensor_sbus, 2));
fprintf('暂时没有对应的 matpower Sbus 计算函数进行比较\n');
end
function ybus = cal_makeybus(case_name)
% 计算 Ybus 矩阵
mpc = loadcase(case_name);
ybus = full(makeYbus(mpc));
end
function jac = cal_makejac(case_name)
% 计算 Jacobian 矩阵
mpc = loadcase(case_name);
jac = full(makeJac(mpc, 1));
end
function field_name = find_field_by_name(results, target_name)
% 在结果结构体中查找包含特定名称的字段
field_names = fieldnames(results);
field_name = '';
for i = 1:length(field_names)
if contains(lower(field_names{i}), lower(target_name))
field_name = field_names{i};
break;
end
end
end
function compare_matrices(matrix1, matrix2, matrix_name)
% 比较两个矩阵的详细函数
fprintf('比较 %s 矩阵:\n', matrix_name);
fprintf(' 张量结果大小: %dx%d\n', size(matrix1, 1), size(matrix1, 2));
fprintf(' Matpower结果大小: %dx%d\n', size(matrix2, 1), size(matrix2, 2));
% 检查矩阵大小是否一致
if ~isequal(size(matrix1), size(matrix2))
fprintf(' ❌ 矩阵大小不一致!\n');
return;
end
% 计算误差
error_matrix = matrix1 - matrix2;
max_error = max(max(abs(error_matrix)));
mean_error = mean(mean(abs(error_matrix)));
relative_error = max_error / max(max(abs(matrix2)));
% 显示比较结果
fprintf(' 最大绝对误差: %.2e\n', max_error);
fprintf(' 平均绝对误差: %.2e\n', mean_error);
fprintf(' 最大相对误差: %.2e\n', relative_error);
% 判断是否相等(使用容差)
tolerance = 1e-10;
if max_error < tolerance
fprintf(' ✅ 矩阵相等 (容差: %.0e)\n', tolerance);
else
fprintf(' ❌ 矩阵不相等 (容差: %.0e)\n', tolerance);
% 显示最大误差的位置
[max_row, max_col] = find(abs(error_matrix) == max_error, 1);
fprintf(' 最大误差位置: (%d, %d)\n', max_row, max_col);
fprintf(' 张量值: %.6f + %.6fi\n', real(matrix1(max_row, max_col)), imag(matrix1(max_row, max_col)));
fprintf(' Matpower值: %.6f + %.6fi\n', real(matrix2(max_row, max_col)), imag(matrix2(max_row, max_col)));
end
% 可选:显示误差矩阵的统计信息
if max_error >= tolerance
fprintf(' 误差矩阵统计:\n');
fprintf(' 实部误差范围: [%.2e, %.2e]\n', min(min(real(error_matrix))), max(max(real(error_matrix))));
fprintf(' 虚部误差范围: [%.2e, %.2e]\n', min(min(imag(error_matrix))), max(max(imag(error_matrix))));
end
end
function update_test_files(case_name)
% 更新所有测试文件以使用指定的算例
% case_name: 算例名称,如 'case14', 'case30', 'case57' 等
if nargin < 1
case_name = 'case14';
end
% 定义数据文件路径
data_file = sprintf('../data/%s.txt', case_name);
% 定义测试文件列表
test_files = {
'test_make_jac.txt'
'test_make_sbus.txt'
'test_make_ybus.txt'
'test_runpf.txt'
};
fprintf('正在更新测试文件以使用算例: %s\n', case_name);
fprintf(' 数据文件: %s\n', data_file);
% 更新每个测试文件
for i = 1:length(test_files)
update_single_test_file(test_files{i}, data_file);
end
fprintf('所有测试文件已更新完成!\n');
end
function update_single_test_file(filename, data_file)
% 更新单个测试文件的算例引用
fprintf(' 更新文件: %s\n', filename);
% 检查文件是否存在
if ~exist(filename, 'file')
fprintf(' 警告: 文件 %s 不存在,跳过更新\n', filename);
return;
end
% 读取现有文件内容
fid = fopen(filename, 'r', 'n', 'UTF-8');
if fid == -1
error('无法读取文件: %s', filename);
end
content = fread(fid, '*char')';
fclose(fid);
lines = strsplit(content, '\n', 'CollapseDelimiters', false);
% 判断是否需要更新第一行
need_update = true;
if ~isempty(lines) && contains(lines{1}, '#include')
% 提取当前include的文件路径
current_include = strtrim(extractAfter(lines{1}, '#include'));
% 判断是否是data目录下的算例文件
if contains(current_include, '../data/case') && contains(current_include, '.txt')
% 是算例文件,检查是否与目标算例相同
if strcmp(current_include, data_file)
fprintf(' 已经是目标算例,无需修改\n');
need_update = false;
else
fprintf(' 当前算例: %s,将更新为: %s\n', current_include, data_file);
end
else
% 不是算例文件,需要在第一行插入算例include
fprintf(' 第一行include不是算例文件,将插入算例include\n');
lines = [{sprintf('#include %s', data_file)}, lines];
need_update = true;
end
else
% 第一行不是include语句,插入算例include
fprintf(' 文件第一行不是include语句,将插入算例include\n');
lines = [{sprintf('#include %s', data_file)}, lines];
end
% 如果需要更新,则写回文件
if need_update
% 更新第一行的 #include 语句(如果是已有的算例include)
if ~isempty(lines) && contains(lines{1}, '#include') && contains(lines{1}, '../data/case') && ~strcmp(lines{1}, sprintf('#include %s', data_file))
lines{1} = sprintf('#include %s', data_file);
end
% 写回文件
fid = fopen(filename, 'w', 'n', 'UTF-8');
if fid == -1
error('无法写入文件: %s', filename);
end
for i = 1:length(lines)
fprintf(fid, '%s', lines{i});
if i < length(lines)
fprintf(fid, '\n');
end
end
fclose(fid);
fprintf(' 文件已更新\n');
end
end
% 使用说明:
% 1. 在 MATLAB 中切换到 rspower/examples 目录
% 2. 运行: compare_results()
% 3. 程序会自动执行测试文件、解析结果,并与 matpower 计算结果比较
% 4. 可选择调用 visualize_comparison() 进行可视化比较
function results = execute_and_parse()
% 执行命令行并解析复数矩阵结果
% 返回一个结构体,包含所有测试文件的解析结果
% 设置路径和文件
tensoreval_path = '..\..\..\eig-rc\target\release\examples\tensoreval.exe';
test_files = {
'test_make_jac.txt'
'test_make_sbus.txt'
'test_make_ybus.txt'
'test_runpf.txt'
};
% 初始化结果结构体
results = struct();
fprintf('开始执行并解析测试文件...\n');
% 逐一执行和解析
for i = 1:length(test_files)
test_file = test_files{i};
fprintf('正在处理: %s\n', test_file);
% 构建并执行命令
cmd = sprintf('%s --complex %s', tensoreval_path, test_file);
[status, output] = system(cmd);
% 解析结果
if status == 0 && ~isempty(strtrim(output))
try
parsed_matrix = parse_complex_matrix_from_output(output);
results.(sprintf('test_%d_%s', i, strrep(test_file, '.txt', ''))) = parsed_matrix;
fprintf(' 成功解析矩阵: %dx%d\n', size(parsed_matrix, 1), size(parsed_matrix, 2));
catch ME
fprintf(' 解析失败: %s\n', ME.message);
results.(sprintf('test_%d_%s', i, strrep(test_file, '.txt', ''))) = [];
end
else
fprintf(' 无输出或执行失败\n');
results.(sprintf('test_%d_%s', i, strrep(test_file, '.txt', ''))) = [];
end
end
fprintf('\n解析完成!\n');
% 显示结果概要
fieldNames = fieldnames(results);
for i = 1:length(fieldNames)
matrix = results.(fieldNames{i});
if ~isempty(matrix)
fprintf('%s: %dx%d 复数矩阵\n', fieldNames{i}, size(matrix, 1), size(matrix, 2));
else
fprintf('%s: 空矩阵\n', fieldNames{i});
end
end
end
function matrix = parse_complex_matrix_from_output(output)
% 从命令输出中解析复数矩阵,自动检测矩阵大小
% 使用正则表达式匹配所有复数对
pattern = 'Complex\s*{\s*re:\s*([-\d.eE]+),\s*im:\s*([-\d.eE]+)\s*}';
[matches] = regexp(output, pattern, 'tokens');
if isempty(matches)
error('未找到复数数据');
end
% 将所有匹配转换为复数
complexNums = zeros(1, length(matches));
for i = 1:length(matches)
re = str2double(matches{i}{1});
im = str2double(matches{i}{2});
complexNums(i) = complex(re, im);
end
% 从输出中提取矩阵维度信息
[rows, cols] = extract_matrix_dimensions(output);
if rows * cols ~= length(complexNums)
error('矩阵维度与元素数量不匹配: %dx%d != %d', rows, cols, length(complexNums));
end
% 重塑为指定维度的矩阵
% 注意:根据输出格式调整reshape的维度顺序
matrix = reshape(complexNums, [cols, rows]).';
end
function [rows, cols] = extract_matrix_dimensions(output)
% 从输出中提取矩阵维度信息
% 优先尝试匹配 shape=[rows, cols] 模式
shape_pattern = 'shape=\[(\d+),\s*(\d+)\]';
shape_matches = regexp(output, shape_pattern, 'tokens');
if ~isempty(shape_matches)
rows = str2double(shape_matches{1}{1});
cols = str2double(shape_matches{1}{2});
fprintf(' 从shape信息中提取矩阵大小: %dx%d\n', rows, cols);
return;
end
% 如果没有找到shape信息,尝试通过解析矩阵结构推断
% 查找矩阵行的模式 [Complex{...}, Complex{...}, ...]
row_pattern = '\[Complex[^\]]+\]';
row_matches = regexp(output, row_pattern, 'match');
if ~isempty(row_matches)
rows = length(row_matches);
% 从第一行计算列数
first_row = row_matches{1};
col_pattern = 'Complex\s*{\s*re:\s*[-\d.eE]+,\s*im:\s*[-\d.eE]+\s*}';
col_matches = regexp(first_row, col_pattern, 'match');
cols = length(col_matches);
fprintf(' 从矩阵结构中推断大小: %dx%d\n', rows, cols);
return;
end
% 如果无法推断,尝试假设为方阵
total_elements = length(regexp(output, 'Complex\s*{\s*re:\s*[-\d.eE]+,\s*im:\s*[-\d.eE]+\s*}', 'match'));
sqrt_total = sqrt(total_elements);
if abs(sqrt_total - round(sqrt_total)) < 1e-10
rows = round(sqrt_total);
cols = round(sqrt_total);
fprintf(' 推断为方阵: %dx%d\n', rows, cols);
return;
end
error('无法确定矩阵维度,总元素数: %d', total_elements);
end
function matrix = parse_complex_matrix_from_file(filename)
% 从文件中解析复数矩阵(兼容原有功能)
% 读取文件内容
fileContent = fileread(filename);
% 使用解析函数
matrix = parse_complex_matrix_from_output(fileContent);
end
% 使用示例函数
function demo_usage()
% 演示如何使用这个函数
fprintf('执行示例...\n');
% 执行并解析所有测试
results = execute_and_parse();
% 访问特定结果
fieldNames = fieldnames(results);
for i = 1:length(fieldNames)
matrix = results.(fieldNames{i});
if ~isempty(matrix)
fprintf('\n%s 矩阵统计:\n', fieldNames{i});
fprintf(' 大小: %dx%d\n', size(matrix, 1), size(matrix, 2));
fprintf(' 最大实部: %.6f\n', max(max(real(matrix))));
fprintf(' 最大虚部: %.6f\n', max(max(imag(matrix))));
fprintf(' 矩阵范数: %.6f\n', norm(matrix, 'fro'));
end
end
end
% 使用说明:
% 1. 在 MATLAB 中切换到 rspower/examples 目录
% 2. 运行: results = execute_and_parse()
% 3. 或者运行演示: demo_usage()
% 4. 结果存储在返回的结构体中,每个字段对应一个测试文件的矩阵
Markdown 格式
0%
您添加了 0 到此讨论。请谨慎行事。
请先完成此评论的编辑!
注册 或者 后发表评论