LQR：Linear Quadratic Regulator 线性二次型调节器

1 前言

如果所研究的系统是线性的，且性能指标为状态变量和控制变量的二次型函数，则最优控制问题称为线性二次型问题。LQR，Linear Quadratic Regulator，即线性二次型调节器， 是求解线性二次型问题常用的求解方法。LQR ，其对象是现代控制理论中以状态空间形式给出的线性系统，而目标函数为对象状态和控制输入的二次型函数。LQR最优设计是指设计出的状态反馈控制器 K要使二次型目标函数J 取最小值，而 K由权矩阵Q 与 R 唯一决定，故此 Q、 R 的选择尤为重要。LQR理论是现代控制理论中发展最早也最为成熟的一种状态空间设计法。特别可贵的是，LQR可得到状态线性反馈的最优控制规律，易于构成闭环最优控制。

2 现代控制理论基础

2.1 状态空间描述

一个线性定常连续系统的状态空间描述包含两个部分：状态方程和输出方程。
1、状态方程（输入维度为 $r$, 状态维度为 $n$, 输出维度为 $m$）。状态方程描述了系统下一状态与当前状态和当前行为的关系，即
$$\dot{x}(t)=Ax(t)+Bu(t)$$ $A$ 是$n\ast n$的常系数矩阵，称为系统矩阵； $B$ 是$n\ast r$的常系数矩阵，称为控制矩阵；$A$和$B$ 都是由系统本身的参数决定的。 $u$是输入信号向量，$x$是状态向量。

2、输出方程
$$y(t)=Cx(t)+Du(t)$$ $C$ 是$m\ast n$的常系数矩阵，称为输出矩阵，它表达了输出变量与状态变量之间的关系；$D$ 是$m\ast r$的常系数矩阵，称为直接转移矩阵，它表示输入变量通过矩阵$D$ 直接转移到输出。（注：大多数时候，$D=0$）

因此，对于线性定常连续系统，其状态空间描述为（或者状态空间表达式）：
$$\dot{x}(t)=Ax(t)+Bu(t)$$$$y(t)=Cx(t)+Du(t)$$ $x(t)\in R^n;u(t)\in R^r;y(t)\in R^m$.

类似地，
对于一般连续时间系统，其状态空间描述为：
$$\dot{x}(t)=f[x(t),u(t),t]$$$$y(t)=g[x(t),u(t),t]$$
对于线性时变连续系统，其空间描述为：
$$\dot{x}(t)=A(t)x(t)+B(t)u(t)$$$$y(t)=C(t)x(t)+D(t)u(t)$$ $x(t)\in R^n;u(t)\in R^r;y(t)\in R^m$.

2.2 线性定常连续系统的状态空间描述框图（开环）

根据线性定常连续系统的状态描述可以得到其一般性的状态空间描述框图，如下图：

2.3 线性定常连续系统的反馈控制

经典控制理论中，我们依据描述对象输入输出行为的传递函数模型来设计控制器，因此只能用系统的可测量输出作为反馈信号。而现代控制理论则是用刻画系统内部特征的状态空间模型来描述对象，除了可测量的输出信号外，还可以用系统的内部状态来作为反馈信号。根据可利用的信息是系统的输出还是状态，相应的反馈控制可分为输出反馈和状态反馈。
现代控制理论中，更多地使用状态反馈，由于状态反馈能提供更丰富的状态信息和可供选择的自由度，因而能使系统更容易获得更为优异的性能。其实，输出反馈是可以看做是部分状态反馈。

2.3.1 全状态反馈控制器

