%% compare2_MPC.m % LPV-MPC 模态切换对比实验：并发协同算法 vs 分段串行算法 % % 修改说明： % 【修改1】初始 A8 由 105%/110% 起点改为 100% 起点； % 【修改2】加入 WFM / A_MSV / A8 的物理速率限制； % 【修改3】整理 Simulink 稳态注入变量，明确 NL/NH 与 NRICDyn 的注入原则； % 【修改4】增加 Simulink/CLM 部件级模型一致性对比： % - 并发协同算法控制序列送入 Simulink； % - 分段串行算法控制序列送入 Simulink； % - 输出 LPV vs CLM 一致性表； % - 基于 CLM 结果重新计算过渡态过程时间。 % 【修改5】增加分段串行策略开关与融合绘图： % - serial_mode_switch=0/1/2/3 对应先切换后跟踪/先跟踪后切换/两者均运行/仅并发协同； % - LPV vs CLM 绘图按开关自动包含启用的分段串行曲线； % - 推力子图按 0.5% 推力误差保持指标标注各方法用时； % - 各子图标注 LPV 与 CLM 的拟合一致性。 % - 推力子图采用色带表示各方法模态切换区间。 % % 【本次修改6】只保留 A_MSV/A8 最终几何目标，删除中间线性参考轨迹； % 【本次修改7】并发与串行推力权重统一为 80，消除权重不一致； % 【本次修改8】任务总完成时间改为 max(几何完成时间, 推力进入目标误差带时间)， % 并输出 LPV/CLM 中文指标说明表，区分总任务时间与推力自身响应时间。 % % 控制量： % u = [WFM; A_MSV; A8] % % 单位： % WFM : lb/s % A_MSV : in^2 % A8 : ratio % % 依赖： % - ../LPVtest_3sched_08MA.m % - sample(sysLPV3,...) % - mpc_build_prediction.m % - mpc_build_qp.m % - Optimization Toolbox: quadprog % - 可选：mpc_plot_compare.m % - Simulink 模型：vceEngineModelCompare_260420.slx 或按 slxName 修改 % 20260502changelog-庞鸿泰 % 改代码用于亚音速点下的MPC实现 direct_plot_only = 0; % 1=直接根据当前工作区数据作图，不运行算法；0=正常运行算法 plot_line_width_scale = 1.25; % 融合图线宽倍率，建议 1.0~1.5 if direct_plot_only clc; warning off; else clearvars -except direct_plot_only plot_line_width_scale; clc; warning off; end close all; cd(fileparts(mfilename('fullpath'))); addpath(fullfile('..')); addpath(fullfile('..','CDFS_VCE_MPC','Stage_work_VCE','vceEngine_ModeSwitch_MPC','MPC_code')); addpath(fullfile('..','CDFS_VCE_MPC','Stage_work_VCE','vceEngine_ModeSwitch_MPC','MPC_code','MPC_code_3u')); if direct_plot_only requiredPlotVars = {'res_concurrent','clm_concurrent','metric_concurrent','metric_clm_concurrent', ... 'serialCases','t_sw_start','t_sw_end','FN_target_N','FN_initial_N', ... 'T4_max_K','SM_LP_min','SM_IP_min','SM_HP_min','FAR_min_phy','far_idx','norm_u','norm_y'}; missingPlotVars = requiredPlotVars(~cellfun(@(v) evalin('base', sprintf('exist(''%s'',''var'')', v)), requiredPlotVars)); if ~isempty(missingPlotVars) error('直接作图模式缺少工作区变量：%s', strjoin(missingPlotVars, ', ')); end figCLM_merged = plot_lpv_clm_compare_merged( ... res_concurrent, clm_concurrent, metric_concurrent, metric_clm_concurrent, ... serialCases, t_sw_start, t_sw_end, ... FN_target_N, FN_initial_N, T4_max_K, SM_LP_min, SM_IP_min, SM_HP_min, ... FAR_min_phy, far_idx, norm_u, norm_y); return; end %% ── 步骤0：加载三调度量 LPV 模型 ─────────────────────────────────────── run( 'LV2_LPV_3sched_0_8MA.m'); Ts = 0.02; norm_x = Engine.DP.Value(indexSVMState)'; % nx×1 norm_u = Engine.DP.Value(indexSVMIn)'; % nu_s×1 norm_y = Engine.DP.Value(indexSVMOut)'; % ny×1 nx = nx_s; % 2 nu_mpc = nu_s; % 3: [WFM, A_MSV, A8] ny = ny_s; far_idx = find(strcmp(outNames, 'FAR'), 1); if isempty(far_idx) error('当前 LPV 输出中未找到 FAR，请先确认 LV2_LPV_3sched_0_8MA.m 已将 FAR 加入 Engine.SVM.OutName。'); end %% ── 步骤1：初始与目标工作点 ───────────────────────────────────────────── grid_XNLC_phy = Engine.PWLM.SchdOutVec{1}; grid_MSV_phy = Engine.PWLM.SchdInVec{4}; grid_A8_phy = Engine.PWLM.SchdInVec{5}; % ========================================================= % 【修改1】A8 初始点改为 100% % ========================================================= A8_INIT_FACTOR = 1; A8_TGT_FACTOR = 1; A8_base_phy = grid_A8_phy(ceil(numel(grid_A8_phy)/2)); A8_CTRL_MIN_FACTOR = 1.00; A8_CTRL_MAX_FACTOR = 1.05; A8_INIT_FACTOR = min(max(A8_INIT_FACTOR, A8_CTRL_MIN_FACTOR), A8_CTRL_MAX_FACTOR); A8_TGT_FACTOR = min(max(A8_TGT_FACTOR, A8_CTRL_MIN_FACTOR), A8_CTRL_MAX_FACTOR); [~, idx_nl_tgt] = min(abs(grid_XNLC_phy - 0.80*norm_x(1))); [~, idx_msv_tgt] = min(abs(grid_MSV_phy - 51.2)); [~, idx_a8_tgt] = min(abs(grid_A8_phy - A8_TGT_FACTOR * A8_base_phy)); [~, idx_nl_init] = min(abs(grid_XNLC_phy - 1.00*norm_x(1))); [~, idx_msv_init] = min(abs(grid_MSV_phy - 2.0)); [~, idx_a8_init] = min(abs(grid_A8_phy - A8_INIT_FACTOR * A8_base_phy)); initInfo = get_lpv_clm_init_point(Engine, [idx_nl_init, idx_msv_init, idx_a8_init], ... indexSVMState, indexSVMIn, indexSVMOut); x0_norm = x_off(:, idx_nl_init, idx_msv_init, idx_a8_init); u0_norm = u_off(:, idx_nl_init, idx_msv_init, idx_a8_init); y0_norm = y_off(:, idx_nl_init, idx_msv_init, idx_a8_init); xr_norm = x_off(:, idx_nl_tgt, idx_msv_tgt, idx_a8_tgt); ur_norm = u_off(:, idx_nl_tgt, idx_msv_tgt, idx_a8_tgt); yr_norm = y_off(:, idx_nl_tgt, idx_msv_tgt, idx_a8_tgt); % 目标点可能落在最近邻填充的缺失格点上，此时强制覆盖 A_MSV/A8 的物理目标。 ur_norm(2) = grid_MSV_phy(idx_msv_tgt) / norm_u(2); ur_norm(3) = grid_A8_phy(idx_a8_tgt) / norm_u(3); fprintf('\n================ 初始/目标工作点 ================\n'); fprintf('A8 基准值 = %.6f, 初始因子 = %.2f, 目标因子 = %.2f, 控制范围 = %.0f%%~%.0f%%\n', ... A8_base_phy, A8_INIT_FACTOR, A8_TGT_FACTOR, 100*A8_CTRL_MIN_FACTOR, 100*A8_CTRL_MAX_FACTOR); fprintf('Initial: NL=%.1f%% A_MSV=%.2f in^2 A8=%.4f WFM=%.4f lb/s\n', ... x0_norm(1)*100, ... u0_norm(2)*norm_u(2), ... u0_norm(3)*norm_u(3), ... u0_norm(1)*norm_u(1)); fprintf('Target : NL=%.1f%% A_MSV=%.2f in^2 A8=%.4f WFM_ref=%.4f lb/s\n', ... xr_norm(1)*100, ... ur_norm(2)*norm_u(2), ... ur_norm(3)*norm_u(3), ... ur_norm(1)*norm_u(1)); fprintf('[诊断] 初始点: T4=%.0f K FN=%.1f FAN_SM=%.1f%% IPC_SM=%.1f%% HPC_SM=%.1f%% FAR=%.5f\n', ... y0_norm(3)*norm_y(3), ... y0_norm(4)*norm_y(4), ... y0_norm(5)*norm_y(5), ... y0_norm(6)*norm_y(6), ... y0_norm(7)*norm_y(7), ... y0_norm(far_idx)*norm_y(far_idx)); fprintf('[诊断] 目标点: T4=%.0f K FN=%.1f FAN_SM=%.1f%% IPC_SM=%.1f%% HPC_SM=%.1f%% FAR=%.5f\n', ... yr_norm(3)*norm_y(3), ... yr_norm(4)*norm_y(4), ... yr_norm(5)*norm_y(5), ... yr_norm(6)*norm_y(6), ... yr_norm(7)*norm_y(7), ... yr_norm(far_idx)*norm_y(far_idx)); close all %% ── 步骤2：统一 MPC 参数与约束 ───────────────────────────────────────── Np = 10; %Nc = 15; Nc = 5; webNp = str2double(getenv('LPV_MPC_NP')); webNc = str2double(getenv('LPV_MPC_NC')); if isfinite(webNp) && webNp >= 1, Np = round(webNp); end if isfinite(webNc) && webNc >= 1, Nc = round(webNc); end serial_mode_switch = 2; plot_LPV_MPC_figs = true; webSerialMode = str2double(getenv('LPV_SERIAL_MODE')); if isfinite(webSerialMode) && ismember(round(webSerialMode), 0:3) serial_mode_switch = round(webSerialMode); end if strcmp(getenv('LPV_DISABLE_PLOTS'), '1'), plot_LPV_MPC_figs = false; end plot_extra_figs = 0; plot_compare_figs = plot_LPV_MPC_figs; runSimulinkCompare = true; if strcmp(getenv('SKIP_CLM_COMPARE'), '1') runSimulinkCompare = false; end run_CLM_compare = runSimulinkCompare; plot_CLM_compare_figs = runSimulinkCompare; plot_CLM_method_figs = runSimulinkCompare && plot_extra_figs ~= 0; inject_simulink_steady = true; % true: 注入 NL/NH/WFM/A_MSV/A8 与 NRICDyn 稳态初值 beautify_figures = true; saveFig = false; saveResults = true; saveDir = pwd; figFormat = 'pdf'; simTime = 20; webSimTime = str2double(getenv('LPV_SIM_TIME')); if isfinite(webSimTime) && webSimTime > 0, simTime = webSimTime; end Nsim = round(simTime / Ts); tVec = (0:Nsim-1) * Ts; t_sw_start = 2.0; t_sw_end = 5.0; settle_tol_rel = 0.005; settle_hold_time = 0.50; % 说明：本版本不再使用 t_geom_end 作为串行算法的硬分段边界。 % t_sw_end 仅保留为并发算法名义窗口、绘图区间和兼容旧变量。 serial_fuel_start = t_sw_start; % 兼容旧绘图/保存变量；事件串行算法从 t_sw_start 开始。 serial_fuel_idx = find(tVec >= serial_fuel_start, 1, 'first'); %#ok if ~ismember(serial_mode_switch, [0 1 2 3]) error('serial_mode_switch 必须为 0、1、2 或 3。'); end serialCases = struct('id',{},'tag',{},'caseName',{},'plotName',{}, ... 'geomStart',{},'geomEnd',{},'fuelStart',{},'color',{}); if serial_mode_switch == 0 || serial_mode_switch == 2 serialCases(end+1) = struct( ... 'id','serial1_three_stage', ... 'tag','serial1_three_stage', ... 'caseName','分段串行算法1-三段事件触发', ... 'plotName','Seq-I', ... 'geomStart',NaN, ... 'geomEnd',NaN, ... 'fuelStart',t_sw_start, ... 'color',[0.78 0.18 0.18]); end if serial_mode_switch == 1 || serial_mode_switch == 2 serialCases(end+1) = struct( ... 'id','serial2_two_stage', ... 'tag','serial2_two_stage', ... 'caseName','分段串行算法2-两段事件触发', ... 'plotName','Seq-II', ... 'geomStart',NaN, ... 'geomEnd',NaN, ... 'fuelStart',t_sw_start, ... 'color',[0.12 0.55 0.32]); end run_serial_cases = ~isempty(serialCases); t_sw_idx = find(tVec >= t_sw_start, 1, 'first'); % 推力目标 FN_target_N = 1610; FN_mid_N = 2474; % 分段串行算法1 a段的中间推力目标 FN_target_norm = FN_target_N / norm_y(4); FN_mid_norm = FN_mid_N / norm_y(4); FN_initial_N = y0_norm(4) * norm_y(4); fprintf('\n[参考] 初始 FN=%.1f, 中间 FN=%.1f, 最终目标 FN=%.1f\n', ... FN_initial_N, FN_mid_N, FN_target_N); fprintf('[参考] 并发算法名义模态切换窗口 %.2f~%.2f s；事件串行算法不再由 t_geom_end 固定切段。\n', ... t_sw_start, t_sw_end); fprintf('[分段串行开关] serial_mode_switch = %d\n', serial_mode_switch); for iSerial = 1:numel(serialCases) if strcmp(serialCases(iSerial).id, 'serial1_three_stage') fprintf([' %s: A段 WFM/A8 跟踪 FN=%.1f；', ... 'B段 WFM/A_MSV/A8 在 %.1fN 附近平稳模态切换；', ... 'C段 WFM/A8 跟踪最终 FN=%.1f。\n'], ... serialCases(iSerial).caseName, FN_mid_N, FN_mid_N, FN_target_N); else fprintf([' %s: A段 WFM/A8 跟踪最终 FN=%.1f；', ... 'B段 WFM/A_MSV/A8 在最终推力附近完成模态切换。\n'], ... serialCases(iSerial).caseName, FN_target_N); end end if ~run_serial_cases fprintf(' 仅运行并发协同算法；跳过分段串行 LPV/CLM 仿真与作图。\n'); end % 输出权重：并发与串行统一为 80，避免“不同权重”干扰算法公平对比。 Qy_concurrent = zeros(ny); Qy_serial = zeros(ny); Qy_concurrent(4,4) = 80; Qy_serial(4,4) = 80; % 终端代价 Pu_terminal = diag([0, 50, 300]); % 并发协同算法保持原设置：A_MSV 与 A8 均为最终几何目标 Pu_terminal_serial_A8 = diag([0, 0, 300]); % 事件串行算法：仅要求 A8 逼近最终目标 Pu_terminal_serial_mode = diag([0, 50, 300]); % 事件串行算法：A_MSV 与 A8 在后续阶段均需逼近最终目标 Px = diag([0, 0]); % 控制增量惩罚 / 绝对偏差惩罚 Ru = diag([1, 8, 5]); % Qu 设为 0：不再要求 A_MSV/A8 跟踪人为中间轨迹；几何目标只通过 Pu_terminal 体现。 Qu = diag([0, 0, 0]); % 控制量上下界 A8_min_phy_ctrl = A8_CTRL_MIN_FACTOR * A8_base_phy; A8_max_phy_ctrl = A8_CTRL_MAX_FACTOR * A8_base_phy; u_min_norm = [0.1; 2.0/norm_u(2); A8_min_phy_ctrl / norm_u(3)]; u_max_norm = [2.0; 51.2/norm_u(2); A8_max_phy_ctrl / norm_u(3)]; % ========================================================= % 【修改2】执行机构速率限制：物理速率 -> 单步物理增量 -> 归一化增量 % ========================================================= %WFM_rate_max = 0.08; % lb/s^2 WFM_rate_max = 0.07; % lb/s^2· AMSV_rate_max = 10; % in^2/s A8_rate_max = 0.1; % ratio/s webWfmRate = str2double(getenv('LPV_WFM_RATE_MAX')); if isfinite(webWfmRate) && webWfmRate > 0, WFM_rate_max = webWfmRate; end du_max_phy = [ ... WFM_rate_max * Ts; ... AMSV_rate_max * Ts; ... A8_rate_max * Ts ... ]; du_max = du_max_phy ./ norm_u(1:3); du_min = -du_max; fprintf('\n[速率限制] 控制量单步最大变化：\n'); fprintf(' WFM : %.4f lb/s per step, 等效速率 %.3f lb/s^2\n', ... du_max_phy(1), WFM_rate_max); fprintf(' A_MSV : %.4f in^2 per step, 等效速率 %.3f in^2/s\n', ... du_max_phy(2), AMSV_rate_max); fprintf(' A8 : %.5f ratio per step, 等效速率 %.3f ratio/s\n', ... du_max_phy(3), A8_rate_max); fprintf('[速率限制] 归一化 du_max = [%.6f, %.6f, %.6f]^T\n', ... du_max(1), du_max(2), du_max(3)); fprintf('[A8控制范围] A8_min = %.6f, A8_max = %.6f ratio，对应 %.0f%%~%.0f%%\n\n', ... A8_min_phy_ctrl, A8_max_phy_ctrl, 100*A8_CTRL_MIN_FACTOR, 100*A8_CTRL_MAX_FACTOR); % 输出硬约束 T4_max_K = 3700; SM_LP_min = 10; SM_IP_min = 10; SM_HP_min = 10; FAR_min_ratio_dp = 0.7; FAR_min_phy = FAR_min_ratio_dp * norm_y(far_idx); T4_max_norm = T4_max_K / norm_y(3); SM_LP_min_norm = SM_LP_min / norm_y(5); SM_IP_min_norm = SM_IP_min / norm_y(6); SM_HP_min_norm = SM_HP_min / norm_y(7); FAR_min_norm = FAR_min_phy / norm_y(far_idx); fprintf('约束: T4 <= %.0f K, FAN/CDFS/HPC SM >= %.1f%%, FAR >= %.5f (%.0f%% DP)\n', T4_max_K, SM_LP_min, FAR_min_phy, 100*FAR_min_ratio_dp); %% ── 步骤2.1：Simulink / CLM 对比配置 ────────────────────────────────── % ========================================================= % 【修改4】Simulink / CLM 部件级模型一致性对比选项 % ========================================================= archiveRoot = fileparts(fileparts(mfilename('fullpath'))); slxFullPath = fullfile(archiveRoot, 'vceEngineModelCompare_0511_speeddown.slx'); [slxDir, slxName, ~] = fileparts(slxFullPath); % 将模型所在目录加入 MATLAB path，保证 load_system(slxName) 与 sim(slxName) 能够找到模型。 slxDirList = { ... slxDir, ... fileparts(mfilename('fullpath')), ... fullfile(fileparts(mfilename('fullpath')), '..'), ... fullfile(fileparts(mfilename('fullpath')), '..', 'CDFS_VCE_MPC', 'Stage_work_VCE', 'vceEngine_ModeSwitch_MPC', 'MPC_code', 'MPC_code_3u') ... }; NL_block_path = [slxName, '/Engine/NL_integrator']; NH_block_path = [slxName, '/Engine/NH_integrator']; %% ── 步骤3：构造两类参考轨迹 ───────────────────────────────────────────── [ur_concurrent, yr_concurrent, xr_concurrent] = build_reference_trajectory( ... tVec, t_sw_start, t_sw_end, t_sw_start, ... x0_norm, xr_norm, u0_norm, ur_norm, y0_norm, FN_target_norm, nu_mpc); for iSerial = 1:numel(serialCases) [serialCases(iSerial).ur, serialCases(iSerial).yr, serialCases(iSerial).xr] = ... build_reference_trajectory( ... tVec, serialCases(iSerial).geomStart, serialCases(iSerial).geomEnd, ... serialCases(iSerial).fuelStart, ... x0_norm, xr_norm, u0_norm, ur_norm, y0_norm, FN_target_norm, nu_mpc); end if run_serial_cases ur_serial = serialCases(1).ur; yr_serial = serialCases(1).yr; xr_serial = serialCases(1).xr; else ur_serial = []; yr_serial = []; xr_serial = []; end %% ── 步骤4：分别运行两种 LPV-MPC 算法 ────────────────────────────────── baseParam = struct(); baseParam.Ts = Ts; baseParam.Nsim = Nsim; baseParam.tVec = tVec; baseParam.t_sw_idx = t_sw_idx; baseParam.Np = Np; baseParam.Nc = Nc; baseParam.nx = nx; baseParam.nu_mpc = nu_mpc; baseParam.ny = ny; baseParam.Ru = Ru; baseParam.Qu = Qu; baseParam.Pu_terminal = Pu_terminal; baseParam.Pu_terminal_serial_A8 = Pu_terminal_serial_A8; baseParam.Pu_terminal_serial_mode = Pu_terminal_serial_mode; baseParam.Pu_zero = zeros(nu_mpc); baseParam.Px = Px; % 【修正】run_lpv_mpc_case / 串行阶段设置均需要 p.Px baseParam.t_start = t_sw_start; baseParam.y0_norm = y0_norm; baseParam.x0_norm = x0_norm; baseParam.u0_norm = u0_norm(1:nu_mpc); baseParam.ur_final_norm = ur_norm(1:nu_mpc); baseParam.xr_final_norm = xr_norm; baseParam.FN_mid_N = FN_mid_N; baseParam.FN_mid_norm = FN_mid_norm; baseParam.FN_target_N = FN_target_N; baseParam.FN_target_norm = FN_target_norm; baseParam.FN_tol_rel = settle_tol_rel; baseParam.stage_hold_time = settle_hold_time; baseParam.AMSV_target_norm = ur_norm(2); baseParam.AMSV_target_phy = ur_norm(2) * norm_u(2); baseParam.AMSV_mode_tol_phy = max(0.02 * abs(ur_norm(2)*norm_u(2) - u0_norm(2)*norm_u(2)), 0.10); baseParam.u_min_norm = u_min_norm; baseParam.u_max_norm = u_max_norm; baseParam.du_max = du_max; baseParam.du_min = du_min; baseParam.T4_max_norm = T4_max_norm; baseParam.SM_LP_min_norm = SM_LP_min_norm; baseParam.SM_IP_min_norm = SM_IP_min_norm; baseParam.SM_HP_min_norm = SM_HP_min_norm; baseParam.FAR_min_norm = FAR_min_norm; baseParam.far_idx = far_idx; baseParam.gMSV = gMSV; baseParam.gA8 = gA8; baseParam.printEvery = round(1/Ts); baseParam.qp_opts = optimoptions('quadprog','Display','off','MaxIterations',500); fprintf('\n================ 开始算法1：并发协同算法 ================\n'); paramConcurrent = baseParam; paramConcurrent.Qy = Qy_concurrent; requiredParamFields = { ... 'Ts','Nsim','tVec','t_sw_idx','Np','Nc','nx','nu_mpc','ny', ... 'Qy','Px','Ru','Qu','Pu_terminal', ... 'u_min_norm','u_max_norm','du_max','du_min', ... 'T4_max_norm','SM_LP_min_norm','SM_IP_min_norm','SM_HP_min_norm','FAR_min_norm','far_idx', ... 'gMSV','gA8','printEvery','qp_opts'}; assert_required_fields(paramConcurrent, 'paramConcurrent', requiredParamFields); res_concurrent = run_lpv_mpc_case('并发协同算法', sysLPV3, ... x0_norm, u0_norm(1:nu_mpc), norm_x, norm_u, norm_y, ... ur_concurrent, yr_concurrent, xr_concurrent, paramConcurrent); paramSerial = baseParam; paramSerial.Qy = Qy_serial; assert_required_fields(paramSerial, 'paramSerial', requiredParamFields); for iSerial = 1:numel(serialCases) fprintf('\n================ 开始算法%d：%s ================\n', ... iSerial + 1, serialCases(iSerial).caseName); paramSerialCase = paramSerial; paramSerialCase.serial_id = serialCases(iSerial).id; serialCases(iSerial).res = run_lpv_mpc_event_serial_case(serialCases(iSerial), sysLPV3, ... x0_norm, u0_norm(1:nu_mpc), norm_x, norm_u, norm_y, ... ur_norm(1:nu_mpc), yr_norm, xr_norm, paramSerialCase); if isfield(serialCases(iSerial).res, 'stageTimes') serialCases(iSerial).geomStart = serialCases(iSerial).res.stageTimes.GeomStartAbsTime_s; serialCases(iSerial).geomEnd = serialCases(iSerial).res.stageTimes.GeomEndAbsTime_s; serialCases(iSerial).fuelStart = t_sw_start; end end if run_serial_cases res_serial = serialCases(1).res; else res_serial = []; end %% ── 步骤5：LPV 层面过渡态过程时间与性能指标 ─────────────────────────── metricOpt = struct(); metricOpt.FN_tol_rel = settle_tol_rel; metricOpt.hold_time = settle_hold_time; metricOpt.Ts = Ts; metricOpt.t_start = t_sw_start; metricOpt.t_geom_end = t_sw_end; metricOpt.geom_check_idx = [2 3]; % 并发协同算法：A_MSV 和 A8 均计入几何最终目标完成判据 metricOpt.FN0 = FN_initial_N; metricOpt.FNtgt = FN_target_N; metricOpt.u_tgt_phy = ur_norm(1:nu_mpc).*norm_u(1:nu_mpc); metricOpt.T4_max_K = T4_max_K; metricOpt.SM_min = [SM_LP_min, SM_IP_min, SM_HP_min]; metricOpt.far_idx = far_idx; metricOpt.FAR_min_phy = FAR_min_phy; metric_concurrent = calc_case_metrics(res_concurrent, metricOpt, t_sw_start); for iSerial = 1:numel(serialCases) metricOptSerial = metricOpt; metricOptSerial.t_geom_end = serialCases(iSerial).geomEnd; metricOptSerial.geom_check_idx = [2 3]; % 事件串行算法：A_MSV 与 A8 均计入最终几何/控制量目标完成判据 serialCases(iSerial).metric = calc_case_metrics( ... serialCases(iSerial).res, metricOptSerial, serialCases(iSerial).fuelStart); end if run_serial_cases metric_serial = serialCases(1).metric; metric_serial_all = [serialCases.metric]; metricTable = struct2table([metric_concurrent; metric_serial_all(:)]); else metric_serial = []; metric_serial_all = []; metricTable = struct2table(metric_concurrent); end % 【本次修改】工作区中文指标表：区分“总任务完成时间”和“推力自身响应时间”。 LPV_metric_table = metricTable; LPV_metric_summary_CN = build_metric_summary_cn('LPV预测模型', metric_concurrent, serialCases, 'metric'); MetricExplain_CN = build_metric_explain_cn(); assignin('base', 'LPV_metric_table', LPV_metric_table); assignin('base', 'LPV_metric_summary_CN', LPV_metric_summary_CN); assignin('base', 'MetricExplain_CN', MetricExplain_CN); fprintf('\n================ LPV 指标中文摘要：任务总完成时间与推力响应时间分开统计 ================\n'); disp(LPV_metric_summary_CN); fprintf('指标含义说明已输出到工作区变量 MetricExplain_CN。\n'); fprintf('\n[LPV metrics] LPV transition metrics are computed for internal comparison only; reported transition times below use CLM responses.\n'); show_lpv_transition_metrics = strcmp(getenv('SHOW_LPV_METRICS'), '1'); if show_lpv_transition_metrics for iSerial = 1:numel(serialCases) metricS = serialCases(iSerial).metric; if ~isnan(metricS.SerialCompleteFromSwitch_s) && ~isnan(metric_concurrent.SerialCompleteFromSwitch_s) save_time = metricS.SerialCompleteFromSwitch_s - metric_concurrent.SerialCompleteFromSwitch_s; save_rate = save_time / max(metricS.SerialCompleteFromSwitch_s, eps) * 100; fprintf('LPV结果：并发协同算法相对%s节省完整过程时间 %.3f s，占串行法 %.2f%%。\n', ... serialCases(iSerial).caseName, save_time, save_rate); else fprintf('警告：LPV 中并发协同算法或%s未在仿真时长内完成推力/几何串行过程，无法计算节省比例。\n', ... serialCases(iSerial).caseName); end end %% ── 步骤6：LPV 层面两种算法对比图 ───────────────────────────────────── end if run_serial_cases res_serial_plot = res_serial; metric_serial_plot = metric_serial; else res_serial_plot = []; metric_serial_plot = []; end if plot_compare_figs figList = plot_algorithm_compare(res_concurrent, res_serial_plot, ... FN_target_N, FN_initial_N, t_sw_start, t_sw_end, serial_fuel_start, ... T4_max_K, SM_LP_min, SM_IP_min, SM_HP_min, FAR_min_phy, FAR_min_ratio_dp, far_idx, metric_concurrent, metric_serial_plot); else figList = gobjects(0); end close all %% ── 步骤7：生成 Simulink 可用控制序列与稳态注入变量 ───────────────────── u_WFM_ts_concurrent = timeseries(res_concurrent.u_phy(1,:), tVec, 'Name','WFM_concurrent'); u_AMSV_ts_concurrent = timeseries(res_concurrent.u_phy(2,:), tVec, 'Name','A_MSV_concurrent'); u_A8_ts_concurrent = timeseries(res_concurrent.u_phy(3,:), tVec, 'Name','A8_concurrent'); for iSerial = 1:numel(serialCases) serialCases(iSerial).u_WFM_ts = timeseries(serialCases(iSerial).res.u_phy(1,:), tVec, ... 'Name',['WFM_', serialCases(iSerial).tag]); serialCases(iSerial).u_AMSV_ts = timeseries(serialCases(iSerial).res.u_phy(2,:), tVec, ... 'Name',['A_MSV_', serialCases(iSerial).tag]); serialCases(iSerial).u_A8_ts = timeseries(serialCases(iSerial).res.u_phy(3,:), tVec, ... 'Name',['A8_', serialCases(iSerial).tag]); end if run_serial_cases u_WFM_ts_serial = serialCases(1).u_WFM_ts; u_AMSV_ts_serial = serialCases(1).u_AMSV_ts; u_A8_ts_serial = serialCases(1).u_A8_ts; else u_WFM_ts_serial = []; u_AMSV_ts_serial = []; u_A8_ts_serial = []; end % 默认指向并发协同算法 u_WFM_ts = u_WFM_ts_concurrent; u_AMSV_ts = u_AMSV_ts_concurrent; u_A8_ts = u_A8_ts_concurrent; % ========================================================= % 【修改3】Simulink 稳态注入变量 % ========================================================= NL_init_phy = initInfo.NL_init_phy; NH_init_phy = initInfo.NH_init_phy; WFM_init_phy = initInfo.WFM_init_phy; AMSV_init_phy = initInfo.AMSV_init_phy; A8_init_phy = initInfo.A8_init_phy; fprintf('\n================ Simulink 稳态注入建议值 ================\n'); fprintf('转速积分器初值：\n'); fprintf(' NL_init_phy = %.6f rpm\n', NL_init_phy); fprintf(' NH_init_phy = %.6f rpm\n', NH_init_phy); fprintf('执行机构初值：\n'); fprintf(' WFM_init_phy = %.6f lb/s\n', WFM_init_phy); fprintf(' AMSV_init_phy = %.6f in^2\n', AMSV_init_phy); fprintf(' A8_init_phy = %.6f ratio\n', A8_init_phy); assignin('base', 'u_WFM_ts', u_WFM_ts); assignin('base', 'u_AMSV_ts', u_AMSV_ts); assignin('base', 'u_A8_ts', u_A8_ts); assignin('base', 'u_WFM_ts_concurrent', u_WFM_ts_concurrent); assignin('base', 'u_AMSV_ts_concurrent', u_AMSV_ts_concurrent); assignin('base', 'u_A8_ts_concurrent', u_A8_ts_concurrent); assignin('base', 'u_WFM_ts_serial', u_WFM_ts_serial); assignin('base', 'u_AMSV_ts_serial', u_AMSV_ts_serial); assignin('base', 'u_A8_ts_serial', u_A8_ts_serial); for iSerial = 1:numel(serialCases) assignin('base', ['u_WFM_ts_', serialCases(iSerial).tag], serialCases(iSerial).u_WFM_ts); assignin('base', ['u_AMSV_ts_', serialCases(iSerial).tag], serialCases(iSerial).u_AMSV_ts); assignin('base', ['u_A8_ts_', serialCases(iSerial).tag], serialCases(iSerial).u_A8_ts); end if inject_simulink_steady assignin('base', 'NL_init_phy', NL_init_phy); assignin('base', 'NH_init_phy', NH_init_phy); assignin('base', 'WFM_init_phy', WFM_init_phy); assignin('base', 'AMSV_init_phy', AMSV_init_phy); assignin('base', 'A8_init_phy', A8_init_phy); else fprintf('[Simulink注入] inject_simulink_steady=false，跳过 NL/NH/WFM/A_MSV/A8 稳态初值写入 base。\n'); end if exist('NRICDyn_init_100A8', 'var') NRICDyn_init = NRICDyn_init_100A8; fprintf('NRICDyn 来源：NRICDyn_init_100A8，即 100%% A8 初始点专用稳态求解器初值。\n'); elseif isfield(Engine, 'WarmStart') && isfield(Engine.WarmStart, 'NRICDyn') NRICDyn_init = Engine.WarmStart.NRICDyn; fprintf(['警告：当前使用 Engine.WarmStart.NRICDyn 作为 NRICDyn_init。\n', ... ' 请确认它确实来自 A8=100%% 初始点；否则建议重新稳态求解并保存为 NRICDyn_init_100A8。\n']); else NRICDyn_init = []; fprintf(['警告：未找到 NRICDyn 初值。\n', ... ' Simulink NR Solver 若使用默认初值，t=0 可能出现初始瞬态。\n']); end if inject_simulink_steady && ~isempty(NRICDyn_init) assignin('base', 'NRICDyn1', NRICDyn_init); if isfield(Engine, 'MWS') && isfield(Engine.MWS, 'Input') Engine.MWS.Input.NRICDyn1 = NRICDyn_init; assignin('base', 'Engine', Engine); end fprintf('NRICDyn1 已写入 base 工作区，元素个数 = %d\n', numel(NRICDyn_init)); fprintf('NRICDyn1 = %s\n', mat2str(NRICDyn_init, 5)); end fprintf('\n已生成控制序列：\n'); fprintf(' 并发协同: u_WFM_ts_concurrent / u_AMSV_ts_concurrent / u_A8_ts_concurrent\n'); for iSerial = 1:numel(serialCases) fprintf(' %s: u_WFM_ts_%s / u_AMSV_ts_%s / u_A8_ts_%s\n', ... serialCases(iSerial).caseName, serialCases(iSerial).tag, ... serialCases(iSerial).tag, serialCases(iSerial).tag); end if run_serial_cases fprintf(' 分段串行默认别名: u_WFM_ts_serial / u_AMSV_ts_serial / u_A8_ts_serial\n'); else fprintf(' 分段串行: 已跳过，串行控制序列别名置为空。\n'); end fprintf('默认 Simulink 输入变量 u_WFM_ts / u_AMSV_ts / u_A8_ts 指向并发协同算法。\n'); %% ── 步骤8：CLM（Simulink部件级模型）一致性对比 ───────────────────────── if run_CLM_compare fprintf('\n================ 步骤8：开始 Simulink / CLM 一致性对比 ================\n'); for iDir = 1:numel(slxDirList) if exist(slxDirList{iDir}, 'dir') addpath(slxDirList{iDir}); end end if ~exist('NRICDyn_init', 'var') NRICDyn_init = []; end clmOpt = struct(); clmOpt.slxName = slxName; clmOpt.simTime = simTime; clmOpt.Ts = Ts; clmOpt.NL_block_path = NL_block_path; clmOpt.NH_block_path = NH_block_path; clmOpt.NL_init_phy = NL_init_phy; clmOpt.NH_init_phy = NH_init_phy; clmOpt.WFM_init_phy = WFM_init_phy; clmOpt.AMSV_init_phy = AMSV_init_phy; clmOpt.A8_init_phy = A8_init_phy; clmOpt.NRICDyn_init = NRICDyn_init; clmOpt.inject_simulink_steady = inject_simulink_steady; clmOpt.Engine = Engine; clm_concurrent = run_clm_compare_case( ... '并发协同算法', ... u_WFM_ts_concurrent, u_AMSV_ts_concurrent, u_A8_ts_concurrent, ... clmOpt); for iSerial = 1:numel(serialCases) serialCases(iSerial).clm = run_clm_compare_case( ... serialCases(iSerial).caseName, ... serialCases(iSerial).u_WFM_ts, ... serialCases(iSerial).u_AMSV_ts, ... serialCases(iSerial).u_A8_ts, ... clmOpt); end if run_serial_cases clm_serial = serialCases(1).clm; else clm_serial = []; end consistency_concurrent = print_lpv_clm_consistency( ... '并发协同算法', res_concurrent, clm_concurrent, t_sw_start); for iSerial = 1:numel(serialCases) serialCases(iSerial).consistency = print_lpv_clm_consistency( ... serialCases(iSerial).caseName, ... serialCases(iSerial).res, serialCases(iSerial).clm, t_sw_start); end metricOptClmConcurrent = metricOpt; metricOptClmConcurrent.t_geom_done_abs = metric_concurrent.GeomSettleAbsTime_s; metric_clm_concurrent = calc_clm_metrics( ... '并发协同算法-CLM', clm_concurrent, metricOptClmConcurrent, t_sw_start); for iSerial = 1:numel(serialCases) metricOptClmSerial = metricOpt; metricOptClmSerial.t_geom_end = serialCases(iSerial).geomEnd; metricOptClmSerial.geom_check_idx = [2 3]; metricOptClmSerial.t_geom_done_abs = serialCases(iSerial).metric.GeomSettleAbsTime_s; serialCases(iSerial).metric_clm = calc_clm_metrics( ... [serialCases(iSerial).caseName, '-CLM'], ... serialCases(iSerial).clm, metricOptClmSerial, serialCases(iSerial).fuelStart); end if run_serial_cases consistency_serial = serialCases(1).consistency; metric_clm_serial = serialCases(1).metric_clm; metric_clm_serial_all = [serialCases.metric_clm]; metricCLMTable = struct2table([metric_clm_concurrent; metric_clm_serial_all(:)]); else consistency_serial = []; metric_clm_serial = []; metric_clm_serial_all = []; metricCLMTable = struct2table(metric_clm_concurrent); end CLM_metric_table = metricCLMTable; CLM_metric_summary_CN = build_metric_summary_cn('CLM部件级模型', metric_clm_concurrent, serialCases, 'metric_clm'); assignin('base', 'CLM_metric_table', CLM_metric_table); assignin('base', 'CLM_metric_summary_CN', CLM_metric_summary_CN); assignin('base', 'MetricExplain_CN', MetricExplain_CN); fprintf('\n================ CLM 过渡态过程时间对比 ================\n'); disp(metricCLMTable); fprintf('\n================ CLM 指标中文摘要：任务总完成时间与推力响应时间分开统计 ================\n'); disp(CLM_metric_summary_CN); print_clm_key_metrics(metric_clm_concurrent, serialCases, settle_tol_rel); for iSerial = 1:numel(serialCases) metricClmS = serialCases(iSerial).metric_clm; if ~isnan(metricClmS.SerialCompleteFromSwitch_s) && ... ~isnan(metric_clm_concurrent.SerialCompleteFromSwitch_s) save_time_clm = metricClmS.SerialCompleteFromSwitch_s - ... metric_clm_concurrent.SerialCompleteFromSwitch_s; save_rate_clm = save_time_clm / ... max(metricClmS.SerialCompleteFromSwitch_s, eps) * 100; fprintf('CLM结果：并发协同算法相对%s节省完整过程时间 %.3f s，占串行法 %.2f%%。\n', ... serialCases(iSerial).caseName, save_time_clm, save_rate_clm); else fprintf('警告：CLM 中并发协同算法或%s未在仿真时长内完成推力/几何串行过程，无法计算节省比例。\n', ... serialCases(iSerial).caseName); end end if plot_CLM_compare_figs figCLM_merged = plot_lpv_clm_compare_merged( ... res_concurrent, clm_concurrent, metric_concurrent, metric_clm_concurrent, ... serialCases, t_sw_start, t_sw_end, ... FN_target_N, FN_initial_N, T4_max_K, SM_LP_min, SM_IP_min, SM_HP_min, ... FAR_min_phy, far_idx, norm_u, norm_y); if plot_CLM_method_figs && run_serial_cases figCLM_method = plot_clm_method_compare( ... clm_concurrent, clm_serial, ... FN_target_N, FN_initial_N, ... t_sw_start, t_sw_end, serial_fuel_start, ... T4_max_K, SM_LP_min, SM_IP_min, SM_HP_min, ... metric_clm_concurrent, metric_clm_serial); end end fprintf('================ 步骤8：Simulink / CLM 一致性对比完成 ================\n\n'); end %% ── 步骤9：可选保存结果 ─────────────────────────────────────────────── if beautify_figures mpc_beautify_figures(findall(groot, 'Type', 'figure')); end %% 保存仿真结果：SP_日期_时分.mat % 当前脚本所在路径 saveDir = fileparts(mfilename('fullpath')); webOutputDir = getenv('LPV_OUTPUT_DIR'); if ~isempty(webOutputDir) saveDir = webOutputDir; if ~exist(saveDir, 'dir'), mkdir(saveDir); end end % 如果在命令行或 Live Script 中运行，mfilename 可能为空，此时退回当前工作路径 if isempty(saveDir) saveDir = pwd; end % 生成时间标签：年月日_时分 timeTag = char(datetime('now', 'Format', 'yyyyMMdd_HHmm')); % 生成保存文件名 saveFileName = ['PS_', timeTag, '.mat']; % 完整保存路径 savePath = fullfile(saveDir, saveFileName); % 保存数据 if saveResults save(savePath, ... 'res_concurrent', 'res_serial', 'serialCases', 'tVec', 'Ts', 'Np', 'Nc', ... 't_sw_start', 't_sw_end', 'serial_fuel_start', 'serial_mode_switch', ... 'FN_target_N', 'T4_max_K', 'SM_LP_min', 'SM_IP_min', 'SM_HP_min', ... 'LPV_metric_table', 'LPV_metric_summary_CN', 'MetricExplain_CN'); fprintf('仿真结果已保存至：%s\n', savePath); save(fullfile(saveDir,'compare2_MPC_result.mat'), ... 'res_concurrent','res_serial','serialCases','metricTable','LPV_metric_table','LPV_metric_summary_CN','MetricExplain_CN', ... 'metric_concurrent','metric_serial','metric_serial_all', ... 'u_WFM_ts_concurrent','u_AMSV_ts_concurrent','u_A8_ts_concurrent', ... 'u_WFM_ts_serial','u_AMSV_ts_serial','u_A8_ts_serial', ... 'u_WFM_ts','u_AMSV_ts','u_A8_ts', ... 'NL_init_phy','NH_init_phy','WFM_init_phy','AMSV_init_phy','A8_init_phy','NRICDyn_init', ... 'tVec','Ts','simTime','t_sw_start','t_sw_end','serial_fuel_start','serial_mode_switch','runSimulinkCompare', ... 'A8_INIT_FACTOR','A8_TGT_FACTOR','A8_CTRL_MIN_FACTOR','A8_CTRL_MAX_FACTOR', ... 'WFM_rate_max','AMSV_rate_max','A8_rate_max'); if exist('clm_concurrent', 'var') && exist('clm_serial', 'var') save(fullfile(saveDir,'compare2_MPC_CLM_result.mat'), ... 'clm_concurrent','clm_serial','serialCases', ... 'consistency_concurrent','consistency_serial', ... 'metric_clm_concurrent','metric_clm_serial','metric_clm_serial_all','metricCLMTable','CLM_metric_table','CLM_metric_summary_CN','MetricExplain_CN'); end end if saveFig if isempty(saveDir) saveDir = pwd; end if ~exist(saveDir,'dir'), mkdir(saveDir); end if plot_compare_figs && exist('figList','var') for iFig = 1:numel(figList) if isgraphics(figList(iFig)) exportgraphics(figList(iFig), fullfile(saveDir, sprintf('compare2_fig%d.%s', iFig, figFormat)), 'ContentType','vector', 'BackgroundColor','white'); end end end if exist('clm_concurrent', 'var') && exist('clm_serial', 'var') if exist('figCLM_merged','var') for iFig = 1:numel(figCLM_merged) if isgraphics(figCLM_merged(iFig)) exportgraphics(figCLM_merged(iFig), ... fullfile(saveDir, sprintf('CLM_LPV_merged_fig%d.%s', iFig, figFormat)), 'ContentType','vector', 'BackgroundColor','white'); end end end if exist('figCLM_concurrent','var') for iFig = 1:numel(figCLM_concurrent) if isgraphics(figCLM_concurrent(iFig)) exportgraphics(figCLM_concurrent(iFig), ... fullfile(saveDir, sprintf('CLM_concurrent_fig%d.%s', iFig, figFormat)), 'ContentType','vector', 'BackgroundColor','white'); end end end if exist('figCLM_serial','var') for iFig = 1:numel(figCLM_serial) if isgraphics(figCLM_serial(iFig)) exportgraphics(figCLM_serial(iFig), ... fullfile(saveDir, sprintf('CLM_serial_fig%d.%s', iFig, figFormat)), 'ContentType','vector', 'BackgroundColor','white'); end end end if exist('figCLM_method','var') && isgraphics(figCLM_method) exportgraphics(figCLM_method, fullfile(saveDir, ['CLM_method_compare.', figFormat]), 'ContentType','vector', 'BackgroundColor','white'); end end fprintf('结果已保存至：%s\n', saveDir); end %% ======================================================================== % 局部函数 % ======================================================================== function [ur_traj, yr_traj, xr_traj] = build_reference_trajectory( ... tVec, t_geom_start, t_geom_end, t_FN_step, ... x0_norm, xr_norm, u0_norm, ur_norm, y0_norm, FN_target_norm, nu_mpc) %#ok % t_geom_end 保留在函数接口中，便于兼容原有调用与串行时序打印。 % 【本次修改】只保留最终几何目标，不再构造 A_MSV/A8 的中间线性参考轨迹。 % 物理含义： % 1) t < t_geom_start 时，几何目标保持初始值； % 2) t >= t_geom_start 后，只告诉 MPC 最终目标 A_MSV_target/A8_target； % 3) 中间过渡路径由 MPC 在速率限制、T4/SM 约束和推力目标下自主生成。 % 注意：Qu 已设置为 diag([0,0,0])，因此该 ur_traj 不作为过程跟踪轨迹， % 主要通过 Pu_terminal 的终端几何目标进入优化问题。 Nsim = numel(tVec); ur_traj = repmat(u0_norm(1:nu_mpc), 1, Nsim); idx_geom = find(tVec >= t_geom_start, 1, 'first'); if ~isempty(idx_geom) ur_traj(2, idx_geom:end) = ur_norm(2); ur_traj(3, idx_geom:end) = ur_norm(3); end yr_traj = repmat(y0_norm, 1, Nsim); idx_FN = find(tVec >= t_FN_step, 1, 'first'); if ~isempty(idx_FN) yr_traj(4, idx_FN:end) = FN_target_norm; end xr_traj = repmat(x0_norm, 1, Nsim); idx_xr = find(tVec >= t_FN_step, 1, 'first'); if ~isempty(idx_xr) xr_traj(:, idx_xr:end) = repmat(xr_norm, 1, Nsim-idx_xr+1); end end function res = run_lpv_mpc_case(caseName, sysLPV3, ... x0_norm, u0_norm_ctrl, norm_x, norm_u, norm_y, ... ur_traj, yr_traj, xr_traj, p) nx = p.nx; nu = p.nu_mpc; ny = p.ny; Nsim = p.Nsim; tVec = p.tVec; x_hist = zeros(nx, Nsim+1); u_hist = zeros(nu, Nsim); y_hist = zeros(ny, Nsim); x_hist(:,1) = x0_norm; u_prev_ctrl = u0_norm_ctrl(:); qp_fail = 0; fprintf(' %5s | %8s %8s | %8s %10s %8s | %8s %10s\n', ... 'Time','NL(%)','NH(%)','A8','A_MSV','WFM','T4','FN'); fprintf(' %s\n', repmat('-',1,92)); for k = 1:Nsim xnlc_s = x_hist(1,k); msv_s = max(min(u_prev_ctrl(2), max(p.gMSV)), min(p.gMSV)); a8_s = max(min(u_prev_ctrl(3), max(p.gA8)), min(p.gA8)); [lpv_ss, lpv_off] = sample(sysLPV3, [], xnlc_s, msv_s, a8_s); A = lpv_ss.A; B = lpv_ss.B; C = lpv_ss.C; D = lpv_ss.D; Xss = lpv_off.x; Uss = lpv_off.u; Yss = lpv_off.y; dx0 = x_hist(:,k) - Xss; [Phi, Gamma, CG_total, Cpred] = mpc_build_prediction(A, B, C, D, p.Np, p.Nc); if k < p.t_sw_idx Qy_t = zeros(ny); Px_t = zeros(nx); Pu_t = zeros(nu); else Qy_t = p.Qy; Px_t = p.Px; Pu_t = p.Pu_terminal; end k_end = min(k+p.Np-1, Nsim); yr_h = yr_traj(:, k:k_end); xr_h = xr_traj(:, k:k_end); if size(yr_h,2) < p.Np yr_h = [yr_h, repmat(yr_h(:,end), 1, p.Np-size(yr_h,2))]; xr_h = [xr_h, repmat(xr_h(:,end), 1, p.Np-size(xr_h,2))]; end ur_hrz = ur_traj(:, k:min(k+p.Nc-1, Nsim)); if size(ur_hrz,2) < p.Nc ur_hrz = [ur_hrz, repmat(ur_hrz(:,end), 1, p.Nc-size(ur_hrz,2))]; end ur_hrz_stk = ur_hrz(:); ur_hrz_end = ur_traj(:, min(k+p.Np-1, Nsim)); Du_prev_offset = repmat(D*(u_prev_ctrl - Uss), p.Np, 1); Yss_stk = repmat(Yss, p.Np, 1); Y_free = Cpred*Phi*dx0 + Du_prev_offset; E_free = yr_h(:) - Yss_stk - Y_free; [H, f, Aineq, bineq] = mpc_build_qp(CG_total, Phi, Cpred, Gamma, ... dx0, Xss, Uss, Yss, u_prev_ctrl, ... Qy_t, p.Ru, p.Qu, Px_t, Pu_t, ... E_free, Y_free, Yss_stk, ... ur_hrz_stk, ur_hrz_end, xr_h, ... p.u_min_norm, p.u_max_norm, p.du_max, p.du_min, ... p.T4_max_norm, p.SM_LP_min_norm, p.SM_IP_min_norm, p.SM_HP_min_norm, ... p.Np, p.Nc, nx, nu, ny); [Aineq, bineq] = append_far_min_constraint(Aineq, bineq, CG_total, Y_free, Yss_stk, ... p.FAR_min_norm, p.far_idx, p.Np, ny); H = (H + H')/2; [DU_opt, ~, exitflag] = quadprog(H, f, Aineq, bineq, [], [], [], [], [], p.qp_opts); if exitflag <= 0 || isempty(DU_opt) qp_fail = qp_fail + 1; du_apply = zeros(nu,1); else du_apply = DU_opt(1:nu); end u_apply = min(max(u_prev_ctrl + du_apply, p.u_min_norm), p.u_max_norm); if isfield(p, 'holdGeomUntilAbs') && tVec(k) < p.holdGeomUntilAbs u_apply(2:3) = u0_norm_ctrl(2:3); end y_hist(:,k) = Yss + C*dx0 + D*(u_apply - Uss); x_hist(:,k+1) = Xss + A*dx0 + B*(u_apply - Uss); u_hist(:,k) = u_apply; u_prev_ctrl = u_apply; if mod(k, p.printEvery) == 0 fprintf(' t=%5.1fs | NL=%7.2f NH=%7.2f | A8=%8.4f A_MSV=%8.3f WFM=%7.4f | T4=%8.1f FN=%9.2f\n', ... tVec(k), ... x_hist(1,k)*100, x_hist(2,k)*100, ... u_apply(3)*norm_u(3), u_apply(2)*norm_u(2), u_apply(1)*norm_u(1), ... y_hist(3,k)*norm_y(3), y_hist(4,k)*norm_y(4)); end end fprintf('[%s] QP 不可行/失败次数：%d / %d\n', caseName, qp_fail, Nsim); res = struct(); res.caseName = caseName; res.tVec = tVec; res.x_norm = x_hist(:,1:Nsim); res.u_norm = u_hist; res.y_norm = y_hist; res.x_phy = x_hist(:,1:Nsim) .* norm_x; res.u_phy = u_hist .* norm_u(1:nu); res.y_phy = y_hist .* norm_y; res.ur_norm = ur_traj; res.yr_norm = yr_traj; res.xr_norm = xr_traj; res.ur_phy = ur_traj .* norm_u(1:nu); res.yr_phy = yr_traj .* norm_y; res.xr_phy = xr_traj .* norm_x; res.qp_fail = qp_fail; end function res = run_lpv_mpc_event_serial_case(caseDef, sysLPV3, ... x0_norm, u0_norm_ctrl, norm_x, norm_u, norm_y, ... ur_final_norm, yr_final_norm, xr_final_norm, p) nx = p.nx; nu = p.nu_mpc; ny = p.ny; Nsim = p.Nsim; tVec = p.tVec; x_hist = zeros(nx, Nsim+1); u_hist = zeros(nu, Nsim); y_hist = zeros(ny, Nsim); ur_hist = zeros(nu, Nsim); yr_hist = zeros(ny, Nsim); xr_hist = zeros(nx, Nsim); stage_id_hist = zeros(1, Nsim); x_hist(:,1) = x0_norm; u_prev_ctrl = u0_norm_ctrl(:); qp_fail = 0; stage = "PRE"; prevStage = stage; stageHoldSteps = max(1, round(p.stage_hold_time / p.Ts)); FN_mid_tol_abs = p.FN_tol_rel * abs(p.FN_mid_N); FN_final_tol_abs = p.FN_tol_rel * abs(p.FN_target_N); AMSV_tol_phy = p.AMSV_mode_tol_phy; holdCounterFNMid = 0; holdCounterFNFinal = 0; holdCounterMode = 0; holdCounterModeAndFN = 0; stageTimes = struct(); stageTimes.StageAStartAbsTime_s = NaN; stageTimes.StageBStartAbsTime_s = NaN; stageTimes.StageCStartAbsTime_s = NaN; stageTimes.DoneAbsTime_s = NaN; stageTimes.GeomStartAbsTime_s = NaN; stageTimes.GeomEndAbsTime_s = NaN; fprintf(' [事件串行] %s\n', caseDef.caseName); if strcmp(caseDef.id, 'serial1_three_stage') fprintf(' A段: WFM/A8自由, A_MSV锁定, FN -> %.1f N\n', p.FN_mid_N); fprintf(' B段: WFM/A_MSV/A8自由, FN保持 %.1f N, A_MSV -> %.3f in^2\n', ... p.FN_mid_N, p.AMSV_target_phy); fprintf(' C段: WFM/A8自由, A_MSV锁定, FN -> %.1f N\n', p.FN_target_N); else fprintf(' A段: WFM/A8自由, A_MSV锁定, FN -> %.1f N\n', p.FN_target_N); fprintf(' B段: WFM/A_MSV/A8自由, FN保持 %.1f N, A_MSV -> %.3f in^2\n', ... p.FN_target_N, p.AMSV_target_phy); end fprintf(' 切换判据: FN误差 <= %.2f%%, A_MSV容差 %.4f in^2, 保持时间 %.2f s\n', ... 100*p.FN_tol_rel, AMSV_tol_phy, p.stage_hold_time); fprintf(' %5s | %8s %8s | %8s %10s %8s | %8s %10s | %8s\n', ... 'Time','NL(%)','NH(%)','A8','A_MSV','WFM','T4','FN','Stage'); fprintf(' %s\n', repmat('-',1,104)); for k = 1:Nsim if strcmp(stage, "PRE") && tVec(k) >= p.t_start stage = "A"; stageTimes.StageAStartAbsTime_s = tVec(k); fprintf(' [Stage] %s enters A at t=%.3f s.\n', caseDef.plotName, tVec(k)); end [y_ref_now, x_ref_now, u_ref_now, Qy_t, Px_t, Pu_t, du_min_stage, du_max_stage, stageID] = ... get_event_serial_stage_setting(stage, caseDef.id, p, x0_norm, u0_norm_ctrl, ... ur_final_norm, yr_final_norm, xr_final_norm); if ~strcmp(stage, prevStage) prevStage = stage; end xnlc_s = x_hist(1,k); msv_s = max(min(u_prev_ctrl(2), max(p.gMSV)), min(p.gMSV)); a8_s = max(min(u_prev_ctrl(3), max(p.gA8)), min(p.gA8)); [lpv_ss, lpv_off] = sample(sysLPV3, [], xnlc_s, msv_s, a8_s); A = lpv_ss.A; B = lpv_ss.B; C = lpv_ss.C; D = lpv_ss.D; Xss = lpv_off.x; Uss = lpv_off.u; Yss = lpv_off.y; dx0 = x_hist(:,k) - Xss; [Phi, Gamma, CG_total, Cpred] = mpc_build_prediction(A, B, C, D, p.Np, p.Nc); yr_h = repmat(y_ref_now, 1, p.Np); xr_h = repmat(x_ref_now, 1, p.Np); ur_hrz = repmat(u_ref_now, 1, p.Nc); ur_hrz_stk = ur_hrz(:); ur_hrz_end = u_ref_now; Du_prev_offset = repmat(D*(u_prev_ctrl - Uss), p.Np, 1); Yss_stk = repmat(Yss, p.Np, 1); Y_free = Cpred*Phi*dx0 + Du_prev_offset; E_free = yr_h(:) - Yss_stk - Y_free; [H, f, Aineq, bineq] = mpc_build_qp(CG_total, Phi, Cpred, Gamma, ... dx0, Xss, Uss, Yss, u_prev_ctrl, ... Qy_t, p.Ru, p.Qu, Px_t, Pu_t, ... E_free, Y_free, Yss_stk, ... ur_hrz_stk, ur_hrz_end, xr_h, ... p.u_min_norm, p.u_max_norm, du_max_stage, du_min_stage, ... p.T4_max_norm, p.SM_LP_min_norm, p.SM_IP_min_norm, p.SM_HP_min_norm, ... p.Np, p.Nc, nx, nu, ny); [Aineq, bineq] = append_far_min_constraint(Aineq, bineq, CG_total, Y_free, Yss_stk, ... p.FAR_min_norm, p.far_idx, p.Np, ny); H = (H + H')/2; [DU_opt, ~, exitflag] = quadprog(H, f, Aineq, bineq, [], [], [], [], [], p.qp_opts); if exitflag <= 0 || isempty(DU_opt) qp_fail = qp_fail + 1; du_apply = zeros(nu,1); else du_apply = DU_opt(1:nu); end u_apply = min(max(u_prev_ctrl + du_apply, p.u_min_norm), p.u_max_norm); y_hist(:,k) = Yss + C*dx0 + D*(u_apply - Uss); x_hist(:,k+1) = Xss + A*dx0 + B*(u_apply - Uss); u_hist(:,k) = u_apply; ur_hist(:,k) = u_ref_now; yr_hist(:,k) = y_ref_now; xr_hist(:,k) = x_ref_now; stage_id_hist(k) = stageID; FN_now = y_hist(4,k) * norm_y(4); AMSV_now = u_apply(2) * norm_u(2); FNMidOK = abs(FN_now - p.FN_mid_N) <= FN_mid_tol_abs; FNFinalOK = abs(FN_now - p.FN_target_N) <= FN_final_tol_abs; ModeOK = abs(AMSV_now - p.AMSV_target_phy) <= AMSV_tol_phy; tNext = min(tVec(k) + p.Ts, tVec(end)); if strcmp(caseDef.id, 'serial1_three_stage') switch stage case "A" holdCounterFNMid = update_hold_counter(holdCounterFNMid, FNMidOK); if holdCounterFNMid >= stageHoldSteps stage = "B"; stageTimes.StageBStartAbsTime_s = tNext; stageTimes.GeomStartAbsTime_s = tNext; holdCounterMode = 0; fprintf(' [Stage] %s A->B at t=%.3f s: FN reached %.1f N band.\n', ... caseDef.plotName, tNext, p.FN_mid_N); end case "B" holdCounterMode = update_hold_counter(holdCounterMode, ModeOK); if holdCounterMode >= stageHoldSteps stage = "C"; stageTimes.StageCStartAbsTime_s = tNext; stageTimes.GeomEndAbsTime_s = tNext; holdCounterFNFinal = 0; fprintf(' [Stage] %s B->C at t=%.3f s: A_MSV reached target mode.\n', ... caseDef.plotName, tNext); end case "C" holdCounterFNFinal = update_hold_counter(holdCounterFNFinal, FNFinalOK); if holdCounterFNFinal >= stageHoldSteps && isnan(stageTimes.DoneAbsTime_s) stage = "DONE"; stageTimes.DoneAbsTime_s = tNext; fprintf(' [Stage] %s DONE at t=%.3f s: final FN reached %.1f N band.\n', ... caseDef.plotName, tNext, p.FN_target_N); end end else switch stage case "A" holdCounterFNFinal = update_hold_counter(holdCounterFNFinal, FNFinalOK); if holdCounterFNFinal >= stageHoldSteps stage = "B"; stageTimes.StageBStartAbsTime_s = tNext; stageTimes.GeomStartAbsTime_s = tNext; holdCounterModeAndFN = 0; fprintf(' [Stage] %s A->B at t=%.3f s: final FN reached %.1f N band.\n', ... caseDef.plotName, tNext, p.FN_target_N); end case "B" holdCounterModeAndFN = update_hold_counter(holdCounterModeAndFN, ModeOK && FNFinalOK); if holdCounterModeAndFN >= stageHoldSteps && isnan(stageTimes.DoneAbsTime_s) stage = "DONE"; stageTimes.DoneAbsTime_s = tNext; stageTimes.GeomEndAbsTime_s = tNext; fprintf(' [Stage] %s DONE at t=%.3f s: mode completed and FN is stable.\n', ... caseDef.plotName, tNext); end end end u_prev_ctrl = u_apply; if mod(k, p.printEvery) == 0 fprintf(' t=%5.1fs | NL=%7.2f NH=%7.2f | A8=%8.4f A_MSV=%8.3f WFM=%7.4f | T4=%8.1f FN=%9.2f | %8s\n', ... tVec(k), ... x_hist(1,k)*100, x_hist(2,k)*100, ... u_apply(3)*norm_u(3), u_apply(2)*norm_u(2), u_apply(1)*norm_u(1), ... y_hist(3,k)*norm_y(3), y_hist(4,k)*norm_y(4), char(stage)); end end fprintf('[%s] QP 不可行/失败次数：%d / %d\n', caseDef.caseName, qp_fail, Nsim); res = struct(); res.caseName = caseDef.caseName; res.tVec = tVec; res.x_norm = x_hist(:,1:Nsim); res.u_norm = u_hist; res.y_norm = y_hist; res.x_phy = x_hist(:,1:Nsim) .* norm_x; res.u_phy = u_hist .* norm_u(1:nu); res.y_phy = y_hist .* norm_y; res.ur_norm = ur_hist; res.yr_norm = yr_hist; res.xr_norm = xr_hist; res.ur_phy = ur_hist .* norm_u(1:nu); res.yr_phy = yr_hist .* norm_y; res.xr_phy = xr_hist .* norm_x; res.stage_id = stage_id_hist; res.mode_factor = res.u_phy(2,:) ./ max(p.AMSV_target_phy, eps); res.stageTimes = stageTimes; StageName = string({'A段开始'; 'B段开始/几何开始'; 'C段开始/几何结束'; '任务完成'}); AbsTime_s = [stageTimes.StageAStartAbsTime_s; ... stageTimes.StageBStartAbsTime_s; ... stageTimes.StageCStartAbsTime_s; ... stageTimes.DoneAbsTime_s]; res.stageTimesTable_CN = table(StageName, AbsTime_s); res.qp_fail = qp_fail; end function [y_ref, x_ref, u_ref, Qy_t, Px_t, Pu_t, du_min_stage, du_max_stage, stageID] = ... get_event_serial_stage_setting(stage, serialID, p, x0_norm, u0_norm, ur_final_norm, yr_final_norm, xr_final_norm) %#ok yr_final_norm 当前保留在接口中，便于后续扩展。 y_ref = p.y0_norm; x_ref = x0_norm; u_ref = u0_norm; Qy_t = zeros(p.ny); Px_t = zeros(p.nx); Pu_t = zeros(p.nu_mpc); du_min_stage = p.du_min; du_max_stage = p.du_max; stageID = 0; if strcmp(stage, "PRE") % 任务启动前：不设置推力/几何目标。控制增量代价会使控制量保持稳态。 return; end Qy_t = p.Qy; Px_t = p.Px; % 串行算法中，A8 虽然始终允许自由调节，但其目标值应为最终 A8。 % 否则优化器只会把 A8 当作“辅助输入”，而不会主动把它推到目标值。 u_ref(3) = ur_final_norm(3); switch serialID case 'serial1_three_stage' switch stage case "A" % A段：WFM/A8自由，A_MSV锁定，FN -> 2474 N，同时要求 A8 逐步到达最终目标。 y_ref(4) = p.FN_mid_norm; du_min_stage(2) = 0; du_max_stage(2) = 0; Pu_t = p.Pu_terminal_serial_A8; stageID = 1; case "B" % B段：WFM/A_MSV/A8自由，FN保持2474 N，A_MSV 与 A8 完成模态/几何切换。 y_ref(4) = p.FN_mid_norm; u_ref(2) = ur_final_norm(2); Pu_t = p.Pu_terminal_serial_mode; stageID = 2; case {"C", "DONE"} % C段：WFM/A8自由，A_MSV锁定在目标模态，FN -> 最终目标；A8 继续保持最终目标。 y_ref(4) = p.FN_target_norm; x_ref = xr_final_norm; u_ref(2) = ur_final_norm(2); du_min_stage(2) = 0; du_max_stage(2) = 0; Pu_t = p.Pu_terminal_serial_A8; stageID = 3; end case 'serial2_two_stage' switch stage case "A" % A段：WFM/A8自由，A_MSV锁定，FN -> 最终目标，同时要求 A8 逐步到达最终目标。 y_ref(4) = p.FN_target_norm; x_ref = xr_final_norm; du_min_stage(2) = 0; du_max_stage(2) = 0; Pu_t = p.Pu_terminal_serial_A8; stageID = 1; case {"B", "DONE"} % B段：WFM/A_MSV/A8自由，FN保持最终目标，A_MSV 与 A8 完成最终几何切换。 y_ref(4) = p.FN_target_norm; x_ref = xr_final_norm; u_ref(2) = ur_final_norm(2); Pu_t = p.Pu_terminal_serial_mode; stageID = 2; end end end function assert_required_fields(s, sName, fields) missing = fields(~isfield(s, fields)); if ~isempty(missing) error('%s 缺少必要字段: %s', sName, strjoin(missing, ', ')); end end function [Aineq, bineq] = append_far_min_constraint(Aineq, bineq, CG_total, Y_free, Yss_stk, FAR_min_norm, far_idx, Np, ny) nVar = size(CG_total, 2); A_far = zeros(Np, nVar); b_far = zeros(Np, 1); for iPred = 1:Np row = (iPred-1)*ny + far_idx; A_far(iPred,:) = -CG_total(row,:); b_far(iPred) = Yss_stk(row) + Y_free(row) - FAR_min_norm; end Aineq = [Aineq; A_far]; bineq = [bineq; b_far]; end function counter = update_hold_counter(counter, isOK) if isOK counter = counter + 1; else counter = 0; end end function m = calc_case_metrics(res, opt, t_fuel_start) t = res.tVec(:); FN = res.y_phy(4,:).'; T4 = res.y_phy(3,:).'; SM = res.y_phy(5:7,:).'; u = res.u_phy; holdSteps = max(1, round(opt.hold_time / opt.Ts)); idxStart = find(t >= opt.t_start, 1, 'first'); FN_tol_abs = opt.FN_tol_rel * abs(opt.FNtgt); tFNSettleAbs = first_hold_time(t, abs(FN - opt.FNtgt) <= FN_tol_abs, idxStart, holdSteps); [tFNStableAbs, ~] = calc_fn_stable_time(t, FN, opt.t_start); deltaFN = opt.FNtgt - opt.FN0; if abs(deltaFN) < eps tFN90Abs = NaN; else FN90 = opt.FN0 + 0.9*deltaFN; if deltaFN > 0 idx90 = find(t >= opt.t_start & FN >= FN90, 1, 'first'); else idx90 = find(t >= opt.t_start & FN <= FN90, 1, 'first'); end if isempty(idx90) tFN90Abs = NaN; else tFN90Abs = t(idx90); end end if isfield(opt, 'geom_check_idx') && ~isempty(opt.geom_check_idx) geomCheckIdx = opt.geom_check_idx(:).'; else geomCheckIdx = [2 3]; end geomOK = true(numel(t),1); for iGeom = geomCheckIdx if iGeom == 2 tol_i = max(0.02*abs(opt.u_tgt_phy(2) - u(2,1)), 0.10); elseif iGeom == 3 tol_i = max(0.02*abs(opt.u_tgt_phy(3) - u(3,1)), 0.002); else tol_i = max(0.02*abs(opt.u_tgt_phy(iGeom) - u(iGeom,1)), 1e-6); end geomOK = geomOK & abs(u(iGeom,:).' - opt.u_tgt_phy(iGeom)) <= tol_i; end tGeomAbs = first_hold_time(t, geomOK, idxStart, holdSteps); tSerialCompleteAbs = max_time_ignore_nan(tGeomAbs, tFNSettleAbs); idxStage1 = t >= opt.t_start & t <= opt.t_geom_end; if any(idxStage1) FNStage1FluctAbs = max(abs(FN(idxStage1) - opt.FN0)); else FNStage1FluctAbs = NaN; end signDelta = sign(deltaFN); if signDelta == 0 signDelta = 1; end idxAfterStart = t >= opt.t_start; overshootAbs = max([0; signDelta*(FN(idxAfterStart) - opt.FNtgt)]); undershootAbs = max([0; -signDelta*(FN(idxAfterStart) - opt.FN0)]); m = struct(); m.Method = string(res.caseName); m.MetricNote_CN = string('任务总完成时间=几何到达最终目标且推力进入目标误差带后的较晚时刻；推力自身响应时间另按从任务开始/从燃油启用分别统计。'); m.GeomSettleAbsTime_s = tGeomAbs; m.GeomSettleFromSwitch_s = subtract_or_nan(tGeomAbs, opt.t_start); m.FN90AbsTime_s = tFN90Abs; m.FN90FromSwitch_s = subtract_or_nan(tFN90Abs, opt.t_start); m.FNSettleAbsTime_s = tFNSettleAbs; m.FNSettleTimeFromSwitch_s = subtract_or_nan(tFNSettleAbs, opt.t_start); m.FNSettleTimeFromFuel_s = subtract_or_nan(tFNSettleAbs, t_fuel_start); m.FNStableAbsTime_s = tFNStableAbs; m.FNStableFromSwitch_s = subtract_or_nan(tFNStableAbs, opt.t_start); m.FNStableFromFuel_s = subtract_or_nan(tFNStableAbs, t_fuel_start); m.SerialCompleteAbsTime_s = tSerialCompleteAbs; m.SerialCompleteFromSwitch_s = subtract_or_nan(tSerialCompleteAbs, opt.t_start); m.SerialCompleteFromFuel_s = subtract_or_nan(tSerialCompleteAbs, t_fuel_start); m.Stage1FNFluctAbs = FNStage1FluctAbs; m.Stage1FNFluctRelPct = FNStage1FluctAbs / max(abs(opt.FN0), eps) * 100; m.FNOvershootAbs = overshootAbs; m.FNOvershootRelPct = overshootAbs / max(abs(deltaFN), eps) * 100; m.FNUndershootAbs = undershootAbs; Yall = res.y_phy.'; [tAllStableAbs, maxFinalRelPct] = calc_all_signal_stable_time( ... t, Yall, opt.t_start, opt.FN_tol_rel, holdSteps); idxTrack = t >= t_fuel_start; if any(idxTrack) errFNTrack = FN(idxTrack) - opt.FNtgt; m.FNTrackMAE_N = mean(abs(errFNTrack)); m.FNTrackRMSE_N = sqrt(mean(errFNTrack.^2)); m.FNTrackMaxAbs_N = max(abs(errFNTrack)); m.FNTrackRMSE_RelPct = m.FNTrackRMSE_N / max(abs(opt.FNtgt), eps) * 100; else m.FNTrackMAE_N = NaN; m.FNTrackRMSE_N = NaN; m.FNTrackMaxAbs_N = NaN; m.FNTrackRMSE_RelPct = NaN; end m.AllParamStableAbsTime_s = tAllStableAbs; m.AllParamStableFromSwitch_s = subtract_or_nan(tAllStableAbs, opt.t_start); m.AllParamMaxFinalRelPct = maxFinalRelPct; m.T4Max = max(T4); m.FAN_SM_Min = min(SM(:,1)); m.IPC_SM_Min = min(SM(:,2)); m.HPC_SM_Min = min(SM(:,3)); if size(res.y_phy,1) >= opt.far_idx m.FAR_Min = min(res.y_phy(opt.far_idx,:)); else m.FAR_Min = NaN; end m.QPFail = res.qp_fail; end function tHit = first_hold_time(t, okFlag, idxStart, holdSteps) tHit = NaN; okFlag = okFlag(:); if isempty(idxStart) || idxStart > numel(okFlag) return; end lastIdx = numel(okFlag) - holdSteps + 1; for i = idxStart:lastIdx if all(okFlag(i:i+holdSteps-1)) tHit = t(i); return; end end end function y = subtract_or_nan(a, b) if isnan(a) y = NaN; else y = a - b; end end function y = max_time_ignore_nan(a, b) vals = [a, b]; vals = vals(isfinite(vals)); if isempty(vals) y = NaN; else y = max(vals); end end function [tStableAbs, maxFinalRelPct] = calc_all_signal_stable_time(t, Y, t_start, tolRel, holdSteps) tStableAbs = NaN; maxFinalRelPct = NaN; if isempty(Y) return; end t = t(:); Y = double(Y); N = min(numel(t), size(Y,1)); t = t(1:N); Y = Y(1:N,:); tailN = min(max(holdSteps, 1), N); tailY = Y(N-tailN+1:N,:); yFinal = zeros(1, size(Y,2)); for iCh = 1:size(Y,2) vals = tailY(:,iCh); vals = vals(isfinite(vals)); if isempty(vals) yFinal(iCh) = NaN; else yFinal(iCh) = mean(vals); end end denom = max(abs(yFinal), 1); relErr = abs(Y - yFinal) ./ denom; okFlag = all(relErr <= tolRel, 2); idxStart = find(t >= t_start, 1, 'first'); tStableAbs = first_hold_time(t, okFlag, idxStart, holdSteps); idxAfterStart = find(t >= t_start); if ~isempty(idxAfterStart) vals = relErr(idxAfterStart,:); vals = vals(isfinite(vals)); if ~isempty(vals) maxFinalRelPct = max(vals) * 100; end end end function tbl = build_metric_summary_cn(dataSource, metricConcurrent, serialCases, metricFieldName) metrics = metricConcurrent; methodCN = {char(metricConcurrent.Method)}; for ii = 1:numel(serialCases) if isfield(serialCases(ii), metricFieldName) metrics(end+1) = serialCases(ii).(metricFieldName); %#ok methodCN{end+1} = serialCases(ii).caseName; %#ok end end n = numel(metrics); DataSource_CN = strings(n,1); Method_CN = strings(n,1); GeomComplete_FromTaskStart_s = nan(n,1); ThrustSettle_FromTaskStart_s = nan(n,1); ThrustSettle_FromFuelStart_s = nan(n,1); TaskComplete_FromTaskStart_s = nan(n,1); TaskComplete_FromFuelStart_s = nan(n,1); FN_RMSE_N = nan(n,1); T4_Max_K = nan(n,1); FAN_SM_Min_pct = nan(n,1); IPC_SM_Min_pct = nan(n,1); HPC_SM_Min_pct = nan(n,1); FAR_Min = nan(n,1); QPFail = nan(n,1); Comment_CN = strings(n,1); for ii = 1:n DataSource_CN(ii) = string(dataSource); Method_CN(ii) = string(methodCN{ii}); GeomComplete_FromTaskStart_s(ii) = metrics(ii).GeomSettleFromSwitch_s; ThrustSettle_FromTaskStart_s(ii) = metrics(ii).FNSettleTimeFromSwitch_s; ThrustSettle_FromFuelStart_s(ii) = metrics(ii).FNSettleTimeFromFuel_s; TaskComplete_FromTaskStart_s(ii) = metrics(ii).SerialCompleteFromSwitch_s; TaskComplete_FromFuelStart_s(ii) = metrics(ii).SerialCompleteFromFuel_s; FN_RMSE_N(ii) = metrics(ii).FNTrackRMSE_N; T4_Max_K(ii) = metrics(ii).T4Max; FAN_SM_Min_pct(ii) = metrics(ii).FAN_SM_Min; IPC_SM_Min_pct(ii) = metrics(ii).IPC_SM_Min; HPC_SM_Min_pct(ii) = metrics(ii).HPC_SM_Min; if isfield(metrics(ii), 'FAR_Min') FAR_Min(ii) = metrics(ii).FAR_Min; end if isfield(metrics(ii), 'QPFail') QPFail(ii) = metrics(ii).QPFail; end Comment_CN(ii) = "TaskComplete 采用 max(几何最终目标到达时间, 推力进入目标误差带时间)；ThrustSettle_FromFuelStart 用于消除 Seq-I 推力目标晚启用造成的天然延迟。"; end tbl = table(DataSource_CN, Method_CN, ... GeomComplete_FromTaskStart_s, ... ThrustSettle_FromTaskStart_s, ThrustSettle_FromFuelStart_s, ... TaskComplete_FromTaskStart_s, TaskComplete_FromFuelStart_s, ... FN_RMSE_N, T4_Max_K, FAN_SM_Min_pct, IPC_SM_Min_pct, HPC_SM_Min_pct, FAR_Min, QPFail, Comment_CN); end function tbl = build_metric_explain_cn() MetricName = string({ ... 'GeomComplete_FromTaskStart_s'; ... 'ThrustSettle_FromTaskStart_s'; ... 'ThrustSettle_FromFuelStart_s'; ... 'TaskComplete_FromTaskStart_s'; ... 'TaskComplete_FromFuelStart_s'; ... 'FN_RMSE_N'; ... 'T4_Max_K'; ... 'FAN_SM_Min_pct / IPC_SM_Min_pct / HPC_SM_Min_pct'; ... 'FAR_Min'; ... 'QPFail'}); ChineseName = string({ ... '几何最终目标完成时间'; ... '推力进入目标误差带时间：从统一任务开始计时'; ... '推力进入目标误差带时间：从该方法燃油/推力目标启用时刻计时'; ... '总任务完成时间：从统一任务开始计时'; ... '总任务完成时间：从该方法燃油/推力目标启用时刻计时'; ... '推力跟踪均方根误差'; ... '过渡过程最高 T4'; ... '三级喘振裕度最小值'; ... '油气比最小值'; ... 'QP 求解失败次数'}); Meaning_CN = string({ ... 'A_MSV 与 A8 同时进入最终目标容差带并保持指定时间后的时刻。'; ... '从 t_sw_start 开始计算，体现整个任务启动后的推力到位速度；Seq-I 会天然包含先切换后跟踪的等待时间。'; ... '从各方法自身推力目标启用时刻计算，主要用于比较推力通道自身响应能力，避免 Seq-I 的天然等待时间干扰。'; ... '采用 max(几何最终目标完成时间, 推力进入目标误差带时间)，不再使用只表示曲线变平的 FNStable 时间。'; ... '与上一项相同，但从各方法自身推力目标启用时刻计时，用于辅助理解串行方法。'; ... '在推力目标启用后的时间段内计算 FN 与 FN_target 的 RMSE。'; ... '用于检查是否触碰或超过 T4 上限。'; ... '用于检查 FAN/CDFS/HPC 喘振裕度是否低于安全下限。'; ... '用于检查 FAR 是否低于设定的设计点油气比下限。'; ... '用于判断某种方法是否依赖大量不可行回退。'}); tbl = table(MetricName, ChineseName, Meaning_CN); end function print_clm_key_metrics(metricC, serialCases, tolRel) metrics = metricC; labels = {'Concurrent'}; for ii = 1:numel(serialCases) if isfield(serialCases(ii), 'metric_clm') metrics(end+1) = serialCases(ii).metric_clm; %#ok labels{end+1} = char(serialCases(ii).plotName); %#ok end end fprintf('\n================ CLM key metrics (%.2f%% band) ================\n', tolRel*100); fprintf(' %-12s %10s %12s %12s %12s %12s %12s\n', ... 'Method', 'FN_settle', 'All_stable', 'Complete', 'FN_RMSE', 'FN_MaxErr', 'T4_max'); fprintf(' %-12s %10s %12s %12s %12s %12s %12s\n', ... '', 'from sw/s', 'from sw/s', 'from sw/s', 'N', 'N', 'K'); for ii = 1:numel(metrics) fprintf(' %-12s %10.3f %12.3f %12.3f %12.4f %12.4f %12.2f\n', ... labels{ii}, ... metrics(ii).FNSettleTimeFromSwitch_s, ... metrics(ii).AllParamStableFromSwitch_s, ... metrics(ii).SerialCompleteFromSwitch_s, ... metrics(ii).FNTrackRMSE_N, ... metrics(ii).FNTrackMaxAbs_N, ... metrics(ii).T4Max); end fprintf(' Constraint minima are in metricCLMTable: FAN_SM_Min / IPC_SM_Min / HPC_SM_Min.\n'); end function figList = plot_algorithm_compare(resC, resS, ... FN_target, FN_initial, t_sw_start, t_sw_end, t_fuel_start, ... T4_max_K, SM_LP_min, SM_IP_min, SM_HP_min, FAR_min_phy, FAR_min_ratio_dp, far_idx, metricC, metricS) t = resC.tVec; simTime = max(t); figList = gobjects(0); hasSerial = ~isempty(resS); clrRef = [0.85 0.33 0.10]; clrC = [0.08 0.38 0.65]; clrS = [0.80 0.25 0.18]; clrGray = [0.45 0.45 0.45]; function add_sw_patch(ax, t1, t2, tMax) return; yl = ylim(ax); patch(ax, [t1 t2 t2 t1], ... [yl(1) yl(1) yl(2) yl(2)], [0.85 0.92 0.98], ... 'EdgeColor','none','FaceAlpha',0.5,'HandleVisibility','off'); xline(ax, t1,'--','Color',[0.5 0.5 0.5],'LineWidth',0.9,... 'Label','切换起','HandleVisibility','off'); xline(ax, t2,'--','Color',[0.5 0.5 0.5],'LineWidth',0.9,... 'Label','切换止','HandleVisibility','off'); if hasSerial && abs(t_fuel_start - t2) > 1e-9 xline(ax, t_fuel_start, ':', 'Color', clrS, 'LineWidth', 1.0, ... 'Label','串行燃油响应','HandleVisibility','off'); end xlim(ax,[0 tMax]); grid(ax,'on'); end fig1 = figure('Color','w','Name','控制量对比：并发协同 vs 分段串行', ... 'Units','centimeters','Position',[2 4 18 20]); uNames = {'WFM (lb/s)', 'A\_MSV (in²)', 'A_8 (ratio)'}; uIdx = [1 2 3]; for k = 1:3 ax = subplot(3,1,k); hold(ax,'on'); plot(ax, t, resC.ur_phy(uIdx(k),:), '--', 'Color', clrRef, ... 'LineWidth',1.5, 'DisplayName','最终目标'); stairs(ax, t, resC.u_phy(uIdx(k),:), '-', 'Color', clrC, ... 'LineWidth',1.8, 'DisplayName','并发协同算法'); if hasSerial stairs(ax, t, resS.u_phy(uIdx(k),:), '-', 'Color', clrS, ... 'LineWidth',1.8, 'DisplayName','分段串行算法'); end add_sw_patch(ax, t_sw_start, t_sw_end, simTime); ylabel(ax, uNames{k}); if k == 1 title(ax,'','Visible','off'); end if k == 3 xlabel(ax,'Time (s)'); end legend(ax,'Location','best','FontSize',8); hold(ax,'off'); end figList(end+1) = fig1; fig2 = figure('Color','w','Name','喘振裕度对比：并发协同 vs 分段串行', ... 'Units','centimeters','Position',[22 4 18 18]); smNames = {'FAN\_SM (%)', 'CDFS\_SM (%)', 'HPC\_SM (%)'}; smIdx = [5 6 7]; smLimits = [SM_LP_min, SM_IP_min, SM_HP_min]; for k = 1:3 ax = subplot(3,1,k); hold(ax,'on'); stairs(ax, t, resC.y_phy(smIdx(k),:), '-', 'Color', clrC, ... 'LineWidth',1.6, 'DisplayName','并发协同算法'); if hasSerial stairs(ax, t, resS.y_phy(smIdx(k),:), '-', 'Color', clrS, ... 'LineWidth',1.6, 'DisplayName','分段串行算法'); end yline(ax, smLimits(k), 'r--', 'LineWidth',1.2, ... 'DisplayName',sprintf('下限 %.0f%%', smLimits(k))); add_sw_patch(ax, t_sw_start, t_sw_end, simTime); ylabel(ax, smNames{k}); if k == 1 title(ax,'','Visible','off'); end if k == 3 xlabel(ax,'Time (s)'); end legend(ax,'Location','best'); hold(ax,'off'); end figList(end+1) = fig2; fig3 = figure('Color','w','Name','状态与推力对比：并发协同 vs 分段串行', ... 'Units','centimeters','Position',[42 4 20 16]); ax31 = subplot(2,2,1); hold(ax31,'on'); stairs(ax31, t, resC.x_phy(1,:), '-', 'Color', clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法'); if hasSerial stairs(ax31, t, resS.x_phy(1,:), '-', 'Color', clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法'); end plot(ax31, t, resC.xr_phy(1,:), '--', 'Color',[0.7 0.7 0.7], 'LineWidth',1.2, 'DisplayName','参考'); add_sw_patch(ax31, t_sw_start, t_sw_end, simTime); ylabel(ax31,'NL (rpm)'); title(ax31,'','Visible','off'); legend(ax31,'Location','best'); hold(ax31,'off'); ax32 = subplot(2,2,2); hold(ax32,'on'); stairs(ax32, t, resC.x_phy(2,:), '-', 'Color', clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法'); if hasSerial stairs(ax32, t, resS.x_phy(2,:), '-', 'Color', clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法'); end plot(ax32, t, resC.xr_phy(2,:), '--', 'Color',[0.7 0.7 0.7], 'LineWidth',1.2, 'DisplayName','参考'); add_sw_patch(ax32, t_sw_start, t_sw_end, simTime); ylabel(ax32,'NH (rpm)'); title(ax32,'','Visible','off'); legend(ax32,'Location','best'); hold(ax32,'off'); ax33 = subplot(2,2,3); hold(ax33,'on'); stairs(ax33, t, resC.y_phy(3,:), '-', 'Color', clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法'); if hasSerial stairs(ax33, t, resS.y_phy(3,:), '-', 'Color', clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法'); end yline(ax33, T4_max_K, 'r--', 'LineWidth',1.2, 'DisplayName',sprintf('上限 %dK',T4_max_K)); add_sw_patch(ax33, t_sw_start, t_sw_end, simTime); ylabel(ax33,'T4 (K)'); xlabel(ax33,'Time (s)'); title(ax33,'','Visible','off'); legend(ax33,'Location','best'); hold(ax33,'off'); ax34 = subplot(2,2,4); hold(ax34,'on'); stairs(ax34, t, resC.y_phy(4,:), '-', 'Color', clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法'); if hasSerial stairs(ax34, t, resS.y_phy(4,:), '-', 'Color', clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法'); end plot(ax34, t, resC.yr_phy(4,:), '--', 'Color',[0.7 0.7 0.7], 'LineWidth',1.2, 'DisplayName','并发参考'); if hasSerial plot(ax34, t, resS.yr_phy(4,:), ':', 'Color',clrGray, 'LineWidth',1.3, 'DisplayName','串行参考'); end yline(ax34, FN_target, '--', 'Color',clrRef, 'LineWidth',1.1, 'DisplayName','目标推力'); yline(ax34, FN_initial, ':', 'Color',clrRef, 'LineWidth',1.1, 'DisplayName','初始推力'); if ~isnan(metricC.FNSettleAbsTime_s) xline(ax34, metricC.FNSettleAbsTime_s, '--', 'Color', clrC, ... 'LineWidth',1.0, ... 'Label',sprintf('并发 %.2fs', metricC.FNSettleTimeFromSwitch_s), ... 'HandleVisibility','off'); end if hasSerial && ~isnan(metricS.FNSettleAbsTime_s) xline(ax34, metricS.FNSettleAbsTime_s, '--', 'Color', clrS, ... 'LineWidth',1.0, ... 'Label',sprintf('串行 %.2fs', metricS.FNSettleTimeFromSwitch_s), ... 'HandleVisibility','off'); end add_sw_patch(ax34, t_sw_start, t_sw_end, simTime); ylabel(ax34,'FN (N)'); xlabel(ax34,'Time (s)'); title(ax34,'','Visible','off'); legend(ax34,'Location','best'); hold(ax34,'off'); figList(end+1) = fig3; fig4 = figure('Color','w','Name','油气比对比：并发协同 vs 分段串行', ... 'Units','centimeters','Position',[62 4 18 8]); ax41 = axes(fig4); hold(ax41,'on'); FAR_dp_phy = FAR_min_phy / FAR_min_ratio_dp; stairs(ax41, t, resC.y_phy(far_idx,:) / FAR_dp_phy * 100, '-', 'Color', clrC, 'LineWidth',1.6, 'DisplayName','Con.'); if hasSerial stairs(ax41, t, resS.y_phy(far_idx,:) / FAR_dp_phy * 100, '-', 'Color', clrS, 'LineWidth',1.6, 'DisplayName','Seq.'); end yline(ax41, 100*FAR_min_ratio_dp, 'r--', 'LineWidth',1.2, ... 'DisplayName',sprintf('%.0f%% DP 下限', 100*FAR_min_ratio_dp)); add_sw_patch(ax41, t_sw_start, t_sw_end, simTime); ylabel(ax41,'FAR (%DP)'); xlabel(ax41,'Time (s)'); title(ax41,'','Visible','off'); legend(ax41,'Location','best'); hold(ax41,'off'); figList(end+1) = fig4; end function figList = plot_lpv_clm_compare_merged(resC, clmC, metricC, metricClmC, ... serialCases, t_sw_start, t_sw_end, FN_target, FN_initial, ... T4_max_K, SM_LP_min, SM_IP_min, SM_HP_min, FAR_min_phy, far_idx, norm_u, norm_y) %#ok metricC 当前保留在接口中，便于后续扩展。 fnEn = 'Times New Roman'; fSz = 9; fSzLeg = 7.2; if evalin('base', 'exist(''plot_line_width_scale'',''var'')') lwScale = evalin('base', 'plot_line_width_scale'); else lwScale = 1.25; end lwMain = 1.35 * lwScale; lwAux = 1.15 * lwScale; clrC = [0.00 0.30 0.60]; clrRef = [0.20 0.20 0.20]; clrLim = [0.60 0.00 0.00]; AMSV_NORM = 51.2; norm_u = norm_u(:); norm_y = norm_y(:); nSerial = numel(serialCases); [uC, yC] = prepare_lpv_data(resC); yClmC = prepare_clm_data(clmC); tC = resC.tVec(:); tClmC = clmC.t_sim(:); tEnd = max([tC; tClmC]); for iCase = 1:nSerial [serialCases(iCase).uPlot, serialCases(iCase).yPlot] = ... prepare_lpv_data(serialCases(iCase).res); serialCases(iCase).yClmPlot = prepare_clm_data(serialCases(iCase).clm); serialCases(iCase).w25Metric = calc_signal_settle_metric( ... serialCases(iCase).clm.t_sim(:), serialCases(iCase).clm.w25(:), ... t_sw_start, 0.005, 0.50); tEnd = max([tEnd; serialCases(iCase).res.tVec(:); serialCases(iCase).clm.t_sim(:)]); end w25MetricC = calc_signal_settle_metric(clmC.t_sim(:), clmC.w25(:), t_sw_start, 0.005, 0.50); figList = gobjects(4,1); figList(1) = figure('Color','w','Name','LPV vs CLM - Fig1 - Control inputs', ... 'Units','centimeters','Position',[2 4 18 5.8],'PaperPositionMode','auto'); tiledlayout(1,3,'TileSpacing','compact','Padding','compact'); uLbl = {'{\it W}_{\rm f} (%)', '\beta_{\rm MSV} (%)', '{\it A}_8 (%DP)'}; for k = 1:3 ax = nexttile(k); plot_control_axis(ax, k, uLbl{k}, k == 1); end figList(2) = figure('Color','w','Name','LPV vs CLM - Fig2 - NL NH FN', ... 'Units','centimeters','Position',[22 4 18 5.8],'PaperPositionMode','auto'); tiledlayout(1,3,'TileSpacing','compact','Padding','compact'); axNL = nexttile(1); plot_response_axis(axNL, 1, '{\it N}_{\rm L} (%)', [], true, false); axNH = nexttile(2); plot_response_axis(axNH, 2, '{\it N}_{\rm H} (%)', [], false, false); axFN = nexttile(3); plot_response_axis(axFN, 4, '{\it F}_{\rm N} (%)', [], false, true); figList(3) = figure('Color','w','Name','LPV vs CLM - Fig3 - Surge margins', ... 'Units','centimeters','Position',[42 4 18 5.8],'PaperPositionMode','auto'); tiledlayout(1,3,'TileSpacing','compact','Padding','compact'); smIdx = [5 6 7]; smLbl = {'{\rm SM}_{\rm FAN} (%)', '{\rm SM}_{\rm CDFS} (%)', '{\rm SM}_{\rm HPC} (%)'}; smLim = [SM_LP_min, SM_IP_min, SM_HP_min]; for k = 1:3 ax = nexttile(k); plot_response_axis(ax, smIdx(k), smLbl{k}, smLim(k), k == 1, false); end figList(4) = figure('Color','w','Name','LPV vs CLM - Fig4 - T4 w25 FAR', ... 'Units','centimeters','Position',[62 4 18 5.8],'PaperPositionMode','auto'); tiledlayout(1,3,'TileSpacing','compact','Padding','compact'); axT4 = nexttile(1); plot_response_axis(axT4, 3, '{\it T}_4 (%)', T4_max_K / norm_y(3) * 100, true, false); axw25 = nexttile(2); plot_w25_axis(axw25); axFAR = nexttile(3); plot_response_axis(axFAR, far_idx, '{\rm FAR} (%DP)', FAR_min_phy / norm_y(far_idx) * 100, false, false); function [uPlot, yPlot] = prepare_lpv_data(res) norm_y_plot = norm_y; norm_y_plot(5:7) = 1; yPlot = bsxfun(@rdivide, res.y_phy, norm_y_plot) * 100; yPlot(5:7,:) = res.y_phy(5:7,:); yPlot(far_idx,:) = res.y_phy(far_idx,:) / norm_y(far_idx) * 100; uPlot = zeros(3, size(res.u_phy,2)); uPlot(1,:) = res.u_phy(1,:) / norm_u(1) * 100; uPlot(2,:) = res.u_phy(2,:) / AMSV_NORM * 100; uPlot(3,:) = res.u_phy(3,:) / norm_u(3) * 100; end function yPlot = prepare_clm_data(clm) yPlot = NaN(numel(norm_y), numel(clm.NL)); yPlot(1:7,:) = [ ... clm.NL(:).' / norm_y(1) * 100; ... clm.NH(:).' / norm_y(2) * 100; ... clm.T4(:).' / norm_y(3) * 100; ... clm.FN(:).' / norm_y(4) * 100; ... clm.SML(:).'; ... clm.SMI(:).'; ... clm.SMH(:).' ... ]; if isfield(clm, 'FAR') && numel(clm.FAR) == numel(clm.NL) yPlot(far_idx,:) = clm.FAR(:).' / norm_y(far_idx) * 100; end end function plot_control_axis(ax, ch, yLbl, showLeg) allDat = uC(ch,:).'; for ii = 1:nSerial allDat = [allDat; serialCases(ii).uPlot(ch,:).']; end [ylo, yhi] = calc_ylim(allDat, []); hold(ax,'on'); stairs(ax, tC, uC(ch,:).', '-', 'Color', clrC, 'LineWidth', lwMain, ... 'HandleVisibility','off'); for ii = 1:nSerial stairs(ax, serialCases(ii).res.tVec(:), serialCases(ii).uPlot(ch,:).', '-', ... 'Color', serialCases(ii).color, 'LineWidth', lwMain, ... 'HandleVisibility','off'); end apply_axis_style(ax, yLbl, ylo, yhi); if showLeg add_method_legend(ax); end hold(ax,'off'); end function plot_response_axis(ax, ch, yLbl, limVal, showLeg, isFN) allDat = [yC(ch,:).'; yClmC(ch,:).']; for ii = 1:nSerial allDat = [allDat; serialCases(ii).yPlot(ch,:).'; serialCases(ii).yClmPlot(ch,:).']; end if isFN allDat = [allDat; FN_target / norm_y(4) * 100; FN_initial / norm_y(4) * 100]; end [ylo, yhi] = calc_ylim(allDat, limVal); hold(ax,'on'); plot(ax, tClmC, yClmC(ch,:).', '-', 'Color', clrC, 'LineWidth', lwMain, ... 'HandleVisibility','off'); plot(ax, tC, yC(ch,:).', '--', 'Color', clrC, 'LineWidth', lwAux, ... 'HandleVisibility','off'); for ii = 1:nSerial plot(ax, serialCases(ii).clm.t_sim(:), serialCases(ii).yClmPlot(ch,:).', '-', ... 'Color', serialCases(ii).color, 'LineWidth', lwMain, ... 'HandleVisibility','off'); plot(ax, serialCases(ii).res.tVec(:), serialCases(ii).yPlot(ch,:).', '--', ... 'Color', serialCases(ii).color, 'LineWidth', lwAux, ... 'HandleVisibility','off'); end if ~isempty(limVal) yline(ax, limVal, ':', 'Color', clrLim, 'LineWidth', 1.0, 'HandleVisibility','off'); end if isFN yline(ax, FN_target / norm_y(4) * 100, ':', 'Color', clrRef, 'LineWidth', 0.9, 'HandleVisibility','off'); yline(ax, FN_initial / norm_y(4) * 100, ':', 'Color', [0.55 0.55 0.55], 'LineWidth', 0.9, 'HandleVisibility','off'); annotate_fn_time(ax); end annotate_consistency(ax, ch); apply_axis_style(ax, yLbl, ylo, yhi); if showLeg add_response_legend(ax); end hold(ax,'off'); end function plot_w25_axis(ax) allDat = clmC.w25(:); for ii = 1:nSerial allDat = [allDat; serialCases(ii).clm.w25(:)]; end [ylo, yhi] = calc_ylim(allDat, []); hold(ax,'on'); plot(ax, clmC.t_sim(:), clmC.w25(:), '-', 'Color', clrC, 'LineWidth', lwMain, 'HandleVisibility','off'); for ii = 1:nSerial plot(ax, serialCases(ii).clm.t_sim(:), serialCases(ii).clm.w25(:), '-', ... 'Color', serialCases(ii).color, 'LineWidth', lwMain, ... 'HandleVisibility','off'); end annotate_w25_time(ax); apply_axis_style(ax, 'w25', ylo, yhi); hold(ax,'off'); end function annotate_fn_time(ax) txt = sprintf('F_N settle: C %.2fs', metricClmC.FNSettleTimeFromSwitch_s); for ii = 1:nSerial txt = sprintf('%s\n%s %.2fs', txt, serialCases(ii).plotName, serialCases(ii).metric_clm.FNSettleTimeFromSwitch_s); end text(ax, 0.02, 0.04, txt, 'Units','normalized', 'HorizontalAlignment','left', 'VerticalAlignment','bottom', ... 'FontName',fnEn,'FontSize',fSzLeg,'Color',[0.12 0.12 0.12], 'BackgroundColor','w', ... 'Margin',2,'EdgeColor',[0.82 0.82 0.82], 'Interpreter','none'); if ~isnan(metricClmC.FNSettleAbsTime_s) xline(ax, metricClmC.FNSettleAbsTime_s, '--', 'Color', clrC, 'LineWidth', 0.9, 'HandleVisibility','off'); end for ii = 1:nSerial if ~isnan(serialCases(ii).metric_clm.FNSettleAbsTime_s) xline(ax, serialCases(ii).metric_clm.FNSettleAbsTime_s, '--', 'Color', serialCases(ii).color, 'LineWidth', 0.9, 'HandleVisibility','off'); end end end function annotate_w25_time(ax) txt = sprintf('w25 settle: C %.2fs', w25MetricC.settleFromSwitch_s); if isfinite(w25MetricC.settleAbsTime_s) xline(ax, w25MetricC.settleAbsTime_s, '--', 'Color', clrC, 'LineWidth', 0.9, 'HandleVisibility','off'); end for ii = 1:nSerial txt = sprintf('%s\n%s %.2fs', txt, serialCases(ii).plotName, serialCases(ii).w25Metric.settleFromSwitch_s); if isfinite(serialCases(ii).w25Metric.settleAbsTime_s) xline(ax, serialCases(ii).w25Metric.settleAbsTime_s, '--', 'Color', serialCases(ii).color, 'LineWidth', 0.9, 'HandleVisibility','off'); end end text(ax, 0.02, 0.04, txt, 'Units','normalized', 'HorizontalAlignment','left', 'VerticalAlignment','bottom', ... 'FontName',fnEn,'FontSize',fSzLeg,'Color',[0.12 0.12 0.12], 'BackgroundColor','w', ... 'Margin',2,'EdgeColor',[0.82 0.82 0.82], 'Interpreter','none'); end function annotate_consistency(ax, ch) s = calc_fit_summary(resC.tVec(:), yC(ch,:).', clmC.t_sim(:), yClmC(ch,:).'); txt = sprintf('NRMSE: C %.1f%%', s.nrmsePct); for ii = 1:nSerial ss = calc_fit_summary(serialCases(ii).res.tVec(:), serialCases(ii).yPlot(ch,:).', ... serialCases(ii).clm.t_sim(:), serialCases(ii).yClmPlot(ch,:).'); txt = sprintf('%s, %s %.1f%%', txt, serialCases(ii).plotName, ss.nrmsePct); end text(ax, 0.02, 0.96, txt, 'Units','normalized', 'HorizontalAlignment','left', 'VerticalAlignment','top', ... 'FontName',fnEn,'FontSize',fSzLeg,'Color',[0.20 0.20 0.20], 'BackgroundColor','w', ... 'Margin',1.5,'EdgeColor','none', 'Interpreter','none'); end function add_method_legend(ax) h = gobjects(nSerial+1,1); names = strings(nSerial+1,1); h(1) = plot(ax, NaN, NaN, '-', 'Color', clrC, 'LineWidth', lwMain); names(1) = "Con."; for jj = 1:nSerial h(jj+1) = plot(ax, NaN, NaN, '-', 'Color', serialCases(jj).color, 'LineWidth', lwMain); names(jj+1) = string(serialCases(jj).plotName); end lg = legend(ax, h, cellstr(names), 'Location','best','FontSize',fSzLeg,'Box','off'); lg.FontName = fnEn; end function add_response_legend(ax) nMethod = nSerial + 1; h = gobjects(nMethod+2,1); names = strings(nMethod+2,1); h(1) = plot(ax, NaN, NaN, '-', 'Color', clrC, 'LineWidth', lwMain); names(1) = "Con."; for jj = 1:nSerial h(jj+1) = plot(ax, NaN, NaN, '-', 'Color', serialCases(jj).color, 'LineWidth', lwMain); names(jj+1) = string(serialCases(jj).plotName); end h(nMethod+1) = plot(ax, NaN, NaN, '-', 'Color', [0 0 0], 'LineWidth', lwMain); names(nMethod+1) = "CLM"; h(nMethod+2) = plot(ax, NaN, NaN, '--', 'Color', [0 0 0], 'LineWidth', lwAux); names(nMethod+2) = "LPV"; lg = legend(ax, h, cellstr(names), 'Location','best','FontSize',fSzLeg,'Box','off'); lg.FontName = fnEn; end function apply_axis_style(ax, yLbl, ylo, yhi) xlim(ax, [0 tEnd]); ylim(ax, [ylo yhi]); grid(ax,'on'); pbaspect(ax, [1 1 1]); set(ax, 'FontName',fnEn,'FontSize',fSz,'Box','on', 'TickDir','in', 'TickLength',[0.012 0.012], ... 'GridAlpha',0.12,'GridLineStyle',':', 'XMinorTick','on','YMinorTick','on','LineWidth',0.75); set(ax, 'LooseInset', max(get(ax,'TightInset'), 0.01)); xlabel(ax, 'Time (s)', 'FontName',fnEn,'FontSize',fSz,'Interpreter','tex'); ylabel(ax, yLbl, 'FontName',fnEn,'FontSize',fSz,'Interpreter','tex'); title(ax, '', 'Visible','off'); legend(ax,'off'); end function [ylo, yhi] = calc_ylim(yData, limVal) yData = yData(isfinite(yData)); if isempty(yData) yData = [0; 1]; end ylo = min(yData); yhi = max(yData); if ~isempty(limVal) && isfinite(limVal) ylo = min(ylo, limVal); yhi = max(yhi, limVal); end yrng = yhi - ylo; if yrng < eps yrng = max(abs(yhi), 1) * 0.08; end ylo = ylo - 0.12 * yrng; yhi = yhi + 0.16 * yrng; end end function [tStable, durStable] = calc_fn_stable_time(t, FN, t_start) t = t(:); FN = FN(:); win = 20; idxStart = find(t >= t_start, 1, 'first'); tStable = NaN; durStable = NaN; if isempty(idxStart) || numel(FN) < win return; end fnRange = max(FN) - min(FN); if fnRange < eps fnRange = max(abs(mean(FN)), 1); end lastIdx = numel(FN) - win + 1; for ii = idxStart:lastIdx seg = FN(ii:ii+win-1); if max(abs(diff(seg))) <= 0.005 * fnRange tStable = t(ii); durStable = tStable - t_start; return; end end end function metric = calc_signal_settle_metric(t, y, t_start, tolRel, holdTime) t = t(:); y = y(:); metric = struct('finalValue', NaN, 'settleAbsTime_s', NaN, 'settleFromSwitch_s', NaN); if isempty(t) || isempty(y) return; end holdSteps = max(1, round(holdTime / max(median(diff(t)), eps))); idxStart = find(t >= t_start, 1, 'first'); if isempty(idxStart) idxStart = 1; end tailN = min(max(holdSteps, 20), numel(y)); yFinal = mean(y(end-tailN+1:end), 'omitnan'); tolAbs = max(tolRel * max(abs(yFinal), eps), 1e-8); okFlag = abs(y - yFinal) <= tolAbs; tSettle = first_hold_time(t, okFlag, idxStart, holdSteps); metric.finalValue = yFinal; metric.settleAbsTime_s = tSettle; if isfinite(tSettle) metric.settleFromSwitch_s = tSettle - t_start; end end function s = calc_fit_summary(tLPV, yLPV, tCLM, yCLM) yLPV_on_CLM = interp1(tLPV(:), yLPV(:), tCLM(:), 'linear', 'extrap'); yCLM = yCLM(:); n = min(numel(yLPV_on_CLM), numel(yCLM)); yLPV_on_CLM = yLPV_on_CLM(1:n); yCLM = yCLM(1:n); err = yCLM - yLPV_on_CLM; s.rmse = sqrt(mean(err.^2)); s.maxAbs = max(abs(err)); span = max(yCLM) - min(yCLM); if span < eps span = max(abs(mean(yCLM)), 1); end s.nrmsePct = s.rmse / span * 100; end function clm = run_clm_compare_case(caseName, u_WFM_case, u_AMSV_case, u_A8_case, opt) fprintf('\n[CLM] 开始运行 Simulink：%s\n', caseName); slxName = opt.slxName; u_WFM_ts = u_WFM_case; u_AMSV_ts = u_AMSV_case; u_A8_ts = u_A8_case; assignin('base', 'u_WFM_ts', u_WFM_ts); assignin('base', 'u_AMSV_ts', u_AMSV_ts); assignin('base', 'u_A8_ts', u_A8_ts); injectSteady = ~isfield(opt, 'inject_simulink_steady') || opt.inject_simulink_steady; if injectSteady assignin('base', 'NL_init_phy', opt.NL_init_phy); assignin('base', 'NH_init_phy', opt.NH_init_phy); assignin('base', 'WFM_init_phy', opt.WFM_init_phy); assignin('base', 'AMSV_init_phy', opt.AMSV_init_phy); assignin('base', 'A8_init_phy', opt.A8_init_phy); end fprintf('[CLM-%s] 注入初值：NL=%.3f rpm, NH=%.3f rpm, WFM=%.5f lb/s, A_MSV=%.3f in^2, A8=%.5f\n', ... caseName, opt.NL_init_phy, opt.NH_init_phy, ... opt.WFM_init_phy, opt.AMSV_init_phy, opt.A8_init_phy); Engine_local = opt.Engine; if injectSteady && isfield(opt, 'NRICDyn_init') && ~isempty(opt.NRICDyn_init) NRICDyn1 = opt.NRICDyn_init; assignin('base', 'NRICDyn1', NRICDyn1); if isfield(Engine_local, 'MWS') && isfield(Engine_local.MWS, 'Input') Engine_local.MWS.Input.NRICDyn1 = NRICDyn1; assignin('base', 'Engine', Engine_local); end fprintf('[CLM-%s] NRICDyn1 已注入，元素个数 = %d\n', caseName, numel(NRICDyn1)); else fprintf('[CLM-%s] 警告：NRICDyn1 未注入，Simulink 初始瞬态可能不完全稳态。\n', caseName); assignin('base', 'Engine', Engine_local); end try load_system(slxName); try if injectSteady set_param(opt.NL_block_path, 'InitialCondition', num2str(opt.NL_init_phy, '%.10f')); set_param(opt.NH_block_path, 'InitialCondition', num2str(opt.NH_init_phy, '%.10f')); end if injectSteady fprintf('[CLM-%s] 积分器初值 set_param 成功。\n', caseName); else fprintf('[CLM-%s] inject_simulink_steady=false，跳过积分器 InitialCondition 注入。\n', caseName); end catch ME_ic fprintf('[CLM-%s] 警告：积分器 set_param 失败，改用工作区变量方式。\n', caseName); fprintf(' 错误信息：%s\n', ME_ic.message); fprintf(' 请确认积分器 IC 是否填写 NL_init_phy / NH_init_phy。\n'); end set_param(slxName, 'StopTime', num2str(opt.simTime)); set_param(slxName, 'FixedStep', num2str(opt.Ts)); set_param(slxName, 'SolverType', 'Fixed-step'); fprintf('[CLM-%s] 开始 sim(%s)，StopTime=%.3f, Ts=%.5f ...\n', ... caseName, slxName, opt.simTime, opt.Ts); simOut = sim(slxName, 'StopTime', num2str(opt.simTime)); try if ~isempty(simOut.ErrorMessage) error('[CLM-%s] Simulink 仿真失败：%s', caseName, simOut.ErrorMessage); end catch end fprintf('[CLM-%s] Simulink 仿真完成。\n', caseName); catch ME fprintf('[CLM-%s] Simulink 启动/运行失败：%s\n', caseName, ME.message); rethrow(ME); end clm = struct(); clm.caseName = caseName; try clm.t_sim = simOut.tout(:); clm.NL = extract_sim_signal(simOut, 'NL'); clm.NH = extract_sim_signal(simOut, 'NH'); clm.T4 = extract_sim_signal(simOut, 'T4'); clm.FN = extract_sim_signal(simOut, 'FN'); clm.w25 = extract_sim_signal(simOut, 'W25'); clm.SML = extract_sim_signal(simOut, 'SML'); clm.SMI = extract_sim_signal(simOut, 'SMI'); clm.SMH = extract_sim_signal(simOut, 'SMH'); try clm.FAR = extract_sim_signal(simOut, 'FAR'); catch clm.FAR = NaN(size(clm.NL)); end Nout = min([ ... numel(clm.t_sim), numel(clm.NL), numel(clm.NH), numel(clm.T4), ... numel(clm.FN), numel(clm.w25), numel(clm.SML), numel(clm.SMI), numel(clm.SMH), numel(clm.FAR)]); clm.t_sim = clm.t_sim(1:Nout); clm.NL = clm.NL(1:Nout); clm.NH = clm.NH(1:Nout); clm.T4 = clm.T4(1:Nout); clm.FN = clm.FN(1:Nout); clm.w25 = clm.w25(1:Nout); clm.SML = clm.SML(1:Nout); clm.SMI = clm.SMI(1:Nout); clm.SMH = clm.SMH(1:Nout); clm.FAR = clm.FAR(1:Nout); clm.y_phy = [ ... clm.NL(:).'; ... clm.NH(:).'; ... clm.T4(:).'; ... clm.FN(:).'; ... clm.SML(:).'; ... clm.SMI(:).'; ... clm.SMH(:).'; ... clm.FAR(:).' ... ]; fprintf('[CLM-%s] 输出提取成功，共 %d 个采样点。\n', caseName, Nout); catch ME_extract fprintf('[CLM-%s] 输出提取失败：%s\n', caseName, ME_extract.message); fprintf('请确认 To Workspace 变量名为 NL/NH/T4/FN/w25/SML/SMI/SMH，FAR 可选。\n'); rethrow(ME_extract); end end function v = extract_sim_signal(simOut, sigName) raw = simOut.get(sigName); if isa(raw, 'timeseries') v = raw.Data; elseif isa(raw, 'Simulink.SimulationData.Signal') v = raw.Values.Data; elseif isstruct(raw) && isfield(raw, 'signals') && isfield(raw.signals, 'values') % 兼容 To Workspace 的 Structure With Time / Structure 格式。 v = raw.signals.values; else v = raw; end v = squeeze(double(v)); if ismatrix(v) && size(v,2) > 1 && size(v,1) >= 1 % 若用户误用 Array 格式导出多列信号，只取第一列，避免直接 v(:) 将多列拼接成伪长序列。 v = v(:,1); end v = v(:); end function tbl = print_lpv_clm_consistency(caseName, resLPV, clm, t_eval_start) fprintf('\n========= LPV vs CLM 一致性速查：%s =========\n', caseName); chNames = {'NL'; 'NH'; 'T4'; 'FN'; 'FAN_SM'; 'IPC_SM'; 'HPC_SM'; 'FAR'}; units = {'rpm'; 'rpm'; 'K'; 'N'; '%'; '%'; '%'; '-'}; nCh = min(numel(chNames), size(resLPV.y_phy,1)); if isfield(clm, 'y_phy') nCh = min(nCh, size(clm.y_phy,1)); end chNames = chNames(1:nCh); units = units(1:nCh); tailN = 50; LPV_TailMean = zeros(nCh,1); CLM_TailMean = zeros(nCh,1); TailDiff = zeros(nCh,1); RMSE_After = zeros(nCh,1); MaxAbs_After = zeros(nCh,1); NRMSE_AfterPct = zeros(nCh,1); t_lpv = resLPV.tVec(:); t_clm = clm.t_sim(:); fprintf(' %-10s %-8s %14s %14s %12s %12s %10s\n', ... '通道', '单位', 'LPV尾均值', 'CLM尾均值', '尾均差', 'MaxAbs', 'NRMSE%'); for ch = 1:nCh y_lpv_raw = resLPV.y_phy(ch,:).'; y_clm_raw = clm.y_phy(ch,:).'; y_lpv_on_clm = interp1(t_lpv, y_lpv_raw, t_clm, 'linear', 'extrap'); N = min(numel(y_lpv_on_clm), numel(y_clm_raw)); y_lpv_on_clm = y_lpv_on_clm(1:N); y_clm_raw = y_clm_raw(1:N); t_use = t_clm(1:N); tailIdx = max(1, N-tailN+1):N; evalIdx = find(t_use >= t_eval_start); if isempty(evalIdx) evalIdx = 1:N; end err = y_clm_raw - y_lpv_on_clm; LPV_TailMean(ch) = mean(y_lpv_on_clm(tailIdx)); CLM_TailMean(ch) = mean(y_clm_raw(tailIdx)); TailDiff(ch) = CLM_TailMean(ch) - LPV_TailMean(ch); RMSE_After(ch) = sqrt(mean(err(evalIdx).^2)); MaxAbs_After(ch) = max(abs(err(evalIdx))); evalSpan = max(y_clm_raw(evalIdx)) - min(y_clm_raw(evalIdx)); if evalSpan < eps evalSpan = max(abs(mean(y_clm_raw(evalIdx))), 1); end NRMSE_AfterPct(ch) = RMSE_After(ch) / evalSpan * 100; fprintf(' %-10s %-8s %14.4f %14.4f %+12.4f %12.4f %10.2f\n', ... chNames{ch}, units{ch}, ... LPV_TailMean(ch), CLM_TailMean(ch), TailDiff(ch), ... MaxAbs_After(ch), NRMSE_AfterPct(ch)); end fprintf('============================================================\n\n'); tbl = table( ... string(chNames), string(units), ... LPV_TailMean, CLM_TailMean, TailDiff, RMSE_After, MaxAbs_After, NRMSE_AfterPct, ... 'VariableNames', {'Channel','Unit','LPV_TailMean','CLM_TailMean','TailDiff','RMSE_AfterSwitch','MaxAbs_AfterSwitch','NRMSE_AfterSwitchPct'}); end function m = calc_clm_metrics(caseName, clm, opt, t_fuel_start) t = clm.t_sim(:); FN = clm.FN(:); T4 = clm.T4(:); SM = [clm.SML(:), clm.SMI(:), clm.SMH(:)]; holdSteps = max(1, round(opt.hold_time / opt.Ts)); idxStart = find(t >= opt.t_start, 1, 'first'); FN_tol_abs = opt.FN_tol_rel * abs(opt.FNtgt); tFNSettleAbs = first_hold_time(t, abs(FN - opt.FNtgt) <= FN_tol_abs, idxStart, holdSteps); [tFNStableAbs, ~] = calc_fn_stable_time(t, FN, opt.t_start); if isfield(opt, 't_geom_done_abs') && isfinite(opt.t_geom_done_abs) tGeomAbs = opt.t_geom_done_abs; elseif isfield(opt, 't_geom_end') && isfinite(opt.t_geom_end) tGeomAbs = opt.t_geom_end; else tGeomAbs = NaN; end tSerialCompleteAbs = max_time_ignore_nan(tGeomAbs, tFNSettleAbs); deltaFN = opt.FNtgt - opt.FN0; if abs(deltaFN) < eps tFN90Abs = NaN; else FN90 = opt.FN0 + 0.9*deltaFN; if deltaFN > 0 idx90 = find(t >= opt.t_start & FN >= FN90, 1, 'first'); else idx90 = find(t >= opt.t_start & FN <= FN90, 1, 'first'); end if isempty(idx90) tFN90Abs = NaN; else tFN90Abs = t(idx90); end end idxStage1 = t >= opt.t_start & t <= opt.t_geom_end; if any(idxStage1) FNStage1FluctAbs = max(abs(FN(idxStage1) - opt.FN0)); else FNStage1FluctAbs = NaN; end signDelta = sign(deltaFN); if signDelta == 0 signDelta = 1; end idxAfterStart = t >= opt.t_start; overshootAbs = max([0; signDelta*(FN(idxAfterStart) - opt.FNtgt)]); undershootAbs = max([0; -signDelta*(FN(idxAfterStart) - opt.FN0)]); m = struct(); m.Method = string(caseName); m.MetricNote_CN = string('任务总完成时间=几何到达最终目标且推力进入目标误差带后的较晚时刻；推力自身响应时间另按从任务开始/从燃油启用分别统计。'); m.GeomSettleAbsTime_s = tGeomAbs; m.GeomSettleFromSwitch_s = subtract_or_nan(tGeomAbs, opt.t_start); m.FN90AbsTime_s = tFN90Abs; m.FN90FromSwitch_s = subtract_or_nan(tFN90Abs, opt.t_start); m.FNSettleAbsTime_s = tFNSettleAbs; m.FNSettleTimeFromSwitch_s = subtract_or_nan(tFNSettleAbs, opt.t_start); m.FNSettleTimeFromFuel_s = subtract_or_nan(tFNSettleAbs, t_fuel_start); m.FNStableAbsTime_s = tFNStableAbs; m.FNStableFromSwitch_s = subtract_or_nan(tFNStableAbs, opt.t_start); m.FNStableFromFuel_s = subtract_or_nan(tFNStableAbs, t_fuel_start); m.SerialCompleteAbsTime_s = tSerialCompleteAbs; m.SerialCompleteFromSwitch_s = subtract_or_nan(tSerialCompleteAbs, opt.t_start); m.SerialCompleteFromFuel_s = subtract_or_nan(tSerialCompleteAbs, t_fuel_start); m.Stage1FNFluctAbs = FNStage1FluctAbs; m.Stage1FNFluctRelPct = FNStage1FluctAbs / max(abs(opt.FN0), eps) * 100; m.FNOvershootAbs = overshootAbs; m.FNOvershootRelPct = overshootAbs / max(abs(deltaFN), eps) * 100; m.FNUndershootAbs = undershootAbs; Yall = [clm.NL(:), clm.NH(:), clm.T4(:), clm.FN(:), clm.SML(:), clm.SMI(:), clm.SMH(:)]; [tAllStableAbs, maxFinalRelPct] = calc_all_signal_stable_time( ... t, Yall, opt.t_start, opt.FN_tol_rel, holdSteps); idxTrack = t >= t_fuel_start; if any(idxTrack) errFNTrack = FN(idxTrack) - opt.FNtgt; m.FNTrackMAE_N = mean(abs(errFNTrack)); m.FNTrackRMSE_N = sqrt(mean(errFNTrack.^2)); m.FNTrackMaxAbs_N = max(abs(errFNTrack)); m.FNTrackRMSE_RelPct = m.FNTrackRMSE_N / max(abs(opt.FNtgt), eps) * 100; else m.FNTrackMAE_N = NaN; m.FNTrackRMSE_N = NaN; m.FNTrackMaxAbs_N = NaN; m.FNTrackRMSE_RelPct = NaN; end m.AllParamStableAbsTime_s = tAllStableAbs; m.AllParamStableFromSwitch_s = subtract_or_nan(tAllStableAbs, opt.t_start); m.AllParamMaxFinalRelPct = maxFinalRelPct; m.T4Max = max(T4); m.FAN_SM_Min = min(SM(:,1)); m.IPC_SM_Min = min(SM(:,2)); m.HPC_SM_Min = min(SM(:,3)); if isfield(clm, 'FAR') m.FAR_Min = min(clm.FAR(:)); else m.FAR_Min = NaN; end m.QPFail = NaN; end function figCLM = plot_clm_method_compare(clmC, clmS, ... FN_target, FN_initial, t_sw_start, t_sw_end, t_fuel_start, ... T4_max_K, SM_LP_min, SM_IP_min, SM_HP_min, metricC, metricS) tC = clmC.t_sim(:); tS = clmS.t_sim(:); simTime = max([tC; tS]); clrC = [0.08 0.38 0.65]; clrS = [0.80 0.25 0.18]; clrRef = [0.85 0.33 0.10]; function add_sw_patch_local(ax) return; yl = ylim(ax); if t_sw_end > t_sw_start patch(ax, [t_sw_start t_sw_end t_sw_end t_sw_start], ... [yl(1) yl(1) yl(2) yl(2)], [0.85 0.92 0.98], ... 'EdgeColor','none','FaceAlpha',0.5,'HandleVisibility','off'); end xline(ax, t_sw_start, '--', 'Color',[0.5 0.5 0.5], ... 'LineWidth',0.9,'Label','切换起','HandleVisibility','off'); xline(ax, t_sw_end, '--', 'Color',[0.5 0.5 0.5], ... 'LineWidth',0.9,'Label','切换止','HandleVisibility','off'); if abs(t_fuel_start - t_sw_end) > 1e-9 xline(ax, t_fuel_start, ':', 'Color',clrS, ... 'LineWidth',1.0,'Label','串行燃油响应','HandleVisibility','off'); end xlim(ax, [0 simTime]); grid(ax,'on'); end figCLM = figure('Color','w','Name','CLM状态与推力对比：并发协同 vs 分段串行', ... 'Units','centimeters','Position',[8 4 22 17]); ax1 = subplot(2,2,1); hold(ax1,'on'); stairs(ax1, tC, clmC.NL, '-', 'Color',clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法-CLM'); stairs(ax1, tS, clmS.NL, '-', 'Color',clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法-CLM'); add_sw_patch_local(ax1); ylabel(ax1,'NL (rpm)'); title(ax1,'','Visible','off'); legend(ax1,'Location','best'); hold(ax1,'off'); ax2 = subplot(2,2,2); hold(ax2,'on'); stairs(ax2, tC, clmC.NH, '-', 'Color',clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法-CLM'); stairs(ax2, tS, clmS.NH, '-', 'Color',clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法-CLM'); add_sw_patch_local(ax2); ylabel(ax2,'NH (rpm)'); title(ax2,'','Visible','off'); legend(ax2,'Location','best'); hold(ax2,'off'); ax3 = subplot(2,2,3); hold(ax3,'on'); stairs(ax3, tC, clmC.T4, '-', 'Color',clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法-CLM'); stairs(ax3, tS, clmS.T4, '-', 'Color',clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法-CLM'); yline(ax3, T4_max_K, 'r--', 'LineWidth',1.2, 'DisplayName',sprintf('上限 %dK', T4_max_K)); add_sw_patch_local(ax3); ylabel(ax3,'T4 (K)'); xlabel(ax3,'Time (s)'); title(ax3,'','Visible','off'); legend(ax3,'Location','best'); hold(ax3,'off'); ax4 = subplot(2,2,4); hold(ax4,'on'); stairs(ax4, tC, clmC.FN, '-', 'Color',clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法-CLM'); stairs(ax4, tS, clmS.FN, '-', 'Color',clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法-CLM'); yline(ax4, FN_target, '--', 'Color',clrRef, 'LineWidth',1.1, 'DisplayName','目标推力'); yline(ax4, FN_initial, ':', 'Color',clrRef, 'LineWidth',1.1, 'DisplayName','初始推力'); if ~isnan(metricC.FNSettleAbsTime_s) xline(ax4, metricC.FNSettleAbsTime_s, '--', 'Color',clrC, ... 'LineWidth',1.0, ... 'Label',sprintf('并发 %.2fs', metricC.FNSettleTimeFromSwitch_s), ... 'HandleVisibility','off'); end if ~isnan(metricS.FNSettleAbsTime_s) xline(ax4, metricS.FNSettleAbsTime_s, '--', 'Color',clrS, ... 'LineWidth',1.0, ... 'Label',sprintf('串行 %.2fs', metricS.FNSettleTimeFromSwitch_s), ... 'HandleVisibility','off'); end add_sw_patch_local(ax4); ylabel(ax4,'FN (N)'); xlabel(ax4,'Time (s)'); title(ax4,'','Visible','off'); legend(ax4,'Location','best'); hold(ax4,'off'); figure('Color','w','Name','CLM喘振裕度对比：并发协同 vs 分段串行', ... 'Units','centimeters','Position',[32 4 18 18]); smName = {'FAN\_SM (%)', 'CDFS\_SM (%)', 'HPC\_SM (%)'}; smLim = [SM_LP_min, SM_IP_min, SM_HP_min]; smC = {clmC.SML, clmC.SMI, clmC.SMH}; smS = {clmS.SML, clmS.SMI, clmS.SMH}; for k = 1:3 ax = subplot(3,1,k); hold(ax,'on'); stairs(ax, tC, smC{k}, '-', 'Color',clrC, 'LineWidth',1.6, 'DisplayName','并发协同算法-CLM'); stairs(ax, tS, smS{k}, '-', 'Color',clrS, 'LineWidth',1.6, 'DisplayName','分段串行算法-CLM'); yline(ax, smLim(k), 'r--', 'LineWidth',1.2, ... 'DisplayName',sprintf('下限 %.0f%%', smLim(k))); add_sw_patch_local(ax); ylabel(ax, smName{k}); if k == 1 title(ax,'','Visible','off'); end if k == 3 xlabel(ax,'Time (s)'); end legend(ax,'Location','best'); hold(ax,'off'); end end