【数据驱动预测控制3】数据驱动预测控制算法及仿真

麦斯威尔逊

2344人浏览 · 2024-12-26 11:42:39

麦斯威尔逊 · 2024-12-26 11:42:39 发布

系列文章目录

【数据驱动预测控制1】数据驱动预测控制的概念
 【数据驱动预测控制2】Willems基本引理
【数据驱动预测控制3】数据驱动预测控制算法及仿真

文章目录

系列文章目录
前言
一、模型预测控制
二、数据驱动预测控制
三、数据驱动预测控制仿真

前言

在上一篇文章中，我们深入探讨了Willems基本引理，这一理论为理解和优化动态系统提供了强有力的工具。现在我们将这一理论应用到模型预测控制（MPC）中，探索其在实际控制系统中的潜力。但在深入讨论之前，让我们先回顾一下传统的模型预测控制。

提示：以下是本篇文章正文内容，下面案例可供参考

一、模型预测控制

考虑一个LTI系统： $x_{t+1} = Ax_t+Bu_t\\ y_t = C x_t+D u_t \tag{1}$
其中 $x_t\in R^{n_x}$ 为系统的状态变量， $y_t\in R^{n_y} 、u_t\in R^{n_u}$ 分别为系统的输入输出， $A$ 、 $B$ 、 $C$ 和 $D$ 为适当维度系统矩阵。我们都知道，MPC拥有三大特点：预测模型、滚动优化、反馈校正。具体来说，MPC通过以下步骤实现控制目标：
预测： 基于当前的系统状态和已知的模型，MPC预测未来一段时间内系统的行为。这通常涉及求解一个优化问题，以找到使系统在未来某个时间窗口内达到期望状态的控制输入序列。
优化： 在预测的基础上，MPC选择一个最优的控制输入序列，通常是在满足一定约束条件下使某个性能指标（如成本函数）最小化。这个性能指标可能包括系统的稳定性、跟踪性能、能耗等。
反馈校正： 在实际执行控制输入时，MPC不断监测系统的实际状态，并与预测状态进行比较。如果存在差异，MPC会重新进行预测和优化，以修正未来的控制输入序列。

当系统矩阵 $A$ 、 $B$ 、 $C$ 和 $D$ 是已知的时候，我们可以将上述步骤整合为如下预测时域长度为 $L$ 的MPC优化问题：
$\begin{aligned} \min_{u,x,y} \quad \sum_{k=0}^{L-1}\Vert \bar{y}_k(t)-r_{t+k}\Vert^2_Q+\Vert \bar{u}_k(t) \Vert ^2_R \end{aligned}\\ \begin{align} s.t\quad x_{k+1} &= Ax_k+Bu_k ,\notag \\ y_k&=Cx_k +Du_k,k\in[0,L-1],\tag{2a}\\ x_0 &= \hat{x}(t), \tag{2b} \end{align}\\ \begin{align} \bar{u}_k(t)\in\mathbb{U},\bar{y}_k(t)\in\mathbb{Y},k\in[0,L-1].\tag{2b} \end{align}$ 其中公式(2a)是预测模型， $r_{t+k}\in \mathbb{R}^{n_y}$ 是 $t + k$ 时刻的期望输出， $x_k$ 和 $y_k$ 表示 $t + k$ 时刻的预测状态和输出， $\hat{x}(t)$ 表示 $t$ 时刻的估计状态。传统MPC算法步骤如下：
算法 MPC
准备: 预测时域L，系统矩阵 $A, B, C, D$
运行:
step1：估计当前系统状态 $x_t$ ，然后求解优化问题(4)。
step2：对系统施加控制输入 $u_t=\bar{u}^*_0(t)$ 。
step3： $t = t + 1$ ，回到step1。

从上述优化问题中可以看出，为了构造预测模型，必须要知道系统矩阵 $A$ 、 $B$ 、 $C$ 和 $D$ 具体的值。而对于复杂的系统，系统矩阵是很难获取的。面对这一困境，学术界和工业界开始探索新的途径，旨在不依赖系统矩阵的具体数值，而是利用系统的输入输出数据来直接构建预测模型。这一思想的核心在于，通过收集和分析系统在实际运行过程中的输入输出数据，我们可以挖掘出系统行为的潜在规律，并据此构建出能够准确预测系统未来状态的模型。