1 设计一个状态反馈控制器，如下图所示：
增加反馈环节K，使得闭环系统能够满足我们期望的系统性能：
$$u(t)=r(t)-Kx(t)\text{\hspace{1em}(1)}$$
已知系统的状态空间描述为：
$$\dot{x}(t)=Ax(t)+Bu(t)$$$$y(t)=Cx(t)+Du(t)\text{\hspace{1em}(2)}$$ $x(t)\in R^n;u(t)\in R^r;y(t)\in R^m$.
将$(1)$式代入$(2)$式得到闭环系统状态空间描述：
$$\dot{x}(t)=Ax(t)+B(r(t)-Kx(t))=(A-BK)x(t)+Br(t)$$$$y(t)=Cx(t)+D(r(t)-Kx(t))=(C-DK)x(t)+Dr(t)\text{\hspace{1em}(3)}$$
2 设计完反馈控制器的架构，下面要保证反馈系统的稳定性。
反馈系统稳定性的充要条件是系统闭环传递函数的所有极点均有负实部，即均在复频域S平面的左侧。
根据稳定性判定的条件，首先求闭环系统的传递函数。
为了书写一致性，我们重写系统状态表达式将$r$替换成$u$：
$$\dot{x}(t)=(A-BK)x(t)+Bu(t)$$$$y(t)=(C-DK)x(t)+Du(t)\text{\hspace{1em}(4)}$$
对系统状态空间表达式进行拉普拉斯变换：
$$sX(s)=(A-BK)X(s)+BU(s)$$$$Y(s)=(C-DK)X(s)+DU(s)\text{\hspace{1em}(5)}$$

求传递函数
$$\frac{Y(s)}{U(s)}=\frac{\ast }{sI-(A-BK)}\text{\hspace{1em}(6)}$$
可以得到系统的闭环传递函数的形式入式$(6)$所示。可见系统传递的极点就是矩阵$A-BK$的特征值。
因此，可以通过配置$K$，使闭环系统达到我们期望的状态。

问题是：当系统变量很多的时候，即使设计好了极点，但是矩阵$K$也不好计算。接下来开始引入LQR帮助求解$K$值。

3 LQR设计控制器的方法

3.1 什么是二次型

在不指定优化标准的前提下，控制领域中的“最优”体现在“输出能够完全跟踪控制，即在每一时刻输出量与控制量完全一致”。实际过程并不是这样完美的过程，每一时刻都会存在误差。退而求其次，追求在整个工作时间的范围内误差最小,与轨迹误差类似，我们研究状态误差。
因此，把整个工作时间内每一时刻状态的误差都累加起来，只要累加值更小，便会更加接近系统性能的期望。
1 首先，假设状态向量$x(t)$的维度为1以及闭环系统稳定。
令 $e(t)=x(t)-x^{\ast }(t)(t)$, $x(t)$为实际状态，$x^{\ast }(t)$为期望状态。
$$J={\int }_0^{\infty }|e(t)|dt={\int }_0^{\infty }|x(t)-x^{\ast }(t)|dt$$
显然，理想状态下$x(t)=x^{\ast }(t)$，$J$为零，达到最小值。非理想情况下$J$越小越好，$J$ 一般被称为代价函数（cost function）。

由于绝对值操作，不太方便数学运算，因此可以改写成平方的形式：
$$J={\int }_0^{\infty }e(t{)}^{2}dt={\int }_0^{\infty }(x(t)-x^{\ast }(t){)}^{2}dt$$
系统稳定之后，系统的状态处于平衡状态，即$x^{\ast }(t)$，不妨取 $x^{\ast }(t)=0$ 作为期望状态。
$$J={\int }_0^{\infty }(x(t){)}^{2}dt$$
2 扩充状态为$n$个，则代价函数为：
$$J={\int }_0^{\infty }{x}_1^2(t)+x_2^2(t)+...+x_n^2(t)dt$$
考虑到不同状态对累计误差所占的权重不同，赋予每个状态权重$w_i$:
$$J={\int }_0^{\infty }(w_1{x}_1{)}^2(t)+(w_2{x}_2{)}^2(t)+...+(w_n{x}_n{)}^2(t)dt$$
写成紧凑表达式为：
$$J={\int }_0^{\infty }(wx{)}^T(wx)dt={\int }_0^{\infty }{x}^T(w^{T}w)xdt$$

