二次型最优控制(一)

本文来自潘神的知乎笔记什么是二次型最优控制

什么是二次型

说到“最”字，我们自然需要一个标准。比如，我们说到“最高”，可以用长度来丈量；说到“最重”，一般用质量来衡量，这些都是比较容易获得标准的。现在如果问，什么是“最美”？我们可能就得思考一会了。同样的道理，现在如果问什么是“最优”控制？——在没有“标准”的情况下讨论这个问题，就没有太大的意义。所以，我们要先来搞一个标准。既然我们搞控制，简单来说就是要“指东打东，指西打西”，翻译一下就是“输出能跟随输入”，那什么是“最优控制”呢？——一个很直观的想法就是：“输入能完全跟随输入，在每一个时间点上都是100%一样”，当然这是理想情况，实际上很难实现。那我们能不能拉开来看呢？——我们要有更大的视野，不要求每一点都一样，那在整个工作时间内是不是可以更接近呢？比如我们可以这样：把每一时刻的状态的误差都累加起来，只要累加的值更小，那是不是就更接近呢？比如我们可以令：
$$J={\int }_0^{\infty }|e|dt={\int }_0^{\infty }\left|x_o(t)-x_i(t)\right|dt$$
其中$x_o(t)$为实际状态，$x_i(t)$为期望状态。这样显然是可以的，如果$x_o(t)$和$x_i(t)$完全一样，$J$则为零，达到最小值，显然$J$ 越小越好——我们一般将$J$称之为代价函数（cost function）。

但是我们知道，绝对值这个操作，数学上不太好处理——因为它的导数不连续，会带来很多麻烦，那怎么办？——聪明的你可能已经想到了，可以用平方啊，它也能表征距离，而且性能非常好。 这样代价函数就变成了：
$$J={\int }_0^{\infty }{e}^{2}dt={\int }_0^{\infty }{\left(x_o(t)-x_i(t)\right)}^{2}dt$$
为简单起见，我们先假设期望状态$x_i(t)$为零，即
$$J={\int }_0^{\infty }{e}^{2}dt={\int }_0^{\infty }{\left(x_o(t)\right)}^{2}dt$$

我们看一下被积分的部分，这只有一个状态，我们可以扩充一下，假如现在有$n$个状态，即
$$J={\int }_0^{\infty }{x}_1^2(t)+x_2^2(t)+\dots +x_n^2(t)dt$$

这样不够紧凑，我们换个写法：

$$J={\int }_0^{\infty }{\mathbf{x}}^T\mathbf{x}dt$$
其中

$$\begin{array}{c}{\mathbf{x}}^T\mathbf{x}=\left[x_1,x_2,\dots x_n\right]\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]\\ =x_1^2+x_2^2+\dots +x_n^2\end{array}$$

这就是一种简单的二次型函数，因为变量的最高次数是2。

当然，我们的工具可以更强大，比如单位矩阵$I$可以扩展为对角矩阵$\lambda$：
$$\mathbf{\lambda }=\left[\begin{array}{cccc}{\lambda }_1& 0& \dots & 0\\ 0& {\lambda }_2& \dots & 0\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & {\lambda }_n\end{array}\right]$$
这相当于给各个变量误差取个权重，然后对角矩阵$\lambda$还可以扩展为更一般的矩阵$Q$：
$$\mathbf{Q}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \dots & a_{1n}\\ a_{21}& a_{22}& \dots & a_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ a_{n1}& a_{n2}& \dots & a_{nn}\end{array}\right]$$
但无论怎么扩展，展开的函数最高次数都是2次，因此我们把这种类型的函数统称为二次型函数。

二阶系统的阻尼比为何选择0.707

一个简单例子，二阶系统，框图如下，现在加入我们要实现最优控制，看看最佳阻尼比应该是什么样。

其开环传递函数为

$$\frac{C(s)}{E(s)}=\frac{1}{s(s+2\zeta )}$$
误差为：
$$E(s)=R(s)-C(s)$$
将上面两个传递函数变成时域方程为：