正是在这样的背景下，数据驱动预测控制应运而生。这种方法摒弃了传统预测控制对系统矩阵的依赖，转而利用丰富的数据资源来优化控制策略。通过先进的数据分析、机器学习和优化算法，数据驱动预测控制能够在不掌握系统精确数学模型的情况下，实现对系统的有效控制。而我们介绍的就是一种以Willems基本引理为基础的数据驱动预测控制算法。

二、数据驱动预测控制

再次回顾一下Willems基本引理：

Willems基本引理： 若 $\{u^d_k\}_{k=0}^{N-1}$ 是满足 $L + n_x$ 阶持续激励性的一段系统输入轨迹， $\{u^d_k，y^d_k\}_{k=0}^{N-1}$ 为系统(1)的一段输入输出轨迹。则结合Hankel矩阵的定义，存在 $g ∈R^{N −L+1}$ 满足如下：
$\begin{bmatrix} H_L(u^d) \\ H_L(y^d) \end{bmatrix}g = \begin{bmatrix} \bar{u}\\ \bar{y}\end{bmatrix} \tag{3}$ 则 $\{\bar{u}_k，\bar{y}_k\}_{k=0}^{L-1}$ 也是系统(1)的一条输入输出轨迹。

根据Willems基本引理我们可以将公式(3)右侧的输入输出轨迹 $\{\bar{u}，\bar{y}\}$ 分为两段：过去时域P和预测时域L。此处过去时域的作用就是传统MPC优化问题中公式(2b)的作用，用来决定当前时刻系统的状态，它们为预测未来的输入输出轨迹提供了必要的出发点和约束条件。过去时域和预测时域的关系，就好比我们的生活一样，过去几天我们过得怎么样，无论是经历了愉快还是难过，都将深刻地影响我们对未来几天的规划。
想象一下，如果你在过去几天里工作顺利、心情愉悦，那么你可能会更加积极地规划未来的活动。这些规划都是基于你对过去经历的正面评价和对未来的乐观预期。
相反，如果你在过去几天里遇到了困难或挫折，那么你可能会更加谨慎地规划未来，专注于解决问题、调整心态。这些规划同样反映了你对过去经历的反思和对未来的审慎预期。
过去时域就像是我们生活的“回忆录”，记录了我们过去的经历、成长和变化。而预测时域则像是我们生活的“蓝图”，描绘了我们对未来的期待、规划和梦想。两者共同构成了我们生活的完整画卷，让我们在回顾过去的同时，也能展望未来，不断前行。因此，就像我们在生活中需要平衡过去和未来的关系一样，在模型预测控制中，我们也需要充分利用过去时域和预测时域的信息，以实现最优的控制和决策。

有了上述的理论支持，接下来可以利用Willems基本引理来替换优化问题(2)中的预测模型，可以得到数据驱动预测控制的优化问题如下:
$\begin{aligned} \min_{u,y} \quad \sum_{k=0}^{L-1}\Vert \bar{y}_k(t)-r_{t+k}\Vert^2_Q+\Vert \bar{u}_k(t)\Vert ^2_R, \end{aligned}\\ \begin{align} s.t\quad \begin{bmatrix} H_{L+P}(u_N^d)\\ H_{L+P}(y_N^d) \end{bmatrix}g(t)= \begin{bmatrix} \bar{u}(t)\\ \bar{y}(t) \end{bmatrix},\tag{4a} \end{align}\\ \begin{align} \begin{bmatrix} \bar{u}_{[-P,-1]}(t)\\ \bar{y}_{[-P,-1]}(t) \end{bmatrix} = \begin{bmatrix} u_{[t-P,t-1]}\\ y_{[t-P,t-1]} \end{bmatrix},\tag{4b} \end{align}\\ \begin{align} \bar{u}_k(t)\in\mathbb{U},\bar{y}_k(t)\in\mathbb{Y},k\in[0,L-1].\tag{4c} \end{align}$
其中 $L$ 代表预测时域步长， $P$ 代表过去时域步长。与优化问题(2)对比，(4a)等价于(2a)，(4b)等价于(2b)，现在的优化问题(4)则完全不需要系统矩阵，而是通过输入输出数据来预测未来输入输出。数据驱动预测控制算法如下所示，其中 $\bar{u}^*_k(t)$ 记为在 $t$ 时刻对 $t + k$ 时刻的最优预测值
算法数据驱动预测控制(DDPC)
准备: 预测时域L，过去时域P，历史输入输出数据 $\{u^d_k，y^d_k\}_{k=0}^{N-1}$ ，其中 $\{u^d_k \}_{k=0}^{N-1}$ 满足 $L+P+n_x$ 阶持续激励。
运行:
step1：采集过去输入输出 ${u_k，y_k\}_{k=-P}^{-1}$ 然后求解优化问题(4)。
step2：对系统施加控制输入 $u_t=\bar{u}^*_0(t)$ 。
step3： $t = t + 1$ ，回到step1。

