【数据驱动预测控制2】Willems基本引理

系列文章目录

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

---

---

前言

 本篇文章我们来介绍一个数据驱动算法中的一个重要引理——Willems基本引理。对于一个线性时不变（LTI）系统，其输入输出关系可以通过系统的状态空间模型来描述。而通过Willems基本引理，如果我们有足够的输入输出数据，我们可以利用这些数据来构建系统模型。
 首先我们考虑一个LTI系统：$$\begin{gathered} x(k+1)=Ax(k)+Bu(k) \\ y(k)=Cx(k)+Du(k)\text{\hspace{1em}(1)} \end{gathered}$$
其中$x(k)\in R^{n_x}$为系统的状态变量，$y(k)\in R^{n_y}、u(k)\in R^{n_u}$分别为系统的输入输出，$A$、$B$、$C$和$D$为适当维度系统矩阵。我们平常要对系统进行仿真或者控制都需要$A$、$B$、$C$和$D$是已知的。而本文介绍的Willems基本引理可以在系统矩阵未知的情况下只通过系统的输入输出轨迹来达到对系统的仿真或控制。

---

一、预备知识

在介绍Willems基本引理之前，我们需要先介绍几个在Willems基本引理中要用到的数学表达和定义。

1.时间序列

 我们记一个时间序列 { u k } k = 0 N − 1 = { u 0 , u 1 , ⋯   , u N − 1 } \{u_k\}_{k=0}^{N-1} = \{u_0,u_1,\cdots,u_{N-1}\} {uk​}k=0N−1​={u0​,u1​,⋯,uN−1​}，定义$$u_{[a,b]}=\left[\begin{array}{c}{u}_a\\ u_{a+1}\\ ⋮\\ u_b\end{array}\right]$$其中$a<b$。

2.汉克矩阵

 首先给定一个时间序列 { u k } k = 0 N − 1 \{u_k\}_{k=0}^{N-1} {uk​}k=0N−1​，我们记由$u_k$构造的汉克矩阵为$$H_L(u)=\left[\begin{array}{cccc}{u}_0& u_1& \cdots & u_{N-L}\\ u_1& u_2& \cdots & u_{N-L+1}\\ ⋮& ⋮& \ddots & ⋮\\ u_{L-1}& u_L& \cdots & u_{N-1}\end{array}\right]$$
汉克矩阵会在Willems基本引理中发挥一个非常重要的作用。

3.系统输入输出轨迹

 系统输入轨迹是指系统施加给定输入信号的时间序列，输入可以是阶跃信号、正弦波、白噪声或其他任何类型的测量信号。系统输出轨迹是指系统在输入作用下，所产生的动态响应，通常也是一个时间序列。为了更加直观的说明，我们用MATLAB代码来演示系统的输入输出轨迹。

%% 初始化
%系统矩阵A,B,C,D
A = [0 0.0734;-6.5529 -0.4997];
B = [-0.07221;-9.6277];
C = [0.1 0.3];
D = 0;
ny = length(C(:,1));%输出向量y的维度
nu = length(B(1,:));%输入向量u的维度
nx = length(A);%状态向量x的维度
N = 200;%轨迹长度
u_tra = zeros(nu,N);%系统输入轨迹
y_tra = nan(ny,N);%对应u_tra的系统输出轨迹
%% 采集系统阶跃响应的轨迹
xk = randn(nx,1);%系统初始状态
for k = 1:N
    uk = 1;%阶跃信号
    yk = C*xk;
    xk = A*xk + B*uk;
    
    u_tra(:,k) = uk;
    y_tra(:,k) = yk;
end
%% 绘图
subplot(2,1,1); plot(1:N, y_tra(1,:)); ylabel('y'); title('trajectory')
subplot(2,1,2); plot(1:N, u_tra(1,:));
ylabel('control input'); xlabel('k')

上述代码所采集到的u_tra和y_tra就是系统的一条长度为$N$的输入输出轨迹。

4.持续激励性

 有了汉克矩阵和输入输出轨迹的概念，我们再来介绍一下持续激励性的定义：
 定义（持续激励性）： 若存在一段输入轨迹 { u k } k = 0 N − 1 \{u_k\}_{k=0}^{N-1} {uk​}k=0N−1​满足：
$rank(H_L(u))=n_{u}L$，则称 { u k } k = 0 N − 1 \{u_k\}_{k=0}^{N-1} {uk​}k=0N−1​是 $L$阶持续激励 的。

二、Willems基本引理

有了上述的几个概念，我们可以正式来看一下Willems基本引理的内容了。
 Willems基本引理： 若 { u k d } k = 0 N − 1 \{u^d_k\}_{k=0}^{N-1} {ukd​}k=0N−1​是满足$L+n_x$阶持续激励性的一段系统输入轨迹， { u k d ， y k d } k = 0 N − 1 \{u^d_k，y^d_k\}_{k=0}^{N-1} {ukd​，ykd​}k=0N−1​为系统(1)的一段输入输出轨迹。则结合Hankel矩阵的定义，存在$g\in R^{N-L+1}$满足如下：