$$\begin{array}{l}\ddot{c}+2\zeta \dot{c}=e\\ e=r-c\end{array}$$
联立以上两式就可以得到
$$\ddot{e}+2\zeta \dot{e}+e=\ddot{r}+2\zeta \dot{r}$$
传递函数只能描述单输入单输出系统，现在我们要把它改成状态方程的形式：

$$\begin{array}{rl}{x}_1& =e\\ x_2& =\dot{e}\end{array}$$

也就是我们要研究误差和误差变化率，为简单起见，我们假设输入为阶跃信号，即$\ddot{r}=0$，$\dot{r}=0$，$r=1$，对应的$e=1$，$\dot{e}=0$。这样状态方程就变成
$$\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}$$
其中
$$A=\left[\begin{array}{cc}0& 1\\ -1& -2\zeta \end{array}\right]$$
完整状态方程为
$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -1& -2\zeta \end{array}\right]\left[\begin{array}{l}{x}_1\\ x_2\end{array}\right]$$

我们定义代价函数为：

$$J={\int }_0^{\infty }\left(e^2+{\dot{e}}^2\right)dt$$

即两个变量之间同等权重，用变量表示

$$J={\int }_0^{\infty }\left(x_1^2+x_2^2\right)dt$$

现在我们想要$J$最小，具体应该怎么做呢？先将代价函数写成紧凑一点的形式：
$$J={\int }_0^{\infty }{\mathbf{x}}^T\mathbf{x}dt$$

这是个积分的问题，所以我们要先找到被积函数 ${\mathbf{x}}^T\mathbf{x}$的原函数，这个数学问题不是我们的重点，假设我们现在已经知道原函数的长相为 ${\mathbf{x}}^T\mathbf{P}\mathbf{x}$，即

$$\frac{d}{dt}\left({\mathbf{x}}^T\mathbf{P}\mathbf{x}\right)=-{\mathbf{x}}^T\mathbf{x}$$

其中 $\mathbf{P}$是一个待定的实对称矩阵，我们把它看成待定系数，现在要研究一下这系数需要满足什么条件，我们把上式的左边微分展开：

$$\frac{d}{dt}\left({\mathbf{x}}^T\mathbf{P}\mathbf{x}\right)={\dot{\mathbf{x}}}^T\mathbf{P}\mathbf{x}+{\mathbf{x}}^T\mathbf{P}\dot{\mathbf{x}}$$

把状态方程 $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}$代入上式，于是我们就得到：

$$\frac{d}{dt}\left({\mathbf{x}}^T\mathbf{P}\mathbf{x}\right)={\mathbf{x}}^T\left({\mathbf{A}}^T\mathbf{P}+\mathbf{P}\mathbf{A}\right)\mathbf{x}$$

如果我们保证 ${\mathbf{A}}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{I}$，由于 $A$ 是已知的，这样我们就能确定 $P$，也就能确定原函数 ${\mathbf{x}}^T\mathbf{P}\mathbf{x}$ ，此时

$$J={\int }_0^{\infty }{\mathbf{x}}^T\mathbf{x}dt=-{\mathbf{x}}^T\mathbf{P}{\mathbf{x}|}_0^{\infty }={\mathbf{x}}^T(0)\mathbf{P}\mathbf{x}(0)-{\mathbf{x}}^T(\infty )\mathbf{P}\mathbf{x}(\infty )$$

我们的系统是稳定的，而且输入为零，这样我们就可以得到 ${\mathbf{x}}_{\infty }=0$ ，所以代价函数可以进一步简化为：

$$J={\mathbf{x}}^T(0)\mathbf{P}\mathbf{x}(0)$$

其中$P$满足${\mathbf{A}}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{I}$(这其实也是李雅普诺夫第二法)

注意：所有的二次型最优控制，形式可能会很复杂，但是处理思路和方式都是这样。

$P$使用待定系数法求解。得到：

$$\mathbf{P}=\left[\begin{array}{cc}\zeta +\frac{1}{2\zeta }& \frac{1}{2}\\ \frac{1}{2}& \frac{1}{2\zeta }\end{array}\right]$$
将$\mathbf{x}(0)={\left[\begin{array}{ll}1& 0\end{array}\right]}^T$带入$J={\mathbf{x}}^T(0)\mathbf{P}\mathbf{x}(0)$
得到代价函数：
$$J=\zeta +\frac{1}{2\zeta }$$