三、数据驱动预测控制仿真

接下来用一个仿真例子对数据驱动预测控制算法更加仔细的了解。被控对象为如下连续LTI系统
$\dot{x}(t) = \begin{bmatrix}-3.8557&1.4467\\0.4690&-4.7081\end{bmatrix}x(t)+\begin{bmatrix}230.6739\\653.5547\end{bmatrix}u(t)\\ y(t)=\begin{bmatrix} 1&0\\0&1\end{bmatrix}x(t)$
对该系统的控制目标为：在满足输入约束 $u_t \in \mathbb{U}=[-2,2]$ 的条件下， $\in [0,2)$ 时间内将输出控制到 $y_t = [0,0]^T$ ； $\in [2,4)$ 时间内将输出控制到 $y_t = [43.3154,56.0351]^T$ ； $\in [4,6)$ 时间内将输出控制到 $y_t = [86.6309,112.0702]^T$ ；6s后将输出控制到 $y_t = [21.6577,28.0176]^T$ 。注意在设计控制器的过程中系统矩阵是未知的。

为了设计控制器，首先采样时间 $T = 0.02 s$ 的开环实验中采集一条长度为N = 50的持续激励输入输出轨迹，并选择预测时域 $L = 7$ 和过去时域 $P = 2$ ，具体代码如下(下面的代码中会用到CasADi求解器，请先确保Matlab中安装了CasADi求解器，未安装的伙伴可以到CasADi官网下载并安装)

import casadi.*
clear; close all; clc;
rng(123);
Ts = 0.02; %采样时间
As = [-3.8557,1.4467;
    0.4690,-4.7081];
Bs = [230.6739;653.5547];
C = [1 0;
    0 1 ];
D = 0;
%% 系统离散化
sys_cont = ss(As, Bs, C, D);
sys_disc = c2d(sys_cont, Ts, 'zoh');  % 使用零阶保持器进行转换   
A = sys_disc.A;
B = sys_disc.B;
nx = size(As,1); nu = size(Bs,2); ny = size(C,1);  % 系统维度

%% 初始化
L = 7;%预测时域
P = 2;%过去时域
N = 50;%采集数据长度
%构造汉克矩阵
Hu = zeros(nu*(P+L),N-(P+L)+1);
Hy = zeros(ny*(P+L),N-(P+L)+1);
yd = nan(ny,P+L);
ud = zeros(nu,P+L);
u_min = -2; u_max = 2;%输入约束
for k = 1:N-(P+L)+1
    for i = 1:P+L
        ud(:,i) = u_min + (u_max-u_min).*rand(nu,1);%输入数据
    end
    xd = randn(nx,1);%随机当前系统状态
    for i=1:P+L
        yd(:,i) = C*xd;%采集输出数据
        xd = A*xd + B*ud(:,i);
    end
    Hy(:,k) = yd(:);
    Hu(:,k) = ud(:);
end
H = [Hu(1:(P+L)*nu,:);Hy(1:(P+L)*ny,:)];%构造汉克矩阵
assert(rank(Hu) == (P+L)*nu,"rank(Hu) is not full rank");%检测是否满足持续激励条件

%% 优化器创建
% 创建优化变量 
y_L_ddpc = SX.sym('y_L_ddpc',ny,L);%预测时域
u_L_ddpc = SX.sym('u_L_ddpc',nu,L);
g = SX.sym('g',N-(P+L)+1,1);

