单自由度轴向磁悬浮轴承Simulink模型
一、系统数学模型
1.1 系统动力学方程
m * d²x/dt² = F_em - F_d - m*g
其中：
- m：转子质量
- x：转子轴向位移（向下为正）
- F_em：电磁力
- F_d：干扰力
- g：重力加速度
1.2 电磁力模型
F_em = (μ₀ * N² * A * i²) / (4 * (x₀ + x)²) - (μ₀ * N² * A * i²) / (4 * (x₀ - x)²)
简化线性化模型：
F_em = k_i * i + k_x * x
其中：
- k_i：电流刚度系数
- k_x：位移刚度系数
二、Simulink模型实现
2.1 主模型（axial_magnetic_bearing.slx）
%% 创建Simulink模型
model_name = 'axial_magnetic_bearing';
new_system(model_name);
open_system(model_name);
%% 设置仿真参数
set_param(model_name, 'Solver', 'ode4', 'FixedStep', '1e-4', ...
'StopTime', '0.5', 'SaveFormat', 'Array');
%% 添加子系统模块
% 1. 控制器子系统
add_block('simulink/User-Defined Functions/Subsystem', ...
[model_name '/Controller'], 'Position', [100, 100, 200, 200]);
% 2. 功率放大器子系统
add_block('simulink/Discrete/Discrete Filter', ...
[model_name '/Power_Amplifier'], 'Position', [300, 150, 350, 200]);
% 3. 电磁铁子系统
add_block('simulink/User-Defined Functions/Subsystem', ...
[model_name '/Electromagnet'], 'Position', [450, 100, 550, 250]);
% 4. 转子动力学子系统
add_block('simulink/Continuous/Integrator', ...
[model_name '/Integrator1'], 'Position', [650, 150, 700, 200]);
add_block('simulink/Continuous/Integrator', ...
[model_name '/Integrator2'], 'Position', [750, 150, 800, 200]);
% 5. 传感器子系统
add_block('simulink/Discrete/Zero-Order Hold', ...
[model_name '/Sensor'], 'Position', [850, 150, 900, 200]);
% 6. 参考输入
add_block('simulink/Sources/Step', ...
[model_name '/Reference'], 'Position', [50, 150, 80, 180]);
% 7. 干扰输入
add_block('simulink/Sources/Step', ...
[model_name '/Disturbance'], 'Position', [600, 50, 630, 80]);
% 8. 示波器
add_block('simulink/Sinks/Scope', ...
[model_name '/Scope'], 'Position', [950, 100, 1000, 300]);
%% 连接模块
% 参考信号到控制器
add_line(model_name, 'Reference/1', 'Controller/1');
% 传感器反馈到控制器
add_line(model_name, 'Sensor/1', 'Controller/2');
% 控制器到功率放大器
add_line(model_name, 'Controller/1', 'Power_Amplifier/1');
% 功率放大器到电磁铁
add_line(model_name, 'Power_Amplifier/1', 'Electromagnet/1');
% 电磁铁到转子
add_line(model_name, 'Electromagnet/1', 'Integrator2/1');
% 转子动力学连接
add_line(model_name, 'Integrator2/1', 'Integrator1/1');
% 干扰力
add_line(model_name, 'Disturbance/1', 'Integrator2/1');
% 位移到传感器
add_line(model_name, 'Integrator1/1', 'Sensor/1');
% 连接到示波器
add_line(model_name, 'Reference/1', 'Scope/1', 'autorouting', 'on');
add_line(model_name, 'Sensor/1', 'Scope/2', 'autorouting', 'on');
add_line(model_name, 'Electromagnet/1', 'Scope/3', 'autorouting', 'on');
save_system(model_name);
2.2 完整的集成模型代码
%% 单自由度轴向磁悬浮轴承Simulink模型
% 文件名: axial_magnetic_bearing_complete.m
function create_axial_magnetic_bearing_model()
% 创建完整的单自由度轴向磁悬浮轴承Simulink模型
%% 清理工作空间
clear all;
close all;
bdclose all;
%% 系统参数定义
define_system_parameters();
%% 创建新模型
model_name = 'axial_magnetic_bearing_complete';
if exist(model_name, 'file')
delete([model_name '.slx']);
end
new_system(model_name);
open_system(model_name);
%% 设置模型参数
set_param(model_name, 'Solver', 'ode4', ... % 固定步长Runge-Kutta
'FixedStep', '1e-5', ... % 固定步长
'StopTime', '0.2', ... % 仿真时间
'SaveFormat', 'Array', ... % 数据格式
'SaveOutput', 'on', ... % 保存输出
'SaveTime', 'on'); % 保存时间
%% 添加主系统模块
add_main_blocks(model_name);
%% 配置子系统
configure_controller_subsystem(model_name);
configure_electromagnet_subsystem(model_name);
configure_rotor_dynamics_subsystem(model_name);
configure_power_amplifier_subsystem(model_name);
%% 连接所有模块
connect_all_blocks(model_name);
%% 添加观察和测量模块
add_measurement_blocks(model_name);
%% 保存模型
save_system(model_name);
disp(['模型已创建: ' model_name '.slx']);
%% 运行仿真
run_simulation(model_name);
%% 绘制结果
plot_simulation_results();
end
%% 函数定义部分
function define_system_parameters()
% 定义系统参数
global params;
% 物理参数
params.m = 2.0; % 转子质量 [kg]
params.g = 9.81; % 重力加速度 [m/s²]
params.x0 = 0.001; % 额定气隙 [m] (1mm)
params.x_max = 0.0005; % 最大允许位移 [m] (0.5mm)
% 电磁参数
params.mu0 = 4*pi*1e-7; % 真空磁导率 [H/m]
params.N = 200; % 线圈匝数
params.A = 1e-3; % 磁极面积 [m²]
params.R_coil = 2.0; % 线圈电阻 [Ω]
params.L_coil = 0.1; % 线圈电感 [H]
% 电磁力参数（线性化系数）
params.ki = 50; % 电流刚度系数 [N/A]
params.kx = -2e5; % 位移刚度系数 [N/m] (负值表示不稳定)
% 控制器参数
params.Kp = 5000; % 比例增益
params.Ki = 1000; % 积分增益
params.Kd = 50; % 微分增益
% 功率放大器参数
params.K_amp = 10; % 放大器增益 [A/V]
params.f_amp = 1000; % 放大器带宽 [Hz]
params.current_max = 10; % 最大电流 [A]
% 传感器参数
params.K_sensor = 1000; % 传感器增益 [V/m]
params.noise_density = 1e-6; % 噪声密度 [V/√Hz]
% 干扰参数
params.Fd_step = 10; % 阶跃干扰力 [N]
params.Fd_step_time = 0.1; % 阶跃干扰时间 [s]
% 参考信号
params.ref_step = 0; % 参考位移 [m] (0表示平衡位置)
params.ref_step_time = 0.05; % 阶跃时间 [s]
% 初始条件
params.x_initial = 0; % 初始位移 [m]
params.v_initial = 0; % 初始速度 [m/s]
% 保存到基础工作空间
assignin('base', 'params', params);
end
function add_main_blocks(model_name)
% 添加主系统模块
% 参考信号（阶跃输入）
add_block('simulink/Sources/Step', ...
[model_name '/Reference']);
set_param([model_name '/Reference'], ...
'After', 'params.ref_step_time', ...
'Before', '0', ...
'SampleTime', '1e-5');
% 控制器
add_block('simulink/User-Defined Functions/Subsystem', ...
[model_name '/Controller'], 'Position', [100, 100, 220, 200]);
% 功率放大器
add_block('simulink/Continuous/Transfer Fcn', ...
[model_name '/Power_Amplifier'], 'Position', [250, 150, 320, 200]);
% 电流限制器
add_block('simulink/Discontinuities/Saturation', ...
[model_name '/Current_Limiter'], 'Position', [340, 150, 370, 200]);
set_param([model_name '/Current_Limiter'], ...
'UpperLimit', 'params.current_max', ...
'LowerLimit', '-params.current_max');
% 电磁铁子系统
add_block('simulink/User-Defined Functions/Subsystem', ...
[model_name '/Electromagnet'], 'Position', [400, 100, 500, 250]);
% 转子动力学
add_block('simulink/User-Defined Functions/Subsystem', ...
[model_name '/Rotor_Dynamics'], 'Position', [550, 100, 650, 250]);
% 干扰力
add_block('simulink/Sources/Step', ...
[model_name '/Disturbance']);
set_param([model_name '/Disturbance'], ...
'After', 'params.Fd_step_time', ...
'Before', '0', ...
'SampleTime', '1e-5');
% 传感器
add_block('simulink/Discontinuities/Saturation', ...
[model_name '/Sensor_Saturation'], 'Position', [750, 150, 780, 200]);
set_param([model_name '/Sensor_Saturation'], ...
'UpperLimit', 'params.x_max', ...
'LowerLimit', '-params.x_max');
% 传感器增益
add_block('simulink/Math Operations/Gain', ...
[model_name '/Sensor_Gain'], 'Position', [800, 150, 830, 200]);
set_param([model_name '/Sensor_Gain'], 'Gain', 'params.K_sensor');
% 传感器噪声
add_block('simulink/Sources/Band-Limited White Noise', ...
[model_name '/Sensor_Noise'], 'Position', [850, 180, 880, 210]);
% 主示波器
add_block('simulink/Sinks/Scope', ...
[model_name '/Main_Scope'], 'Position', [900, 50, 950, 350]);
set_param([model_name '/Main_Scope'], 'NumInputPorts', '6');
% 工作空间输出
add_block('simulink/Sinks/To Workspace', ...
[model_name '/To_Workspace']);
set_param([model_name '/To_Workspace'], ...
'VariableName', 'simout', ...
'SaveFormat', 'Array');
end
function configure_controller_subsystem(model_name)
% 配置控制器子系统
controller_path = [model_name '/Controller'];
open_system(controller_path);
% 清除默认内容
delete_line(controller_path, 'In1/1', 'Out1/1');
% 添加PID控制器
add_block('simulink/Continuous/PID Controller', ...
[controller_path '/PID']);
set_param([controller_path '/PID'], ...
'P', 'params.Kp', ...
'I', 'params.Ki', ...
'D', 'params.Kd', ...
'N', '100');
% 添加误差计算
add_block('simulink/Math Operations/Sum', ...
[controller_path '/Error_Sum'], 'Position', [50, 100, 80, 130]);
set_param([controller_path '/Error_Sum'], ...
'IconShape', 'round', ...
'Inputs', '|+-');
% 连接
add_line(controller_path, 'In1/1', 'Error_Sum/1');
add_line(controller_path, 'In2/1', 'Error_Sum/2');
add_line(controller_path, 'Error_Sum/1', 'PID/1');
add_line(controller_path, 'PID/1', 'Out1/1');
% 关闭子系统
close_system(controller_path, 1); % 1表示不保存
end
function configure_electromagnet_subsystem(model_name)
% 配置电磁铁子系统
electromagnet_path = [model_name '/Electromagnet'];
open_system(electromagnet_path);
% 清除默认内容
delete_line(electromagnet_path, 'In1/1', 'Out1/1');
% 添加电磁力计算（非线性模型）
add_block('simulink/User-Defined Functions/Fcn', ...
[electromagnet_path '/Electromagnetic_Force']);
set_param([electromagnet_path '/Electromagnetic_Force'], ...
'Expr', 'params.mu0*params.N^2*params.A*u(1)^2/(4*(params.x0+u(2))^2) - params.mu0*params.N^2*params.A*u(1)^2/(4*(params.x0-u(2))^2)');
% 添加线性化模型（用于比较）
add_block('simulink/Math Operations/Gain', ...
[electromagnet_path '/Ki_Gain'], 'Position', [100, 50, 130, 80]);
set_param([electromagnet_path '/Ki_Gain'], 'Gain', 'params.ki');
add_block('simulink/Math Operations/Gain', ...
[electromagnet_path '/Kx_Gain'], 'Position', [100, 120, 130, 150]);
set_param([electromagnet_path '/Kx_Gain'], 'Gain', 'params.kx');
add_block('simulink/Math Operations/Sum', ...
[electromagnet_path '/Force_Sum'], 'Position', [200, 90, 230, 130]);
set_param([electromagnet_path '/Force_Sum'], 'Inputs', '++');
% 添加多路开关选择模型
add_block('simulink/Signal Routing/Multiport Switch', ...
[electromagnet_path '/Model_Selector'], 'Position', [150, 200, 180, 230]);
% 添加常数（选择器）
add_block('simulink/Sources/Constant', ...
[electromagnet_path '/Model_Choice'], 'Position', [50, 200, 80, 230]);
set_param([electromagnet_path '/Model_Choice'], 'Value', '1'); % 1=线性模型，2=非线性模型
% 连接
add_line(electromagnet_path, 'In1/1', 'Ki_Gain/1');
add_line(electromagnet_path, 'In2/1', 'Kx_Gain/1');
add_line(electromagnet_path, 'Ki_Gain/1', 'Force_Sum/1');
add_line(electromagnet_path, 'Kx_Gain/1', 'Force_Sum/2');
add_line(electromagnet_path, 'Force_Sum/1', 'Model_Selector/1');
add_line(electromagnet_path, 'In1/1', 'Electromagnetic_Force/1');
add_line(electromagnet_path, 'In2/1', 'Electromagnetic_Force/2');
add_line(electromagnet_path, 'Electromagnetic_Force/1', 'Model_Selector/2');
add_line(electromagnet_path, 'Model_Choice/1', 'Model_Selector/3');
add_line(electromagnet_path, 'Model_Selector/1', 'Out1/1');
% 关闭子系统
close_system(electromagnet_path, 1);
end
function configure_rotor_dynamics_subsystem(model_name)
% 配置转子动力学子系统
rotor_path = [model_name '/Rotor_Dynamics'];
open_system(rotor_path);
% 清除默认内容
delete_line(rotor_path, 'In1/1', 'Out1/1');
% 添加积分器（速度）
add_block('simulink/Continuous/Integrator', ...
[rotor_path '/Integrator_v'], 'Position', [200, 100, 230, 130]);
set_param([rotor_path '/Integrator_v'], ...
'InitialCondition', 'params.v_initial');
% 添加积分器（位移）
add_block('simulink/Continuous/Integrator', ...
[rotor_path '/Integrator_x'], 'Position', [300, 100, 330, 130]);
set_param([rotor_path '/Integrator_x'], ...
'InitialCondition', 'params.x_initial');
% 添加增益（1/m）
add_block('simulink/Math Operations/Gain', ...
[rotor_path '/Mass_Gain'], 'Position', [100, 100, 130, 130]);
set_param([rotor_path '/Mass_Gain'], 'Gain', '1/params.m');
% 添加求和点（总力）
add_block('simulink/Math Operations/Sum', ...
[rotor_path '/Force_Sum'], 'Position', [50, 100, 80, 140]);
set_param([rotor_path '/Force_Sum'], 'Inputs', '++-');
% 添加重力
add_block('simulink/Sources/Constant', ...
[rotor_path '/Gravity'], 'Position', [20, 150, 50, 180]);
set_param([rotor_path '/Gravity'], 'Value', 'params.m*params.g');
% 连接
add_line(rotor_path, 'In1/1', 'Force_Sum/1'); % 电磁力
add_line(rotor_path, 'In2/1', 'Force_Sum/2'); % 干扰力
add_line(rotor_path, 'Gravity/1', 'Force_Sum/3'); % 重力
add_line(rotor_path, 'Force_Sum/1', 'Mass_Gain/1');
add_line(rotor_path, 'Mass_Gain/1', 'Integrator_v/1');
add_line(rotor_path, 'Integrator_v/1', 'Integrator_x/1');
add_line(rotor_path, 'Integrator_x/1', 'Out1/1');
% 关闭子系统
close_system(rotor_path, 1);
end
function configure_power_amplifier_subsystem(model_name)
% 配置功率放大器子系统
set_param([model_name '/Power_Amplifier'], ...
'Numerator', '[params.K_amp]', ...
'Denominator', '[1/(2*pi*params.f_amp) 1]');
end
function connect_all_blocks(model_name)
% 连接所有模块
% 参考信号到控制器
add_line(model_name, 'Reference/1', 'Controller/1');
% 传感器反馈到控制器
add_line(model_name, 'Sensor_Gain/1', 'Controller/2');
% 控制器到功率放大器
add_line(model_name, 'Controller/1', 'Power_Amplifier/1');
% 功率放大器到电流限制器
add_line(model_name, 'Power_Amplifier/1', 'Current_Limiter/1');
% 电流限制器到电磁铁（电流输入）
add_line(model_name, 'Current_Limiter/1', 'Electromagnet/1');
% 位移反馈到电磁铁（位移输入）
add_line(model_name, 'Rotor_Dynamics/1', 'Electromagnet/2');
% 电磁铁到转子动力学（电磁力）
add_line(model_name, 'Electromagnet/1', 'Rotor_Dynamics/1');
% 干扰力到转子动力学
add_line(model_name, 'Disturbance/1', 'Rotor_Dynamics/2');
% 转子位移到传感器
add_line(model_name, 'Rotor_Dynamics/1', 'Sensor_Saturation/1');
% 传感器饱和后加增益
add_line(model_name, 'Sensor_Saturation/1', 'Sensor_Gain/1');
% 添加传感器噪声
add_line(model_name, 'Sensor_Gain/1', 'Sensor_Noise/1');
% 连接到主示波器
add_line(model_name, 'Reference/1', 'Main_Scope/1');
add_line(model_name, 'Rotor_Dynamics/1', 'Main_Scope/2');
add_line(model_name, 'Sensor_Saturation/1', 'Main_Scope/3');
add_line(model_name, 'Current_Limiter/1', 'Main_Scope/4');
add_line(model_name, 'Electromagnet/1', 'Main_Scope/5');
add_line(model_name, 'Disturbance/1', 'Main_Scope/6');
% 连接到工作空间
add_line(model_name, 'Rotor_Dynamics/1', 'To_Workspace/1');
end
function add_measurement_blocks(model_name)
% 添加额外的测量模块
% 添加XY图显示相轨迹
add_block('simulink/Sinks/XY Graph', ...
[model_name '/Phase_Plot'], 'Position', [900, 400, 950, 450]);
% 添加位移-速度连接
add_block('simulink/Signal Routing/Demux', ...
[model_name '/Demux'], 'Position', [850, 400, 870, 450]);
% 连接相轨迹
add_line(model_name, 'Rotor_Dynamics/1', 'Demux/1');
add_line(model_name, 'Demux/1', 'Phase_Plot/1');
add_line(model_name, 'Demux/2', 'Phase_Plot/2');
% 添加性能指标计算
add_block('simulink/User-Defined Functions/Subsystem', ...
[model_name '/Performance'], 'Position', [800, 300, 900, 450]);
% 打开性能计算子系统
performance_path = [model_name '/Performance'];
open_system(performance_path);
% 添加均方根误差计算
add_block('simulink/User-Defined Functions/Fcn', ...
[performance_path '/RMSE'], 'Position', [100, 50, 200, 80]);
set_param([performance_path '/RMSE'], 'Expr', 'sqrt(mean(u.*u))');
% 添加过冲计算
add_block('simulink/Logic and Bit Operations/Relational Operator', ...
[performance_path '/Overshoot_Detect'], 'Position', [50, 100, 80, 130]);
add_block('simulink/Signal Routing/Max', ...
[performance_path '/Max_Value'], 'Position', [150, 100, 180, 130]);
% 连接
add_line(performance_path, 'In1/1', 'RMSE/1');
add_line(performance_path, 'In1/1', 'Overshoot_Detect/1');
add_line(performance_path, 'In1/1', 'Max_Value/1');
add_block('simulink/Sinks/Display', ...
[performance_path '/Display_RMSE'], 'Position', [250, 50, 300, 80]);
add_line(performance_path, 'RMSE/1', 'Display_RMSE/1');
close_system(performance_path, 1);
end
function run_simulation(model_name)
% 运行仿真
disp('开始仿真...');
simout = sim(model_name, 'ReturnWorkspaceOutputs', 'on');
disp('仿真完成。');
% 保存仿真结果
assignin('base', 'simulation_results', simout);
end
function plot_simulation_results()
% 绘制仿真结果
global params;
% 检查是否有仿真结果
if ~exist('simout', 'var')
error('请先运行仿真！');
end
% 提取数据
t = simout(:,1);
x_ref = simout(:,2); % 参考位移
x_actual = simout(:,3); % 实际位移
x_sensor = simout(:,4); % 传感器输出
i_coil = simout(:,5); % 线圈电流
F_em = simout(:,6); % 电磁力
F_dist = simout(:,7); % 干扰力
% 计算误差
error = x_actual - x_ref;
% 创建图形
figure('Position', [100, 100, 1200, 800]);
% 1. 位移响应
subplot(3, 3, 1);
plot(t, x_ref*1000, 'r--', 'LineWidth', 2); hold on;
plot(t, x_actual*1000, 'b-', 'LineWidth', 1.5);
xlabel('时间 [s]');
ylabel('位移 [mm]');
title('转子位移响应');
legend('参考位移', '实际位移', 'Location', 'best');
grid on;
ylim([-params.x_max*1500, params.x_max*1500]);
% 2. 跟踪误差
subplot(3, 3, 2);
plot(t, error*1e6, 'r-', 'LineWidth', 1.5);
xlabel('时间 [s]');
ylabel('误差 [μm]');
title('跟踪误差');
grid on;
% 3. 控制电流
subplot(3, 3, 3);
plot(t, i_coil, 'b-', 'LineWidth', 1.5);
xlabel('时间 [s]');
ylabel('电流 [A]');
title('控制电流');
grid on;
ylim([-params.current_max*1.1, params.current_max*1.1]);
% 4. 电磁力
subplot(3, 3, 4);
plot(t, F_em, 'g-', 'LineWidth', 1.5); hold on;
plot(t, params.m*params.g*ones(size(t)), 'k--', 'LineWidth', 1);
xlabel('时间 [s]');
ylabel('电磁力 [N]');
title('电磁力');
legend('电磁力', '重力', 'Location', 'best');
grid on;
% 5. 干扰力
subplot(3, 3, 5);
plot(t, F_dist, 'm-', 'LineWidth', 1.5);
xlabel('时间 [s]');
ylabel('干扰力 [N]');
title('干扰力');
grid on;
% 6. 相轨迹
subplot(3, 3, 6);
plot(x_actual*1000, gradient(x_actual)./gradient(t), 'b-', 'LineWidth', 1.5);
xlabel('位移 [mm]');
ylabel('速度 [m/s]');
title('相轨迹');
grid on;
% 7. 功率谱密度（位移）
subplot(3, 3, 7);
[pxx, f] = pwelch(error, [], [], [], 1/(t(2)-t(1)));
semilogy(f, pxx, 'b-', 'LineWidth', 1.5);
xlabel('频率 [Hz]');
ylabel('PSD [m²/Hz]');
title('误差功率谱密度');
grid on;
xlim([0, 500]);
% 8. 控制性能指标
subplot(3, 3, 8);
hold on;
% 计算指标
rmse = sqrt(mean(error.^2));
max_overshoot = max(abs(error)) * 1e6;
settling_time = find(abs(error) < 1e-6, 1) * (t(2)-t(1));
if isempty(settling_time)
settling_time = t(end);
end
% 显示指标
text(0.1, 0.8, sprintf('RMSE: %.2f μm', rmse*1e6), 'FontSize', 10);
text(0.1, 0.6, sprintf('最大过冲: %.2f μm', max_overshoot), 'FontSize', 10);
text(0.1, 0.4, sprintf('调节时间: %.3f s', settling_time), 'FontSize', 10);
text(0.1, 0.2, sprintf('最大电流: %.2f A', max(abs(i_coil))), 'FontSize', 10);
axis off;
title('性能指标');
% 9. 频域分析（Bode图概念展示）
subplot(3, 3, 9);
% 计算频率响应
window_size = min(1024, length(error));
[Pxy, F] = cpsd(x_ref, x_actual, hamming(window_size), [], [], 1/(t(2)-t(1)));
magnitude = abs(Pxy);
plot(F, 20*log10(magnitude+eps), 'b-', 'LineWidth', 1.5);
xlabel('频率 [Hz]');
ylabel('幅值 [dB]');
title('频率响应');
grid on;
xlim([0, 200]);
sgtitle('单自由度轴向磁悬浮轴承仿真结果', 'FontSize', 14, 'FontWeight', 'bold');
% 保存图形
saveas(gcf, 'axial_bearing_simulation_results.png');
end
2.3 控制器设计模块
%% PID控制器调优函数
function tune_pid_controller()
% PID控制器参数调优
% 定义系统传递函数（线性化模型）
% 开环传递函数：G(s) = ki / (m*s² - kx)
global params;
% 定义开环传递函数
num = params.ki;
den = [params.m, 0, -params.kx];
sys_open = tf(num, den);
% 显示系统极点
poles = pole(sys_open);
disp('开环系统极点：');
disp(poles);
disp('注意：有一个正实部极点，系统开环不稳定！');
% PID控制器设计
% 使用根轨迹法或频率响应法
% 方法1：使用根轨迹设计
figure('Position', [100, 100, 1000, 400]);
subplot(1, 2, 1);
rlocus(sys_open);
title('开环系统根轨迹');
grid on;
% 方法2：尝试不同的PID参数
Kp_test = [1000, 5000, 10000];
Ki_test = [100, 1000, 5000];
Kd_test = [10, 50, 100];
% 创建对比图
figure('Position', [100, 100, 1200, 800]);
for i = 1:length(Kp_test)
% 创建PID控制器
C = pid(Kp_test(i), Ki_test(i), Kd_test(i));
% 闭环系统
sys_cl = feedback(C * sys_open, 1);
% 阶跃响应
subplot(2, 3, i);
step(sys_cl, 0.1);
title(sprintf('PID: Kp=%.0f, Ki=%.0f, Kd=%.0f', ...
Kp_test(i), Ki_test(i), Kd_test(i)));
grid on;
% 计算性能指标
info = stepinfo(sys_cl);
fprintf('PID参数 [%.0f, %.0f, %.0f]: 调节时间=%.3fs, 超调=%.1f%%\n', ...
Kp_test(i), Ki_test(i), Kd_test(i), info.SettlingTime, info.Overshoot);
end
% 推荐参数
disp('推荐PID参数：');
disp(sprintf('Kp = %.0f', params.Kp));
disp(sprintf('Ki = %.0f', params.Ki));
disp(sprintf('Kd = %.0f', params.Kd));
end
%% 鲁棒性分析
function robustness_analysis()
% 系统鲁棒性分析
global params;
% 定义标称系统
sys_nominal = tf(params.ki, [params.m, 0, -params.kx]);
% 参数变化范围
m_variation = params.m * [0.8, 1.0, 1.2]; % ±20%质量变化
ki_variation = params.ki * [0.9, 1.0, 1.1]; % ±10%电流刚度变化
kx_variation = params.kx * [0.9, 1.0, 1.1]; % ±10%位移刚度变化
% 创建PID控制器
C = pid(params.Kp, params.Ki, params.Kd);
% 分析质量变化的影响
figure('Position', [100, 100, 1200, 400]);
subplot(1, 3, 1);
hold on;
for i = 1:length(m_variation)
sys_var = tf(params.ki, [m_variation(i), 0, -params.kx]);
sys_cl = feedback(C * sys_var, 1);
step(sys_cl, 0.1);
end
title('质量变化影响 (±20%)');
legend('m=1.6kg', 'm=2.0kg', 'm=2.4kg');
grid on;
subplot(1, 3, 2);
hold on;
for i = 1:length(ki_variation)
sys_var = tf(ki_variation(i), [params.m, 0, -params.kx]);
sys_cl = feedback(C * sys_var, 1);
step(sys_cl, 0.1);
end
title('电流刚度变化影响 (±10%)');
legend('ki=45', 'ki=50', 'ki=55');
grid on;
subplot(1, 3, 3);
hold on;
for i = 1:length(kx_variation)
sys_var = tf(params.ki, [params.m, 0, -kx_variation(i)]);
sys_cl = feedback(C * sys_var, 1);
step(sys_cl, 0.1);
end
title('位移刚度变化影响 (±10%)');
legend('kx=-1.8e5', 'kx=-2.0e5', 'kx=-2.2e5');
grid on;
sgtitle('系统鲁棒性分析');
end
2.4 高级控制算法
%% 滑模控制器
function smc_controller_design()
% 滑模控制器设计
global params;
% 系统状态空间模型
A = [0, 1;
params.kx/params.m, 0];
B = [0;
params.ki/params.m];
C = [1, 0];
D = 0;
sys_ss = ss(A, B, C, D);
% 滑模面设计：s = λ*e + de/dt
lambda = 100; % 滑模面参数
% 控制律：u = u_eq + u_sw
% u_eq: 等效控制
% u_sw: 切换控制
% 设计参数
eta = 10; % 切换增益
phi = 0.1; % 边界层厚度
% 创建滑模控制器Simulink模型
smc_model = 'sliding_mode_controller';
if exist(smc_model, 'file')
delete([smc_model '.slx']);
end
new_system(smc_model);
open_system(smc_model);
% 添加模块
add_block('simulink/Continuous/Integrator', ...
[smc_model '/Integrator1']);
add_block('simulink/Continuous/Integrator', ...
[smc_model '/Integrator2']);
add_block('simulink/Math Operations/Gain', ...
[smc_model '/Gain_lambda']);
set_param([smc_model '/Gain_lambda'], 'Gain', 'lambda');
add_block('simulink/Math Operations/Sum', ...
[smc_model '/Sum1']);
set_param([smc_model '/Sum1'], 'Inputs', '++');
add_block('simulink/Math Operations/Sum', ...
[smc_model '/Sum2']);
set_param([smc_model '/Sum2'], 'Inputs', '++');
add_block('simulink/Math Operations/Product', ...
[smc_model '/Product1']);
add_block('simulink/Discontinuities/Saturation', ...
[smc_model '/Saturation']);
add_block('simulink/User-Defined Functions/Fcn', ...
[smc_model '/Switching_Control']);
set_param([smc_model '/Switching_Control'], ...
'Expr', '-eta*sat(s/phi)');
% 连接模块
add_line(smc_model, 'In1/1', 'Sum1/1'); % 参考位移
add_line(smc_model, 'Integrator2/1', 'Sum1/2'); % 实际位移
add_line(smc_model, 'Sum1/1', 'Gain_lambda/1'); % 误差
add_line(smc_model, 'Gain_lambda/1', 'Sum2/1'); % λ*e
add_line(smc_model, 'Integrator1/1', 'Sum2/2'); % 速度（负反馈）
add_line(smc_model, 'Sum2/1', 'Product1/1'); % 滑模面s
add_line(smc_model, 'Product1/1', 'Switching_Control/1');
add_line(smc_model, 'Switching_Control/1', 'Saturation/1');
add_line(smc_model, 'Saturation/1', 'Out1/1');
save_system(smc_model);
close_system(smc_model);
disp('滑模控制器模型已创建。');
end
%% 自适应控制器
function adaptive_controller()
% 模型参考自适应控制器
global params;
% 参考模型（期望的动态特性）
omega_n = 100; % 自然频率 [rad/s]
zeta = 0.7; % 阻尼比
A_ref = [0, 1;
-omega_n^2, -2*zeta*omega_n];
B_ref = [0;
omega_n^2];
C_ref = [1, 0];
D_ref = 0;
ref_model = ss(A_ref, B_ref, C_ref, D_ref);
% 显示参考模型特性
disp('参考模型特性：');
disp(['自然频率: ', num2str(omega_n), ' rad/s']);
disp(['阻尼比: ', num2str(zeta)]);
step(ref_model);
title('参考模型阶跃响应');
grid on;
% 自适应律参数
Gamma = diag([1e6, 1e6]); % 自适应增益矩阵
% 创建MRAC Simulink模型
mrac_model = 'mrac_controller';
if exist(mrac_model, 'file')
delete([mrac_model '.slx']);
end
new_system(mrac_model);
open_system(mrac_model);
% 添加参考模型
add_block('simulink/Continuous/State-Space', ...
[mrac_model '/Reference_Model']);
set_param([mrac_model '/Reference_Model'], ...
'A', mat2str(A_ref), ...
'B', mat2str(B_ref), ...
'C', mat2str(C_ref), ...
'D', mat2str(D_ref));
% 添加自适应控制器模块
add_block('simulink/User-Defined Functions/Subsystem', ...
[mrac_model '/Adaptive_Controller']);
% 保存模型
save_system(mrac_model);
close_system(mrac_model);
disp('模型参考自适应控制器框架已创建。');
end
参考代码 单自由度轴向磁悬浮轴承的Simulink模型 www.youwenfan.com/contentcnp/96473.html
三、模型使用说明
3.1 运行仿真
%% 主运行脚本
% 文件: run_axial_bearing_simulation.m
clear all;
close all;
clc;
% 1. 创建模型
create_axial_magnetic_bearing_model();
% 2. 打开模型
open_system('axial_magnetic_bearing_complete');
% 3. 运行控制器调优（可选）
% tune_pid_controller();
% 4. 运行鲁棒性分析（可选）
% robustness_analysis();
% 5. 运行滑模控制器设计（可选）
% smc_controller_design();
% 6. 运行自适应控制器（可选）
% adaptive_controller();
disp('仿真准备就绪！');
disp('在Simulink中点击运行按钮开始仿真。');
3.2 参数修改指南
%% 参数修改示例
function modify_parameters()
% 修改系统参数并重新运行仿真
global params;
% 修改转子质量
params.m = 1.5; % 从2.0kg改为1.5kg
% 修改控制器参数
params.Kp = 8000; % 增加比例增益
params.Ki = 2000; % 增加积分增益
params.Kd = 80; % 增加微分增益
% 修改干扰大小
params.Fd_step = 20; % 增大干扰力
% 更新工作空间参数
assignin('base', 'params', params);
disp('参数已更新。请重新运行仿真。');
end
四、故障诊断和调试
%% 故障诊断函数
function fault_diagnosis()
% 诊断常见问题
global params;
disp('=== 磁悬浮轴承故障诊断 ===');
% 检查参数合理性
if params.ki <= 0
warning('电流刚度系数应为正值！');
end
if params.kx >= 0
warning('位移刚度系数应为负值（开环不稳定）！');
end
if params.Kp <= 0
warning('比例增益应为正值！');
end
% 检查稳定性条件
% 线性化闭环系统特征方程：m*s² + (Kd*ki)*s + (Kp*ki - kx) = 0
a2 = params.m;
a1 = params.Kd * params.ki;
a0 = params.Kp * params.ki - params.kx;
% Routh-Hurwitz判据
if a2 <= 0 || a1 <= 0 || a0 <= 0
warning('系统可能不稳定！请检查控制器参数。');
else
disp('稳定性条件检查通过。');
end
% 计算自然频率和阻尼比
omega_n = sqrt(a0 / a2);
zeta = a1 / (2 * sqrt(a0 * a2));
disp(['闭环自然频率: ', num2str(omega_n), ' rad/s']);
disp(['闭环阻尼比: ', num2str(zeta)]);
if zeta < 0.4
warning('阻尼比过小，系统响应可能振荡强烈！');
elseif zeta > 1.2
warning('阻尼比过大，系统响应可能过慢！');
end
% 检查电流限制
required_current = (params.m * params.g) / params.ki;
if required_current > params.current_max * 0.8
warning(['所需平衡电流(', num2str(required_current), ...
'A)接近最大电流限制！']);
end
disp('=== 诊断完成 ===');
end
五、扩展功能
5.1 实时监控界面
%% 实时监控GUI
function create_monitoring_gui()
% 创建实时监控界面
fig = uifigure('Name', '磁悬浮轴承监控系统', ...
'Position', [100, 100, 1000, 600]);
% 创建面板
control_panel = uipanel(fig, 'Title', '控制参数', ...
'Position', [20, 400, 300, 180]);
status_panel = uipanel(fig, 'Title', '系统状态', ...
'Position', [340, 400, 300, 180]);
plot_panel = uipanel(fig, 'Title', '实时曲线', ...
'Position', [660, 400, 320, 180]);
main_plot = uipanel(fig, 'Position', [20, 20, 960, 360]);
% 控制参数滑块
uilabel(control_panel, 'Text', '比例增益 Kp:', ...
'Position', [10, 140, 100, 20]);
Kp_slider = uislider(control_panel, ...
'Position', [120, 140, 150, 3], ...
'Limits', [1000, 20000], ...
'Value', 5000);
uilabel(control_panel, 'Text', '积分增益 Ki:', ...
'Position', [10, 100, 100, 20]);
Ki_slider = uislider(control_panel, ...
'Position', [120, 100, 150, 3], ...
'Limits', [100, 10000], ...
'Value', 1000);
% 状态显示
uilabel(status_panel, 'Text', '位移:', ...
'Position', [10, 140, 80, 20]);
displacement_label = uilabel(status_panel, ...
'Text', '0.000 mm', ...
'Position', [100, 140, 100, 20]);
uilabel(status_panel, 'Text', '电流:', ...
'Position', [10, 100, 80, 20]);
current_label = uilabel(status_panel, ...
'Text', '0.00 A', ...
'Position', [100, 100, 100, 20]);
% 按钮
start_btn = uibutton(fig, 'push', ...
'Text', '开始仿真', ...
'Position', [20, 10, 100, 30], ...
'ButtonPushedFcn', @(btn,event) start_simulation());
stop_btn = uibutton(fig, 'push', ...
'Text', '停止仿真', ...
'Position', [140, 10, 100, 30]);
% 回调函数
function start_simulation()
% 更新参数
params.Kp = Kp_slider.Value;
params.Ki = Ki_slider.Value;
assignin('base', 'params', params);
% 运行仿真
sim('axial_magnetic_bearing_complete');
% 更新显示
plot_simulation_results();
end
disp('监控界面已创建。');
end
这个完整的Simulink模型提供了单自由度轴向磁悬浮轴承的全面仿真环境，包括：
- 完整的物理模型：非线性电磁力、转子动力学
- 控制系统：PID控制器，可扩展为高级控制算法
- 传感器模型：包含噪声和饱和特性
- 干扰模型：可配置的干扰力
- 分析工具：性能指标计算、稳定性分析
- 可视化界面：实时监控和结果展示
浙公网安备 33010602011771号