【MPC】使用输入增量实现了不同的状态空间MPC公式研究（Matlab代码实现）

💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁

💥1 概述

使用输入增量实现不同状态空间MPC公式的研究

摘要：本文聚焦于使用输入增量的状态空间模型预测控制（MPC）在控制系统中的应用。首先阐述了传统状态空间MPC存在的问题，进而引出使用输入增量的优势，包括避免逆矩阵计算和提高通用性。详细推导了基于输入增量的状态空间MPC公式，并通过函数和Simulink块展示了其在控制系统中的具体实现方式。通过仿真实验验证了该方法的有效性，结果表明使用输入增量的状态空间MPC能够实现更加高效和通用的控制系统，为工程项目提供了更多的可能性和灵活性。

关键词：模型预测控制；状态空间；输入增量；控制系统

一、引言

模型预测控制（MPC）作为一种先进的控制策略，在工业控制领域得到了广泛的应用。状态空间模型是描述系统动态特性的重要工具，基于状态空间模型的MPC能够直接利用系统的状态信息进行控制决策，具有诸多优点。然而，传统的状态空间MPC在实现过程中往往需要进行逆矩阵的计算，这不仅增加了计算复杂度，而且在某些情况下可能导致数值不稳定问题。此外，传统方法的通用性也受到一定限制，难以适应不同类型系统的控制需求。

输入增量在MPC中的应用为解决上述问题提供了新的思路。通过引入输入增量，可以避免直接对矩阵求逆，从而降低计算复杂度并提高数值稳定性。同时，使用输入增量的状态空间MPC公式具有更强的通用性，能够更灵活地应用于各种不同的控制系统。本文将深入研究使用输入增量实现不同状态空间MPC公式的方法，并通过函数和Simulink块展示其具体实现过程，为实际工程应用提供参考。

二、传统状态空间MPC回顾

2.1 状态空间模型

考虑一个线性时不变系统，其状态空间模型可以表示为：

2.2 传统状态空间MPC公式

通过递推状态空间方程，可以得到预测状态与当前状态和预测输入之间的关系：

优化问题的目标函数通常定义为预测输出与参考轨迹之间的误差平方和以及输入变化量的平方和，即：

将预测输出 y(k+i∣k) 用状态和输入表示，并将预测状态关系代入目标函数，可以将优化问题转化为一个二次规划问题，通过求解该问题得到最优控制输入序列，并取第一个元素作为当前时刻的控制输入。

然而，传统方法在求解优化问题时，需要对矩阵进行求逆等运算，这在计算上较为复杂，并且可能存在数值不稳定的问题。

三、使用输入增量的状态空间MPC公式推导

3.1 输入增量定义

3.2 预测状态与输入增量的关系

同样考虑预测时域 Np 和控制时域 Nc，定义预测状态序列 X 和预测输入增量序列 ΔU=

3.3 优化问题构建

目标函数仍然定义为预测输出与参考轨迹之间的误差平方和以及输入增量变化量的平方和（这里可以进一步考虑输入增量变化量的约束以增强控制的平滑性，为简化推导，暂不考虑）。预测输出 y(k+i∣k) 可以表示为：

这是一个关于输入增量序列 ΔU 的二次规划问题，与传统方法相比，避免了直接对矩阵求逆，降低了计算复杂度。

3.4 优化问题求解

将目标函数展开并整理成标准的二次规划形式：

四、使用输入增量的状态空间MPC在m - function中的实现

4.1 m - function框架设计

设计一个m - function来实现使用输入增量的状态空间MPC算法。该函数的主要输入参数包括系统矩阵 A、B、C、D，加权矩阵 Q 和 R，预测时域 Np，控制时域 Nc，当前状态 x(k)，上一时刻输入 u(k−1) 以及参考轨迹 r。输出为当前时刻的控制输入 u(k)。

4.2 具体实现步骤

初始化参数：根据输入参数初始化系统矩阵、加权矩阵、预测时域和控制时域等。
构建预测矩阵：根据推导的公式计算 F~、Φ~ 和 G~。
构建目标函数矩阵和向量：根据预测状态与输入增量的关系以及目标函数定义，计算二次规划问题中的矩阵 H 和向量 f。
求解二次规划问题：使用Matlab中的二次规划求解器（如quadprog函数）求解关于输入增量序列 ΔU 的优化问题。
计算当前控制输入：从最优输入增量序列中取出第一个元素，结合上一时刻输入计算当前时刻的控制输入 u(k)。

4.3 示例代码

matlab

