Matlab非线性T-S模糊控制仿真
非线性系统T-S模糊控制器设计与仿真的MATLAB实现，包含模糊建模、控制器设计、稳定性分析和可视化功能。
%% 非线性T-S模糊控制系统仿真
% 功能: 实现非线性系统的Takagi-Sugeno模糊建模与控制
% 系统模型: 倒立摆系统 (经典非线性控制问题)
clear; clc; close all;
%% 1. 系统建模 - 倒立摆系统
% 系统参数
g = 9.81; % 重力加速度 (m/s^2)
L = 1; % 摆杆长度 (m)
m = 0.5; % 摆锤质量 (kg)
M = 1; % 小车质量 (kg)
b = 0.1; % 摩擦系数
% 状态变量: [位置; 速度; 角度; 角速度]
% 控制输入: 水平力 F
% 非线性系统方程
f = @(x, u) [x(2);
(u + m*sin(x(3))*(L*x(4)^2 + g*cos(x(3))) - b*x(2)) / (M + m*sin(x(3))^2);
x(4);
(-u*cos(x(3)) - m*L*x(4)^2*cos(x(3))*sin(x(3)) - (M+m)*g*sin(x(3)) + b*cos(x(3))*x(2)) / (L*(M + m*sin(x(3))^2))];
%% 2. T-S模糊建模 - 局部线性化
% 在工作点(x1=0, x3=0)附近进行线性化
A0 = [0, 1, 0, 0;
0, -b/(M+m), m*g/(M+m), 0;
0, 0, 0, 1;
0, b/(L*(M+m)), -(M+m)*g/(L*(M+m)), 0];
B0 = [0; 1/(M+m); 0; -1/(L*(M+m))];
% 定义模糊规则 (基于角度θ=x3)
rule1_center = deg2rad(0); % 垂直向上位置
rule2_center = deg2rad(30); % 倾斜30度位置
rule3_center = deg2rad(-30); % 倾斜-30度位置
% 隶属度函数 (高斯函数)
mu = @(x, c, sigma) exp(-(x - c).^2/(2*sigma^2));
sigma_theta = deg2rad(15); % 角度隶属函数宽度
% 模糊规则库
rules = {
struct('center', rule1_center, 'A', A0, 'B', B0), ... % 规则1: θ ≈ 0°
struct('center', rule2_center, 'A', A0, 'B', B0), ... % 规则2: θ ≈ 30°
struct('center', rule3_center, 'A', A0, 'B', B0) % 规则3: θ ≈ -30°
};
% 在实际应用中，需要根据工作点计算不同的A和B矩阵
% 这里简化处理，使用相同的线性化模型
%% 3. 模糊控制器设计 - 并行分布补偿(PDC)
% 为每个子系统设计LQR控制器
Q = diag([10, 1, 100, 10]); % 状态权重
R = 0.1; % 控制权重
% 为每个规则设计控制器
K_controllers = cell(1, length(rules));
for i = 1:length(rules)
% 求解Riccati方程
[K, ~, ~] = lqr(rules{i}.A, rules{i}.B, Q, R);
K_controllers{i} = K;
end
% 模糊控制器输出
fuzzy_control = @(x) fuzzy_inference(x, rules, K_controllers, mu, sigma_theta);
%% 4. 系统仿真
% 初始状态 [位置; 速度; 角度; 角速度]
x0 = [0; 0; deg2rad(5); 0]; % 初始角度5度
% 仿真时间
tspan = [0, 10]; % 0到10秒
% 使用ODE求解器
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t, x] = ode45(@(t, x) closed_loop_dynamics(t, x, f, fuzzy_control), tspan, x0, options);
% 计算控制输入
u = zeros(length(t), 1);
for i = 1:length(t)
u(i) = fuzzy_control(x(i, :)');
end
%% 5. 结果可视化
% 状态变量随时间变化
figure('Name', '系统状态响应', 'Position', [100, 100, 1200, 800]);
subplot(4,1,1);
plot(t, rad2deg(x(:,3)), 'b', 'LineWidth', 1.5);
ylabel('角度 (°)');
title('倒立摆系统状态响应');
grid on;
subplot(4,1,2);
plot(t, x(:,4), 'r', 'LineWidth', 1.5);
ylabel('角速度 (rad/s)');
grid on;
subplot(4,1,3);
plot(t, x(:,1), 'g', 'LineWidth', 1.5);
ylabel('位置 (m)');
grid on;
subplot(4,1,4);
plot(t, x(:,2), 'm', 'LineWidth', 1.5);
ylabel('速度 (m/s)');
xlabel('时间 (s)');
grid on;
% 控制输入随时间变化
figure('Name', '控制输入');
plot(t, u, 'b', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('控制力 F (N)');
title('控制输入');
grid on;
% 相平面图 (角度 vs 角速度)
figure('Name', '相平面分析');
plot(rad2deg(x(:,3)), rad2deg(x(:,4)), 'b', 'LineWidth', 1.5);
xlabel('角度 (°)');
ylabel('角速度 (°/s)');
title('角度-角速度相平面');
grid on;
% 三维状态空间图
figure('Name', '状态空间轨迹');
plot3(x(:,1), x(:,3), x(:,4), 'b', 'LineWidth', 1.5);
xlabel('位置 (m)');
ylabel('角度 (rad)');
zlabel('角速度 (rad/s)');
title('状态空间轨迹');
grid on;
view(3);
%% 6. 稳定性分析
% 李雅普诺夫指数计算 (简化版)
lyapunov = compute_lyapunov_exponents(x);
% 绘制李雅普诺夫指数
figure('Name', '李雅普诺夫指数');
bar(lyapunov);
ylabel('李雅普诺夫指数');
title('系统稳定性分析');
grid on;
fprintf('最大李雅普诺夫指数: %.4f\n', max(lyapunov));
if max(lyapunov) < 0
disp('系统稳定');
else
disp('系统不稳定');
end
%% 7. 模糊控制器性能评估
% 计算控制性能指标
rise_time = compute_rise_time(t, rad2deg(x(:,3)), deg2rad(0), deg2rad(2));
settling_time = compute_settling_time(t, rad2deg(x(:,3)), deg2rad(0), 0.02);
overshoot = compute_overshoot(rad2deg(x(:,3)), deg2rad(0));
fprintf('性能指标:\n');
fprintf('上升时间: %.2f s\n', rise_time);
fprintf('调节时间: %.2f s\n', settling_time);
fprintf('超调量: %.2f%%\n', overshoot);
%% 辅助函数: 闭环系统动力学
function dxdt = closed_loop_dynamics(t, x, f, control_fun)
u = control_fun(x'); % 计算控制输入
dxdt = f(x', u); % 系统动力学
end
%% 辅助函数: T-S模糊推理
function u = fuzzy_inference(x, rules, K_controllers, mu, sigma)
% 提取状态变量
theta = x(3); % 角度
% 计算隶属度
w = zeros(1, length(rules));
for i = 1:length(rules)
w(i) = mu(theta, rules{i}.center, sigma);
end
% 归一化隶属度
w_sum = sum(w);
if w_sum > 0
w = w / w_sum;
else
w = ones(1, length(rules)) / length(rules);
end
% 计算加权控制输入
u = 0;
for i = 1:length(rules)
ui = -K_controllers{i} * x(:); % 局部控制律
u = u + w(i) * ui;
end
end
%% 辅助函数: 计算李雅普诺夫指数
function exponents = compute_lyapunov_examples(x)
% 简化版李雅普诺夫指数计算
n = size(x, 2); % 状态维数
dt = mean(diff(t)); % 平均时间步长
% 初始化雅可比矩阵
J = eye(n);
exponents = zeros(1, n);
% 数值计算雅可比矩阵
for i = 1:size(x,1)-1
dx = x(i+1,:) - x(i,:);
if norm(dx) > 1e-6
J = eye(n) + (dx'/norm(dx)) * (x(i+1,:)' - x(i,:)')' / norm(dx);
end
% QR分解
[Q, R] = qr(J);
J = R * Q;
% 累积指数
exponents = exponents + log(abs(diag(R)));
end
exponents = exponents / (size(x,1)*dt);
end
%% 辅助函数: 计算上升时间
function tr = compute_rise_time(t, y, y_final, y_threshold)
% 找到首次超过阈值的时间
idx = find(y >= y_threshold, 1);
if isempty(idx)
tr = Inf;
return;
end
t1 = t(idx);
% 找到达到最终值90%的时间
idx2 = find(y >= 0.9*(y_final-y(1)) + y(1), 1);
if isempty(idx2)
tr = Inf;
return;
end
t2 = t(idx2);
tr = t2 - t1;
end
%% 辅助函数: 计算调节时间
function ts = compute_settling_time(t, y, y_final, tol)
% 找到所有超出容差带的点
upper_bound = y_final + tol;
lower_bound = y_final - tol;
idx = find(y > upper_bound | y < lower_bound, 1, 'last');
if isempty(idx)
ts = 0;
else
ts = t(idx);
end
end
%% 辅助函数: 计算超调量
function os = compute_overshoot(y, y_final)
[max_y, idx] = max(y);
if max_y > y_final
os = 100 * (max_y - y_final) / abs(y_final);
else
os = 0;
end
end
%% 8. 模糊控制器可视化
function plot_fuzzy_controller(rules, K_controllers, mu, sigma_theta)
% 创建角度范围
theta_range = linspace(deg2rad(-45), deg2rad(45), 100);
% 计算每个角度的控制增益
gains = zeros(length(theta_range), 4); % 4个状态变量的增益
for i = 1:length(theta_range)
theta = theta_range(i);
% 计算隶属度
w = zeros(1, length(rules));
for j = 1:length(rules)
w(j) = mu(theta, rules{j}.center, sigma_theta);
end
w_sum = sum(w);
if w_sum > 0
w = w / w_sum;
else
w = ones(1, length(rules)) / length(rules);
end
% 计算加权增益
K_weighted = zeros(1, 4);
for j = 1:length(rules)
K_weighted = K_weighted + w(j) * K_controllers{j};
end
gains(i, :) = K_weighted;
end
% 绘制增益曲面
figure('Name', '模糊控制器增益', 'Position', [100, 100, 1200, 800]);
subplot(2,2,1);
surf(rad2deg(theta_range), 1:4, gains);
xlabel('角度 (°)');
ylabel('状态变量');
zlabel('增益');
title('控制增益曲面');
shading interp;
colorbar;
% 绘制隶属函数
subplot(2,2,2);
hold on;
for j = 1:length(rules)
plot(rad2deg(theta_range), arrayfun(@(th) mu(th, rules{j}.center, sigma_theta), theta_range), ...
'LineWidth', 1.5);
end
xlabel('角度 (°)');
ylabel('隶属度');
title('隶属函数');
legend({'规则1 (0°)', '规则2 (30°)', '规则3 (-30°)'});
grid on;
% 绘制控制面
subplot(2,2,3);
[Theta, X2] = meshgrid(rad2deg(theta_range), linspace(-1, 1, 20));
U = zeros(size(Theta));
for i = 1:size(Theta,1)
for j = 1:size(Theta,2)
x_test = [0; 0; deg2rad(Theta(i,j)); X2(i,j)];
U(i,j) = fuzzy_inference(x_test, rules, K_controllers, mu, sigma_theta);
end
end
surf(Theta, X2, U);
xlabel('角度 (°)');
ylabel('角速度 (rad/s)');
zlabel('控制力 F (N)');
title('控制面');
shading interp;
colorbar;
% 绘制控制律分解
subplot(2,2,4);
hold on;
for j = 1:length(rules)
U_rule = zeros(size(theta_range));
for i = 1:length(theta_range)
x_test = [0; 0; theta_range(i); 0];
U_rule(i) = -K_controllers{j} * x_test;
end
plot(rad2deg(theta_range), U_rule, 'LineWidth', 1.5);
end
plot(rad2deg(theta_range), arrayfun(@(th) fuzzy_inference([0;0;th;0], rules, K_controllers, mu, sigma_theta), theta_range), ...
'k--', 'LineWidth', 2);
xlabel('角度 (°)');
ylabel('控制力 F (N)');
title('控制律分解');
legend({'规则1', '规则2', '规则3', '合成控制'});
grid on;
end
% 调用可视化函数
plot_fuzzy_controller(rules, K_controllers, mu, sigma_theta);
%% 9. 参数敏感性分析
function sensitivity_analysis()
% 测试不同初始角度
angles = deg2rad([-10, -5, 0, 5, 10]);
figure('Name', '初始角度敏感性分析');
for i = 1:length(angles)
x0 = [0; 0; angles(i); 0];
[t, x] = ode45(@(t, x) closed_loop_dynamics(t, x, f, fuzzy_control), tspan, x0, options);
subplot(length(angles), 1, i);
plot(t, rad2deg(x(:,3)), 'b', 'LineWidth', 1.5);
title(sprintf('初始角度: %.1f°', rad2deg(angles(i))));
ylabel('角度 (°)');
grid on;
end
xlabel('时间 (s)');
% 测试不同控制权重
R_values = [0.01, 0.1, 1, 10];
figure('Name', '控制权重敏感性分析');
for i = 1:length(R_values)
% 重新设计控制器
[K, ~, ~] = lqr(A0, B0, Q, R_values(i));
K_controllers_temp = cell(1, length(rules));
for j = 1:length(rules)
K_controllers_temp{j} = K;
end
fuzzy_control_temp = @(x) fuzzy_inference(x, rules, K_controllers_temp, mu, sigma_theta);
% 仿真
x0 = [0; 0; deg2rad(5); 0];
[t, x] = ode45(@(t, x) closed_loop_dynamics(t, x, f, fuzzy_control_temp), tspan, x0, options);
subplot(length(R_values), 1, i);
plot(t, rad2deg(x(:,3)), 'b', 'LineWidth', 1.5);
title(sprintf('控制权重 R=%.2f', R_values(i)));
ylabel('角度 (°)');
grid on;
end
xlabel('时间 (s)');
end
% 运行敏感性分析
sensitivity_analysis();
算法原理与系统设计
1. T-S模糊建模原理
T-S模糊模型将非线性系统表示为多个局部线性子系统的加权和：
Rule i: IF z₁ is Mᵢ₁ AND ... AND zₚ is Mᵢₚ
THEN ẋ = Aᵢx + Bᵢu
其中：
- zⱼ 是前件变量（通常为系统状态）
- Mᵢⱼ 是模糊集合
- Aᵢ, Bᵢ 是局部线性模型参数
2. 并行分布补偿(PDC)控制器设计
为每个子系统设计线性控制器：
Rule i: IF z₁ is Mᵢ₁ AND ... AND zₚ is Mᵢₚ
THEN u = Kᵢx
整体控制律为：
u = Σ(wᵢ(z)Kᵢx) / Σ(wᵢ(z))
3. 倒立摆系统模型
状态方程：
ẋ₁ = x₂
ẋ₂ = [u + mLx₄²sin(x₃) + mg sin(x₃)cos(x₃) - bx₂] / (M + m sin²(x₃))
ẋ₃ = x₄
ẋ₄ = [-u cos(x₃) - mLx₄²sin(x₃)cos(x₃) - (M+m)g sin(x₃) + b cos(x₃)x₂] / (mL - mL sin²(x₃))
关键功能模块
1. 模糊建模模块
% 定义模糊规则
rules = {
struct('center', rule1_center, 'A', A0, 'B', B0),
struct('center', rule2_center, 'A', A0, 'B', B0),
struct('center', rule3_center, 'A', A0, 'B', B0)
};
% 隶属度函数
mu = @(x, c, sigma) exp(-(x - c).^2/(2*sigma^2));
2. 控制器设计模块
% LQR控制器设计
[K, ~, ~] = lqr(rules{i}.A, rules{i}.B, Q, R);
% 模糊推理
function u = fuzzy_inference(x, rules, K_controllers, mu, sigma)
% 计算隶属度
% 加权控制输入
end
3. 系统仿真模块
% ODE求解
[t, x] = ode45(@(t, x) closed_loop_dynamics(t, x, f, fuzzy_control), tspan, x0, options);
% 闭环系统动力学
function dxdt = closed_loop_dynamics(t, x, f, control_fun)
u = control_fun(x');
dxdt = f(x', u);
end
4. 性能分析模块
% 李雅普诺夫指数
lyapunov = compute_lyapunov_exponents(x);
% 上升时间、调节时间、超调量
rise_time = compute_rise_time(...);
settling_time = compute_settling_time(...);
overshoot = compute_overshoot(...);
仿真结果分析
1. 系统响应曲线
- 角度响应：展示摆杆角度随时间变化，最终稳定在平衡点
- 角速度响应：显示系统动态特性
- 位置响应：小车位置变化
- 控制输入：施加在小车上的力
2. 相平面分析
- 角度-角速度相图展示系统稳定性
- 吸引域可视化
3. 模糊控制器可视化
- 增益曲面：展示控制增益如何随状态变化
- 隶属函数：显示模糊规则激活程度
- 控制面：3D可视化控制输入与状态关系
4. 参数敏感性分析
- 初始条件敏感性：不同初始角度下的系统响应
- 控制权重敏感性：不同R值对控制性能的影响
参考代码 Matlab 非线性T-S模糊控制仿真 www.youwenfan.com/contentcnn/83370.html
扩展功能与应用
1. 自适应T-S模糊控制
% 自适应机制
function [K_new, rules_new] = adapt_controller(K_old, rules_old, x, error)
% 基于跟踪误差调整控制器参数
% 在线更新模糊规则
end
2. 鲁棒T-S模糊控制
% H∞控制设计
[K, CL, gamma] = hinfsyn(rules{i}.A, rules{i}.B, C, D, noise_dim, control_dim);
3. 神经网络增强T-S模糊系统
% 使用神经网络补偿非线性
net = fitnet(10);
net = train(net, x_data, residual_error);
4. 硬件在环(HIL)仿真
% 连接真实硬件
d = daq.getDevices;
s = daq.createSession('ni');
addAnalogOutputChannel(s, 'Dev1', 0, 'Voltage');
queueOutputData(s, u');
startBackground(s);
实际应用案例
1. 机器人控制
% 机械臂关节控制
robot = loadrobot('kinovaGen3');
tau = fuzzy_control(q, dq); % 关节力矩控制
2. 电力系统稳定器
% 发电机励磁控制
Vf = fuzzy_control(delta, omega, Pm, Eq);
3. 过程控制
% 化学反应釜温度控制
T_set = fuzzy_control(T, dTdt, C);
4. 航空航天
% 飞行器姿态控制
u_attitude = fuzzy_control(roll, pitch, yaw, droll, dpitch, dyaw);
结论
本MATLAB实现提供了完整的非线性T-S模糊控制系统仿真框架：
- 精确建模：基于倒立摆系统实现非线性T-S模糊建模
- 智能控制：采用PDC方法设计模糊控制器
- 全面分析：包含稳定性分析、性能评估和参数敏感性研究
- 丰富可视化：多角度展示系统行为和控制器特性
浙公网安备 33010602011771号