多变量隐式广义预测控制的MATLAB实现
多变量隐式广义预测控制（Generalized Predictive Control, GPC）是一种先进的控制策略，特别适用于多输入多输出（MIMO）系统。
理论基础
广义预测控制基于受控自回归积分滑动平均（CARIMA）模型：
\[A(z^{-1})y(t) = B(z^{-1})u(t-1) + C(z^{-1})\frac{\xi(t)}{\Delta}
\]
其中：
- \(A(z^{-1})\), \(B(z^{-1})\), \(C(z^{-1})\) 是多项式矩阵
- \(\Delta = 1 - z^{-1}\) 是差分算子
- \(\xi(t)\) 是不相关随机序列
对于多变量系统，我们需要处理多项式矩阵而不是标量多项式。
MATLAB实现
主控制器类
classdef MultivariableGPC < handle
% MULTIVARIABLEGPC 多变量隐式广义预测控制器
properties
% 系统参数
na % A多项式阶数
nb % B多项式阶数
nc % C多项式阶数 (通常设为1)
nu % 输入变量数
ny % 输出变量数
% GPC参数
N1 % 最小预测时域
N2 % 最大预测时域
Nu % 控制时域
lambda % 控制权重矩阵
alpha % 输出权重矩阵
% 系统模型
A_poly % A多项式矩阵系数
B_poly % B多项式矩阵系数
C_poly % C多项式矩阵系数
% 数据存储
past_inputs
past_outputs
past_errors
% 约束
u_min % 输入最小值
u_max % 输入最大值
du_min % 输入变化率最小值
du_max % 输入变化率最大值
% 优化选项
use_constraints
solver_options
end
methods
function obj = MultivariableGPC(na, nb, nc, nu, ny, N1, N2, Nu, lambda, alpha)
% MULTIVARIABLEGPC 构造函数
obj.na = na;
obj.nb = nb;
obj.nc = nc;
obj.nu = nu;
obj.ny = ny;
obj.N1 = N1;
obj.N2 = N2;
obj.Nu = Nu;
obj.lambda = lambda;
obj.alpha = alpha;
% 初始化多项式矩阵
obj.A_poly = zeros(ny, ny, na+1);
obj.B_poly = zeros(ny, nu, nb+1);
obj.C_poly = zeros(ny, ny, nc+1);
% 设置默认单位矩阵
for i = 1:ny
obj.C_poly(i,i,1) = 1;
end
% 初始化数据存储
max_history = max([na, nb, nc]);
obj.past_inputs = zeros(nu, max_history);
obj.past_outputs = zeros(ny, max_history);
obj.past_errors = zeros(ny, max_history);
% 默认无约束
obj.use_constraints = false;
obj.u_min = -inf(nu, 1);
obj.u_max = inf(nu, 1);
obj.du_min = -inf(nu, 1);
obj.du_max = inf(nu, 1);
% 设置优化选项
obj.solver_options = optimoptions('quadprog', ...
'Display', 'none', ...
'Algorithm', 'interior-point-convex');
end
function set_model(obj, A_coeffs, B_coeffs, C_coeffs)
% SET_MODEL 设置系统模型系数
obj.A_poly = A_coeffs;
obj.B_poly = B_coeffs;
obj.C_poly = C_coeffs;
end
function set_constraints(obj, u_min, u_max, du_min, du_max)
% SET_CONSTRAINTS 设置输入约束
obj.u_min = u_min;
obj.u_max = u_max;
obj.du_min = du_min;
obj.du_max = du_max;
obj.use_constraints = true;
end
function update(obj, y, u)
% UPDATE 更新控制器状态
obj.past_outputs = circshift(obj.past_outputs, [0, -1]);
obj.past_outputs(:, end) = y;
obj.past_inputs = circshift(obj.past_inputs, [0, -1]);
obj.past_inputs(:, end) = u;
% 计算预测误差
y_pred = obj.predict_one_step();
error = y - y_pred;
obj.past_errors = circshift(obj.past_errors, [0, -1]);
obj.past_errors(:, end) = error;
end
function y_pred = predict_one_step(obj)
% PREDICT_ONE_STEP 一步预测
y_pred = zeros(obj.ny, 1);
% A多项式贡献
for i = 1:obj.na
if size(obj.past_outputs, 2) >= i
y_pred = y_pred - squeeze(obj.A_poly(:, :, i+1)) * obj.past_outputs(:, end-i+1);
end
end
% B多项式贡献
for i = 1:obj.nb
if size(obj.past_inputs, 2) >= i
y_pred = y_pred + squeeze(obj.B_poly(:, :, i+1)) * obj.past_inputs(:, end-i+1);
end
end
% C多项式贡献 (误差项)
for i = 1:obj.nc
if size(obj.past_errors, 2) >= i
y_pred = y_pred + squeeze(obj.C_poly(:, :, i+1)) * obj.past_errors(:, end-i+1);
end
end
end
function [G, F] = calculate_step_response_matrices(obj)
% CALCULATE_STEP_RESPONSE_MATRICES 计算阶跃响应矩阵
np = obj.N2 - obj.N1 + 1;
G = zeros(obj.ny * np, obj.nu * obj.Nu);
F = zeros(obj.ny * np, 1);
% 创建扩展的历史数据
extended_inputs = [zeros(obj.nu, obj.N2), obj.past_inputs];
extended_outputs = [zeros(obj.ny, obj.N2), obj.past_outputs];
extended_errors = [zeros(obj.ny, obj.N2), obj.past_errors];
% 计算自由响应F
for k = 1:np
N = obj.N1 + k - 1;
f = zeros(obj.ny, 1);
% A多项式贡献
for i = 1:obj.na
if N - i > 0 && N - i <= size(extended_outputs, 2)
f = f - squeeze(obj.A_poly(:, :, i+1)) * extended_outputs(:, end - (N-i) + 1);
end
end
% B多项式贡献 (过去输入)
for i = 1:obj.nb
if N - i > 0 && N - i <= size(extended_inputs, 2)
f = f + squeeze(obj.B_poly(:, :, i+1)) * extended_inputs(:, end - (N-i) + 1);
end
end
% C多项式贡献
for i = 1:obj.nc
if N - i > 0 && N - i <= size(extended_errors, 2)
f = f + squeeze(obj.C_poly(:, :, i+1)) * extended_errors(:, end - (N-i) + 1);
end
end
F((k-1)*obj.ny+1:k*obj.ny) = f;
end
% 计算阶跃响应矩阵G
for j = 1:obj.Nu
for k = 1:np
N = obj.N1 + k - 1;
if j <= N
% 计算第j步输入对第N步输出的影响
g = zeros(obj.ny, obj.nu);
for i = 1:min(obj.nb, N-j+1)
if i+1 <= size(obj.B_poly, 3)
g = g + squeeze(obj.B_poly(:, :, i+1));
end
end
row_start = (k-1)*obj.ny + 1;
row_end = k*obj.ny;
col_start = (j-1)*obj.nu + 1;
col_end = j*obj.nu;
G(row_start:row_end, col_start:col_end) = g;
end
end
end
end
function [du_opt, cost] = compute_control(obj, y_ref)
% COMPUTE_CONTROL 计算最优控制增量
% 计算阶跃响应矩阵
[G, F] = obj.calculate_step_response_matrices();
% 构建参考轨迹
np = obj.N2 - obj.N1 + 1;
Y_ref = zeros(obj.ny * np, 1);
for i = 1:np
Y_ref((i-1)*obj.ny+1:i*obj.ny) = y_ref;
end
% 构建优化问题
H = 2 * (G' * kron(eye(np), obj.alpha) * G + kron(eye(obj.Nu), obj.lambda));
f = 2 * G' * kron(eye(np), obj.alpha) * (F - Y_ref);
if ~obj.use_constraints
% 无约束优化
du_opt = -H \ f;
cost = 0.5 * du_opt' * H * du_opt + f' * du_opt;
else
% 带约束优化
A = [];
b = [];
% 输入变化率约束
if any(isfinite(obj.du_min)) || any(isfinite(obj.du_max))
A_du = [eye(obj.nu * obj.Nu); -eye(obj.nu * obj.Nu)];
b_du = [repmat(obj.du_max, obj.Nu, 1); -repmat(obj.du_min, obj.Nu, 1)];
A = [A; A_du];
b = [b; b_du];
end
% 输入幅值约束（需要预测未来的输入）
if any(isfinite(obj.u_min)) || any(isfinite(obj.u_max))
% 构建累积输入矩阵
M = tril(ones(obj.Nu * obj.nu));
current_u = obj.past_inputs(:, end);
A_u = [M; -M];
b_u = [repmat(obj.u_max - current_u, obj.Nu, 1);
repmat(-obj.u_min + current_u, obj.Nu, 1)];
A = [A; A_u];
b = [b; b_u];
end
% 求解QP问题
[du_opt, cost] = quadprog(H, f, A, b, [], [], [], [], [], obj.solver_options);
end
% 提取第一个控制增量
du_opt = du_opt(1:obj.nu);
end
function u = control(obj, y, y_ref)
% CONTROL 主控制函数
u_prev = obj.past_inputs(:, end);
% 更新控制器状态
obj.update(y, u_prev);
% 计算最优控制增量
du_opt = obj.compute_control(y_ref);
% 应用控制增量
u = u_prev + du_opt;
% 确保输入在约束范围内
if obj.use_constraints
u = max(min(u, obj.u_max), obj.u_min);
end
end
end
end
系统仿真示例
function multivariable_gpc_simulation()
% MULTIVARIABLE_GPC_SIMULATION 多变量GPC仿真示例
% 系统参数
na = 2; nb = 2; nc = 1;
nu = 2; ny = 2;
% GPC参数
N1 = 1; N2 = 10; Nu = 3;
lambda = diag([0.1, 0.1]); % 控制权重
alpha = diag([1, 1]); % 输出权重
% 创建GPC控制器
gpc = MultivariableGPC(na, nb, nc, nu, ny, N1, N2, Nu, lambda, alpha);
% 设置系统模型 (示例: 双输入双输出系统)
A_coeffs = zeros(ny, ny, na+1);
A_coeffs(:,:,1) = eye(ny); % A0
A_coeffs(:,:,2) = [-0.8, 0.1; 0.2, -0.9]; % A1
A_coeffs(:,:,3) = [-0.1, 0.05; 0.1, -0.2]; % A2
B_coeffs = zeros(ny, nu, nb+1);
B_coeffs(:,:,1) = zeros(ny, nu); % B0 (通常为0)
B_coeffs(:,:,2) = [0.5, 0.2; 0.1, 0.6]; % B1
B_coeffs(:,:,3) = [0.2, 0.1; 0.05, 0.3]; % B2
C_coeffs = zeros(ny, ny, nc+1);
C_coeffs(:,:,1) = eye(ny); % C0
gpc.set_model(A_coeffs, B_coeffs, C_coeffs);
% 设置约束
u_min = [-1; -1];
u_max = [1; 1];
du_min = [-0.2; -0.2];
du_max = [0.2; 0.2];
gpc.set_constraints(u_min, u_max, du_min, du_max);
% 仿真参数
T = 100; % 仿真时间
t = 0:T-1;
% 参考轨迹
y_ref = zeros(ny, T);
y_ref(1, 20:end) = 1; % 输出1在t=20时阶跃变化
y_ref(2, 50:end) = 0.5; % 输出2在t=50时阶跃变化
% 初始化变量
y = zeros(ny, T);
u = zeros(nu, T);
% 初始条件
y(:, 1:max([na, nb, nc])) = 0;
u(:, 1:max([na, nb, nc])) = 0;
% 主仿真循环
for k = max([na, nb, nc])+1:T-1
% 获取当前输出
y_current = y(:, k);
% 获取参考值
y_ref_current = y_ref(:, k);
% 计算控制输入
u_current = gpc.control(y_current, y_ref_current);
% 存储控制输入
u(:, k) = u_current;
% 系统仿真 (使用真实系统模型)
y_next = zeros(ny, 1);
% A多项式贡献
for i = 1:na
y_next = y_next - squeeze(A_coeffs(:, :, i+1)) * y(:, k-i+1);
end
% B多项式贡献
for i = 1:nb
y_next = y_next + squeeze(B_coeffs(:, :, i+1)) * u(:, k-i+1);
end
% 添加噪声
noise = 0.01 * randn(ny, 1);
y_next = y_next + noise;
% 存储输出
y(:, k+1) = y_next;
end
% 绘制结果
figure('Position', [100, 100, 1200, 800]);
% 输出响应
subplot(2, 1, 1);
plot(t, y(1, :), 'b-', 'LineWidth', 1.5); hold on;
plot(t, y(2, :), 'r-', 'LineWidth', 1.5);
plot(t, y_ref(1, :), 'b--', 'LineWidth', 1.5);
plot(t, y_ref(2, :), 'r--', 'LineWidth', 1.5);
xlabel('时间');
ylabel('输出');
legend('输出1', '输出2', '参考1', '参考2');
title('多变量GPC控制 - 输出响应');
grid on;
% 控制输入
subplot(2, 1, 2);
plot(t, u(1, :), 'b-', 'LineWidth', 1.5); hold on;
plot(t, u(2, :), 'r-', 'LineWidth', 1.5);
xlabel('时间');
ylabel('控制输入');
legend('输入1', '输入2');
title('控制输入');
grid on;
% 计算性能指标
mse1 = mean((y(1, :) - y_ref(1, :)).^2);
mse2 = mean((y(2, :) - y_ref(2, :)).^2);
control_effort = mean(sum(abs(u), 1));
fprintf('性能指标:\n');
fprintf('输出1 MSE: %.4f\n', mse1);
fprintf('输出2 MSE: %.4f\n', mse2);
fprintf('控制能量: %.4f\n', control_effort);
end
系统辨识辅助函数
function [A_coeffs, B_coeffs] = identify_mimo_system(u, y, na, nb)
% IDENTIFY_MIMO_SYSTEM MIMO系统辨识函数
% 使用最小二乘法估计多项式矩阵系数
[ny, T] = size(y);
nu = size(u, 1);
% 构建回归矩阵
max_lag = max(na, nb);
Phi = [];
for k = max_lag+1:T
% 输出项
phi_row = [];
for i = 1:na
phi_row = [phi_row, -y(:, k-i)'];
end
% 输入项
for i = 1:nb
phi_row = [phi_row, u(:, k-i)'];
end
Phi = [Phi; phi_row];
end
% 构建输出矩阵
Y = y(:, max_lag+1:T)';
% 最小二乘估计
Theta = Phi \ Y;
% 提取系数
A_coeffs = zeros(ny, ny, na+1);
A_coeffs(:,:,1) = eye(ny); % A0
for i = 1:na
start_idx = (i-1)*ny + 1;
end_idx = i*ny;
A_coeffs(:,:,i+1) = Theta(start_idx:end_idx, :)';
end
B_coeffs = zeros(ny, nu, nb+1);
B_coeffs(:,:,1) = zeros(ny, nu); % B0
for i = 1:nb
start_idx = na*ny + (i-1)*nu + 1;
end_idx = na*ny + i*nu;
B_coeffs(:,:,i+1) = Theta(start_idx:end_idx, :)';
end
end
参考代码 多变量隐式广义预测控制 www.youwenfan.com/contentcnf/23922.html
使用说明
-
初始化控制器：
gpc = MultivariableGPC(na, nb, nc, nu, ny, N1, N2, Nu, lambda, alpha); -
设置系统模型：
gpc.set_model(A_coeffs, B_coeffs, C_coeffs); -
设置约束（可选）：
gpc.set_constraints(u_min, u_max, du_min, du_max); -
控制循环：
u = gpc.control(y_measured, y_reference);
关键特性
- 多变量支持：完整支持MIMO系统
- 约束处理：支持输入幅值和变化率约束
- 隐式实现：无需Diophantine方程显式解
- 数值鲁棒性：使用QP求解器处理约束优化
- 实时性能：优化计算效率
参数整定建议
- 预测时域：N2应覆盖系统主要动态，通常选择系统上升时间的1.5-2倍
- 控制时域：Nu通常选择3-5，过大会增加计算负担
- 权重矩阵：
- λ控制输入权重，较大值得到更保守的控制
- α控制输出跟踪性能，较大值加强跟踪
浙公网安备 33010602011771号