广义预测控制（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.mclear 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('+ '); endendfprintf('\n');fprintf('B(z⁻¹) = ');for i = 1:length(B)fprintf('%.4f z^{-%d} ', B(i), i-1);if i < length(B), fprintf('+ '); endendfprintf('\n');% 计算系统滞后d = find(B ~= 0, 1) - 1; % B多项式中第一个非零项的阶次fprintf('系统滞后: d = %d\n', d);% 调整N1为滞后+1N1 = 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 > 1y(k) = simulate_system(A, B, u, y, k-1) + noise_std * randn;elsey(k) = 0;end% 当前参考值if k + N2 <= n_samplesR = 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+1R = [R; R(end)*ones(N2-N1+1-length(R), 1)];endend% ========== 2.1 Diophantine方程求解 ==========% 计算阶次na = length(A) - 1;nb = length(B) - 1;nc = length(C) - 1;% 扩展多项式 E和FE = 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+1j = N1 + i - 1;E_j = E{j};% 计算G_j = E_j * BG_j = conv(E_j, B);% 填充G矩阵for m = 1:Nuif m <= length(G_j)G(i, m) = G_j(m);endendend% 计算自由响应ff = zeros(N2-N1+1, 1);% 过去输入和输出对未来的影响for i = 1:N2-N1+1j = 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)% 计算Δuif k-m >= 1du_past = u(k-m+1) - u(k-m);elsedu_past = 0;end% 计算E_j * B * ΔuEB = conv(E{j}, B);if m <= length(EB)f_i = f_i + EB(m) * du_past;endendf(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;elseu(k) = 0;end% ========== 2.4 预测输出 ==========% 计算预测输出y_pred_vec = G * u_opt(1:Nu, max(k,1)) + f;y_pred(:, k) = y_pred_vec;endfprintf('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 >= 1y_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 >= 1y_k = y_k - A(i) * y(k-i+2);end
end% 除以A(1)
if A(1) ~= 0y_k = y_k / A(1);
end
endfunction [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+1M(i, k) = A_delta(i-k+1);endend% F的贡献（延迟j步）if i >= jM(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;
endfunction 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');
endfunction 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);
elsesettling_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');
endfunction 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 == 1hold 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:4subplot(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;
endsgtitle('GPC鲁棒性测试', 'FontSize', 14, 'FontWeight', 'bold');
endfunction [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 > 1y(k) = simulate_system(A, B, u, y, k-1);endif k > 5 % 开始控制% 计算参考轨迹R = zeros(N2-N1+1, 1);for i = 1:N2-N1+1idx = 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_samplesu(k) = u(k-1) + du_opt(1);u(k) = max(min(u(k), 2), -2); % 饱和限制endend
end
end%% 第三部分：高级GPC功能function gpc_with_online_identification()
% 带有在线辨识的GPCfprintf('\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 >= 1y_k = y_k + B_true(i) * u(k-i+1);endendfor i = 2:min(na+1, k)if k-i+1 >= 1y_k = y_k - A_true(i) * y(k-i+1);endendy(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)]'; % 假设滞后为1end% ========== 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_estG_j = conv(E_j, B_est);i = j - N1 + 1;for m = 1:min(Nu, length(G_j))G(i, m) = G_j(m);endend% 计算自由响应f = zeros(N2-N1+1, 1);for i = 1:N2-N1+1j = 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);endf(i) = f_i;end% 参考轨迹R = zeros(N2-N1+1, 1);for i = 1:N2-N1+1idx = 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_samplesu(k) = u(k-1) + du_opt(1);u(k) = max(min(u(k), 2), -2);endelse% 初始阶段使用简单控制if k > 1u(k) = 0.1 * (reference(k) - y(k));endend
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 choicecase 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 = '自定义系统';
endfprintf('选择的系统: %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;
endfprintf('预测时域 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);
endfunction 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 < 1000fprintf('稳定性: 良好\n');
elseif cond_number < 10000fprintf('稳定性: 中等\n');
elsefprintf('稳定性: 需要注意\n');fprintf('建议: 增加λ或减少Nu\n');
end% 极点分析（概念性）
fprintf('\n极点分析:\n');
fprintf('GPC的稳定性取决于闭环极点的位置\n');
fprintf('适当的选择N2, Nu, λ可以保证稳定性\n');
endfunction 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算法实现
- 在线参数辨识
- 鲁棒性分析
- 性能评估工具
- 代码生成功能
用户可以根据具体应用需求调整参数和修改系统模型。