$$\left[\begin{array}{c}{H}_L(u^d)\\ H_L(y^d)\end{array}\right]g=\left[\begin{array}{c}\bar{u}\\ \bar{y}\end{array}\right]\text{\hspace{1em}(2)}$$则 { u ˉ k ， y ˉ k } k = 0 L − 1 \{\bar{u}_k，\bar{y}_k\}_{k=0}^{L-1} {uˉk​，yˉ​k​}k=0L−1​也是系统(1)的一条输入输出轨迹。
  { u k d } k = 0 N − 1 \{u^d_k\}_{k=0}^{N-1} {ukd​}k=0N−1​要满足$L+n_x$阶持续激励性的原因是为了完整的表达系统的所有动态过程。
接下来演示如何利用公式(2)把系统(1)的状态空间表达式替换掉，我们以系统仿真的例子说明。

三、基于Willems基本引理的系统仿真

 我们记 { u k d ， y k d } k = 0 N − 1 \{u^d_k，y^d_k\}_{k=0}^{N-1} {ukd​，ykd​}k=0N−1​是系统(1)一条长度为$N$的先验测量输入输出轨迹。另外令预测时域为$L$(注意这里的$L$和上文的$L$概念不同)，过去时域为$P$，$u_k(t)$记为在$t$时刻对$t+k$时刻的预测值，$k\in [-P,L-1]$。由于这里我们是做系统仿真，因此令预测时域$L=1$，$L>1$的情况将用于预测控制的算法中。接下来我们可以将公式(2)改写成：
$$\left[\begin{array}{c}{H}_{P+1}(u^d)\\ H_{P+1}(y^d)\end{array}\right]g=\left[\begin{array}{c}{H}_P(u_{[0,N-2]}^d)\\ H_1(u_{[P-1,N-1]}^d)\\ H_P(y_{[0,N-2]}^d)\\ H_1(y_{[P-1,N-1]}^d)\end{array}\right]g=\left[\begin{array}{c}{\bar{u}}_{[-P,-1]}(t)\\ {\bar{u}}_0(t)\\ {\bar{y}}_{[-P,-1]}(t)\\ {\bar{y}}_0(t)\end{array}\right]\text{\hspace{1em}(3)}$$注意上述的轨迹 { u k d } k = 0 N − 1 \{u^d_k\}_{k=0}^{N-1} {ukd​}k=0N−1​需要满足$P+1+n_x$阶持续激励。
接下来的问题是，在当前时刻$t$，我们对系统施加${\bar{u}}_0(t)$的控制作用下，系统产生的输出${\bar{y}}_0(t)$为多少？
回归到公式(3)，在公式(3)中我们只有$g$和${\bar{y}}_0(t)$是未知的。首先求$g$的表达式，可以得到
$$\left[\begin{array}{c}{H}_P(u_{[0,N-2]}^d)\\ H_1(u_{[P-1,N-1]}^d)\\ H_P(y_{[0,N-2]}^d)\end{array}\right]g=\left[\begin{array}{c}{\bar{u}}_{[-P,-1]}(t)\\ {\bar{u}}_0(t)\\ {\bar{y}}_{[-P,-1]}(t)\end{array}\right]\text{\hspace{1em}(4)}$$
通过 { u k d } k = 0 N − 1 \{u^d_k\}_{k=0}^{N-1} {ukd​}k=0N−1​是$P+1+n_x$阶持续激励，可以证明矩阵$H:=\left[\begin{array}{c}{H}_P(u_{[0,N-2]}^d)\\ H_1(u_{[P-1,N-1]}^d)\\ H_P(y_{[0,N-2]}^d)\end{array}\right]$具有行满秩。因此我们记$H^{+}$为矩阵$H$的左逆矩阵，即$H^{+}H=I$，然后可以得到
$$g=H^{+}\left[\begin{array}{c}{\bar{u}}_{[-P,-1]}(t)\\ {\bar{u}}_0(t)\\ {\bar{y}}_{[-P,-1]}(t)\end{array}\right]\text{\hspace{1em}(5)}$$接着把公式(5)代入到公式(3)中，我们可以得到${\bar{y}}_0(t)$的表达式
$${\bar{y}}_0(t)=H_1(y_{[P-1,N-1]}^d)g=H_1(y_{[P-1,N-1]}^d)H^{+}\left[\begin{array}{c}{\bar{u}}_{[-P,-1]}(t)\\ {\bar{u}}_0(t)\\ {\bar{y}}_{[-P,-1]}(t)\end{array}\right]\text{\hspace{1em}(6)}$$
因此我们可以利用公式(6)来代替系统(1)的状态空间表达式。
备注（过去时域P）： 过去时域$P$的长度的是有条件的，过去时域要足够的长以至于能够完整表达系统的当前状态，一般情况下过去时域都选择为状态向量$x$的维度，即$P=n_x$。$P$太小会导致预测时域$L$预测出来的值不唯一从而导致预测值不准确；$P$太大会导致计算量增加，但不会影响预测准确度。
备注（轨迹长度N）： 为了使 { u k d } k = 0 N − 1 \{u^d_k\}_{k=0}^{N-1} {ukd​}k=0N−1​满足$P+L+n_x$阶持续激励，轨迹长度$N$要满足$N\ge (n_u+1)(P+L)-1$。
备注（$u_d$的持续激励性）： 为了得到一个满足持续激励性的输入轨迹 { u k d } k = 0 N − 1 \{u^d_k\}_{k=0}^{N-1} {ukd​}k=0N−1​，我们一般都是施加随机控制信号，下面的仿真代码也将演示如何得到一个满足持续激励性条件的输入轨迹 { u k d } k = 0 N − 1 \{u^d_k\}_{k=0}^{N-1} {ukd​}k=0N−1​。

