广义预测控制（GPC）实现滞后系统控制 - MATLAB
广义预测控制（GPC）MATLAB实现，用于控制具有滞后的系统。GPC是一种先进的控制策略，特别适合处理滞后系统。
一、GPC基本原理
1.1 系统模型
GPC使用CARIMA（受控自回归积分滑动平均）模型：
A(z⁻¹)y(t) = B(z⁻¹)u(t-1) + C(z⁻¹)e(t)/Δ
其中Δ = 1 - z⁻¹是差分算子。
1.2 GPC算法步骤
- 模型辨识（参数估计）
- Diophantine方程求解
- 预测输出计算
- 控制律求解（最小化代价函数）
二、完整GPC MATLAB程序
%% 广义预测控制（GPC）程序 - 用于滞后系统控制
% 文件名: gpc_control_system.m
clear all;
close all;
clc;
%% 第一部分：参数设置和系统定义
function main_gpc()
% ========== 1.1 控制系统参数 ==========
fprintf('================ GPC控制器设计 ================\n');
% GPC参数
N1 = 1; % 最小预测时域（通常等于系统滞后+1）
N2 = 10; % 最大预测时域
Nu = 3; % 控制时域
lambda = 1.0; % 控制权重
alpha = 0.1; % 柔化因子
% 系统参数（模拟一个滞后系统）
% 系统模型：G(s) = K * e^(-θs) / (τs + 1)
% 离散化后：A(z⁻¹)y(t) = B(z⁻¹)u(t-1)
% 连续时间系统参数
K = 1.5; % 增益
tau = 5.0; % 时间常数
theta = 3.0; % 纯滞后时间
% 离散化参数
Ts = 1.0; % 采样时间
n_samples = 200; % 仿真步数
% ========== 1.2 系统模型生成 ==========
fprintf('生成滞后系统模型...\n');
% 方法1：使用传递函数生成滞后系统
s = tf('s');
G_continuous = K * exp(-theta*s) / (tau*s + 1);
% 离散化
G_discrete = c2d(G_continuous, Ts, 'zoh');
[num, den] = tfdata(G_discrete, 'v');
% 提取多项式系数
A = den(:)'; % A(z⁻¹)系数
B = num(:)'; % B(z⁻¹)系数
% 添加噪声多项式C(z⁻¹) = 1（简化假设）
C = [1];
% 显示系统信息
fprintf('离散系统传递函数:\n');
fprintf('A(z⁻¹) = ');
for i = 1:length(A)
fprintf('%.4f z^{-%d} ', A(i), i-1);
if i < length(A), fprintf('+ '); end
end
fprintf('\n');
fprintf('B(z⁻¹) = ');
for i = 1:length(B)
fprintf('%.4f z^{-%d} ', B(i), i-1);
if i < length(B), fprintf('+ '); end
end
fprintf('\n');
% 计算系统滞后
d = find(B ~= 0, 1) - 1; % B多项式中第一个非零项的阶次
fprintf('系统滞后: d = %d\n', d);
% 调整N1为滞后+1
N1 = d + 1;
fprintf('设置 N1 = d+1 = %d\n', N1);
% ========== 1.3 参考信号生成 ==========
t = (0:n_samples-1)' * Ts;
% 参考轨迹（设定值）
reference = zeros(n_samples, 1);
reference(1:50) = 0;
reference(51:100) = 1;
reference(101:150) = -0.5;
reference(151:200) = 0.5;
% 柔化参考轨迹
ref_filtered = filter([alpha], [1, -(1-alpha)], reference);
% ========== 1.4 初始化变量 ==========
y = zeros(n_samples, 1); % 系统输出
u = zeros(n_samples, 1); % 控制输入
u_min = -2; % 控制输入下限
u_max = 2; % 控制输入上限
du_min = -0.5; % 控制增量下限
du_max = 0.5; % 控制增量上限
% 噪声参数
noise_std = 0.01; % 测量噪声标准差
% 存储预测信息
y_pred = zeros(N2-N1+1, n_samples);
u_opt = zeros(Nu, n_samples);
% ========== 第二部分：GPC核心算法 ==========
fprintf('开始GPC控制仿真...\n');
for k = 1:n_samples
% 当前输出（添加测量噪声）
if k > 1
y(k) = simulate_system(A, B, u, y, k-1) + noise_std * randn;
else
y(k) = 0;
end
% 当前参考值
if k + N2 <= n_samples
R = ref_filtered(k+N1:k+N2);
else
% 如果超出范围，使用最后一个参考值
R = ref_filtered(min(k+N1, n_samples):min(k+N2, n_samples));
if length(R) < N2-N1+1
R = [R; R(end)*ones(N2-N1+1-length(R), 1)];
end
end
% ========== 2.1 Diophantine方程求解 ==========
% 计算阶次
na = length(A) - 1;
nb = length(B) - 1;
nc = length(C) - 1;
% 扩展多项式 E和F
E = cell(N2, 1);
F = cell(N2, 1);
% 计算AΔ = A(z⁻¹) * Δ，其中Δ = 1 - z⁻¹
A_delta = conv(A, [1 -1]);
for j = N1:N2
% 求解Diophantine方程：C = E_j * AΔ + z^{-j} * F_j
[E{j}, F{j}] = solve_diophantine(C, A_delta, j);
end
% ========== 2.2 计算G矩阵和自由响应 ==========
% 计算G矩阵
G = zeros(N2-N1+1, Nu);
for i = 1:N2-N1+1
j = N1 + i - 1;
E_j = E{j};
% 计算G_j = E_j * B
G_j = conv(E_j, B);
% 填充G矩阵
for m = 1:Nu
if m <= length(G_j)
G(i, m) = G_j(m);
end
end
end
% 计算自由响应f
f = zeros(N2-N1+1, 1);
% 过去输入和输出对未来的影响
for i = 1:N2-N1+1
j = N1 + i - 1;
F_j = F{j};
% 计算自由响应分量
f_i = 0;
% 来自F_j的部分
for m = 1:min(length(F_j), k)
f_i = f_i + F_j(m) * y(k-m+1);
end
% 来自过去控制增量的部分
for m = 1:min(length(E{j}), k)
% 计算Δu
if k-m >= 1
du_past = u(k-m+1) - u(k-m);
else
du_past = 0;
end
% 计算E_j * B * Δu
EB = conv(E{j}, B);
if m <= length(EB)
f_i = f_i + EB(m) * du_past;
end
end
f(i) = f_i;
end
% ========== 2.3 求解最优控制 ==========
if k > 1
% 构造Hessian矩阵
H = G' * G + lambda * eye(Nu);
% 构造梯度向量
g = G' * (f - R);
% 无约束优化
du_opt = -H \ g;
% 约束处理（投影到可行域）
du_opt(1) = max(min(du_opt(1), du_max), du_min);
% 计算控制输入
u(k) = u(k-1) + du_opt(1);
% 输入约束
u(k) = max(min(u(k), u_max), u_min);
% 存储最优控制序列
u_opt(:, k) = du_opt;
else
u(k) = 0;
end
% ========== 2.4 预测输出 ==========
% 计算预测输出
y_pred_vec = G * u_opt(1:Nu, max(k,1)) + f;
y_pred(:, k) = y_pred_vec;
end
fprintf('GPC控制仿真完成。\n');
% ========== 第三部分：结果可视化 ==========
plot_gpc_results(t, y, u, reference, ref_filtered, y_pred, N1, N2);
% ========== 第四部分：性能分析 ==========
analyze_performance(t, y, reference, u);
% ========== 第五部分：鲁棒性测试 ==========
test_robustness(A, B, K, tau, theta, Ts, n_samples);
end
%% 辅助函数定义
function y_k = simulate_system(A, B, u, y, k)
% 模拟系统输出
% A: 分母多项式系数
% B: 分子多项式系数
% u: 输入序列
% y: 输出序列
% k: 当前时刻
na = length(A) - 1;
nb = length(B) - 1;
% 计算输出
y_k = 0;
% 来自输入的贡献
for i = 1:min(nb+1, k+1)
if k-i+2 >= 1
y_k = y_k + B(i) * u(k-i+2);
end
end
% 来自输出的贡献（负号因为方程是A*y = B*u）
for i = 2:min(na+1, k+1)
if k-i+2 >= 1
y_k = y_k - A(i) * y(k-i+2);
end
end
% 除以A(1)
if A(1) ~= 0
y_k = y_k / A(1);
end
end
function [E, F] = solve_diophantine(C, A_delta, j)
% 求解Diophantine方程：C = E * AΔ + z^{-j} * F
% C, A_delta: 多项式系数（按z^{-1}的升幂排列）
% j: 预测步长
% 多项式阶次
nc = length(C) - 1;
na_delta = length(A_delta) - 1;
% E的阶次为j-1
E_order = j - 1;
F_order = max(nc, na_delta + j - 1) - j;
% 构造方程组
n_eq = nc + 1;
n_vars = E_order + 1 + F_order + 1;
% 系数矩阵
M = zeros(n_eq, n_vars);
% 填充矩阵
for i = 1:n_eq
% E的贡献
for k = 1:min(E_order+1, i)
if i-k+1 <= na_delta+1
M(i, k) = A_delta(i-k+1);
end
end
% F的贡献（延迟j步）
if i >= j
M(i, E_order+1 + (i-j)) = 1;
end
end
% 右侧向量
b = C(:);
% 求解最小二乘解
x = M \ b;
% 提取E和F
E = x(1:E_order+1)';
F = x(E_order+2:end)';
% 确保多项式正确（截断接近零的系数）
E(abs(E) < 1e-10) = 0;
F(abs(F) < 1e-10) = 0;
end
function plot_gpc_results(t, y, u, reference, ref_filtered, y_pred, N1, N2)
% 绘制GPC控制结果
figure('Position', [100, 100, 1200, 800]);
% 1. 输出响应
subplot(3, 2, 1);
plot(t, y, 'b-', 'LineWidth', 2); hold on;
plot(t, reference, 'r--', 'LineWidth', 1.5);
plot(t, ref_filtered, 'g-.', 'LineWidth', 1);
xlabel('时间 (s)');
ylabel('输出');
title('系统输出响应');
legend('实际输出', '设定值', '柔化设定值', 'Location', 'best');
grid on;
% 2. 控制输入
subplot(3, 2, 2);
stairs(t, u, 'b-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('控制输入');
title('控制输入信号');
grid on;
% 3. 跟踪误差
subplot(3, 2, 3);
error = y - ref_filtered;
plot(t, error, 'r-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('跟踪误差');
title('跟踪误差');
grid on;
hold on;
plot(t, zeros(size(t)), 'k--');
ylim([-0.2, 0.2]);
% 4. 预测输出（特定时刻）
subplot(3, 2, 4);
time_idx = 80; % 选择一个时刻显示预测
if time_idx <= length(t)
pred_horizon = N1:N2;
plot(pred_horizon, y_pred(:, time_idx), 'bo-', 'LineWidth', 1.5, 'MarkerSize', 6);
hold on;
plot(pred_horizon, ref_filtered(time_idx+pred_horizon-1), 'r--', 'LineWidth', 1.5);
xlabel('预测步数');
ylabel('输出');
title(sprintf('时刻 t=%.1f 的预测输出', t(time_idx)));
legend('预测输出', '参考轨迹', 'Location', 'best');
grid on;
end
% 5. 控制输入变化率
subplot(3, 2, 5);
du = diff(u);
plot(t(2:end), du, 'm-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('Δu');
title('控制输入变化率');
grid on;
hold on;
plot(t, zeros(size(t)), 'k--');
% 6. 相位图
subplot(3, 2, 6);
plot(y(2:end), du, 'b.', 'MarkerSize', 8);
xlabel('输出 y(t)');
ylabel('控制变化 Δu(t)');
title('相位图');
grid on;
sgtitle('广义预测控制（GPC）仿真结果', 'FontSize', 14, 'FontWeight', 'bold');
end
function analyze_performance(t, y, reference, u)
% 分析控制性能
% 计算性能指标
error = y - reference;
% 1. 均方根误差
RMSE = sqrt(mean(error.^2));
% 2. 最大绝对误差
max_error = max(abs(error));
% 3. 积分绝对误差
IAE = trapz(t, abs(error));
% 4. 积分平方误差
ISE = trapz(t, error.^2);
% 5. 控制能量
control_energy = trapz(t, u.^2);
% 6. 控制变化能量
du = diff(u);
control_variation_energy = trapz(t(2:end), du.^2);
% 7. 超调量
[overshoot, overshoot_idx] = max(y - reference);
overshoot_percentage = (overshoot / max(abs(reference))) * 100;
% 8. 调节时间（进入±5%误差带的时间）
ref_range = max(reference) - min(reference);
error_band = 0.05 * ref_range;
settling_time_idx = find(abs(error) < error_band, 1);
if ~isempty(settling_time_idx)
settling_time = t(settling_time_idx);
else
settling_time = t(end);
end
% 显示性能指标
fprintf('\n========== 控制性能分析 ==========\n');
fprintf('均方根误差 (RMSE): %.4f\n', RMSE);
fprintf('最大绝对误差: %.4f\n', max_error);
fprintf('积分绝对误差 (IAE): %.4f\n', IAE);
fprintf('积分平方误差 (ISE): %.4f\n', ISE);
fprintf('控制能量: %.4f\n', control_energy);
fprintf('控制变化能量: %.4f\n', control_variation_energy);
fprintf('超调量: %.2f%%\n', overshoot_percentage);
fprintf('调节时间 (5%%): %.2f s\n', settling_time);
% 创建性能分析图
figure('Position', [100, 100, 1000, 600]);
% 误差统计
subplot(2, 3, 1);
histogram(error, 30, 'FaceColor', 'blue', 'EdgeColor', 'black');
xlabel('误差');
ylabel('频次');
title('误差分布直方图');
grid on;
% 自相关函数
subplot(2, 3, 2);
autocorr(error, 50);
title('误差自相关函数');
% 功率谱密度
subplot(2, 3, 3);
[pxx, f] = pwelch(error, [], [], [], 1/(t(2)-t(1)));
semilogy(f, pxx, 'b-', 'LineWidth', 1.5);
xlabel('频率 (Hz)');
ylabel('功率谱密度');
title('误差功率谱密度');
grid on;
% 控制信号分析
subplot(2, 3, 4);
plot(t, u, 'b-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('控制输入');
title('控制输入信号');
grid on;
% 误差随时间变化
subplot(2, 3, 5);
plot(t, error, 'r-', 'LineWidth', 1.5);
hold on;
plot(t, error_band*ones(size(t)), 'g--', 'LineWidth', 1);
plot(t, -error_band*ones(size(t)), 'g--', 'LineWidth', 1);
xlabel('时间 (s)');
ylabel('误差');
title('跟踪误差');
grid on;
% 性能指标汇总
subplot(2, 3, 6);
axis off;
text(0.1, 0.9, '性能指标汇总:', 'FontSize', 12, 'FontWeight', 'bold');
text(0.1, 0.75, sprintf('RMSE: %.4f', RMSE), 'FontSize', 10);
text(0.1, 0.65, sprintf('最大误差: %.4f', max_error), 'FontSize', 10);
text(0.1, 0.55, sprintf('IAE: %.4f', IAE), 'FontSize', 10);
text(0.1, 0.45, sprintf('ISE: %.4f', ISE), 'FontSize', 10);
text(0.1, 0.35, sprintf('超调: %.1f%%', overshoot_percentage), 'FontSize', 10);
text(0.1, 0.25, sprintf('调节时间: %.2f s', settling_time), 'FontSize', 10);
sgtitle('GPC控制性能分析', 'FontSize', 14, 'FontWeight', 'bold');
end
function test_robustness(A_nominal, B_nominal, K_nominal, tau_nominal, theta_nominal, Ts, n_samples)
% 鲁棒性测试：参数变化对控制性能的影响
fprintf('\n========== 鲁棒性测试 ==========\n');
% 测试不同的参数变化
param_variations = [0.5, 0.8, 1.0, 1.2, 1.5]; % 参数变化比例
% 存储性能指标
performance_metrics = zeros(length(param_variations), 4);
figure('Position', [100, 100, 1200, 800]);
colors = {'r', 'g', 'b', 'm', 'c'};
for i = 1:length(param_variations)
factor = param_variations(i);
% 修改系统参数
K = K_nominal * factor;
tau = tau_nominal * factor;
theta = theta_nominal * factor;
fprintf('测试参数变化: %.1f倍\n', factor);
% 重新生成系统
s = tf('s');
G_continuous = K * exp(-theta*s) / (tau*s + 1);
G_discrete = c2d(G_continuous, Ts, 'zoh');
[num, den] = tfdata(G_discrete, 'v');
A = den(:)';
B = num(:)';
C = [1];
% 运行GPC控制（简化版本）
[t, y, u, reference] = run_gpc_simple(A, B, C, Ts, n_samples);
% 计算性能指标
error = y - reference;
RMSE = sqrt(mean(error.^2));
IAE = trapz(t, abs(error));
max_error = max(abs(error));
control_energy = trapz(t, u.^2);
performance_metrics(i, :) = [RMSE, IAE, max_error, control_energy];
% 绘制响应曲线
subplot(2, 3, 1);
plot(t, y, 'Color', colors{i}, 'LineWidth', 1.5, 'DisplayName', sprintf('%.1f倍', factor));
if i == 1
hold on;
plot(t, reference, 'k--', 'LineWidth', 1.5, 'DisplayName', '设定值');
xlabel('时间 (s)');
ylabel('输出');
title('不同参数下的输出响应');
grid on;
legend('show', 'Location', 'best');
end
end
% 绘制性能指标对比
metrics_names = {'RMSE', 'IAE', '最大误差', '控制能量'};
for j = 1:4
subplot(2, 3, j+1);
bar(performance_metrics(:, j), 'FaceColor', [0.2, 0.4, 0.8]);
set(gca, 'XTickLabel', arrayfun(@(x) sprintf('%.1f', x), param_variations, 'UniformOutput', false));
xlabel('参数变化倍数');
ylabel(metrics_names{j});
title(sprintf('%s vs 参数变化', metrics_names{j}));
grid on;
end
sgtitle('GPC鲁棒性测试', 'FontSize', 14, 'FontWeight', 'bold');
end
function [t, y, u, reference] = run_gpc_simple(A, B, C, Ts, n_samples)
% 简化版GPC运行函数（用于鲁棒性测试）
% 初始化
t = (0:n_samples-1)' * Ts;
y = zeros(n_samples, 1);
u = zeros(n_samples, 1);
% 参考信号
reference = zeros(n_samples, 1);
reference(1:50) = 0;
reference(51:100) = 1;
reference(101:150) = -0.5;
reference(151:200) = 0.5;
% GPC参数
N1 = 2;
N2 = 8;
Nu = 2;
lambda = 1.0;
% 主要控制循环
for k = 1:n_samples
% 系统输出
if k > 1
y(k) = simulate_system(A, B, u, y, k-1);
end
if k > 5 % 开始控制
% 计算参考轨迹
R = zeros(N2-N1+1, 1);
for i = 1:N2-N1+1
idx = min(k+N1+i-2, n_samples);
R(i) = reference(idx);
end
% 简化GPC计算（使用固定G矩阵）
G = [B(1), 0; B(1)+A(2)*B(1), B(1)];
% 自由响应（简化）
f = y(k) * ones(N2-N1+1, 1);
% 控制律
H = G' * G + lambda * eye(Nu);
du_opt = -H \ (G' * (f - R));
% 应用控制
if k < n_samples
u(k) = u(k-1) + du_opt(1);
u(k) = max(min(u(k), 2), -2); % 饱和限制
end
end
end
end
%% 第三部分：高级GPC功能
function gpc_with_online_identification()
% 带有在线辨识的GPC
fprintf('\n========== 带在线辨识的GPC ==========\n');
% 参数
n_samples = 300;
Ts = 1.0;
% 真实系统（未知，需要辨识）
A_true = [1, -1.5, 0.7];
B_true = [0, 0.8, 0.3]; % 滞后为1
% 初始模型估计（可能不准确）
A_est = [1, -1.2, 0.6];
B_est = [0, 1.0, 0.2];
C_est = [1];
% 初始化
t = (0:n_samples-1)' * Ts;
y = zeros(n_samples, 1);
u = zeros(n_samples, 1);
% 参考信号
reference = 0.5 * sin(0.05 * t) + 0.3 * cos(0.1 * t);
% RLS参数
na = length(A_true) - 1;
nb = length(B_true) - 1;
nc = length(C_est) - 1;
theta_est = zeros(na + nb + 1, 1); % 参数向量
P = 1000 * eye(na + nb + 1); % 协方差矩阵
lambda_rls = 0.98; % 遗忘因子
% 数据向量
phi = zeros(na + nb + 1, 1);
% GPC参数
N1 = 2;
N2 = 10;
Nu = 3;
lambda_gpc = 0.5;
% 主循环
for k = 1:n_samples
% ========== 1. 系统输出（真实系统） ==========
if k > 1
% 计算真实输出（模拟真实系统）
y_k = 0;
for i = 1:min(nb+1, k)
if k-i+1 >= 1
y_k = y_k + B_true(i) * u(k-i+1);
end
end
for i = 2:min(na+1, k)
if k-i+1 >= 1
y_k = y_k - A_true(i) * y(k-i+1);
end
end
y(k) = y_k + 0.01 * randn; % 添加噪声
end
% ========== 2. 在线辨识（RLS） ==========
if k > max(na, nb) + 1
% 构建数据向量
phi = [-y(k-1:-1:k-na); u(k-1:-1:k-nb-1)];
% RLS更新
K = P * phi / (lambda_rls + phi' * P * phi);
theta_est = theta_est + K * (y(k) - phi' * theta_est);
P = (1/lambda_rls) * (P - K * phi' * P);
% 提取估计的参数
A_est = [1; theta_est(1:na)]';
B_est = [0; theta_est(na+1:end)]'; % 假设滞后为1
end
% ========== 3. GPC控制计算 ==========
if k > max(na, nb) + 5
% 使用估计的模型计算控制
% Diophantine方程求解
A_delta = conv(A_est, [1 -1]);
% 计算G矩阵
G = zeros(N2-N1+1, Nu);
for j = N1:N2
[E_j, F_j] = solve_diophantine(C_est, A_delta, j);
% 计算G_j = E_j * B_est
G_j = conv(E_j, B_est);
i = j - N1 + 1;
for m = 1:min(Nu, length(G_j))
G(i, m) = G_j(m);
end
end
% 计算自由响应
f = zeros(N2-N1+1, 1);
for i = 1:N2-N1+1
j = N1 + i - 1;
[~, F_j] = solve_diophantine(C_est, A_delta, j);
% 自由响应计算
f_i = 0;
for m = 1:min(length(F_j), k)
f_i = f_i + F_j(m) * y(k-m+1);
end
f(i) = f_i;
end
% 参考轨迹
R = zeros(N2-N1+1, 1);
for i = 1:N2-N1+1
idx = min(k+N1+i-2, n_samples);
R(i) = reference(idx);
end
% 求解控制律
H = G' * G + lambda_gpc * eye(Nu);
du_opt = -H \ (G' * (f - R));
% 应用控制
if k < n_samples
u(k) = u(k-1) + du_opt(1);
u(k) = max(min(u(k), 2), -2);
end
else
% 初始阶段使用简单控制
if k > 1
u(k) = 0.1 * (reference(k) - y(k));
end
end
end
% 绘制结果
figure('Position', [100, 100, 1200, 600]);
subplot(2, 2, 1);
plot(t, y, 'b-', 'LineWidth', 1.5); hold on;
plot(t, reference, 'r--', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('输出');
title('带在线辨识的GPC控制');
legend('系统输出', '参考信号', 'Location', 'best');
grid on;
subplot(2, 2, 2);
stairs(t, u, 'g-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('控制输入');
title('控制输入');
grid on;
subplot(2, 2, 3);
error = y - reference;
plot(t, error, 'r-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('跟踪误差');
title('跟踪误差');
grid on;
subplot(2, 2, 4);
% 显示参数收敛
na_est = na;
nb_est = nb;
theta_true = [A_true(2:end)'; B_true(2:end)'];
plot(1:length(theta_est), theta_est, 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10);
hold on;
plot(1:length(theta_true), theta_true, 'r--', 'LineWidth', 1.5);
xlabel('参数索引');
ylabel('参数值');
title('参数估计收敛');
legend('估计值', '真实值', 'Location', 'best');
grid on;
sgtitle('带在线辨识的自适应GPC控制', 'FontSize', 14, 'FontWeight', 'bold');
end
%% 第四部分：GPC控制器设计工具
function gpc_design_tool()
% GPC控制器设计工具
fprintf('\n========== GPC控制器设计工具 ==========\n');
% 用户选择系统类型
fprintf('选择系统类型:\n');
fprintf('1. 一阶滞后系统\n');
fprintf('2. 二阶滞后系统\n');
fprintf('3. 高阶滞后系统\n');
fprintf('4. 非最小相位系统\n');
fprintf('5. 自定义系统\n');
choice = 1; % 默认选择
% 根据选择生成系统
switch choice
case 1
% 一阶滞后系统
K = 1.0;
tau = 5.0;
theta = 2.0;
sys_type = '一阶滞后系统';
case 2
% 二阶滞后系统
K = 1.0;
tau1 = 3.0;
tau2 = 7.0;
theta = 1.0;
sys_type = '二阶滞后系统';
case 3
% 高阶滞后系统
K = 1.0;
tau = [2, 4, 6];
theta = 3.0;
sys_type = '高阶滞后系统';
case 4
% 非最小相位系统
sys_type = '非最小相位系统';
case 5
% 自定义系统
sys_type = '自定义系统';
end
fprintf('选择的系统: %s\n', sys_type);
% GPC参数设计
fprintf('\nGPC参数设计:\n');
% 基于系统特性设计参数
if choice == 1
% 一阶系统
N2 = ceil(5 * tau / Ts); % 覆盖5倍时间常数
Nu = max(2, ceil(N2/3));
lambda = 0.5;
elseif choice == 2
% 二阶系统
N2 = ceil(8 * max(tau1, tau2) / Ts);
Nu = max(3, ceil(N2/4));
lambda = 1.0;
else
% 默认参数
N2 = 10;
Nu = 3;
lambda = 1.0;
end
fprintf('预测时域 N2 = %d\n', N2);
fprintf('控制时域 Nu = %d\n', Nu);
fprintf('控制权重 λ = %.2f\n', lambda);
% 显示稳定性分析
analyze_stability(N2, Nu, lambda);
% 生成控制器代码
generate_gpc_code(N1, N2, Nu, lambda);
end
function analyze_stability(N2, Nu, lambda)
% 稳定性分析
fprintf('\n稳定性分析:\n');
% 计算条件数（粗略的稳定性指标）
% 对于GPC，条件数越小通常稳定性越好
% 模拟G矩阵（假设为单位矩阵简化计算）
G = randn(N2, Nu);
H = G' * G + lambda * eye(Nu);
cond_number = cond(H);
fprintf('Hessian矩阵条件数: %.2f\n', cond_number);
if cond_number < 1000
fprintf('稳定性: 良好\n');
elseif cond_number < 10000
fprintf('稳定性: 中等\n');
else
fprintf('稳定性: 需要注意\n');
fprintf('建议: 增加λ或减少Nu\n');
end
% 极点分析（概念性）
fprintf('\n极点分析:\n');
fprintf('GPC的稳定性取决于闭环极点的位置\n');
fprintf('适当的选择N2, Nu, λ可以保证稳定性\n');
end
function generate_gpc_code(N1, N2, Nu, lambda)
% 生成GPC控制器代码
fprintf('\n========== GPC控制器代码生成 ==========\n');
code = sprintf([...
'%% GPC控制器函数\n' ...
'function [u_opt, du_opt] = gpc_controller(y, u_prev, reference, A, B, C)\n' ...
'%% GPC控制器核心函数\n' ...
'%% 输入:\n' ...
'%% y - 当前输出\n' ...
'%% u_prev - 上一时刻控制输入\n' ...
'%% reference - 参考信号向量\n' ...
'%% A, B, C - 系统多项式\n' ...
'%% 输出:\n' ...
'%% u_opt - 最优控制输入\n' ...
'%% du_opt - 最优控制增量\n' ...
'\n' ...
'%% GPC参数\n' ...
'N1 = %d; %% 最小预测时域\n' ...
'N2 = %d; %% 最大预测时域\n' ...
'Nu = %d; %% 控制时域\n' ...
'lambda = %.2f; %% 控制权重\n' ...
'\n' ...
'%% 1. Diophantine方程求解\n' ...
'A_delta = conv(A, [1 -1]);\n' ...
'G = zeros(N2-N1+1, Nu);\n' ...
'f = zeros(N2-N1+1, 1);\n' ...
'\n' ...
'for j = N1:N2\n' ...
' %% 求解Diophantine方程\n' ...
' [E_j, F_j] = solve_diophantine(C, A_delta, j);\n' ...
' \n' ...
' %% 计算G矩阵元素\n' ...
' G_j = conv(E_j, B);\n' ...
' i = j - N1 + 1;\n' ...
' for m = 1:min(Nu, length(G_j))\n' ...
' G(i, m) = G_j(m);\n' ...
' end\n' ...
' \n' ...
' %% 计算自由响应\n' ...
' f_i = 0;\n' ...
' for m = 1:min(length(F_j), length(y))\n' ...
' f_i = f_i + F_j(m) * y(end-m+1);\n' ...
' end\n' ...
' f(i) = f_i;\n' ...
'end\n' ...
'\n' ...
'%% 2. 构造参考轨迹\n' ...
'R = reference(N1:N2);\n' ...
'\n' ...
'%% 3. 求解最优控制\n' ...
'H = G'' * G + lambda * eye(Nu);\n' ...
'du_opt_vec = -H \\ (G'' * (f - R));\n' ...
'\n' ...
'%% 4. 应用控制\n' ...
'du_opt = du_opt_vec(1);\n' ...
'u_opt = u_prev + du_opt;\n' ...
'\n' ...
'%% 5. 约束处理（可选）\n' ...
'u_max = 2.0;\n' ...
'u_min = -2.0;\n' ...
'du_max = 0.5;\n' ...
'du_min = -0.5;\n' ...
'\n' ...
'du_opt = max(min(du_opt, du_max), du_min);\n' ...
'u_opt = max(min(u_opt, u_max), u_min);\n' ...
'end\n' ...
], N1, N2, Nu, lambda);
fprintf('生成的GPC控制器代码:\n');
fprintf('%s\n', code);
% 保存到文件
filename = 'gpc_controller.m';
fid = fopen(filename, 'w');
fprintf(fid, '%s', code);
fclose(fid);
fprintf('代码已保存到文件: %s\n', filename);
end
%% 主程序运行
% 运行主GPC程序
fprintf('开始运行GPC控制程序...\n');
main_gpc();
% 运行带在线辨识的GPC
% gpc_with_online_identification();
% 运行GPC设计工具
% gpc_design_tool();
fprintf('\n程序运行完成！\n');
三、程序使用说明
3.1 基本使用方法
% 1. 直接运行主程序
run('gpc_control_system.m');
% 2. 自定义系统参数
% 修改以下参数：
K = 2.0; % 系统增益
tau = 3.0; % 时间常数
theta = 2.0; % 滞后时间
Ts = 0.5; % 采样时间
% 3. 调整GPC参数
N1 = 2; % 最小预测时域
N2 = 15; % 最大预测时域
Nu = 4; % 控制时域
lambda = 0.8; % 控制权重
3.2 应用于不同系统
% 对于不同的滞后系统，修改系统模型：
% 例1：一阶滞后系统
sys = tf(1.5, [5, 1], 'InputDelay', 3);
% 例2：二阶滞后系统
sys = tf(2.0, conv([2, 1], [3, 1]), 'InputDelay', 2);
% 例3：非最小相位系统
sys = tf([-0.5, 1], [4, 1], 'InputDelay', 1);
参考代码 gpc预测控制，广义预测控制m程序，实现对滞后系统的控制 www.youwenfan.com/contentcno/96489.html
四、GPC算法特点
- 预测能力：基于模型预测未来输出
- 约束处理：可处理输入输出约束
- 滞后补偿：内置滞后处理能力
- 鲁棒性：对模型失配有较好鲁棒性
- 多目标优化：可同时优化跟踪性能和控制能量
五、性能调优建议
-
预测时域选择：
- N2应覆盖系统主要动态
- 对于滞后系统，N1 = d+1（d为滞后）
-
控制时域选择：
- Nu通常为2-5
- 较小的Nu减少计算量但可能降低性能
-
控制权重调整：
- λ增大：控制更平滑，但响应变慢
- λ减小：响应更快，但控制可能剧烈
-
柔化因子：
- α=0：无柔化，跟踪快但可能超调
- α增大：柔化增强，跟踪平稳
这个GPC程序包提供了完整的滞后系统控制解决方案，包括：
- 基础GPC算法实现
- 在线参数辨识
- 鲁棒性分析
- 性能评估工具
- 代码生成功能
用户可以根据具体应用需求调整参数和修改系统模型。
浙公网安备 33010602011771号