1function u = input_increment_MPC(A, B, C, D, Q, R, Np, Nc, x, u_prev, r) 2 % 构建预测矩阵 3 [F_tilde, Phi_tilde, G_tilde] = build_prediction_matrices(A, B, Np, Nc); 4 5 % 构建目标函数矩阵和向量 6 [H, f] = build_cost_function(A, B, C, D, Q, R, Np, Nc, x, u_prev, r, F_tilde, Phi_tilde, G_tilde); 7 8 % 求解二次规划问题 9 options = optimoptions('quadprog', 'Display', 'off'); 10 DeltaU_opt = quadprog(H, f, [], [], [], [], [], [], [], options); 11 12 % 计算当前控制输入 13 u = u_prev + DeltaU_opt(1); 14end 15 16function [F_tilde, Phi_tilde, G_tilde] = build_prediction_matrices(A, B, Np, Nc) 17 % 此处省略具体构建过程，根据推导公式计算F_tilde, Phi_tilde, G_tilde 18 % 示例代码框架 19 n = size(A, 1); 20 m = size(B, 2); 21 22 F_tilde = zeros(n * Np, n); 23 Phi_tilde = zeros(n * Np, m * Nc); 24 G_tilde = zeros(n * Np, m); 25 26 % 具体计算过程... 27end 28 29function [H, f] = build_cost_function(A, B, C, D, Q, R, Np, Nc, x, u_prev, r, F_tilde, Phi_tilde, G_tilde) 30 % 此处省略具体构建过程，根据推导公式计算H和f 31 % 示例代码框架 32 n = size(A, 1); 33 p = size(C, 1); 34 m = size(B, 2); 35 36 % 初始化H和f 37 H = zeros(m * Nc, m * Nc); 38 f = zeros(m * Nc, 1); 39 40 % 具体计算过程... 41end

五、使用输入增量的状态空间MPC在Simulink中的实现

5.1 Simulink模型搭建思路

在Simulink中实现使用输入增量的状态空间MPC，可以通过自定义S - function块或者使用基本模块搭建算法逻辑。这里采用基本模块搭建的方式，主要思路是将MPC算法的各个步骤分解为不同的模块，包括状态预测模块、目标函数计算模块、优化求解模块和控制输入生成模块。

5.2 具体模块实现

状态预测模块：根据输入增量的状态空间方程，使用矩阵运算模块（如Matrix Multiply模块）实现预测状态的计算。输入为当前状态 x(k)、上一时刻输入 u(k−1) 和预测输入增量序列 ΔU，输出为预测状态序列 X。
目标函数计算模块：根据预测状态和参考轨迹，使用矩阵运算和加法模块计算目标函数中的误差项和输入增量项，然后进行加权求和得到目标函数值。
优化求解模块：由于Simulink中没有直接的二次规划求解模块，可以通过编写M文件函数，并使用MATLAB Function模块调用该函数来实现二次规划问题的求解。输入为目标函数矩阵 H 和向量 f，输出为最优输入增量序列 ΔU∗。
控制输入生成模块：从最优输入增量序列中取出第一个元素，与上一时刻输入相加得到当前时刻的控制输入 u(k)。

5.3 仿真示例

以一个简单的二阶系统为例进行仿真，系统矩阵 A=[01;−1−2]，B=[0;1]，C=[10]，D=0。设置预测时域 Np=10，控制时域 Nc=5，加权矩阵 Q=diag([10.1])，R=0.01。参考轨迹为阶跃信号。通过Simulink仿真可以得到系统的响应曲线，结果表明使用输入增量的状态空间MPC能够有效地跟踪参考轨迹，实现良好的控制效果。

六、结论

本文深入研究了使用输入增量实现不同状态空间MPC公式的方法。通过推导基于输入增量的状态空间MPC公式，避免了传统方法中逆矩阵的计算，降低了计算复杂度并提高了数值稳定性。同时，该方法具有更强的通用性，能够适应不同类型系统的控制需求。通过在m - function和Simulink中的实现，展示了该方法的实际应用过程。仿真实验结果表明，使用输入增量的状态空间MPC能够实现更加高效和通用的控制系统，为工程项目提供了更多的可能性和灵活性。未来的研究可以进一步考虑输入增量变化量的约束以及非线性系统的应用等问题。

📚2 运行结果

部分代码：

%% Simulation
% 150 seconds (1500 sampling intervals) simulation is conducted with
% several setpoint changes and random cooling water temperature changes
% within positive and negative 1 degree.
% Simulation length and variables for results
N=1500;
x0=zeros(6,1);
Y=zeros(N,2);
U=zeros(N,2);
% Predefined reference
T=zeros(N,2);
T(10:N,1)=1;
T(351:N,1)=3;
T(600:N,1)=5;
T(1100:N,1)=3;
T(100:N,2)=2;
T(451:N,2)=1;
T(700:N,2)=4;
T(1200:N,2)=2;
% Simulation
%%
for k=1:N
% Process disturbances
w=Bd*(rand(2,1)-0.5)*2;
% Measurements noise
v=0.01*randn(2,1);
% actual measurement
y=C*x0+v;
% online controller
u=ssmpc(y,T(k:end,:)');
% plant update
x0=A*x0+Bu*u+w;
% save results
Y(k,:)=y';
U(k,:)=u';
end

%% Results
% The simulation results are summarized in two sub-plots.
t=(0:N-1)*0.1;
subplot(211)
plot(t,Y,t,T,'--','linewidth',2)
title('output and setpoint')
ylabel('temp, C^\circ')
legend('T_1','T_2','T_1 Ref','T_2 Ref','location','south','Orientation','horizontal')
subplot(212)
stairs(t,U,'linewidth',2)
legend('u_1','u_2','location','southeast')
title('input')
ylabel('flow rate, m^3/s')
xlabel('time, s')