接下来我们用matlab代码来演示公式(6)和状态空间表达式的等价性。

%% 初始化
%系统矩阵A,B,C,D
rng(123);
A = [0 0.0734;-6.5529 -0.4997];
B = [-0.07221;-9.6277];
C = [0.1 0.3];
D = 0;
ny = length(C(:,1));%输出向量y的维度
nu = length(B(1,:));%输入向量u的维度
nx = length(A);%状态向量x的维度
N = 20;%预先采集轨迹长度
L = 1;
P = nx;
ud = zeros(nu,N);%系统输入轨迹
yd = nan(ny,N);%对应ud的系统输出轨迹
%% 采集系统轨迹
xk = randn(nx,1);%系统初始状态
for k = 1:N
    uk = rand(1,nu);%产生持续激励信号
    yk = C*xk;
    xk = A*xk + B*uk;
    
    ud(:,k) = uk;
    yd(:,k) = yk;
end

Hu = zeros(nu*(L+P),N-L-P+1);
Hy = zeros(ny*(L+P),N-L-P+1);
for i = 0:N-L-P
   Hu(:,i+1) = ud(i*nu+1:i*nu+(L+P)*nu);
   Hy(:,i+1) = yd(i*ny+1:i*ny+(L+P)*p);
end
assert(rank(Hu) == (L+P)*nu,"rank(Hu) is not full rank")%检验ud是否满足持续激励条件

%% 数据驱动方法仿真
k_steps = 200;
udd = zeros(nu,k_steps);%系统输入
ydd = nan(ny,k_steps);%系统输出
up = zeros(P*nu,1);%过去输入
yp = zeros(P*ny,1);%过去输出
H_inv = ([Hu;Hy(1:2,:)]' * [Hu;Hy(1:2,:)]) \ [Hu;Hy(1:2,:)]';%求左逆矩阵
for k = 1:k_steps
    uk = sin(k*0.2);%正弦输入信号
    yk = Hy(P+1,:)*H_inv*[up;uk;yp];%求取系统输出
    
    yp = circshift(yp,-ny,1);
    yp(end-ny+1:end,:) = yk;%更新过去输出
    up = circshift(up,-nu,1);
    up(end-nu+1:end,:) = uk;%更新过去输入
    
    udd(:,k) = uk;
    ydd(:,k) = yk;
end

%% 状态空间方法仿真
uss = zeros(nu,k_steps);%系统输入
yss = nan(ny,k_steps);%系统输出
xk = zeros(nx,1);
for k = 1:k_steps
    uk = sin(k*0.2);%正弦输入信号
    yk = C*xk;
    xk = A*xk + B*uk;
    
    uss(:,k) = uk;
    yss(:,k) = yk;
end
subplot(2,1,1); plot(1:k_steps, yss(1,:),'b'); ylabel('output');
hold on; plot(1:k_steps, ydd(1,:),'r--'); 
hold off;legend('数据驱动','状态空间');
subplot(2,1,2); plot(1:k_steps, uss(1,:),'b');ylabel('input')
hold on;subplot(2,1,2); plot(1:k_steps, udd(1,:),'r--');
hold off;legend('数据驱动','状态空间');

---

运行结果：

从上面的仿真结果我们可以看出利用Willems基本引理和传统的状态空间仿真结果是一样的，而且在Willems基本引理仿真过程中没有用到系统矩阵，只利用了系统的输入输出数据。
下一篇文章我们将把Willems基本引理应用到预测控制中，即一种数据驱动预测控制算法。