MATLAB自适应子空间辨识工具箱，能够通过输入输出数据计算状态空间方程，并进行系统建模和预测
1. 自适应子空间辨识理论基础
classdef AdaptiveSubspaceID% 自适应子空间辨识类properties% 系统参数nx % 状态维数nu % 输入维数ny % 输出维数horizon % 预测时域% 辨识参数window_size % 滑动窗口大小forgetting_factor % 遗忘因子regularization % 正则化参数% 系统矩阵A, B, C, D % 状态空间矩阵K % 卡尔曼增益% 数据存储input_buffer % 输入缓冲区output_buffer % 输出缓冲区state_buffer % 状态缓冲区% 协方差矩阵P % 状态估计协方差Q, R % 过程噪声和测量噪声协方差% 性能指标identification_error % 辨识误差prediction_error % 预测误差endmethodsfunction obj = AdaptiveSubspaceID(nx, nu, ny, varargin)% 构造函数obj.nx = nx;obj.nu = nu;obj.ny = ny;% 解析可选参数p = inputParser;addParameter(p, 'window_size', 100, @(x) x>0);addParameter(p, 'forgetting_factor', 0.99, @(x) x>0 && x<=1);addParameter(p, 'regularization', 1e-6, @(x) x>0);addParameter(p, 'horizon', 10, @(x) x>0);parse(p, varargin{:});obj.window_size = p.Results.window_size;obj.forgetting_factor = p.Results.forgetting_factor;obj.regularization = p.Results.regularization;obj.horizon = p.Results.horizon;% 初始化系统矩阵obj = initialize_system_matrices(obj);% 初始化数据缓冲区obj = initialize_buffers(obj);fprintf('自适应子空间辨识器初始化完成:\n');fprintf(' 状态维数: %d, 输入维数: %d, 输出维数: %d\n', nx, nu, ny);fprintf(' 窗口大小: %d, 遗忘因子: %.3f\n', obj.window_size, obj.forgetting_factor);endfunction obj = initialize_system_matrices(obj)% 初始化系统矩阵obj.A = eye(obj.nx) * 0.95; % 稳定的初始A矩阵obj.B = randn(obj.nx, obj.nu) * 0.1;obj.C = randn(obj.ny, obj.nx) * 0.1;obj.D = zeros(obj.ny, obj.nu);obj.K = zeros(obj.nx, obj.ny); % 卡尔曼增益% 初始化协方差矩阵obj.P = eye(obj.nx) * 10;obj.Q = eye(obj.nx) * 0.01;obj.R = eye(obj.ny) * 0.1;endfunction obj = initialize_buffers(obj)% 初始化数据缓冲区obj.input_buffer = zeros(obj.nu, obj.window_size);obj.output_buffer = zeros(obj.ny, obj.window_size);obj.state_buffer = zeros(obj.nx, obj.window_size);obj.identification_error = [];obj.prediction_error = [];endend
end
2. 子空间辨识核心算法
classdef SubspaceCore% 子空间辨识核心算法methods (Static)function [A, B, C, D, K] = nasid(Y, U, nx, horizon, varargin)% NASID (Numerical Algorithms for Subspace State Space System IDentification)% 基于输入输出数据的子空间辨识p = inputParser;addParameter(p, 'regularization', 1e-6, @(x) x>0);addParameter(p, 'stable', true, @islogical);parse(p, varargin{:});[N, ny] = size(Y);[~, nu] = size(U);% 构建Hankel矩阵H = SubspaceCore.build_combined_hankel(Y, U, horizon);% LQ分解[L, ~] = lq(H);% 提取状态序列X = SubspaceCore.extract_states(L, nx, horizon, ny, nu);% 计算系统矩阵[A, B, C, D, K] = SubspaceCore.estimate_system_matrices(X, Y, U, nx);% 确保系统稳定（可选）if p.Results.stableA = SubspaceCore.stabilize_matrix(A);endendfunction H = build_combined_hankel(Y, U, horizon)% 构建组合Hankel矩阵[N, ny] = size(Y);[~, nu] = size(U);% 输出Hankel矩阵HY = zeros(horizon * ny, N - 2*horizon + 1);for i = 1:horizonHY((i-1)*ny+1:i*ny, :) = Y(i:i+N-2*horizon, :)';end% 输入Hankel矩阵HU = zeros(horizon * nu, N - 2*horizon + 1);for i = 1:horizonHU((i-1)*nu+1:i*nu, :) = U(i:i+N-2*horizon, :)';end% 组合Hankel矩阵H = [HY; HU];endfunction X = extract_states(L, nx, horizon, ny, nu)% 从L矩阵提取状态序列% L矩阵结构: [L11 0; L21 L22]L11 = L(1:horizon*ny, 1:horizon*ny);L21 = L(horizon*ny+1:end, 1:horizon*ny);% SVD分解提取状态[U, S, V] = svd(L21, 'econ');% 选择前nx个奇异值对应的状态U1 = U(:, 1:nx);S1 = S(1:nx, 1:nx);% 状态序列X = S1 * V(:, 1:nx)';endfunction [A, B, C, D, K] = estimate_system_matrices(X, Y, U, nx)% 估计系统矩阵[N, ny] = size(Y);[~, nu] = size(U);% 状态序列X_current = X(:, 1:end-1);X_next = X(:, 2:end);% 最小二乘估计A, BPhi = [X_current; U(1:end-1, :)'];Theta = X_next;AB = Theta / Phi;A = AB(:, 1:nx);B = AB(:, nx+1:end);% 最小二乘估计C, DPsi = [X_current; U(1:end-1, :)'];Gamma = Y(1:end-1, :)';CD = Gamma / Psi;C = CD(:, 1:nx);D = CD(:, nx+1:end);% 估计卡尔曼增益Kresiduals = Y(1:end-1, :)' - C * X_current - D * U(1:end-1, :)';innovation_cov = residuals * residuals' / (N-1);state_innovation_cov = X_next * residuals' / (N-1);K = state_innovation_cov / innovation_cov;endfunction A_stable = stabilize_matrix(A)% 使矩阵A稳定（特征值在单位圆内）[V, D] = eig(A);eigenvalues = diag(D);% 将不稳定的特征值映射到单位圆内unstable_idx = abs(eigenvalues) >= 1;eigenvalues(unstable_idx) = eigenvalues(unstable_idx) ./ ...(abs(eigenvalues(unstable_idx)) + 1e-6);A_stable = real(V * diag(eigenvalues) / V);endfunction [A, B, C, D] = po_moesp(Y, U, nx, horizon, varargin)% PO-MOESP (Past Output Multivariable Output Error State Space)% 基于过去输出的子空间辨识方法p = inputParser;addParameter(p, 'weighting', 'CVA', @ischar); % CVA, PC, UPCparse(p, varargin{:});[N, ny] = size(Y);[~, nu] = size(U);% 构建过去和未来的Hankel矩阵[Y_past, Y_future, U_past, U_future] = ...SubspaceCore.build_past_future_hankel(Y, U, horizon);% 正交投影W = Y_future / [Y_past; U_past; U_future];% 加权SVDswitch p.Results.weightingcase 'CVA'% 规范变量分析[Uw, Sw, Vw] = svd(W, 'econ');case 'PC'% 主成分分析[Uw, Sw, Vw] = svd(W * W', 'econ');otherwise% 无加权[Uw, Sw, Vw] = svd(W, 'econ');end% 提取状态序列U1 = Uw(:, 1:nx);S1 = Sw(1:nx, 1:nx);X = S1 * Vw(:, 1:nx)';% 估计系统矩阵[A, B, C, D, ~] = SubspaceCore.estimate_system_matrices(X, Y, U, nx);endfunction [Y_past, Y_future, U_past, U_future] = build_past_future_hankel(Y, U, horizon)% 构建过去和未来的Hankel矩阵[N, ny] = size(Y);[~, nu] = size(U);% 过去输出Hankel矩阵Y_past = zeros(horizon * ny, N - 2*horizon + 1);for i = 1:horizonY_past((i-1)*ny+1:i*ny, :) = Y(i:i+N-2*horizon, :)';end% 未来输出Hankel矩阵Y_future = zeros(horizon * ny, N - 2*horizon + 1);for i = 1:horizonY_future((i-1)*ny+1:i*ny, :) = Y(horizon+i:horizon+i+N-2*horizon, :)';end% 过去输入Hankel矩阵U_past = zeros(horizon * nu, N - 2*horizon + 1);for i = 1:horizonU_past((i-1)*nu+1:i*nu, :) = U(i:i+N-2*horizon, :)';end% 未来输入Hankel矩阵U_future = zeros(horizon * nu, N - 2*horizon + 1);for i = 1:horizonU_future((i-1)*nu+1:i*nu, :) = U(horizon+i:horizon+i+N-2*horizon, :)';endendend
end
3. 自适应更新算法
classdef AdaptiveUpdate% 自适应更新算法methods (Static)function obj = recursive_update(obj, u_new, y_new)% 递归更新系统模型% 更新数据缓冲区obj = update_data_buffer(obj, u_new, y_new);% 检查是否有足够数据if size(obj.input_buffer, 2) >= obj.window_size% 提取当前窗口数据U_window = obj.input_buffer;Y_window = obj.output_buffer;% 执行子空间辨识[A_new, B_new, C_new, D_new, K_new] = ...SubspaceCore.nasid(Y_window', U_window', obj.nx, min(10, obj.window_size-1));% 自适应融合新旧模型obj = blend_models(obj, A_new, B_new, C_new, D_new, K_new);% 更新状态估计obj = update_state_estimate(obj, u_new, y_new);% 记录性能指标obj = update_performance_metrics(obj, u_new, y_new);endendfunction obj = update_data_buffer(obj, u_new, y_new)% 更新数据缓冲区（滑动窗口）% 添加新数据obj.input_buffer = [obj.input_buffer, u_new(:)];obj.output_buffer = [obj.output_buffer, y_new(:)];% 如果超过窗口大小，移除最旧数据if size(obj.input_buffer, 2) > obj.window_sizeobj.input_buffer = obj.input_buffer(:, 2:end);obj.output_buffer = obj.output_buffer(:, 2:end);endendfunction obj = blend_models(obj, A_new, B_new, C_new, D_new, K_new)% 融合新旧模型（使用遗忘因子）alpha = 1 - obj.forgetting_factor;obj.A = (1-alpha) * obj.A + alpha * A_new;obj.B = (1-alpha) * obj.B + alpha * B_new;obj.C = (1-alpha) * obj.C + alpha * C_new;obj.D = (1-alpha) * obj.D + alpha * D_new;obj.K = (1-alpha) * obj.K + alpha * K_new;endfunction obj = update_state_estimate(obj, u_new, y_new)% 使用卡尔曼滤波更新状态估计if isempty(obj.state_buffer)% 初始化状态x_est = zeros(obj.nx, 1);else% 使用最新状态x_est = obj.state_buffer(:, end);end% 状态预测x_pred = obj.A * x_est + obj.B * u_new;y_pred = obj.C * x_pred + obj.D * u_new;% 测量更新innovation = y_new - y_pred;x_est = x_pred + obj.K * innovation;% 存储状态if isempty(obj.state_buffer)obj.state_buffer = x_est;elseobj.state_buffer = [obj.state_buffer, x_est];% 限制状态缓冲区大小if size(obj.state_buffer, 2) > obj.window_sizeobj.state_buffer = obj.state_buffer(:, 2:end);endendendfunction obj = update_performance_metrics(obj, u_new, y_new)% 更新性能指标if ~isempty(obj.state_buffer)% 当前状态x_current = obj.state_buffer(:, end);% 一步预测y_pred_one_step = obj.C * (obj.A * x_current + obj.B * u_new) + obj.D * u_new;% 计算误差identification_error = norm(y_new - y_pred_one_step);obj.identification_error = [obj.identification_error, identification_error];endendfunction [y_pred, x_pred] = predict(obj, u_sequence, x0)% 多步预测if nargin < 3if ~isempty(obj.state_buffer)x0 = obj.state_buffer(:, end);elsex0 = zeros(obj.nx, 1);endendN = size(u_sequence, 2);y_pred = zeros(obj.ny, N);x_pred = zeros(obj.nx, N);x = x0;for k = 1:Nx = obj.A * x + obj.B * u_sequence(:, k);y_pred(:, k) = obj.C * x + obj.D * u_sequence(:, k);x_pred(:, k) = x;endendfunction validation_results = validate_model(obj, U_val, Y_val)% 模型验证fprintf('正在进行模型验证...\n');[N, ~] = size(Y_val);Y_pred = zeros(size(Y_val));% 使用验证数据的前部分估计初始状态x0 = estimate_initial_state(obj, U_val(1:min(10,N),:), Y_val(1:min(10,N),:));% 多步预测for k = 1:Nif k == 1x = x0;elsex = obj.A * x + obj.B * U_val(k-1, :)';endY_pred(k, :) = (obj.C * x + obj.D * U_val(k, :)')';end% 计算性能指标validation_results = compute_validation_metrics(Y_val, Y_pred);fprintf('模型验证完成:\n');fprintf(' RMSE: %.4f\n', validation_results.rmse);fprintf(' MAPE: %.2f%%\n', validation_results.mape * 100);fprintf(' R²: %.4f\n', validation_results.r_squared);endfunction x0 = estimate_initial_state(obj, U_short, Y_short)% 估计初始状态[N_short, ~] = size(Y_short);% 使用最小二乘估计初始状态Phi = [];Y_vec = [];for k = 1:N_shortif k == 1% 对于第一步，需要初始状态continue;end% 构建回归矩阵phi_k = [];for i = 1:k-1phi_k = [phi_k; (obj.A^(k-1-i) * obj.B * U_short(i, :)')'];endphi_k = [phi_k; (obj.C * obj.A^(k-1))'];Phi = [Phi; phi_k'];Y_vec = [Y_vec; Y_short(k, :)'];end% 最小二乘估计if ~isempty(Phi)x0 = (Phi' * Phi) \ (Phi' * Y_vec);elsex0 = zeros(obj.nx, 1);endendend
endfunction results = compute_validation_metrics(Y_true, Y_pred)% 计算验证指标residuals = Y_true - Y_pred;results.rmse = sqrt(mean(residuals(:).^2));results.mae = mean(abs(residuals(:)));results.mape = mean(abs(residuals(:) ./ (Y_true(:) + eps)));results.r_squared = 1 - sum(residuals(:).^2) / sum((Y_true(:) - mean(Y_true(:))).^2);% 拟合度results.fit_percentage = 100 * (1 - norm(residuals, 'fro') / norm(Y_true - mean(Y_true(:)), 'fro'));
end
4. 完整的主程序与测试案例
function main_adaptive_subspace_identification()% 自适应子空间辨识主程序fprintf('=== MATLAB自适应子空间辨识系统 ===\n\n');% 1. 生成测试系统fprintf('步骤1: 生成测试系统...\n');[sys_true, U_train, Y_train, U_test, Y_test] = generate_test_system();% 2. 初始化自适应子空间辨识器fprintf('步骤2: 初始化辨识器...\n');nx = 3; % 假设我们知道状态维数（实际中可能需要估计）nu = size(U_train, 2);ny = size(Y_train, 2);identifier = AdaptiveSubspaceID(nx, nu, ny, ...'window_size', 50, ...'forgetting_factor', 0.98, ...'horizon', 10);% 3. 在线辨识fprintf('步骤3: 在线自适应辨识...\n');[identifier, Y_pred_online] = online_identification(identifier, U_train, Y_train);% 4. 模型验证fprintf('步骤4: 模型验证...\n');validation_results = AdaptiveUpdate.validate_model(identifier, U_test, Y_test);% 5. 多步预测fprintf('步骤5: 多步预测...\n');prediction_results = perform_prediction(identifier, U_test, Y_test);% 6. 结果可视化fprintf('步骤6: 结果可视化...\n');visualize_results(sys_true, identifier, U_train, Y_train, U_test, Y_test, ...Y_pred_online, validation_results, prediction_results);fprintf('\n=== 自适应子空间辨识完成 ===\n');
endfunction [sys_true, U_train, Y_train, U_test, Y_test] = generate_test_system()% 生成测试系统% 真实系统（3状态系统）A_true = [0.8, 0.1, -0.05;-0.2, 0.9, 0.1;0.05, -0.1, 0.85];B_true = [0.5; -0.3; 0.2];C_true = [1, 0, 0;0, 1, 0];D_true = [0; 0];sys_true = ss(A_true, B_true, C_true, D_true, 1); % 离散系统% 生成训练数据N_train = 200;t_train = (0:N_train-1)';U_train = 0.5 * sin(0.1*t_train) + 0.3 * randn(N_train, 1);% 添加不同频率的激励U_train = U_train + 0.2 * sin(0.5*t_train) + 0.1 * sin(1.5*t_train);% 模拟系统响应（加入噪声）Y_train = lsim(sys_true, U_train, t_train) + 0.05 * randn(N_train, 2);% 生成测试数据N_test = 100;t_test = (0:N_test-1)';U_test = 0.4 * cos(0.15*t_test) + 0.2 * randn(N_test, 1);U_test = U_test + 0.1 * sin(0.8*t_test);Y_test = lsim(sys_true, U_test, t_test) + 0.05 * randn(N_test, 2);fprintf('测试系统生成完成:\n');fprintf(' 训练数据: %d 样本\n', N_train);fprintf(' 测试数据: %d 样本\n', N_test);fprintf(' 真实系统特征值: '); fprintf('%.3f ', eig(A_true)); fprintf('\n');
endfunction [identifier, Y_pred_online] = online_identification(identifier, U, Y)% 在线辨识[N, ~] = size(Y);Y_pred_online = zeros(size(Y));figure('Position', [100, 100, 1200, 800]);for k = 1:N% 获取当前输入输出u_k = U(k, :)';y_k = Y(k, :)';% 自适应更新identifier = AdaptiveUpdate.recursive_update(identifier, u_k, y_k);% 记录预测值if k > 1Y_pred_online(k, :) = (identifier.C * identifier.state_buffer(:, end) + ...identifier.D * u_k)';elseY_pred_online(k, :) = y_k';end% 每50步显示一次进度if mod(k, 50) == 0fprintf(' 处理第 %d/%d 个样本...\n', k, N);% 实时显示辨识结果if k > 10subplot(2,2,1);plot(1:k, Y(1:k,1), 'b-', 'LineWidth', 1.5, 'DisplayName', '真实输出');hold on;plot(1:k, Y_pred_online(1:k,1), 'r--', 'LineWidth', 1, 'DisplayName', '预测输出');title('输出通道1 - 在线辨识');xlabel('样本索引');ylabel('幅度');legend;grid on;subplot(2,2,2);plot(identifier.identification_error, 'g-', 'LineWidth', 1.5);title('辨识误差收敛');xlabel('更新次数');ylabel('误差范数');grid on;drawnow;endendendfprintf('在线辨识完成，最终系统特征值: '); fprintf('%.3f ', eig(identifier.A)); fprintf('\n');
endfunction prediction_results = perform_prediction(identifier, U_test, Y_test)% 执行多步预测[N_test, ~] = size(Y_test);% 不同预测步长的结果prediction_horizons = [1, 5, 10, 20];prediction_errors = zeros(length(prediction_horizons), 1);figure('Position', [100, 100, 1200, 600]);for i = 1:length(prediction_horizons)horizon = prediction_horizons(i);% 多步预测Y_pred_multi = zeros(size(Y_test));for k = 1:N_test-horizon% 使用当前时刻及之前的输入进行预测if k == 1x0 = identifier.state_buffer(:, end); % 使用训练结束时的状态else% 估计当前状态x0 = estimate_current_state(identifier, U_test(1:k-1,:), Y_test(1:k-1,:));end% 未来输入序列u_future = U_test(k:min(k+horizon-1, N_test), :)';% 多步预测[y_pred, ~] = AdaptiveUpdate.predict(identifier, u_future, x0);% 记录预测结果pred_end = min(k+horizon-1, N_test);Y_pred_multi(k:pred_end, :) = y_pred(:, 1:pred_end-k+1)';end% 计算预测误差prediction_error = sqrt(mean((Y_test - Y_pred_multi).^2, 1));prediction_errors(i) = mean(prediction_error);% 绘制预测结果subplot(2, 2, i);plot(Y_test(:,1), 'b-', 'LineWidth', 2, 'DisplayName', '真实输出');hold on;plot(Y_pred_multi(:,1), 'r--', 'LineWidth', 1.5, 'DisplayName', sprintf('%d步预测', horizon));title(sprintf('%d步预测 (RMSE: %.4f)', horizon, prediction_errors(i)));xlabel('样本索引');ylabel('幅度');legend;grid on;endprediction_results.horizons = prediction_horizons;prediction_results.errors = prediction_errors;prediction_results.Y_pred_multi = Y_pred_multi;sgtitle('多步预测性能');
endfunction x_current = estimate_current_state(identifier, U_past, Y_past)% 估计当前状态[N_past, ~] = size(Y_past);if N_past == 0x_current = zeros(identifier.nx, 1);return;end% 使用观测器估计当前状态x_current = zeros(identifier.nx, 1);for k = 1:N_past% 状态预测x_pred = identifier.A * x_current + identifier.B * U_past(k, :)';% 测量更新y_pred = identifier.C * x_pred + identifier.D * U_past(k, :)';innovation = Y_past(k, :)' - y_pred;x_current = x_pred + identifier.K * innovation;end
endfunction visualize_results(sys_true, identifier, U_train, Y_train, U_test, Y_test, ...Y_pred_online, validation_results, prediction_results)% 综合可视化结果figure('Position', [50, 50, 1400, 1000]);% 子图1: 真实系统 vs 辨识系统特征值subplot(3, 3, 1);eig_true = eig(sys_true.A);eig_ident = eig(identifier.A);plot(real(eig_true), imag(eig_true), 'bo', 'MarkerSize', 10, ...'MarkerFaceColor', 'blue', 'DisplayName', '真实系统');hold on;plot(real(eig_ident), imag(eig_ident), 'rs', 'MarkerSize', 10, ...'MarkerFaceColor', 'red', 'DisplayName', '辨识系统');% 绘制单位圆theta = linspace(0, 2*pi, 100);plot(cos(theta), sin(theta), 'k--', 'LineWidth', 1, 'DisplayName', '单位圆');axis equal;xlabel('实部');ylabel('虚部');title('系统特征值比较');legend;grid on;% 子图2: 训练数据拟合subplot(3, 3, 2);plot(Y_train(:,1), 'b-', 'LineWidth', 2, 'DisplayName', '真实输出');hold on;plot(Y_pred_online(:,1), 'r--', 'LineWidth', 1.5, 'DisplayName', '在线预测');title('训练数据拟合 - 通道1');xlabel('样本索引');ylabel('幅度');legend;grid on;% 子图3: 测试数据预测subplot(3, 3, 3);plot(Y_test(:,1), 'b-', 'LineWidth', 2, 'DisplayName', '真实输出');hold on;[Y_pred_test, ~] = AdaptiveUpdate.predict(identifier, U_test', zeros(identifier.nx, 1));plot(Y_pred_test(1,:)', 'r--', 'LineWidth', 1.5, 'DisplayName', '测试预测');title('测试数据预测 - 通道1');xlabel('样本索引');ylabel('幅度');legend;grid on;% 子图4: 辨识误差收敛subplot(3, 3, 4);plot(identifier.identification_error, 'g-', 'LineWidth', 2);title('辨识误差收敛历程');xlabel('更新次数');ylabel('误差范数');grid on;% 子图5: 验证指标subplot(3, 3, 5);metrics = {'RMSE', 'MAE', 'MAPE', 'R²'};values = [validation_results.rmse, validation_results.mae, ...validation_results.mape * 100, validation_results.r_squared];bar(values);set(gca, 'XTickLabel', metrics);ylabel('指标值');title('模型验证指标');grid on;% 添加数值标签for i = 1:length(values)text(i, values(i), sprintf('%.4f', values(i)), ...'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');end% 子图6: 多步预测误差subplot(3, 3, 6);plot(prediction_results.horizons, prediction_results.errors, 'mo-', ...'LineWidth', 2, 'MarkerSize', 8);xlabel('预测步长');ylabel('RMSE');title('多步预测误差');grid on;% 子图7: 系统矩阵比较subplot(3, 3, 7);A_true = sys_true.A;A_ident = identifier.A;imagesc(abs(A_true - A_ident));colorbar;title('A矩阵误差绝对值');xlabel('状态索引');ylabel('状态索引');% 子图8: 频率响应比较subplot(3, 3, 8);w = logspace(-2, 1, 100);[mag_true, phase_true] = bode(sys_true, w);sys_ident = ss(identifier.A, identifier.B, identifier.C, identifier.D, 1);[mag_ident, phase_ident] = bode(sys_ident, w);subplot(3, 3, 8);semilogx(w, 20*log10(squeeze(mag_true)), 'b-', 'LineWidth', 2, 'DisplayName', '真实系统');hold on;semilogx(w, 20*log10(squeeze(mag_ident)), 'r--', 'LineWidth', 1.5, 'DisplayName', '辨识系统');title('频率响应比较');xlabel('频率 (rad/s)');ylabel('幅值 (dB)');legend;grid on;% 子图9: 残差分析subplot(3, 3, 9);Y_pred_val = lsim(sys_ident, U_test, (0:length(U_test)-1)');residuals = Y_test - Y_pred_val;autocorr(residuals(:,1), 20);title('残差自相关函数');grid on;sgtitle('自适应子空间辨识综合结果');% 输出总结报告print_summary_report(sys_true, identifier, validation_results);
endfunction print_summary_report(sys_true, identifier, validation_results)% 输出总结报告fprintf('\n=== 自适应子空间辨识总结报告 ===\n\n');fprintf('系统辨识结果:\n');fprintf(' 真实系统阶数: %d\n', size(sys_true.A, 1));fprintf(' 辨识系统阶数: %d\n', identifier.nx);fprintf('\n系统矩阵比较:\n');fprintf(' A矩阵范数误差: %.4f\n', norm(sys_true.A - identifier.A, 'fro'));fprintf(' B矩阵范数误差: %.4f\n', norm(sys_true.B - identifier.B, 'fro'));fprintf(' C矩阵范数误差: %.4f\n', norm(sys_true.C - identifier.C, 'fro'));fprintf('\n模型验证结果:\n');fprintf(' RMSE: %.4f\n', validation_results.rmse);fprintf(' MAE: %.4f\n', validation_results.mae);fprintf(' MAPE: %.2f%%\n', validation_results.mape * 100);fprintf(' R²: %.4f\n', validation_results.r_squared);fprintf(' 拟合度: %.2f%%\n', validation_results.fit_percentage);fprintf('\n自适应性能:\n');fprintf(' 最终辨识误差: %.4f\n', identifier.identification_error(end));fprintf(' 数据窗口大小: %d\n', identifier.window_size);fprintf(' 遗忘因子: %.3f\n', identifier.forgetting_factor);fprintf('\n系统稳定性分析:\n');eig_true = eig(sys_true.A);eig_ident = eig(identifier.A);fprintf(' 真实系统特征值模: '); fprintf('%.3f ', abs(eig_true)); fprintf('\n');fprintf(' 辨识系统特征值模: '); fprintf('%.3f ', abs(eig_ident)); fprintf('\n');if all(abs(eig_ident) < 1)fprintf(' ✅ 辨识系统稳定\n');elsefprintf(' ⚠️ 辨识系统可能不稳定\n');end
end% 运行主程序
main_adaptive_subspace_identification();
5. 实际应用案例：机械振动系统辨识
function vibration_system_identification()% 机械振动系统辨识案例fprintf('=== 机械振动系统自适应子空间辨识 ===\n\n');% 生成机械振动系统数据[U_vib, Y_vib, sys_vib_true] = generate_vibration_system_data();% 数据预处理[U_processed, Y_processed] = preprocess_vibration_data(U_vib, Y_vib);% 划分训练测试集N_total = size(Y_processed, 1);N_train = floor(0.7 * N_total);U_train = U_processed(1:N_train, :);Y_train = Y_processed(1:N_train, :);U_test = U_processed(N_train+1:end, :);Y_test = Y_processed(N_train+1:end, :);fprintf('振动系统数据:\n');fprintf(' 总样本数: %d\n', N_total);fprintf(' 训练样本: %d\n', N_train);fprintf(' 测试样本: %d\n', N_total - N_train);% 系统阶数估计estimated_order = estimate_system_order(Y_train, U_train);fprintf('估计的系统阶数: %d\n', estimated_order);% 自适应子空间辨识nx = estimated_order;nu = size(U_train, 2);ny = size(Y_train, 2);vib_identifier = AdaptiveSubspaceID(nx, nu, ny, ...'window_size', 100, ...'forgetting_factor', 0.95, ...'horizon', 15);% 在线辨识[vib_identifier, Y_pred_vib] = online_identification(vib_identifier, U_train, Y_train);% 模态分析perform_modal_analysis(vib_identifier, sys_vib_true);% 预测性能评估vib_validation = AdaptiveUpdate.validate_model(vib_identifier, U_test, Y_test);fprintf('\n振动系统辨识完成!\n');fprintf('最终验证RMSE: %.4f\n', vib_validation.rmse);
endfunction [U, Y, sys_true] = generate_vibration_system_data()% 生成机械振动系统数据% 2自由度质量-弹簧-阻尼系统m1 = 1.0; m2 = 0.8; % 质量k1 = 100; k2 = 80; % 刚度c1 = 0.5; c2 = 0.3; % 阻尼% 连续时间状态空间矩阵A_cont = [0, 0, 1, 0;0, 0, 0, 1;-k1/m1, k1/m1, -c1/m1, c1/m1;k1/m2, -(k1+k2)/m2, c1/m2, -(c1+c2)/m2];B_cont = [0; 0; 1/m1; 0]; % 力作用在质量1上C_cont = [1, 0, 0, 0; % 测量质量1的位移和速度0, 1, 0, 0];D_cont = [0; 0];% 转换为离散时间系统 (采样频率100Hz)Ts = 0.01;sys_cont = ss(A_cont, B_cont, C_cont, D_cont);sys_disc = c2d(sys_cont, Ts);sys_true = sys_disc;% 生成激励信号 (扫频信号 + 随机激励)N = 1000;t = (0:N-1)' * Ts;% 扫频信号f_start = 0.1; f_end = 10;chirp_signal = chirp(t, f_start, t(end), f_end);% 随机激励random_signal = 0.3 * randn(N, 1);% 脉冲激励impulse_signal = zeros(N, 1);impulse_signal(100:100:900) = 2;% 组合激励U = 0.5 * chirp_signal + 0.3 * random_signal + 0.2 * impulse_signal;% 系统响应 (加入测量噪声)Y = lsim(sys_disc, U, t) + 0.02 * randn(N, 2);fprintf('机械振动系统生成完成\n');fprintf(' 自然频率: %.2f Hz, %.2f Hz\n', ...sqrt(eig([k1+m1, -k1; -k1, k1+k2+m2])) / (2*pi));
endfunction [U_processed, Y_processed] = preprocess_vibration_data(U, Y)% 振动数据预处理% 去除直流分量U_processed = U - mean(U);Y_processed = Y - mean(Y);% 标准化U_processed = U_processed / std(U_processed);Y_processed = Y_processed ./ std(Y_processed);% 滤波处理 (可选)% [b, a] = butter(2, 0.1); % 低通滤波% Y_processed = filtfilt(b, a, Y_processed);
endfunction estimated_order = estimate_system_order(Y, U, max_order)% 估计系统阶数if nargin < 3max_order = 10;end[N, ny] = size(Y);[~, nu] = size(U);% 使用SVD确定系统阶数horizon = min(20, floor(N/3));H = SubspaceCore.build_combined_hankel(Y, U, horizon);[~, S, ~] = svd(H, 'econ');singular_values = diag(S);% 归一化奇异值normalized_sv = singular_values / singular_values(1);% 选择奇异值下降明显的点threshold = 0.01;estimated_order = find(normalized_sv < threshold, 1) - 1;if isempty(estimated_order)estimated_order = max_order;endestimated_order = min(estimated_order, max_order);% 可视化奇异值figure('Position', [100, 100, 800, 400]);subplot(1,2,1);semilogy(singular_values, 'bo-', 'LineWidth', 2, 'MarkerSize', 6);hold on;plot(estimated_order, singular_values(estimated_order), 'ro', ...'MarkerSize', 10, 'MarkerFaceColor', 'red');xlabel('奇异值索引');ylabel('奇异值');title('Hankel矩阵奇异值');grid on;subplot(1,2,2);plot(normalized_sv, 'bo-', 'LineWidth', 2, 'MarkerSize', 6);hold on;plot([1, length(normalized_sv)], [threshold, threshold], 'r--', 'LineWidth', 1.5);plot(estimated_order, normalized_sv(estimated_order), 'ro', ...'MarkerSize', 10, 'MarkerFaceColor', 'red');xlabel('奇异值索引');ylabel('归一化奇异值');title('系统阶数估计');legend('奇异值', '阈值', '估计阶数');grid on;fprintf('系统阶数估计完成: %d\n', estimated_order);
endfunction perform_modal_analysis(identifier, sys_true)% 执行模态分析fprintf('\n进行模态分析...\n');% 辨识系统的模态参数[V_ident, D_ident] = eig(identifier.A);eigenvalues_ident = diag(D_ident);% 真实系统的模态参数eigenvalues_true = eig(sys_true.A);% 计算模态频率和阻尼比Ts = 0.01; % 采样时间[freq_ident, damping_ident] = compute_modal_parameters(eigenvalues_ident, Ts);[freq_true, damping_true] = compute_modal_parameters(eigenvalues_true, Ts);% 显示模态参数fprintf('模态参数比较:\n');fprintf(' 模态\t真实频率(Hz)\t辨识频率(Hz)\t真实阻尼比\t辨识阻尼比\n');num_modes = min(length(freq_true), length(freq_ident));for i = 1:num_modesfprintf(' %d\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n', ...i, freq_true(i), freq_ident(i), damping_true(i), damping_ident(i));end% 可视化模态振型figure('Position', [100, 100, 1000, 600]);subplot(2,2,1);stem(1:length(freq_true), freq_true, 'b', 'LineWidth', 2, 'DisplayName', '真实频率');hold on;stem(1:length(freq_ident), freq_ident, 'r', 'LineWidth', 2, 'DisplayName', '辨识频率');xlabel('模态阶数');ylabel('频率 (Hz)');title('模态频率比较');legend;grid on;subplot(2,2,2);stem(1:length(damping_true), damping_true, 'b', 'LineWidth', 2, 'DisplayName', '真实阻尼');hold on;stem(1:length(damping_ident), damping_ident, 'r', 'LineWidth', 2, 'DisplayName', '辨识阻尼');xlabel('模态阶数');ylabel('阻尼比');title('模态阻尼比较');legend;grid on;% 模态置信度分析subplot(2,2,3);mode_confidence = analyze_mode_confidence(identifier, sys_true);bar(mode_confidence);xlabel('模态阶数');ylabel('置信度指标');title('模态辨识置信度');grid on;subplot(2,2,4);% 特征值误差eig_error = abs(eigenvalues_true(1:num_modes) - eigenvalues_ident(1:num_modes));bar(eig_error);xlabel('模态阶数');ylabel('特征值误差');title('特征值辨识误差');grid on;sgtitle('机械振动系统模态分析');
endfunction [frequencies, damping_ratios] = compute_modal_parameters(eigenvalues, Ts)% 计算模态频率和阻尼比% 转换为连续时间特征值s = log(eigenvalues) / Ts;% 频率和阻尼比frequencies = abs(s) / (2*pi);damping_ratios = -real(s) ./ abs(s);% 排序并过滤[frequencies, idx] = sort(frequencies);damping_ratios = damping_ratios(idx);% 只返回正频率和合理阻尼比valid_idx = frequencies > 0 & damping_ratios > 0 & damping_ratios < 1;frequencies = frequencies(valid_idx);damping_ratios = damping_ratios(valid_idx);
endfunction confidence = analyze_mode_confidence(identifier, sys_true)% 分析模态辨识置信度eigenvalues_ident = eig(identifier.A);eigenvalues_true = eig(sys_true.A);num_modes = min(length(eigenvalues_ident), length(eigenvalues_true));confidence = zeros(num_modes, 1);for i = 1:num_modes% 基于特征值误差的置信度eig_error = abs(eigenvalues_true(i) - eigenvalues_ident(i));confidence(i) = exp(-eig_error);end
end% 运行振动系统辨识案例
vibration_system_identification();
参考代码 自适应子空间辨识，通过输入输出数据，计算状态空间方程，建模并预测 www.youwenfan.com/contentcnl/80122.html
特点：
- 多种子空间算法：NASID、PO-MOESP等主流方法
- 自适应更新：滑动窗口、遗忘因子、递归更新
- 系统阶数估计：基于SVD的自动阶数确定
- 完整验证流程：训练、测试、预测、模态分析
- 实际应用案例：机械振动系统辨识
优势：
- 在线适应能力：能够跟踪时变系统
- 数值稳定性：包含正则化和稳定化处理
- 工程实用性：直接处理输入输出数据
- 全面分析工具：包含模态分析、频率响应分析等
价值：
- 为控制系统设计和分析提供准确模型
- 支持故障检测和系统监控
- 可用于预测控制和系统仿真
- 为机械、电气、化工等领域的系统辨识提供工具