现在我们要计算 $J$的极小值，先对变量微分

$$\frac{dJ}{d\zeta }=1-\frac{1}{2{\zeta }^2}$$
当$\frac{dJ}{d\zeta }=0$时，可以计算得到：
$$\zeta =\frac{\sqrt{2}}{2}\approx 0.707$$
也就是说，当阻尼比取0.707时，$J={\int }_0^{\infty }\left(e^2+{\dot{e}}^2\right)dt$这个损失函数取极小值——也就是说从整体来看，误差和误差变化率累计最小，这就是为什么我们一般把二阶系统的阻尼比设置在0.707左右的原因。

代价函数矩阵需要由$I$变成$Q$

前面选择的代价函数为

$$J={\int }_0^{\infty }\left(x_1^2+x_2^2\right)dt$$

我们认为状态变量的权重都一样，没有谁更重要。实际上，并不都是如此，比如还是前面的例子，有时候我们可能更关心误差，对误差的变化率容忍度高一些，这时候我们就可以给它们分配一个权重：

$$J={\int }_0^{\infty }\left(x_1^2+\mu x_2^2\right)dt\mu >0$$

此时代价函数可以演变为

$$\begin{array}{rl}J& ={\int }_0^{\infty }\left(e^2+\mu {\dot{e}}^2\right)dt={\int }_0^{\infty }\left(x_1^2+\mu x_2^2\right)dt\\ & ={\int }_0^{\infty }\left[x_1{x}_2\right]\left[\begin{array}{cc}1& 0\\ 0& \mu \end{array}\right]\left[\begin{array}{l}{x}_1\\ x_2\end{array}\right]dt\\ & ={\int }_0^{\infty }{\mathbf{x}}^T\mathbf{Q}\mathbf{x}dt\end{array}$$
其中
$$\mathbf{Q}=\left[\begin{array}{ll}1& 0\\ 0& \mu \end{array}\right]$$
同样前面的处理方式一样，代价函数的表达式仍然为
$$J=\mathbf{x}(0{)}^T\mathbf{P}\mathbf{x}(0)$$

不过此时$P$需要满足不再是 ${\mathbf{A}}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{I}$，而是 ${\mathbf{A}}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{Q}$
还是前面的状态方程：

$$\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\end{array}\right]=\underset{A}{\left[\begin{array}{cc}0& 1\\ -1& -2\zeta \end{array}\right]}\left[\begin{array}{l}{x}_1\\ x_2\end{array}\right]$$

此时可计算满足 ${\mathbf{A}}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{Q}$ 的$P$为

$$\mathbf{P}=\left[\begin{array}{cc}\zeta +\frac{1+\mu }{4\zeta }& \frac{1}{2}\\ \frac{1}{2}& \frac{1+\mu }{4\zeta }\end{array}\right]$$

将初始条件$\mathbf{x}(0)={\left[\begin{array}{ll}1& 0\end{array}\right]}^T$带入$J=\mathbf{x}(0{)}^T\mathbf{P}\mathbf{x}(0)$，即可得到：

$$J=\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{cc}\zeta +\frac{1+\mu }{4\zeta }& \frac{1}{2}\\ \frac{1}{2}& \frac{1+\mu }{4\zeta }\end{array}\right]\left[\begin{array}{l}1\\ 0\end{array}\right]=\zeta +\frac{1+\mu }{4\zeta }$$

可见，此时代价函数取决于阻尼比 $\zeta$和输入的权重 $\mu$ ，为了计算 $J$的极值，我们同样要计算$J$ 对 $\zeta$ 的微分
$$\frac{\partial J}{\partial \zeta }=1-\frac{1+\mu }{4{\zeta }^2}$$
当$\frac{\partial J}{\partial \zeta }=0$时，$\zeta =\frac{\sqrt{1+\mu }}{2}$
显然，当$\mu =1$时，$\zeta =\frac{\sqrt{2}}{2}$，这与我们前面的计算结果一致。如果我们现在增加 $\mu$ ，也就意味着误差变化率 $\dot{e}$的权重更大，此时也就意味着我们需要增加阻尼比 $\zeta$。