y_P_ddpc = SX.sym('y_P_ddpc',ny,P);%过去时域
u_P_ddpc = SX.sym('u_P_ddpc',nu,P);
rk = SX.sym('rk',ny,1);%目标输出
%权重矩阵
R = 1;
Q = [5 0;0 1];
% 创建目标函数  
obj = 0;  
for k = 1:L  
    obj = obj + (y_L_ddpc(:,k)-rk)'*Q*(y_L_ddpc(:,k)-rk) + (u_L_ddpc(:,k))'*R*(u_L_ddpc(:,k));
end
% 创建约束(4b)
b = vertcat( u_P_ddpc(:),u_L_ddpc(:),y_P_ddpc(:), y_L_ddpc(:));
cons_ddpc = H* g - b;
% 创建变量界限(4c)
lbx_ddpc = [repmat(-2,nu*L,1);-inf(L*ny,1);-inf(N-(P+L)+1,1)];%下界
ubx_ddpc = [repmat(2,nu*L,1);inf(L*ny,1);inf(N-(P+L)+1,1)];  %上界
%生成求解器
nlp_dpc = struct('x', [u_L_ddpc(:); y_L_ddpc(:);g(:)], 'f', obj, 'g', cons_ddpc,'p',[u_P_ddpc(:); y_P_ddpc(:);rk(:)]);
S_dpc = nlpsol('S_dpc', 'ipopt', nlp_dpc);

%% 采集测试数据
n_exc = P;%采集过去时域长度的输入输出，确保能够进行优化问题求解
u_test = zeros(nu,n_exc);
y_test = zeros(ny,n_exc);
xk = randn(nx,1);%随机初始状态
for i=1:n_exc
    y_test(:,i) = C*xk ;
    xk = A*xk+B*u_test(:,i);
end
y_P= y_test(:,end-P+1:end);
u_P = u_test(:,end-P+1:end);
%过去时域序列
y_P = y_P(:);
u_P = u_P(:);
%% 闭环预测优化求解
time_sim = 10;%仿真时间10s
k_steps = time_sim/Ts;
y_history = zeros(ny,k_steps);
reference = zeros(ny,k_steps);
u_history = zeros(nu,k_steps);
rk = [0;0];
for k = 1:k_steps
    if(k*Ts == 2)
        rk=[43.3154;56.0351];
    elseif(k*Ts == 4)
        rk = [86.6309;112.0702];
    elseif(k*Ts == 6)
        rk = [21.6577;28.0176];
    end
    %求解二次规划问题
    r = S_dpc('p',[u_P;y_P;rk],'lbx',lbx_ddpc,'ubx',ubx_ddpc,'lbg',0,'ubg',0);
    x_pot=full(r.x);
    u_L_ddpc = x_pot(1:L*nu);
    u_L_ddpc = reshape(u_L_ddpc,nu,L);
    uk = u_L_ddpc(:,1);
    %将uk施加到控制系统
    yk = C*xk ;
    xk = A*xk + B*uk ;
    y_history(:,k) = yk;
    u_history(:,k) = uk;
    reference(:,k) = rk;
    %更新过去时域序列
    y_P = circshift(y_P,-ny,1);
    y_P(end-ny+1:end,:) = yk;
    u_P = circshift(u_P,-nu,1);
    u_P(end-nu+1:end,:) = uk;

end

%% 绘制图像
figure(1)
subplot(211)
hold on
plot(Ts:Ts:time_sim,y_history(1,:),'b-','LineWidth',1.5);
plot(Ts:Ts:time_sim,reference(1,:),'black--','LineWidth',1.5);
hold off
grid on
ylabel("Output  y_1",'FontSize',12);
xlabel('time(s)','FontSize',12,'FontName','Times New Roman');
legend("DDPC","reference");
subplot(212)

hold on
plot(Ts:Ts:time_sim,y_history(2,:),'b-','LineWidth',1.5);
plot(Ts:Ts:time_sim,reference(2,:),'black--','LineWidth',1.5);
hold off
grid on
ylabel("Output  y_2",'FontSize',12);
xlabel('time(s)','FontSize',12,'FontName','Times New Roman');
legend("DDPC","reference");