类似的$J$的函数称为二次型函数，变量的最高次数是2。
所有展开的函数最高次数为2的，这种类型的函数统称为二次型函数。

3.2 二次型最优控制

上一节介绍了二次型，这一节继续上一节的分析。
1、令$Q=w^{T}w$，则原来的代价方程改写为：
$$J={\int }_0^{\infty }{x}^{T}Qxdt$$
显然 Q是半正定矩阵。
思考：系统的状态发生的原因是在上一个状态时，有外界干扰或者系统输入发生改变引起的。外界干扰因素我们暂不考虑，说明代价函数不仅仅和系统的状态有关，而且和系统的输入有关。
因此，参考误差的处理方式，将控制量$u$也引入到代价函数里面，而且形式保持一致，具体如下：
$$J={\int }_0^{\infty }(x^{T}Qx+u^{T}Ru)dt\text{\hspace{1em}(7)}$$其中，Q和R均为自己设计的半正定矩阵。
式$(7)$也称为能量函数，二次型最优控制的目的就是使能量函数最小。

那如何设计反馈矩阵$K$使得能量函数$J$最小呢？

3.3 LQR 调节器设计

LQR控制器是计算反馈矩阵$K$主要分为三步：

1. QR矩阵选取；
2. 求解Riccati 方程得到矩阵P；
3. 计算$K=R^{-1}{B}^{T}P$；

3.3.1 QR矩阵选取

矩阵Q,R的选取，一般来说，Q值选得大意味着，要使得$J$小，那$x(t)$需要更小，也就是意味着闭环系统的矩阵(A-BK)的特征值处于S平面左边更远的地方，这样状态$x(t)$就以更快的速度衰减到0。另一方面，大的R表示更加关注输入变量$u(t)$,$u(t)$的减小，意味着状态衰减将变慢。

3.3.2 求解Riccati 方程得到矩阵$P$

我们求解的前提是假定系统处于稳定状态，此时的状态反馈为$u(t)=-Kx(t)$
将状态方程代入到代价函数中
状态方程：$\dot{x}(t)=(A-BK)x(t)$
代价函数：
$$\begin{gathered} J={\int }_0^{\infty }(x^{T}Qx+u^{T}Ru)dt \\ ={\int }_0^{\infty }{x}^T(t)(Q+K^{T}RK)x(t)dt \end{gathered}$$
我们要找到这个积分的原函数，为了找到$K$，假设存在一个常量矩阵$P$，使得
$$\frac{d}{dt}(x(t{)}^{T}Px(t))=-x^T(t)(Q+K^{T}RK)x(t)\text{\hspace{1em}(8)}$$
对式$(8)$两边取微分
$${\dot{x}}^T(t)Px(t)+x^T(t)P\dot{x}(t)+x^T(t)Qx(t)+x^T(t)K^{T}RKx(t)=0\text{\hspace{1em}(9)}$$
将状态方程代入$(9)$:
$$x^T(t)((A-BK{)}^{T}P+P(A-BK)+Q+K^{T}RK)x(t)=0$$
这个式子要成立的话，括号里的项必须为恒为零。
$$(A-BK{)}^{T}P+P(A-BK)+Q+K^{T}RK=0$$
进一步化简：
$$A^{T}P+PA+Q+K^{T}RK-K^T{B}^{T}P-PBK=0$$
取$K=R^{-1}{B}^{T}P$，代入上式得：
$$A^{T}P+PA+Q-PB{R}^{-1}{B}^{T}P=0\text{\hspace{1em}(10)}$$
$K$的二次项没有了，可$K$的取值和$P$有关，而P是我们假设的一个量，P只要使得的$(10)$式成立就行了。而$(10)$式在现代控制理论中极其重要，它就是有名的Riccati 方程。

3.3.3 计算$K=R^{-1}{B}^{T}P$

上个步骤中，可以利用工具箱求得矩阵P，最后根据公式计算出反馈矩阵$K$即可。